© ESO 2022

1. Introduction

In their late stages of evolution, prior to reaching the white dwarf cooling track, most low- and intermediate-mass stars go through the asymptotic giant branch (AGB) phase (e.g., Iben & Renzini 1983; Catelan & Smith 2015, and references therein). When such stars undergo high-amplitude brightness variations, they are generally classified as long-period variables (LPVs), given the long periods, in excess of about 80 d, that are typically observed. More massive stars may also undergo LPV pulsations while on the red supergiant and/or super-AGB phase, and thus at least some LPVs may also be core-collapse supernova progenitors (Levesque 2017; O’Grady et al. 2020, and references therein). LPV-like pulsations have also been observed in other, more exotic objects whose physical nature is still not fully understood (O’Grady et al. 2020, and references therein).

Based on their light-curve morphology, the range of their amplitudes and their periodicity, LPVs are commonly subgrouped into four categories: Miras, semi-regular variables (SRVs), irregulars (Irrs), and OSARGs (OGLE small-amplitude red giants). A recent review was provided by Catelan & Smith (2015), to which we refer the reader for further details and references. In our work, we mainly focus on Miras presenting periods of 80–1400 days and amplitudes greater than approximately 0.4 mag in the K_s-band. Their initial masses have a lower limit of about 0.8–1 M_⊙ and an upper limit of 6–9 M_⊙, with ages between ∼10⁸ and ∼10¹⁰ yr (Iben & Renzini 1983; Catelan & Smith 2015).

Since Glass & Evans (1981) and Feast et al. (1989) first presented the Miras period-luminosity relation, this class of variable stars has been used as distance indicators, and thus as tracers of the structure of the different Milky Way components (e.g., Catchpole et al. 2016; Grady et al. 2020; Urago et al. 2020). Recently, Gaia’s Data Release 2 (DR2) provided a catalog of LPV candidates based on Gaia data collected over a time span of 22 months (Mowlavi et al. 2018). The Gaia DR2 LPV catalog includes both Miras and SRVs with variability amplitudes larger than 0.2 mag in the Gaia G-band. Accurate parallaxes are also included, which can be used to obtain their distances (but see also Whitelock 2012). On the other hand, Miras also display a period-age relation (Feast & Whitelock 2000b; Feast 2007), which, together with the fact that they span a broad range in ages and follow a period-luminosity relation, can make them useful probes not only of Galactic structure but also of the Galaxy’s evolution. Since ages thus determined depend on the period, they are insensitive to systematic errors that may affect the reddening or distances. Even relative age estimates for an abundant variable population can supply constraints on the existence of age gradients. Using Miras as age indicators provides an opportunity not only of mapping the Galactic disk and bulge as a function of age but also the halo structure and even farther away, since Miras are intrinsically bright and detectable out to large distances.

Miras are classified into two subclasses based on their surface C/O abundance ratio: O-rich Miras, presenting H₂O, TiO, and SiO on their surfaces, and C-rich Miras, presenting C₂ and CN, with the C being dredged up from the stellar interior to the surface. In the Galactic bulge, the Mira population is dominated by the O-rich subclass, with only a few C-rich Miras being known (Whitelock et al. 2006; Matsunaga et al. 2017; Grady et al. 2020). In the lower-metallicity environments of the Magellanic Clouds, on the other hand, the proportion of C-rich Miras is larger (e.g., Whitelock et al. 2003).

It has long been known empirically that O-rich Mira kinematics is correlated with their periods in the Galactic disk. Thus, the relationship between age and velocity dispersion calibrated for stars in the solar neighborhood can be used to derive a period-age relationship for nearby O-rich Miras (Feast & Whitelock 2000b; Feast 2007). Recently, Grady et al. (2019) introduced a period-age relation for Miras based on a tentative association between some of these stars and star clusters. Their resulting relation presented a large offset in comparison to other expressions available in the literature (Feast 2007; López-Corredoira 2017) dealing with ages versus kinematics. That notwithstanding, all these works consistently find that longer-period Miras are younger, as expected since the period of a Mira increases with increasing luminosity, and hence mass.

To study long-term variability phenomena in the bulge, a long-baseline near-infrared (IR) survey is desired. The IR regime allows one to overcome the high extinction levels along these sight lines, and the long baseline affords multiepoch photometry covering long variability timescales. Prior near- and mid-IR surveys, such as the Two-Mass All-Sky Survey (2MASS; Skrutskie et al. 2006), DENIS (Deep Near Infrared Survey of the Southern Sky; Epchtein et al. 1999), GLIMPSE (Galactic Legacy Infrared Mid-Plane Survey Extraordinaire; Benjamin & GLIMPSE Team 2005), and UKIDSS-GPS (United Kingdom Infrared telescope Deep Sky Galactic Plane Survey; Lucas et al. 2008), could alleviate the problem of high extinction in the inner part of the Milky Way (MW; see also, e.g., Groenewegen & Blommaert 2005; Jiménez-Esteban & Engels 2015), but still the absence of multiepoch observations hindered detailed analyses of variable stars, including LPVs. On the other hand, optical long-baseline surveys of the bulge region, such as OGLE (Optical Gravitational Lensing Experiment; Udalski et al. 2015), MACHO (Massive Compact Halo Objects survey; Alcock et al. 1996), EROS (Expérience pour la Recherche d’Objets Sombres; Aubourg et al. 1993), and Gaia (Global Astrometric Interferometer for Astrophysics; Mowlavi et al. 2018), provided a great demonstration of the power of multiepoch observations, with a plethora of important results for numerous different types of variable stars. On the other hand, the inner part of the bulge could not be penetrated due to the severe extinction in that region. Thus, even in the Gaia era, IR observations remain crucial to properly understand the stellar populations of the innermost MW. This is especially true in the case of faint stars or variables close to minimum light.

The VISTA Variables in the Vía Láctea (VVV) ESO Public Survey (Minniti et al. 2010) was the first systematic, wide-field near-IR variability survey of the MW bulge and inner disk. A total of 1929 h of observations were obtained between 2010–2016, yielding on average about 100 datapoints per field, and a total sky coverage of order 520 deg². Other near-IR variability studies of MW bulge fields, by comparison, have sometimes covered sky areas not exceeding 1 deg² (e.g., Matsunaga et al. 2009). In addition to its extended sky coverage, the VVV images offer significantly improved resolution and higher depth by up to several magnitudes, as compared to 2MASS (e.g., Saito et al. 2012).

In this paper, we explore the VVV dataset in order to extract LPVs in regions presenting overlap with the OGLE footprint. We provide a first catalog of LPV candidates in this field, including periods and (when possible) distances. In future papers of this series, we will expand our search toward other regions covered by the VVV survey, with the ultimate goal of building a comprehensive catalog of LPVs covering the entire VVV footprint and studying their spatial distribution and ages.

2. Observations and data reduction

Observations were taken with the 4 m Visible and Infrared Survey Telescope for Astronomy (VISTA) at ESO’s Cerro Paranal Observatory, northern Chile, using the wide-field VISTA InfraRed Camera (VIRCAM; Dalton et al. 2006; Emerson & Sutherland 2010) in the time span of 2010–2016. The data consist of multiepoch K_s images and a couple of observations in each of the ZYJH bandpasses per field. Each exposure (commonly referred to as a “pawprint image”) is taken with a mosaic of 16 CCDs covering 2048 × 2048 pixels each at 0.34″ resolution, with exposure time 5 s. The seeing value is typically around 1″. Multiple pawprint exposures of a given field are later combined to generate so-called “tile images”, which (by design) remove the gaps that are present in between individual CCDs in the pawprint images. Each such tile covers 1.64 deg² (Minniti et al. 2010). Disk tile names are labeled by a d followed by an ID number ranging between 001 and 152; bulge tiles are instead denoted by a b, followed by an ID number in the range 201–396 (e.g., Catelan et al. 2011; Saito et al. 2012).

Due to the very high stellar density in our studied fields, we consider point-spread function (PSF) photometry to be superior to aperture photometry, as adopted in some previous variability studies that also used VVV data (e.g., Saito et al. 2013; Teixeira et al. 2018; Dékány et al. 2019; Ferreira Lopes et al. 2020; Herpich et al. 2021). Accordingly, PSF photometry was obtained for all sources by fitting appropriate models on individual pawprint images, using the DAOPHOT/ALLSTAR software package (Stetson 1987). The seeing of the images used in this work was on average of order 0.8 ± 0.2″, and images with seeing poorer than 1.7″ were discarded from further analysis. The individual images were aligned using the DAOMASTER routine, which computes the astrometric transformation equation coefficients from the DAOPHOT/ALLSTAR results. All available images were combined using the MONTAGE2 routine to make the master image of the area. More details can be found in Contreras Ramos et al. (2017).

In this paper, we focused on VVV tiles b302, b305–b309, b319, b320, b364, and b365, which are covered by our PSF photometry and which present maximum overlap with the OGLE survey, in the sense that a minimum of 50 OGLE Miras are already known in each tile (Table 1). This choice allows us to check and refine our search procedures, which will be extended to other VVV tiles in the near future. The OGLE data were taken with the 1.3 m Warsaw telescope at Las Campanas Observatory, northern Chile, in the course of the OGLE-III¹ campaigns (Udalski 2003; Soszyński et al. 2013b), and are available online². Our resulting LPV catalog covering these ten tiles is presented in this paper; analysis of the remaining VVV tiles is currently in progress and will be the subject of forthcoming papers.

Figure 1 provides a schematic map of the innermost region of the MW, with the locations of the centers of VVV tiles indicated by their respective IDs and the region for which our PSF photometry is completed indicated by the large rectangle and grid at −10.0 < ℓ(deg) <  + 10.0, −2.5 < b(deg) <  + 2.8. In this plot, red dots represent stars classified as Miras by the OGLE-III project. As can be readily inferred from this plot, tiles b302, b305–b309, b319, b320, b364, and b365 (highlighted in blue) are those with the largest number of previously known Miras, which is the reason why we specifically select them for this first paper in this series (see Table 1). According to our PSF photometry, these tiles comprise an initial sample of 61 627 100 individual sources with at least 20 VVV K_s-band observations each, which is the minimum number we adopt to ensure we only use sufficiently well-sampled light curves in our analysis.

Fig. 1. Footprint of the VVV survey in the central bulge of the MW, expressed in Galactic coordinates. Each VVV tile is marked by its respective ID in this diagram, those analyzed in this work being highlighted in blue and using a larger font size. Red dots represent stars classified as Miras in the OGLE-III project. The gray grids present a region of the VVV footprint where our PSF photometry is already completed.

3. Variable star search

In our approach, and in order to reduce the impact of statistical fluctuations that may impact some indices more than others, we use a combination of both correlated and noncorrelated statistical indices (see Ferreira Lopes et al. 2015a, 2018; Ferreira Lopes & Cross 2016, 2017) to carry out our search for bona-fide variables. Specifically, our candidates were extracted according to the following procedure:

1. Outlier removal: for each star, measurements outside five standard deviations from the median value were excluded from our analysis. In addition, measurements with unusually large associated error bars were also discarded. To accomplish the latter, we first computed the mean error bar and the corresponding standard deviation, using all available data points for the star. In those cases where the individual photometric error was higher than five standard deviations from the mean error bar, the measurement was discarded.

2. The reduced χ² statistic, ${\chi }_{\nu }^2$, was computed with respect to a constant model (Bellm 2021). This quantity is defined as

$$\begin{array}{c}\hfill {\chi }_{\nu }^2=\frac{1}{N-1}{\displaystyle \sum _{i=1}^N\frac{{(m_i-\overline{m})}^2}{{\sigma }_i^2},}\end{array}$$(1)

where $\overline{m}$ represents the weighted mean magnitude, σ_i is the error associated with magnitude m_i, and N is the number of valid observations. This is a commonly used noncorrelated index that is frequently employed to select Miras and other high-amplitude variables from photometric databases (e.g., Carpenter et al. 2001; Braga et al. 2019).

3. In parallel, the correlated index $K_{(fi)}^{(2)}$ was computed. It is defined as follows:

$$\begin{array}{c}\hfill K_{(fi)}^{(2)}=\frac{N^{+}}{N_C},\end{array}$$

where N⁺ is the number of positive correlations and N_C is the total number of correlations, for an assumed ΔT value – the latter quantity corresponding to the time interval between observations. In our case, the correlation δ_i, j between a pair of measurements i, j is defined as follows:

$$\begin{array}{c}\hfill {\delta }_{i,j}=\left(\frac{m_i-\overline{m}}{{\sigma }_i}\right)\times \left(\frac{m_j-\overline{m}}{{\sigma }_j}\right)={\mathrm{\Lambda }}_{i,j}\frac{|m_i-\overline{m}|}{{\sigma }_i}\times \frac{|m_j-\overline{m}|}{{\sigma }_j},\end{array}$$(2)

where $\overline{m}$ denotes the mean magnitude value, σ the error bar, and Λ_i, j reflects the sign of the multiplication, i.e., +1 or −1. According to this definition, δ_i, j describes the degree to which two measurements vary together in a time series. The $K_{(\mathrm{fi})}^{(2)}$ index is computed as the mean value over all (N_C) i, j combinations, i.e., $K_{(fi)}^{(2)}=\sum {\mathrm{\Lambda }}_{i,j}/N_C=N^{+}/N_C$. The $K_{(\mathrm{fi})}^{(2)}$ has a weak dependence on outliers, and an accurate estimation is obtained when only observations carried out close in time are considered, as done in this work (for more details, see Sect. 3 in Ferreira Lopes & Cross 2016). While correlated indices in general, and the one used here in particular, are at least three times more efficient than noncorrelated indices (Ferreira Lopes et al. 2015a, 2018; Ferreira Lopes & Cross 2016, 2017), they only work if certain conditions are satisfied (Ferreira Lopes & Cross 2016). In particular, the ratio between ΔT and the variability period P should be smaller than 0.1 (ΔT/P ≪ 0.1; see Fig. 2 in Ferreira Lopes & Cross 2016), in order to ensure accurate correlated values. To satisfy this criterion, in our work ΔT was increased until the number of correlated measurements reached 10, and ΔT = 5 d was set as an upper limit. This ensured ΔT/P ≪ 0.1, since Miras have periods typically in excess of ∼80 days. Indeed, ΔT can be quite large, if one considers the long periods typically found for LPVs. Such an approach allows us to compute the $K_{(\mathrm{fi})}^{(2)}$ index for about 97% of the initial sample. Ferreira Lopes & Cross (2016) also defined an empirical relation to determine the expected value for the noise (Fap_K(fi), their Eq. (16)), based on which unreliable index values can be identified. Indeed, those stars having periods shorter than ΔT can mimic noise according to those authors. Therefore, our approach improves the detection of LPVs and also reduces the number of stars having short periods.

4. All sources having $K_{(\mathrm{fi})}^{(2)}/Fap$ and ${\chi }_{\nu }^2$ inside the region depicted by the dashed line in Fig. 2 were then selected for further analysis. The borders of this selection region were adjusted manually in order to provide an optimum balance between the number of bona-fide LPV candidates and spurious sources (for more details, see Appendix 6). This leads to a preliminary sample having 365 696 candidates. The minimum cutoff in ${\chi }_{\nu }^2$ is smaller than typically used to select high-amplitude Miras, but was adopted in order to ensure detection of smaller-amplitude LPVs, as in the case of SRVs, and Miras whose amplitudes might be spuriously reduced in our photometry on account of saturation (see Sect. 3.2). On the other hand, the $K_{(\mathrm{fi})}^{(2)}/Fap>1$ line is expected to enclose all nonsaturated LPVs.

Fig. 2. Distribution of VVV-PSF sources in the $K_{(\mathrm{fi})}^{(2)}$ vs. ${\chi }_{\nu }^2$ plane. The Miras identified by Matsunaga et al. (2009, labeled MiraMat in the plot) are shown as blue circles. The region marked by dashed lines corresponds to the one from which we extracted our initial Mira candidates. The latter, upon further analysis and inspection, eventually led to our final selection of Miras, shown in this plot as orange squares. The maximum number of targets per pixel is 61 697 in this plot, corresponding to the dark orange-red region in the main panel.

5. The periods of variable star candidates were estimated using three different methods, namely phase dispersion minimization (PDM; Stellingwerf 1978), Lomb-Scargle (Lomb 1976; Scargle 1982), and generalized Lomb-Scargle (Zechmeister & Kürster 2009). A frequency interval between 0.5 cycles/day (corresponding to a period of 2 days) and 1/range(t), where range(t) is the total time span in days, was adopted. In most cases, the periods found by the different methods were in good agreement. In the end, upon visual inspection of the results, the Lomb-Scargle solutions were finally adopted.

6. The periods thus estimated were used to compute the signal-to-noise ratio (S/N), defined as the ratio of the amplitude (ΔK, measured peak-to-peak from the model) and the standard deviation of the residuals. The model was obtained using a Fourier series having three harmonics:

$$\begin{array}{c}\hfill m(t)=m_0+{\displaystyle \sum _{j=1}^3{a}_{j}cos(\frac{2\pi jt}{P}+{\theta }_j),}\end{array}$$(3)

where a_j and θ_j are the Fourier amplitude and phase coefficients, respectively, and m₀ is the star’s mean magnitude (e.g., Debosscher et al. 2007; Ferreira Lopes et al. 2015b,c; Matsunaga et al. 2017).

7. The period and its respective power, together with the S/N estimated according to the procedure described in the previous step, were used to select our targets, which were required to satisfy one of the following criteria:

• (a)

  Power > 20 and S/N > 1;

• (b)

  9< Power < 20 and S/N > 2;

• (c)

  5< Power < 9 and S/N > 3;

• (d)

  Power > 10, S/N > 1, and ⟨K_s⟩< 12.

We note that criterion (d) includes only sources that did not pass criteria (a) or (b). In addition, a minimum period of 80 days was also imposed. The power cutoff value is motivated by Horne & Baliunas (1986), whose Eq. (13) indicates that, for a number of observations exceeding ∼20 epochs, a Fap of ∼20%, ∼2%, ∼0.2%, ∼0.0% results for the power used in the (a), (b), (c), and (d) criteria, respectively (see horizontal dashed line in Fig. 3). Compared with criterion (a), criteria (b) and (c) relax the constraint on the periodogram power, provided a sufficiently large S/N is still present in the light curve. Finally, compared to criterion (a), lower limits in the periodogram power and S/N were adopted for sources having many saturated observations (⟨K_s⟩< 12) but still clearly defined variability, as saturated measurements can reduce the S/N in particular to values smaller than two. We emphasize that these criteria are sufficiently lax that bona-fide Miras, including fainter ones with large photometric errors, are unlikely to be missed at this stage, even if at the expense of a large number of contaminants. As a result, 12 476 sources were selected; we call these our VVV-PSF LPV candidates.

Fig. 3. Logarithm of the power as a function of the period for the selected VVV-PSF sources (see region delimited by black dashed line in Fig. 2). Symbols and colors are the same as in Fig. 2) as well. The maximum number of targets per pixel is 3146 in this plot, corresponding to the dark orange-red region in the main panel.

8. Given that, as already stated, the aforementioned criteria imply the presence of a significant number of contaminants, visual inspection was then performed, in order to remove the most obvious sources that are unlikely to be bona-fide LPVs. As a consequence, 11 333 of these 12 476 sources were considered to have variability signals that were either unclear or deemed unlikely to correspond to LPV stars. The remaining 1143 stars were finally assigned to two subgroups, namely “Mira candidates” (130 stars) and “LPV⁺ candidates” (1013 sources). The latter correspond mainly to those cases where period determination was particularly problematic (for more details, see Sect. 3.2). Table 2 provides a summary with the number of selected targets according to each of the criteria mentioned above.

3.1. Photometry recalibration

Our PSF photometry was originally calibrated to the VISTA photometric system by CASU (González-Fernández et al. 2018). CASU provides ZYJHK_s aperture photometry catalogs directly calibrated from 2MASS point source data. We calibrated our instrumental PSF photometry in each chip through the comparison of common sets of stars with respect to CASU, which would place our photometry in the 2MASS photometric scale. However, as pointed out by Hajdu et al. (2020), the standard CASU calibration procedure has some issues biasing the photometric scale, especially in regions of high stellar density. Therefore, we decided to recalibrate our PSF CASU-system photometry directly from 2MASS by carefully choosing the set of standards used in the procedure. We closely followed the general recommendations laid out in Hajdu et al. (2020), but since we introduced a few modifications and additions, we briefly describe our procedure here.

For each chip catalog, we select stars with JHK_s magnitudes (according to our PSF CASU-tied photometry) in the range of 12.0–15 mag. This range was adopted because we want to find matches in an equivalent (but shorter) range of 2MASS magnitudes, keeping in mind the fact that our zero points might still be wrong by up to ±0.2 mag or so at this stage (Hajdu et al. 2020). Indeed, the presence of these potential errors in the zero points are the motivation behind this recalibration. We remove all sources with companions closer than 3″ (corresponding to the mode value of the 2MASS PSF), unless the companions are all fainter than 8 times the faint limit of the initial selection. This eliminates sources which might be blended with a significant (in terms of brightness) contaminant in 2MASS, and completes the definition of our initial catalog. Next, we cross-match this VVV selection with the 2MASS point source catalog, adopting a stringent radius of 0.4″. We keep all sources with 2MASS photometric quality flag “AAA”, with color J − K_s < 2, and with an extinction-corrected color of 0.0 < J − K_s < 1.0, according to the 2MASS photometry (see below for a discussion of the adopted extinction values).

A few more cuts are needed to define the calibration sample. They require to transform the retrieved 2MASS photometry to the VISTA photometric system. To this end, we adopt the version 1.5 set of transformation equations defined by CASU:

$$\begin{array}{c}\hfill J_{2\mathrm{M}-\mathrm{VVV}}=J_{2\mathrm{M}}-0.031(J_{2\mathrm{M}}-K_{2\mathrm{M}}),\end{array}$$(4)

$$\begin{array}{c}\hfill H_{2\mathrm{M}-\mathrm{VVV}}=H_{2\mathrm{M}}+0.015(J_{2\mathrm{M}}-K_{2\mathrm{M}}),\end{array}$$(5)

$$\begin{array}{c}\hfill K_{s,2\mathrm{M}-\mathrm{VVV}}=K_{2\mathrm{M}}-0.006(J_{2\mathrm{M}}-K_{s,2\mathrm{M}})+0.005E{(B-V)}_{\mathrm{corrected}},\end{array}$$(6)

where 2M stands for 2MASS. As it can be seen, we need a color excess term, E(B − V)_corrected, to compute the K_s-band transformation. Following González-Fernández et al. (2018), we retrieve E(B − V) at the coordinates of each star from the extinction map of Schlegel et al. (1998), and in all cases where the extinction-corrected (J − K_s) color of the star is ≤0.25, we recompute the extinction E(B − V)_corrected by imposing (J − K)₀ = 0.5. We note that, since the extinction term in Eq. (6) is very small, the impact of this choice upon the derived magnitudes is exceedingly small in K_s and, according to Eqs. (4) and (5), strictly zero in the cases of J and H (see also Hajdu et al. 2020).

To avoid VISTA nonlinearity and/or saturation issues for bright targets, at this stage we remove stars with J, H, K_s < 12.5, according to their 2MASS magnitudes converted to the VVV system (see also Contreras Ramos et al. 2017). This completes the selection of the calibration sample. However, if the several applied cuts yield a sample of less than 25 stars, we remove the color cut on 2MASS photometry; if this is still not enough, we remove the cut on the reddening-corrected color as well.

Finally, the photometric zero-points in each of the J, H, and K_s bands are computed from the median of the residuals of the respective mag_{2M − VVV} − mag_VVV distribution. The zero-points are added to the input original PSF magnitudes to cast them into the corrected 2MASS system. On the other hand, if the calibration sample had less than 3 stars, no zero-points were computed for that chip’s catalog. We arrived at that situation only in a handful of cases, and therefore, for these chips, we interpolated zero-points from the respective (ℓ,b) zero-point maps, as the spatial variations of the zero-points are sufficiently smooth.

3.2. LPV selection

Steps 1 to 4 in the procedure described in the beginning of this section were applied in order to ensure the recovery of the majority of the LPVs contained in our initial sample. About 95% of the OGLE and Matsunaga et al. (2009) Miras follow these criteria. As a result, 0.6% of the entire VVV initial sample was selected as candidates (see Fig. 2). Steps 5–8 were then applied to evaluate the reliability of the variability signal. We note that, since LPVs are intrinsically very bright stars, many fall in the nonlinear regime of the VIRCAM detectors, and thus turned out to be saturated in the VVV data. We opted not to reject all of these saturated sources (i.e., with K_s ≲ 12 mag; Contreras Ramos et al. 2017). While our PSF photometry (Sect. 2) allowed us to recover and measure these saturated LPVs, the photometry for these sources (and hence their reported mean magnitudes and amplitudes) should of course be treated with extreme caution. On the other hand, we have found that their variability periods can often still be properly recovered (an example is shown in Fig. 4), which is the main reason we include them in our catalog. Saturated sources with unclear periodicity signals, on the other hand, were removed from our sample. We note that the ages of even the saturated stars can still be estimated, even when reliably measured mean magnitudes are not available, as long as their periods can still be computed on the basis of the available measurements (see Sect. 6.4 below). We encourage follow-up studies of these sources, so that their photometric properties can also be properly established, and thus estimates of their distances obtained. We note that numerous LPVs brighter than K_s = 10 are known to exist in the direction of the Milky Way bulge; in particular, Qin et al. (2018) show, using 2MASS data, that they cover the range K_s ≈ 4 − 10, and peak around K_s ≈ 7.5 (see also Groenewegen & Blommaert 2005, who arrived at similar conclusions using 2MASS and DENIS near-IR measurements). The vast majority of these very bright Miras are inaccessible to VVV, and thus remain outside of the scope of this series of papers.

Fig. 4. Example of raw (left) and phased (right) light curves, for a source cross-matched between OGLE (top panel) and VVV (bottom panel). The OGLE and VVV IDs, as well as the period found by OGLE, are displayed in the top panel.

An amplitude criterion was then applied in order to distinguish Miras. In the visual, Miras are defined as LPVs whose amplitudes exceed 2.5 mag in V. A similar threshold is however not well-established in the near-IR, although the majority of Miras have been shown to have K-band amplitudes exceeding 0.4 mag (Whitelock et al. 2000). Matsunaga et al. (2009), in turn, used amplitude cutoffs of 0.4 mag in any of J, H, or K_s (depending on data availability) in their working definition of Miras. In this work, we similarly used a cutoff of 0.4 mag in the K_s band to distinguish between our Mira and SRV candidates. Finally, a visual inspection of the resulting light curves was carried out, which led to the rejection of a number of non-LPV contaminants.

Our thus derived catalog of LPVs is comprised of two subgroups, namely Mira candidates (130 sources) and LPV⁺ candidates (1013 sources). The latter group includes all variable sources satisfying the aforementioned criteria and whose variability is thus likely to be real, but for which the available data did not allow us to obtain fully convincing periods, some of which are likely to be multiperiodic. The LPV nature of some of these stars – i.e., whether they are Miras, SRVs, SRVs with long secondary periods (Soszyński et al. 2021, and references therein), or other types of objects – is thus not conclusively settled in this study. In fact, a number of SRVs may potentially also be present among the Mira candidates; conversely, a number of Miras may be present in the LPV⁺ sample, in addition to the OGLE-III Miras that are present there because of their unreliable VVV-PSF photometry. In this paper, even though we report all of these sources, analysis of the derived period distribution, ages, and distances is restricted to the Mira candidates sample (130 stars) only.

The light curves of all our Mira candidates, plus a representative sample of LPV⁺ candidates, are shown in Appendix A. Their photometric properties are summarized in Appendix B.

4. Cross-identification and tests

We cross-matched our initial sample of VVV-PSF LPV candidates, which as we have seen is comprised of 12 476 sources (see Table 2), with the published OGLE-III Mira catalog (Soszyński et al. 2013b), resulting in 202 sources in common. We next performed a visual inspection of their respective I- and K_s-band light curves. As a result, we find that the vast majority of these sources are saturated in the VVV photometry and present poorly defined light curves. The only OGLE source in our sample of Mira candidates is OGLE-BLG-LPV-075099, whose light curve is displayed in Fig. 4.

The fact that the majority of OGLE sources are saturated in VVV is not an unexpected result, for two main reasons. First, the VVV survey was specifically designed to go deeper in regions of high interstellar extinction than afforded by previous surveys, both in the visible (such as OGLE) and near-IR (such as 2MASS). Thus, sources in the VVV footprint that are bright in the visible can easily saturate the VIRCAM detectors. On the other hand, faint sources detected in the near-IR using VVV data will typically not be present in visual catalogs, as confirmed by our analysis.

Whenever available, OGLE periods were used, as they are based on much longer and robust datasets, where the periodicity is much more clearly defined than in our study (see Fig. 4). In any case, for the one source in common, our period (255.8 d) and the one from OGLE (254.1 d) differ by only 1.7 d. Figure 4 shows one cross-matched source of our initial sample with the published OGLE-III Mira catalog. While this source is saturated in the VVV K_s-band data, the cyclical nature of its variability can still be observed.

In addition, our sample was cross-matched with the Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD; Wenger et al. 2000) and the American Association of Variable Star Observers (AAVSO) Variable Star Index (VSX; Watson et al. 2006). In total, 43 matches were found, with those sources being classified as Mira candidates (2 sources), asymptotic giant branch stars (AGB, 13 sources), and infrared sources (28 sources). We also cross-matched with the Gaia DR2 LPV catalog (Mowlavi et al. 2018). Positive matches mainly correspond to saturated sources in our VVV-PSF photometry. However, there are two Gaia DR2 matches that are not saturated in our photometry, and whose periods also fall in the right range for the computation of photometric distances; these will be further discussed in Sect. 6.3 below. Lastly, we cross-matched our sources with known clusters located in our studied fields, which includes four globular clusters (Terzan 9, Djorg 2, Terzan 10, and NGC 6544) and a number of open clusters (Dias et al. 2002). No cluster members could be confirmed in our analysis.

The significant number of outliers in the partially saturated VVV-PSF light curves (see Fig. 4, bottom panels) for the OGLE cross-matched sources and the different passbands employed hinder the analysis of the criteria we used to select Miras. Therefore, the catalog of Miras with K_s-band light curves obtained by Matsunaga et al. (2009, hereafter MiraMat) was also added to check our approach. This catalog contains 936 Miras with reliable photometry. Next, variability indices, periods and feature parameters were computed using the same algorithms used to analyze the VVV-PSF data (Sect. 3). We find that, among the MiraMat sources flagged in their study as having good light curves and for which the authors could provide periods, about 98% are selected, and their estimated periods are in agreement with those of Matsunaga et al. (2009).

Histograms of selected properties of MiraMat and VVV-PSF Miras are presented in Fig. 5. Additional comparisons are provided in Figs. 2 and 3. Based on these, the following main conclusions can be reached:

• The VVV-PSF Mira magnitudes are distributed in the range of 8–16 mag in K_s. The MiraMat sample does not reach as deep, covering 9–14 mag in K_s.

• The number of correlated measurements (N_C) found for the MiraMat sample is at least twice that found in the VVV-PSF Miras, which shows the different observation strategies of these two surveys, keeping in mind that the total numbers of observations (N) are not markedly different.

• According to Fig. 2, about 50% of the MiraMat sample stars assume values along the upper limit of the $K_{(\mathrm{fi})}^{(2)}/Fap$ distribution. As expected, higher N_C values lead to a better discrimination of Miras from noise.

• We caution that the analysis of the ${\chi }_{\nu }^2$ values for the MiraMat sample is biased, since the magnitude errors are unavailable (Fig. 2). The majority of MiraMat sources are bright stars that are not saturated in their data, which can imply error bars smaller than adopted in our analysis (σ_i, MiraMat = 0.01), and hence larger ${\chi }_{\nu }^2$ values according to Eq. (1). On the other hand, the faint sources could have lower ${\chi }_{\nu }^2$ values if their corresponding observational error bars are larger than assumed in our calculations.

• Our VVV-PSF Mira candidates cover a wider range in ${\chi }_{\nu }^2$ than do the Miras in the MiraMat sample. In particular, the former sample has proportionally more stars with lower ${\chi }_{\nu }^2$ values than the latter.

Fig. 5. Comparison of VVV-PSF candidates with the sample of Miras from Matsunaga et al. (2009). Histograms are shown of mean Ks-band magnitudes (upper left), amplitudes ΔKs (upper right), number of observations N (bottom left), and number of correlated measurements NC (bottom right). The VVV-PSF Miras, MiraMat sample, and VVV-PSF LPV candidates are shown in blue, yellow, and gray, respectively.

In summary, the fact that ∼98% of MiraMat sources are successfully recovered and that our periods are in agreement with those from Matsunaga et al. (2009) lends support to our using our proposed methodology (Sect. 3). On the other hand, we caution that, while application of our method to the MiraMat sample demonstrates that it can be very successful with high-quality photometric data such as provided in their study, the method may not necessarily perform equally well when applied to different datasets, having different cadences, time coverage, total number of measurements, and photometric error distributions. To account for such differences, and depending on the desired purity of the sample to be selected prior to the visual inspection stage, the selection criteria laid out in Sect. 3 may need to be tweaked accordingly.

5. Light curve features

Mean K_s magnitudes and amplitudes (ΔK_s) were estimated based on the models given by Eq. (3). Indeed, in the K_s case, the relative difference between a straight average over the available magnitudes (circles in Fig. 6) and mean magnitudes computed on the basis of detailed light curve fits (colored bands in Fig. 6) is on average ∼1.8%. The largest differences are found when the available observations do not cover uniformly the whole light curve – for instance, when most data are placed close to maximum or minimum light. In the case of the star shown on the bottom right panel of Fig. 6, for instance, the relative lack of measurements close to minimum light biases the mean K_s value, which is significantly brighter than found from the detailed light curve fit.

Fig. 6. Examples of folded light curves in J (blue), H (green), and Ks (yellow). The solid and dashed lines represent the mean magnitudes considering a straight average over the available observations and best-fit template model (if available), respectively. The color bands around data points correspond to a 2σ range around the mean. The VVV-ID and periods are included on top of each panel. In the bottom right panel, no data are available in the J and H bands, and so only the Ks photometry shown. We note that distances could be computed for both stars in the upper row, but not for those in the bottom row, due to saturation (bottom left panel) or a period that places the star outside the distance indicator regime (bottom right panel).

VVV-PSF measurements in J and H are not available for all sources in our sample (Table B.1). In those cases, we retrieved 2MASS magnitudes and examined both 2MASS and VVV images, in the hope of replacing the missing VVV values with 2MASS ones. Unfortunately, however, we concluded that for no such star was a reliable 2MASS magnitude available, especially due to the faintness of the stars and/or the presence of nearby companions, unresolved in the 2MASS images.

The VVV-PSF sample includes typically around three observations in the J and H bands, as compared to over 70 in the K_s band. Therefore, we can precisely determine mean magnitudes ⟨K_s⟩ and amplitudes ΔK_s using the available data, but a more careful approach is needed in the case of J and H, since the available observations are generally insufficient to completely define the light curves. In order to obtain reliable J and H mean magnitudes and amplitudes, we resort to describe the J and H light curves by means of suitably scaled templates obtained initially using the K_s-band data. Our approach is based on the empirical relations for amplitudes and phase lags given by Yuan et al. (2018), and the results are illustrated in Fig. 6.

More in detail, we initially obtain a K_s-band template for each star, fitting the coefficients of the model given by Eq. (3) to our K_s data. From the thus derived K_s templates, we then obtain the corresponding templates in J and H by using the empirical relations for amplitude ratios and phase lags provided by Yuan et al. (2018). The latter authors assume a linear dependence of these two quantities on the log-period; the corresponding coefficients are provided in their Table 2. Errors are unfortunately not available for these coefficients. However, their Fig. 2 shows that the scatter around the values inferred from such linear expressions is of order ≈ ± 0.5 in the amplitude ratios and ≈0.05 in the phase lags (1 − σ level). By allowing the amplitude ratios and phase lags to vary within these ranges, we finally adopted, for each individual star and bandpass, the combination of amplitude ratio and phase lag that minimized the residuals of the fits to our data. The thus-fitted templates allow us to obtain reliable mean J and H magnitudes, along with the corresponding amplitudes, for all our stars with at least one J or H VVV-PSF measurement. Examples of our resulting template fits are shown in Fig. 6, which also indicates the resulting mean values.

We note that, in our procedure, the coefficients corresponding to both C-rich and O-rich Mira provided by Yuan et al. (2018) were initially used. However, since the Mira population found in the Galactic bulge is mainly comprised of O-rich Miras (Soszyński et al. 2013a), we finally decided to adopt the solutions corresponding to the latter. However, as we shall see later, a significant fraction of the Miras for which we could compute distances are not located in the bulge proper, and thus the assumption of an O-rich composition is not necessarily valid for those in particular. In this sense, Guo et al. (2022) used VVV-PSF photometry from VIRAC2 (Smith et al. 2018, and in prep.) to perform a search for young stellar objects (YSOs) in the VVV footprint, independently finding a large number of Mira candidates. In those tiles in common with our analysis, they were able to recover most of our candidate Miras; in addition, based on WISE (Wide-field Infrared Survey Explorer; Wright et al. 2010) mid-IR photometry, and using diagnostic color-color plots as in Lian et al. (2014), they were able to confirm that the majority of these Miras are indeed O-rich, even though there are a few that may be C-rich. Follow-up spectroscopy would be of significant interest to verify their color-based classifications. C-rich stars in the bulge are of particular interest as they may either belong to a young population or be the result of mergers (Feast et al. 2013; Matsunaga et al. 2017). A full discussion of this topic is beyond the scope of this paper, and will be presented in Guo et al. (2022).

6. Results

Appendix 6 reports the main measured properties for the stars in our final catalog, containing, as we have seen, 130 Mira candidates and 1013 LPV⁺ candidates. Columns (1) through (9) in Table B.1 list, respectively, our Mira sources’ IDs, coordinates, mean magnitudes, and amplitudes in each of J, H, and K_s, in addition to the period (in days), which is given in Col. (10). In Cols. (1) through (8) of Table B.2, in turn, one can find a list of our LPV⁺ sources’ IDs, OGLE IDs, coordinates, mean magnitudes in each of J, H, and K_s, and amplitudes in K_s, respectively.

6.1. Period distribution

Our main goal in this work is to detect Mira stars. The period histograms of our VVV-PSF Miras, along with the Miras of the OGLE-III and MiraMat surveys, as well as Gaia LPVs, are shown in Fig. 7. Because of selection biases, the period histograms of these surveys are somewhat different. The derived period distributions depend on the observed Galactic regions (affecting the relative proportion of young versus old Miras), wavelength regime (bright and/or close versus faint and/or distant ones), and time span of the observation (short- versus long-period variables). For instance, the OGLE survey, as an optical experiment, is unable to detect Miras located in high-extinction areas, while being successful to cover Miras either situated quite close to us, or intrinsically young and bright. On the other hand, the Gaia survey, being an all-sky experiment, covers a comparatively broader sample of Miras, in terms of period range. We caution, however, that the Gaia LPV sample includes non-Miras as well, and so a more refined classification of the Gaia LPVs will be needed before stronger conclusions can be drawn.

Fig. 7. Period histograms of Gaia LPVs (black), OGLE-III Miras (green), VVV-PSF Miras (orange), and MiraMat Miras (blue) are presented.

In the case of the VVV survey, the brightest Miras are mainly saturated, therefore we do not identify their contribution in our sample unless they are located sufficiently far. Even though we need to cover more tiles of the VVV survey to allow more detailed comparisons with other available surveys that cover wider fields, the bimodality of the VVV-PSF Miras histogram is readily apparent. More specifically, the VVV histogram reveals two notable peaks, around ∼200 d and ∼500 d. Other studies, in addition to those that are shown in the figure, have variously found a single peak similar to our long-period one (Whitelock et al. 1991; Matsunaga et al. 2005) or a bimodality that resembles ours to some extent (Glass et al. 2001; Matsunaga et al. 2009), though with a short-period peak located at a longer period value than we find (see also Matsunaga et al. 2009). On the other hand, also using VVV data but with much smaller samples of Miras and SRVs than in the present study, Gran (2013) and Molina et al. (2019) also found signs of bimodality in their period distributions, with minima located in the range 300 ≲ P(d)≲400.

The bimodality that is seen in our sample thus indicates the presence of both young (long-period) and old (short-period) LPVs. Molina et al. (2019) suggested that the demise of VVV LPVs in the 200–400 d period range could be due to the fact that Miras located in the bulge proper would be saturated for periods of 200 d or longer, as they would have ⟨K_s⟩≈11 or less. Indeed, we see a minimum in the magnitude distribution of our Miras and Mira candidates at ⟨K_s⟩≈11 (see Fig. 5), and have also verified that the bulk of the OGLE Miras that have cross-matches with our sources have ⟨K_s⟩≈11 and very poor VVV-PSF light curves. In addition, a selection bias may plausibly also be present, caused by the seasonal nature of VVV observations: we note that the minimum in our histogram is found close to 365 d (Fig. 7), where most other non-VVV surveys have found significant numbers of LPVs. Whatever the exact reason(s) for this dip may be, more observations will be required in order to properly detect and establish the multiband light curve properties of intermediate-period LPVs located toward the Galactic bulge.

6.2. Period-amplitude relation

Figure 8 indicates the JHK_s amplitudes of VVV-PSF Miras as a function of their periods. The same parameters found in the Matsunaga et al. (2009) catalog are plotted as reference. The ΔJ, ΔH, and ΔK_s amplitudes are in agreement with Matsunaga’s results (see also Gran 2013). As can be seen, the period-amplitude distribution of Miras gets tighter in the K_s-band, since the amplitudes are reduced in this bandpass.

Fig. 8. ΔJ, ΔH, and ΔKs amplitudes as function of logarithmic period, from top to bottom, respectively. The VVV-PSF Miras (orange circles) and MiraMat Miras (blue circles) are displayed. The dashed line marks the short-period limit of the Conroy et al. (2015) relation (Eq. (7)). The lines show our fits for the combined sample of VVV-PSF and MiraMat Miras.

Conroy et al. (2015) introduced a period-amplitude relation in the I-band for Miras (log P > 2.2):

$$\begin{array}{c}\hfill \mathrm{\Delta }I={10}^{({\alpha }_{I}logP+{\beta }_I)},\end{array}$$(7)

where P is the period in days, ΔI represents the amplitude in the I band, and α_I and β_I are constants (see Table 3). In the same fashion, we also establish α and β for each of the JHK_s wavebands. Our results are shown in Table 3, using the combined sample of VVV-PSF and MiraMat Miras.

Our amplitude-period relation in the K_s band can be compared with the similar ones that were introduced recently by Molina et al. (2019), which however are valid only for log(P/d) > 2.6. As a result, we find that our expression provides ΔK_s values that are in better agreement with their Eq. (2) in the period range 2.6 < log(P/d) < 3.18, which is also the range of applicability of our expression (see Fig. 8). This is not unexpected, as Eq. (2) in Molina et al. (2019) was also obtained by combining VVV and MiraMat amplitudes, and predicts amplitudes that are smaller than the ones based on VVV-only amplitudes in the long-period regime (which is also apparent from our Fig. 8). For further discussions on the implications of these results, we refer the reader to Molina et al. (2019) and Ita et al. (2021), and the references cited therein.

6.3. Distance indicators and extinction

Thirty-one VVV-PSF Miras fall in the period range 100 ≲ P(d)≲350, over which the near-IR period-luminosity relation of Miras is known to be linear, and they are consequently expected to be reliable distance indicators (e.g., Feast et al. 2002; Matsunaga et al. 2009). Among these, 13 have ⟨K_s⟩< 12 mag, and are thus in the regime of saturation according to Contreras Ramos et al. (2017, see also Sect. 3.1). We exclude these stars from further analysis, and we are left with a final sample of 18 sources.

In what follows, we refer to stars satisfying these criteria as distance-indicator Miras. Follow-up photometry of the remaining 13 Miras in this crucial period range is strongly encouraged, as reliable distances could be obtained for them, thus providing useful insight into the structure of the inner Milky Way, and indeed beyond.

The absolute magnitudes of the distance-indicator Miras in our catalog were obtained by means of the period-luminosity relation derived by Matsunaga et al. (2009) in the Infrared Survey Facility (IRSF)/SIRIUS photometric system, suitably transformed into the VISTA system. This was accomplished by first converting the IRSF relations to their 2MASS counterparts following Carpenter et al. (2001), and the latter to the VISTA system according to the prescriptions given in González-Fernández et al. (2018). The thus transformed equations are as follows:

$$\begin{array}{c}\hfill M_J=-5.82-2.88(logP-2.3),\end{array}$$(8)

$$\begin{array}{c}\hfill M_H=-6.53-3.18(logP-2.3),\end{array}$$(9)

$$\begin{array}{c}\hfill M_{K_s}=-6.89-3.55(logP-2.3),\end{array}$$(10)

where M_J, H, Ks is the absolute magnitude and the period P is in days. Naturally, stars in the inner regions of the Milky Way can be subject to large amounts of foreground extinction, even in the near-IR, and this needs to be properly accounted for in order for accurate distances to be derived.

While dust maps obtained with tracers such as red clump and red giant stars are available (Gonzalez et al. 2012; Surot et al. 2020), a comparison between the observed and expected absolute magnitudes of Miras in different passbands also provides a route toward measuring foreground extinction on a star-by-star basis. Thus, for instance, following Matsunaga et al. (2009), extinction in the K_s band can be obtained from the observed J − K_s or H − K_s colors using the following expressions:

$$\begin{array}{c}\hfill A_{K_s}^{J{K}_s}=\frac{(⟨J⟩-⟨K_s⟩)-(M_J-M_{K_s})}{r_J-1},\end{array}$$(11)

$$\begin{array}{c}\hfill A_{K_s}^{H{K}_s}=\frac{(⟨H⟩-⟨K_s⟩)-(M_H-M_{K_s})}{r_H-1},\end{array}$$(12)

where r_J ≡ A_J/A_Ks and r_H ≡ A_H/A_Ks.

Both r_J and r_H in these expressions depend on the adopted extinction law. The extinction law toward the Galactic bulge is complex and nonstandard, and several authors have provided their own calibrations. Distances derived using Eqs. (11) and (12) will naturally depend on the adopted extinction law, and the differences, along with calibration issues such as those detected by Hajdu et al. (2020) in the case of VVV (which were corrected for in this study as described in Sect. 3.1), may affect the astrophysical interpretation of the results (e.g., Matsunaga et al. 2016; Dékány et al. 2019). In this work, we have performed calculations using different prescriptions, including a (standard) Cardelli et al. (1989) law. For ease of reference, the r_J and r_H values corresponding to the different laws examined in this paper are summarized in Table 4.

In our work, when deriving distances and extinctions, intensity-mean magnitudes in J, H, and K_s from our VVV-PSF photometry are used when available. Absolute magnitudes are in turn obtained from the periods using the Mira period-luminosity relations (Eqs. (8)–(10)). We call the thus derived, star-by-star extinction values intrinsic extinctions in what follows.

In addition, we used both the Bulge Extinction and Metallicity Calculator³ (BEAM; Gonzalez et al. 2012) and the latest Surot et al. (2020) extinction maps⁴ to obtain E(J − K_s) reddening values, based on which extinction values $A_{K_s}^{\mathrm{BEAM}}$ and $A_{K_s}^{\mathrm{Surot}}$ were obtained, assuming an extinction law as proposed by Surot et al. (2020). We note that the recent Surot maps have higher spatial resolution than provided by BEAM. The latter is designed to provide extinction values at spatial resolutions not better than 2′. Therefore, $A_{K_s}^{J{K}_s}$ and $A_{K_s}^{H{K}_s}$ should in principle provide better estimations, if reliable VVV-PSF photometry is available, as $A_{K_s}^{J{K}_s}$ and $A_{K_s}^{H{K}_s}$ are both obtained on a star-by-star basis. The Surot maps, on the other hand, have sub-arcmin spatial resolution, and thus its extinction values should be more directly comparable to those that we derive.

The final extinctions for our distance indicator Miras are provided in Table 5, whereas Fig. 9 shows how the values computed using different extinction laws and reddening maps compare with one another other. There is a clear tendency for our color-based values to be higher than their respective BEAM and Surot ones. This can plausibly be ascribed to depth effects, as the BEAM and Surot maps are two-dimensional and based primarily on red clump bulge stars, thus being incapable of tracing any dust that may be present on the far side of the bulge and beyond, where most of our candidates appear to be located (see below).

Fig. 9. For nonsaturated VVV-PSF stars, the intrinsic extinction $A_{K_s}^{J{K}_s}$ (top panel) and $A_{K_s}^{H{K}_s}$ (bottom panel) is shown as a function of $A_{K_s}^{\mathrm{Surot}}$, according to the indicated extinction laws (see Table 4). The same $A_{K_s}^{\mathrm{BEAM}}$ values are shown in both panels, obtained from E(J − Ks) reddenings using a Surot et al. (2020) extinction law, which is also used to obtain $A_{K_s}^{\mathrm{Surot}}$ values. The dashed gray line indicates a perfect correlation.

True distance moduli are computed as (m − M)₀ = (⟨K_s⟩−A_Ks)−M_Ks, using the different prescriptions for A_Ks that were just described, whereupon distances D (in pc) were finally obtained using (m − M)₀ = 5(log D − 1). The latter distances are shown in Table 6. In that table, we refer to distances derived using $A_{K_s}^{J{K}_s}$ as D^JK_s, whereas those obtained using $A_{K_s}^{H{K}_s}$ are termed D^HK_s instead. Also, D^BEAM and D^Surot are the distances obtained adopting extinctions based on the BEAM and Surot maps, respectively, in the manner described in the previous paragraph.

The resulting D^JK_s and D^HK_s distance distributions of the VVV-PSF distance indicator Miras are shown in Fig. 10, plotted as a function of the Miras’ ages (see Sect. 6.4). As can be seen from this figure and Table 6, most of our distance-indicator Miras appear to be very distant, though their exact distances remain uncertain, with significant differences found for some stars based on the method used and extinction law assumed. As a rule, however, our D^JK_s and D^HK_s distances place them much closer than implied by D^BEAM or D^Surot. This too can be understood in terms of depth effects, to which our color-based method is sensitive but which are not properly described by the BEAM and Surot maps. The end result is that the BEAM and Surot extinctions are likely underestimated, and the corresponding distances likely overestimated, compared with those computed using the stars’ colors. We note, in particular, and in accordance with Table 5 and Fig. 9, that a significant fraction of our stars may be subject to quite extreme levels of extinction, at a level that is not captured by either the BEAM or Surot maps. Such high extinction levels are consistent with the fact that some of these sources could not be reliably detected in the J band. One possible problem with this explanation is the fact that the four stars that appear to be closest to us according to their D^HK_s values (i.e., those with VVV IDs 319302072957, 319509072340, 320201006873, and 320201014922), are also the ones with the highest $A_{K_s}^{H{K}_s}$ extinction values (see Tables 5 and 6), which would appear to be contradictory, particularly since their corresponding $A_{K_s}^{\mathrm{BEAM}}$ and $A_{K_s}^{\mathrm{Surot}}$ values, which measure the integrated dust extinction up to the distance of the Milky Way bulge, are significantly smaller. Note, in addition, that there is one extinction law, namely AG17 (see Table 4), that places all these stars sufficiently far away that they could still be on the far side of the bulge. Other possible explanations may include poor measurements (these stars are all quite faint in H, and indeed not detected in J) or their being C-rich stars. Further analysis will be required to properly understand their properties.

Fig. 10. Distribution of distances (in kpc) versus ages (in Gyr) for the distance indicator Mira stars in our sample. DJKs (top panel) and DHKs (bottom panel) values are computed using the different indicated extinction laws (see Table 4). DBEAM and DSurot values are shown in both panels, and were computed from the E(J − Ks) values provided by the BEAM and Surot maps assuming a Surot et al. (2020) extinction law. These maps are based primarily on bulge red clump stars, hence are not expected to accurately measure the reddening toward objects located beyond the Galactic bulge; accordingly, the corresponding distances should be considered as upper limits only. Due to lack of J-band measurements, several of the stars that are shown in the bottom panel with their respective DHKs values do not have DJKs data shown in the top panel (see also Table 6. Ages were computed following Eq. (13). We note that the age range is restricted in this plot, as only Miras that have periods in the range [100, 350] d (i.e., distance indicators) are included, corresponding to an age range of [5.4, 13.6] Gyr according to Eq. (13).

As a further check of the colors of the stars in our sample, in Fig. 11 we show the color-color distribution of the sources for which we could obtain both J − H and H − Ks measurements, compared with the loci occupied by the Mira and SRV sources previously studied by Kerschbaum et al. (2001) and Matsunaga et al. (2009) (light blue and red polygons, respectively). As we can see from this plot, the majority of the stars in our samples for which colors could be measured fall reasonably close to the loci defined by the distributions of Miras and SRVs studied by these authors. However, as in Matsunaga et al. (2009), there is a significant fraction of stars in our sample that fall to the right of the main distribution, i.e., toward redder H − K_s colors. Some could be symbiotic Miras (Whitelock 1987, 2003). On the other hand, uncertainties in the zero-point corrections that were introduced in Sect. 3.1 to alleviate the calibration issues that were pointed out by Hajdu et al. (2020) could affect the derived colors, particularly in the higher-density regions, as could the fact that for many stars the number of available J and H measurements is extremely small. However, there are several other possible explanations for these anomalous colors. As pointed out by Matsunaga et al. (2009), unresolved blends, some of which may be difficult to detect even upon close examination of the images, may affect the measured colors, by moving them either to the blue (if the contaminant is a foreground blue source) or to the red (if the contaminant is redder than the candidate Mira or SRV). In addition, we have verified that some of these stars fall within the regime of saturation in the VVV images, so that, as previously pointed out (Sect. 3.2), their magnitudes and colors must be treated with extreme caution. Also, for a few of them the available J and/or H measurements were all obtained close to the maximum of the light curve, thus rendering their corresponding mean magnitudes particularly uncertain: in fact, this is exactly what happens in the case of the two outliers in the bottom right of the diagram, whose VVV IDs are 320211006775 and 364207072166. Lastly, some stars were detected close to the limit of the VVV-PSF photometry, and thus their measurements can be quite unreliable. In particular, we have checked that these are indeed suitable explanations for the three sources that deviate the most to the red in the color-color diagram, at around H − K_s ≈ 4, and the one that has a J − H color close to zero. Further analysis of all of these sources with peculiar colors would be very valuable to place their near-IR magnitudes and colors, and indeed their very classification as Mira candidates, on a firmer footing.

Fig. 11. Color-color (J − H vs. H − Ks) plot for the nonsaturated candidate Miras in our sample for which measurements in the three bandpasses could be obtained (crosses). The two stars in our distance indicator sample are marked with blue circles. The shaded blue and red regions are schematic depictions of the loci occupied by Miras and SRVs in the Kerschbaum et al. (2001) and Matsunaga et al. (2009) samples, respectively. The arrow at the top of the figure indicates the reddening vector for an absolute extinction of 1 mag in Ks, assuming a D19 extinction law (see Table 4).

Be as it may, our sample clearly differs from the one previously studied by Matsunaga et al. (2009), which presents a sharp peak around 8.2 kpc, which is essentially the distance to the Galactic center (Bland-Hawthorn & Gerhard 2016; GRAVITY Collaboration 2019, and references therein). The VVV-PSF Miras for which distances could be reliably computed extend much farther out, with many stars seemingly located well beyond the Galactic bulge. This can be understood if we recall not only that the VISTA telescope that was used in the VVV survey has a 4-m aperture, as compared to IRSF’s 1.4-m aperture, but also that the exposure times adopted in the VVV survey were sufficiently long that unprecedented depth could be reached. Conversely, the more nearby Miras are missing in our catalog because they are saturated in the VVV-PSF images, hence reliable distances could not be computed for them. Lastly, the volume surveyed on the far side of the bulge is vastly greater than that on its near side.

Gaia parallaxes can in principle also be used to obtain distances to LPVs. However, such stars may often not be point sources for the Gaia survey, so that their corresponding Gaia distances should be taken with due caution, considering that these are large stars that can have complex, structured photospheres, and angular diameters comparable to (or even larger than) their measured parallaxes (Whitelock 2012). Keeping this caveat in mind, we have cross-matched our catalog against Gaia’s DR2 and early DR3 (eDR3; Gaia Collaboration 2021). As pointed out in Sect. 4, most of the matches that we found correspond to saturated sources in our VVV-PSF photometry. There are, however, two exceptions, namely stars with VVV IDs 307408087155 and 309203069542. Both of these sources have periods that place them within the realm of distance indicators, and thus we were able to compute distances for all of them. Indeed, according to Table 6, these sources are located at distances of roughly 144 kpc and 150 kpc, respectively. These are rather extreme distances which, if correct, would make derivation of their parallaxes extremely challenging, even for Gaia (see also Mauron et al. 2021, who find negative or insignificant parallaxes for their outer halo LPVs, whose distances they report to be in the range between 50 and 119 kpc). On the other hand, since the only available extinction values for these stars are the ones obtained on the basis of the BEAM and Surot maps, their distances could be significantly overestimated (see Fig. 10 and Table 6). Accordingly, we carefully checked Gaia’s eDR3 catalog, and found the following results:

– VVV 307408087155: There are actually two Gaia’s eDR3 sources within 1″ of this star. The closest one, Gaia eDR3 4062695911426083840, has a parallax value of 1.37 ± 0.31 mas, thus placing it at a distance of only about 739 pc from us. It is also relatively bright, with a reported G = 16.7. This source, which also looks bright in VVV images, is offset from the position of VVV 307408087155 by about 0.6″. We deem it very unlikely that our source, which is so faint as to even go undetected in J and H, could correspond to this bright, nearby star. There is another source in the vicinity that is present in Gaia’s eDR3, namely Gaia eDR3 4062695915664384256, that is 0.8″ away from the position of the VVV source. This star has a reported G = 19.3, and no parallax measurements reported. We believe that it is sufficiently far away in the sky that it is also unlikely to correspond to our detected VVV source.

– VVV 309203069542: There is only one source matching the VVV coordinates within a radius or 1″, namely Gaia eDR3 4064127338032497024. At G = 20.1, this is a very faint source in Gaia’s eDR3, which assigns to it a negative parallax, thus making it impossible to compute a reliable astrometric distance.

The majority of our distance-indicator Miras appear to be located well beyond the Galactic bulge: some of them may be associated to distant MW arms, background halo stars, heretofore unknown globular clusters or dwarf spheroidal galaxies, or even the Sagittarius dwarf spheroidal galaxy (Ibata et al. 1994), whose main body is located at a distance of about 26 kpc (McConnachie 2012), and whose tidal arms extend across much of the sky (e.g., Navarrete et al. 2017; Ramos et al. 2020, and references therein).

We note that the distances of some of our stars place them far out in the MW halo, some being even close to (one even more distant than) the predicted halo edge, which Deason et al. (2020) place at a distance of 292 ± 61 kpc. One should keep in mind, however, that for many of the stars with the largest distances, the latter were derived on the basis of BEAM and Surot extinctions, and are thus likely overestimated. On the other hand, using light curves from the Catalina Sky Surveys (Drake et al. 2017), Mauron et al. (2021) were similarly able to identify ten LPVs in the outer halo, with distances in the range ∼50–120 kpc, and thus not reaching as far out into the outermost halo as we do in this study. On the other hand, Mauron et al. (2021) were able to use data from Gaia’s DR3 (Gaia Collaboration 2021) to show that the tangential velocities of their 10 distant LPVs are consistent with expectations based on other bona-fide halo objects, such as globular clusters and dwarf galaxies.

Given the fact that the period range of our distance indicator Miras is capable of tracing intermediate- and old-age populations only (see Sect. 6.4), even the closest among our stars with measured distances are unlikely members of spiral arms, but may more likely belong to the intermediate-age components of the MW thick disk (e.g., Robin et al. 2003; Giammaria et al. 2021), dwarf satellite galaxies, or even halo (see also Grady et al. 2019, 2020). Some of the closer stars, if their D^HK_s distances are to be trusted, could in principle belong to the MW bulge, the latter being a complex system with a radius of ∼2.5 kpc (Zoccali & Valenti 2016; Barbuy et al. 2018; Zoccali 2019, and references therein), though this may perhaps seem unlikely from a statistical perspective, given the bulge’s fairly restricted volume of space compared to that covered by our sample, and also the fact that the BEAM and Surot maps place them much farther away. On the other hand, BEAM- and Surot-based extinction values for individual stars, and particularly the more distant ones, are likely underestimated, since, as previously discussed, the BEAM and Surot maps do not count the extinction component that is due to dust lanes located beyond the Galactic bulge. If so, BEAM- and Surot-based distances for our more distant stars may have been correspondingly overestimated, and thus their true distances may potentially be much smaller, given the very significant differences between the Surot/BEAM extinction values and the “intrinsic” ones computed for a few of our stars (see Tables 5 and 6). On the other hand, color-based extinctions can also be inaccurate, due to insufficient and/or incorrect data and/or measurements in the J and/or H bands. Clearly, further analysis will be needed in order to place the distances of our Miras at a firmer footing.

To close, we note that VVV data have also afforded the detection of faraway structures beyond the MW bulge based on other tracers, such as Cepheids, red clump stars, and microlensing events (e.g., Dékány et al. 2015a,b, 2019; Minniti et al. 2018, 2020; Navarro et al. 2020; Saito et al. 2020).

6.4. Period-age relation

As pointed out in Sect. 1, different prescriptions that are available in the literature for the relationship between periods and ages of Miras provide different results. This can be clearly seen from Fig. 12, where the analytical expressions provided by López-Corredoira (2017) and Grady et al. (2019) are compared, and data inferred from a variety of studies and compilations (Feast & Whitelock 1987, 2014; Feast 2007; López-Corredoira 2017) are also overplotted. As can be clearly seen, the Grady et al. (2019) relation does not match the available data well. The López-Corredoira (2017) relation, on the other hand, provides a reasonable description of the data at the intermediate ages where his relation was calibrated, but its quadratic nature leads to unphysical results when extrapolated toward both the short and long period ends. To compensate for these shortcomings, here we propose the following expression:

$$\begin{array}{c}\hfill \mathrm{Age}=\frac{a}{1+{\left(\frac{P}{c}\right)}^b},\end{array}$$(13)

Fig. 12. Period-age relations for Mira variables, shown in a linear scale (upper panel) and using a logarithmic scale for the ages (bottom panel). Blue crosses, blue squares, yellow circles, and magenta diamonds represent datapoints extracted from Feast & Whitelock (1987, 2014), Feast (2007), and López-Corredoira (2017), respectively. The blue arrows indicate that only an upper limit on the age is available. The green dashed line represents the relation suggested by Grady et al. (2019) whereas the light blue solid line the relation by López-Corredoira (2017). The light blue dotted line segments illustrate the expected behavior of the López-Corredoira (2017) relation, when it is extrapolated beyond its formal range of validity. The red line is the relation proposed in this work (Eq. (13)), and the light red band represents a dispersion of Δlog t = ±0.15 around it, which basically encompasses all available measurements in the literature, but may be an underestimate at the long-period/short-age end in particular, where empirical data are not yet available.

with a = 13.8 Gyr, b = 3.6, and c = 310 days. In this expression, the age is in Gyr, and the period P is in days. As Fig. 12 shows, this relation has the clear advantage of describing correctly the expected dependence of periods on ages at both the young and old age extremes, while at the same time providing a reasonable description of the intermediate-age regime as well. Note, in particular, that Eq. (13) tends asymptotically to a value of 13.8 Gyr at the short-period end, which is the current best estimate for the age of the Universe (Bennett et al. 2013; Planck Collaboration XIII 2016; see also Catelan 2018 and Valcin et al. 2020, for recent reviews and additional references).

Several previous authors have used the Mira period-age relation to study how the ages of Miras change as a function of position across the different Milky Way components, including halo, disk, bulge, and bar (e.g., Catchpole et al. 2016; López-Corredoira 2017; Grady et al. 2019, 2020). We adopt our improved relation with a similar purpose in what follows. However, we do caution that more age measurements, particularly at the young and old regimes where data are currently lacking, and at different metallicities, would be of great importance to place it on a firmer footing, as well as to constrain its possible dependence on chemical composition. Indeed, as also shown in Fig. 12, the still very limited data that are available for calibration purposes imply an uncertainty of Δlog t ∼ 0.15 in the ages inferred from this expression, which may actually be an underestimate, particularly in the long-period/short-age regime. Similarly, theoretical models indicate that individual LPVs may undergo significant period changes in the course of their evolution (e.g., Marigo et al. 2017; Trabucchi et al. 2021), which would imply the existence of an additional, empirically unconstrained, intrinsic spread around the values implied by Eq. (13). Figure 4 in Marigo et al. (2017), where evolutionary models and isochrones that include the thermally pulsing AGB phase were presented, shows that the period evolution, for a given pulsation mode, may be comprised of a relatively smooth (and small) long-term trend, with larger superimposed variations taking place on much shorter timescales. While the long-term component would likely dominate in a random sample of Miras of a given age, a larger sample and detailed Monte Carlo modeling would be required to properly evaluate the most likely period values and corresponding dispersions. Such modeling would be of significant interest, but is outside the scope of our paper.

It should also be noted that the limited data used to obtain this relation corresponds to disk stars in the solar neighborhood. Thus, it is unclear whether Eq. (13) can be applied in the case of stars belonging to different populations and with different metallicities. Observations of Miras in globular clusters suggest that the metallicity has a large impact on the period, at least among low-mass Miras (Frogel & Whitelock 1998; Feast & Whitelock 2000a). There are known Miras in moderately metal-rich globular clusters (Sloan et al. 2010) whose periods, according to Eq. (13), would imply ages of only a few Gyr, which is likely several Gyr younger than the actual ages of these clusters (e.g., Leaman et al. 2013; VandenBerg et al. 2013; Cohen et al. 2021). It has been suggested that the presence of long-period, C-rich Miras in globular clusters could be ascribed to the mergers of two old stars (Feast et al. 2013; Matsunaga et al. 2017). Naturally, for any such stars that are merger products, whether in clusters or in the field, the age implied by Eq. (13) would be an underestimate, compared with the ages of the original binary pairs that led to their formation. This has indeed been suggested as a possible explanation for the observed properties of C-rich stars in the Milky Way bulge, whose periods, when taken at face value, would imply that they have ages of only a few Gyr (Matsunaga et al. 2017).

A histogram showing the age distribution of our Mira candidates is presented in Fig. 13. Bimodality is apparent in this plot, and is a direct consequence of the bimodal distribution in periods for the stars in our sample (Fig. 7), the minimum close to 365 d corresponding to an age of 4.9 Gyr according to Eq. (13).

Fig. 13. Histogram of ages for the Mira stars in our sample. The ages are given in Gyr, and were computed based on the stars’ periods according to Eq. (13). The bimodal age distribution is a direct consequence of the bimodal distribution in periods (Fig. 7): the minimum close to 365 d corresponds to an age of 4.9 Gyr, according to Eq. (13).

The distribution of distances (listed in Table 6) as a function of age is displayed in Fig. 10 for our distance indicator Miras. D^JK_s (top panel) and D^HK_s (bottom panel) values are computed using the different indicated extinction laws (Table 4), based on VVV-PSF magnitudes. In addition, D^BEAM and D^Surot values are shown in both panels. Since only distance indicator Miras are included in this plot, the age range is significantly limited, compared with Fig. 13. Indeed, distance indicator Miras have periods in the range [100, 350] d, which corresponds to an age range of [5.4, 13.6] Gyr, according to Eq. (13).

In this very limited sample, and over this limited age range, there is no clear trend between distance and age, although a tendency of ages increasing with distance may be implied by the D^HK_s values only. Indeed, there is significant dispersion in ages at any given distance. This could in principle imply that Miras belonging to different stellar populations are present in our dataset. On the other hand, and as the bottom plot of Fig. 10 illustrates particularly well, any firm conclusions can only be derived once the uncertainty in the extinction values toward individual stars can be greatly reduced. Note, in this sense, that the seemingly reduced dispersion that is seen in the upper panel of Fig. 10 is, at least in part, an artifact of the fact that J-band magnitudes are not available for several stars, whose D^JK_s values are accordingly missing from that plot.

7. Conclusions

Near-IR time-series data from the VVV survey were used to perform a systematic search for LPVs toward selected fields in the MW bulge. These fields were chosen in order to maximize overlap with the OGLE-III LPV catalog. Based on PSF photometry performed on the VVV images, we developed a new technique that uses both correlated and noncorrelated statistical indices in order to select Mira candidates, which was tested against the previously known OGLE stars in those fields and was found to produce excellent results. As a consequence, a sample of 130 Mira candidates was obtained. An additional 1013 LPV candidates are also reported (our LPV⁺ candidates) for which the available data proved insufficient to obtain reliable periods, and thus provide a better assessment of their LPV status (Mira, SRV, or other). In both sets, and especially in the latter, other types of LPVs, in particular SRVs, may also be present. We provide a catalog with the full list of Miras and Mira candidates detected in this work. Our catalog includes coordinates, mean magnitudes, and amplitudes (obtained from a template-fitting procedure) in each of J, H, and K_s, periods, and ages. The latter were derived on the basis of a reassessment of previously proposed period-age relations in the literature, with a new analytical expression presented in this work.

Among the stars in our catalog, a significant number of distance indicator Miras are present – i.e., Miras with periods in the range [100, 350] d, over which the PL relation of Miras is known to be well behaved. Unfortunately, saturation affects the photometry of many of the stars in our sample, so that only a relatively small subset possesses reliable data that can be used for distance determination purposes. Conversely, extinction toward some of these stars can be so extreme, and their distances so large, that few or no reliable measurements in J and/or H exist. When the latter are available, however, extinction values were obtained on a star-by-star basis by comparing the intrinsic colors expected from the Mira PL relations in J, H, and K_s with the actually observed colors (when available), and assuming different extinction laws that have been proposed in the literature. Distances were also obtained for the same stars using both the BEAM and Surot reddening maps, assuming a Surot et al. (2020) extinction law. We find that the distances obtained using reddening maps are systematically larger than those computed using our color-based method, which may plausibly be due to the fact that those maps fail to describe the presence of dust on the far side of the Milky Way bulge and beyond.

We find that the stars in our sample have a bimodal distribution in periods, implying also a bimodal distribution in ages, with the minimum located close to 365 d, corresponding to an age of approximately 4.9 Gyr. This result is likely a consequence of different sources of bias affecting the VVV data. For the small subsample of distance indicator Miras which are not saturated in the VVV data, we do not find clear evidence of an age-distance relation, obtaining instead a significant spread in ages at any given distance. Unlike previous studies that lacked the sensitivity and depth of VVV observations, we find a significant number of very distant Miras, some of which may belong to the halo or dwarf spheroidal galaxies (including Sagittarius). A few of our stars may even be located close to the edge of the Milky Way halo, although this conclusion hinges heavily on reddening maps that do not properly count the extinction caused by dust lanes on the other side of the Milky Way bulge, likely leading to overestimated distances. Accordingly, any definite conclusions regarding the spatial and age distribution of our distance indicator Miras will require a significant reduction in the uncertainty associated with extinction along the line of sight, including, most crucially, three-dimensional (i.e., distance-dependent) effects. In future papers of this series, we will apply the same techniques developed in this study to other MW fields covered by VVV.

Acknowledgments

We gratefully acknowledge the constructive feedback provided by the referees, whose comments and suggestions have led to a significantly improved manuscript. We warmly thank Felipe Gran, Phil Lucas, Claudio A. Navarro Molina, and Zhen Guo for helpful discussions. Support for this project is provided by ANID’s Millennium Science Initiative through grant ICN12_12009, awarded to the Millennium Institute of Astrophysics (MAS); by ANID’s Basal projects AFB-170002 and FB210003; and by FONDECYT grant #1171273. F.N is grateful for financial support by Proyecto Gemini CONICYT grants #32130013 and #32140036, and by P. Universidad Católica de Chile’s Vicerrectoría de Investigación (VRI). C.E.F.L. acknowledges PCI/CNPQ/MCTIC post-doctoral support. M.Z. thanks Fondecyt Regular #1191505. A.R.A acknowledges support from FONDECYT through grant #3180203. The authors thank MCTIC/FINEP (CT-INFRA grant 0112052700) and the Embrace Space Weather Program for the computing facilities at INPE. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Some of the plots in this article were prepared using the python (van Rossum 1995) package matplotlib (Hunter 2007).

References

Appendix A: Light curves

The VVV-PSF photometry for our sample of Mira and LPV⁺ candidates is provided in Table A.1, available at the Centre de Données astronomiques de Strasbourg (CDS). A portion is shown here just for guidance in regard to its contents.

VVV-PSF K_s-band light curves are shown in Figs. A.1 and A.2 for our full sample of Mira candidates and for a representative sample of LPV⁺ candidates, respectively.

Fig. A.1. VVV-PSF light curves in the Ks band for the sample of Mira candidates. For each star, both the raw time series (left panels) and phased light curves (right panels) are shown. The VVV ID is displayed on top of the left panel, whereas our derived period is given at the top of the right panel.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.2. Representative VVV-PSF light curves in the Ks band for the sample of LPV+ candidates. The VVV ID is displayed on top of each panel.

Appendix B: Catalog

This appendix reports the main measured photometric parameters for our stars. Columns 1 through 9 in Table B.1 list our Mira sources’ IDs, coordinates, mean magnitudes, and amplitudes in each of J, H, and K_s, respectively, in addition to their periods, which can be found in column 10. Columns 1 through 8 in Table B.2 provide our LPV⁺ sources’ IDs, OGLE IDs, coordinates, mean magnitudes in each of J, H, and K_s, and amplitudes in K_s, respectively.

All Tables

All Figures

Fig. 1. Footprint of the VVV survey in the central bulge of the MW, expressed in Galactic coordinates. Each VVV tile is marked by its respective ID in this diagram, those analyzed in this work being highlighted in blue and using a larger font size. Red dots represent stars classified as Miras in the OGLE-III project. The gray grids present a region of the VVV footprint where our PSF photometry is already completed.

In the text

Fig. 2. Distribution of VVV-PSF sources in the $K_{(\mathrm{fi})}^{(2)}$ vs. ${\chi }_{\nu }^2$ plane. The Miras identified by Matsunaga et al. (2009, labeled MiraMat in the plot) are shown as blue circles. The region marked by dashed lines corresponds to the one from which we extracted our initial Mira candidates. The latter, upon further analysis and inspection, eventually led to our final selection of Miras, shown in this plot as orange squares. The maximum number of targets per pixel is 61 697 in this plot, corresponding to the dark orange-red region in the main panel.

In the text

	Fig. 3. Logarithm of the power as a function of the period for the selected VVV-PSF sources (see region delimited by black dashed line in Fig. 2). Symbols and colors are the same as in Fig. 2) as well. The maximum number of targets per pixel is 3146 in this plot, corresponding to the dark orange-red region in the main panel.
In the text

	Fig. 4. Example of raw (left) and phased (right) light curves, for a source cross-matched between OGLE (top panel) and VVV (bottom panel). The OGLE and VVV IDs, as well as the period found by OGLE, are displayed in the top panel.
In the text

	Fig. 5. Comparison of VVV-PSF candidates with the sample of Miras from Matsunaga et al. (2009). Histograms are shown of mean Ks-band magnitudes (upper left), amplitudes ΔKs (upper right), number of observations N (bottom left), and number of correlated measurements NC (bottom right). The VVV-PSF Miras, MiraMat sample, and VVV-PSF LPV candidates are shown in blue, yellow, and gray, respectively.
In the text

Fig. 6. Examples of folded light curves in J (blue), H (green), and Ks (yellow). The solid and dashed lines represent the mean magnitudes considering a straight average over the available observations and best-fit template model (if available), respectively. The color bands around data points correspond to a 2σ range around the mean. The VVV-ID and periods are included on top of each panel. In the bottom right panel, no data are available in the J and H bands, and so only the Ks photometry shown. We note that distances could be computed for both stars in the upper row, but not for those in the bottom row, due to saturation (bottom left panel) or a period that places the star outside the distance indicator regime (bottom right panel).

In the text

	Fig. 7. Period histograms of Gaia LPVs (black), OGLE-III Miras (green), VVV-PSF Miras (orange), and MiraMat Miras (blue) are presented.
In the text

	Fig. 8. ΔJ, ΔH, and ΔKs amplitudes as function of logarithmic period, from top to bottom, respectively. The VVV-PSF Miras (orange circles) and MiraMat Miras (blue circles) are displayed. The dashed line marks the short-period limit of the Conroy et al. (2015) relation (Eq. (7)). The lines show our fits for the combined sample of VVV-PSF and MiraMat Miras.
In the text

Fig. 9. For nonsaturated VVV-PSF stars, the intrinsic extinction $A_{K_s}^{J{K}_s}$ (top panel) and $A_{K_s}^{H{K}_s}$ (bottom panel) is shown as a function of $A_{K_s}^{\mathrm{Surot}}$, according to the indicated extinction laws (see Table 4). The same $A_{K_s}^{\mathrm{BEAM}}$ values are shown in both panels, obtained from E(J − Ks) reddenings using a Surot et al. (2020) extinction law, which is also used to obtain $A_{K_s}^{\mathrm{Surot}}$ values. The dashed gray line indicates a perfect correlation.

In the text

Fig. 10. Distribution of distances (in kpc) versus ages (in Gyr) for the distance indicator Mira stars in our sample. DJKs (top panel) and DHKs (bottom panel) values are computed using the different indicated extinction laws (see Table 4). DBEAM and DSurot values are shown in both panels, and were computed from the E(J − Ks) values provided by the BEAM and Surot maps assuming a Surot et al. (2020) extinction law. These maps are based primarily on bulge red clump stars, hence are not expected to accurately measure the reddening toward objects located beyond the Galactic bulge; accordingly, the corresponding distances should be considered as upper limits only. Due to lack of J-band measurements, several of the stars that are shown in the bottom panel with their respective DHKs values do not have DJKs data shown in the top panel (see also Table 6. Ages were computed following Eq. (13). We note that the age range is restricted in this plot, as only Miras that have periods in the range [100, 350] d (i.e., distance indicators) are included, corresponding to an age range of [5.4, 13.6] Gyr according to Eq. (13).

In the text

Fig. 11. Color-color (J − H vs. H − Ks) plot for the nonsaturated candidate Miras in our sample for which measurements in the three bandpasses could be obtained (crosses). The two stars in our distance indicator sample are marked with blue circles. The shaded blue and red regions are schematic depictions of the loci occupied by Miras and SRVs in the Kerschbaum et al. (2001) and Matsunaga et al. (2009) samples, respectively. The arrow at the top of the figure indicates the reddening vector for an absolute extinction of 1 mag in Ks, assuming a D19 extinction law (see Table 4).

In the text

Fig. 12. Period-age relations for Mira variables, shown in a linear scale (upper panel) and using a logarithmic scale for the ages (bottom panel). Blue crosses, blue squares, yellow circles, and magenta diamonds represent datapoints extracted from Feast & Whitelock (1987, 2014), Feast (2007), and López-Corredoira (2017), respectively. The blue arrows indicate that only an upper limit on the age is available. The green dashed line represents the relation suggested by Grady et al. (2019) whereas the light blue solid line the relation by López-Corredoira (2017). The light blue dotted line segments illustrate the expected behavior of the López-Corredoira (2017) relation, when it is extrapolated beyond its formal range of validity. The red line is the relation proposed in this work (Eq. (13)), and the light red band represents a dispersion of Δlog t = ±0.15 around it, which basically encompasses all available measurements in the literature, but may be an underestimate at the long-period/short-age end in particular, where empirical data are not yet available.

In the text

	Fig. 13. Histogram of ages for the Mira stars in our sample. The ages are given in Gyr, and were computed based on the stars’ periods according to Eq. (13). The bimodal age distribution is a direct consequence of the bimodal distribution in periods (Fig. 7): the minimum close to 365 d corresponds to an age of 4.9 Gyr, according to Eq. (13).
In the text

	Fig. A.1. VVV-PSF light curves in the Ks band for the sample of Mira candidates. For each star, both the raw time series (left panels) and phased light curves (right panels) are shown. The VVV ID is displayed on top of the left panel, whereas our derived period is given at the top of the right panel.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.2. Representative VVV-PSF light curves in the Ks band for the sample of LPV+ candidates. The VVV ID is displayed on top of each panel.
In the text