© ESO 2018

1. Introduction

In the Galactic bulge, the Red Clump (RC) stars, which are core-helium burning low-mass stars with ages from intermediate (1 ≤ age < 10 Gyr) to old (age ≥ 10 Gyr) and average to high metallicities ([Fe/H] ≳ −1.5dex, Cole 1998; Hill et al. 2011), are used extensively to study their kinematics, chemical abundances, and spatial distributions (McWilliam & Zoccali 2010; Saito et al. 2011; Gonzalez et al. 2015; Zoccali et al. 2017). High-resolution spectroscopic studies of RCs suggest that only stars that are more metal-rich than [Fe/H] ∼ −0.5 dex trace an X-shaped structure that appears to be in the form of a peanutboxy bulge (Ness et al. 2013; Zoccali & Valenti 2016). Recently, mid-infrared data obtained with the Wide-field Infrared Survey Explorer (WISE) also showed a clear large-scale X structure (Ness & Lang 2016). In contrast, stars that are more metal-poor than [Fe/H] ∼− 0.5 dex not only display a centrally concentrated axisymmetric spatial distribution, but also reveal different kinematics (Ness et al. 2013; Valenti et al. 2016; Zoccali et al. 2017).

However, the spatial distribution of the old and more metalpoor population of stars in the bulge, which is traced by RR Lyræ (RRLs), is still under discussion. Pietrukowicz et al. (2015) found an ellipsoidal distribution (axis ratios of 1:0.49 ± 0.02:0.39 ± 0.02), elongated along the same direction as the bar traced by the metal-rich red giants, while Dékány et al. (2013) and Kunder et al. (2016) found a spheroidal distribution.

Type II Cepheids (T2Cs) are old (> 10 Gyr), low-mass post-horizontal branch, asymptotic giant branch (AGB) and post-AGB stars. Like the RRLs, T2Cs trace old stellar populations, but they have longer periods (1–80 days), are brighter by 1–3 mag, and their amplitudes can be up to twice as large as those of the RRLs. In his pioneering work that led to the separation of Population I and Population II stars, Baade (1944) showed that T2Cs are distance indicators that obey a period-luminosity (PL) relation different from that of Classical Cepheids. T2Cs have been widely used in the literature as distance indicators, both in the optical (Harris 1985; Nemec et al. 1994) and in the near-infrared (NIR) bands (Matsunaga et al. 2006, 2011; Bhardwaj et al. 2017a, b), although not as frequently as RRLs and Classical Cepheids. A key feature of the PL relations is that their intrinsic dispersion becomes smaller from the optical to the NIR (Di Criscienzo et al. 2007). This means that NIR PL relations are not only more accurate because the reddening is less severe, but they are also intrinsically more precise.

Near-infrared photometry of T2Cs in the bulge obtained with the Son of ISAAC (SOFI) telescope was used by Groenewegen et al. (2008) to estimate the distance of the Galactic center (R₀ = 7.99 ± 0.09 kpc) using a sample of 38 T2Cs. Bhardwaj et al. (2017b) matched photometry of the VISTA Variables in the Vía Láctea (VVV) survey (Minniti et al. 2010; Saito et al. 2012) with the OGLE III version of the catalog of T2Cs (Soszyński et al. 2011, 335 T2Cs) and obtained individual distances. They estimated R₀ = 8.34 ± 0.03 kpc and ruled out a barred structure. All recent estimates of R₀ based on other diagnostics and recent reviews agree that the official IAU value of R₀ = 8.5 kpc is overestimated (R₀ = 8.2 ± 0.1 kpc, Bland-Hawthorn & Gerhard 2016; R₀ = 8.3 ± 0.2 ± 0.4 kpc, de Grijs & Bono 2016).

Recently, the Optical Gravitational Lensing Experiment (OGLE) IV survey (Udalski et al. 2015) has generated the largest homogeneous sample of T2Cs known to date, amounting to 924 objects projected toward the Galactic bulge (Soszyński et al. 2017). This is almost three times the size of the previous sample. The VVV survey has collected NIR K_s-band time series toward the Galactic bulge in a sky area that covers almost the entire OGLE survey area and provides an optimal framework to characterize the structure of the old population of the Galactic bulge with stellar tracers such as RRLs and T2Cs, for which the optical photometry and accurate periods are available from the OGLE survey. The increase of the sample size with respect to previous works is a unique opportunity to achieve new insight into the old stellar population in the bulge, especially for a detailed comparison with RRLs, which has always been hampered by the small sample size of T2Cs. Furthermore, we have the unprecedented opportunity to combine the T2C NIR catalog with the proper motion measurements from VVV itself (Contreras Ramos et al. 2017; Smith et al. 2018) and Gaia Data Release 2 (DR2; Gaia Collaboration 2016, 2018) to constrain the kinematic properties of the old stellar population.

The paper is organized as follows: in Sect. 2 we present our photometric and astrometric databases. We analyze the light curves and derive their properties in Sect. 3. Section 4 is dedicated to estimating individual distances of T2Cs and their overall distribution, while in Sect. 5 we discuss the kinematic properties of our targets. We discuss and summarize our results in Sect. 6.

2. Data

Light curves. We used the aperture photometry of VVV DR4 data (Minniti et al. 2010; Saito et al. 2012) that is publicly available through VISTA Science Archive (VSA)¹. We compared point spread function (PSF) with aperture photometry for a sample of our targets and found the differences to be negligible, but aperture photometry has the advantage of being available in the entire VVV survey area. As a first step, we matched the OGLE IV catalog of T2Cs (Soszyński et al. 2017) with the source detection catalog of VVV, adopting a matching radius of 2″. This allowed us to retrieve 894 of 924 targets within the VVV survey area. A posteriori, we checked that all the good matches are within 1″̣3 of the OGLE coordinates. Of those that were not retrieved, 25 are outside the VVV area, and for five of them we could not find a good match, even with a larger searching radius of 10″. Of these 894 T2Cs, according to the classification of Soszyński et al. (2017), 369 are BL Herculis (BLHs), 343 are W Virginis (WVs), 28 are peculiar W Virginis (pWVs), and 154 are RV Tauri (RVTs). We discuss the different types of T2Cs in more detail in Sect. 4.2.

We collected both the single-epoch JH-band photometry and the K_s-band time series. ZY photometry was neglected because it is not useful for our goals and we cannot even estimate mean magnitudes in these bands, since light-curve templates (Bhardwaj et al. 2017a) are available only for the JHK_s bands. The number of valid K_s-band phase points per variable (those with good photometric solution) ranges from four to 185, with a median of 51. Only four variables have fewer than ten phase points, and 90% of the variables have more than 47 observations, with a good phase coverage over the whole range of periods of our targets (1–85 days). The number of available K_s epochs varies across the VVV survey area, with the majority of the VVV pointings (so-called tiles) having between 50 and 100 epochs. Given the overlap between the adjacent tiles, which amounts to about 1′ on a side, a small fraction of our variables (34) were observed in two tiles, with up to 185 data points.

Reddening. We adopted the two Galactic bulge E(J − K_s) reddening maps of Gonzalez et al. (2012) and Schultheis et al. (2014), henceforth, G12 and S14. The reasons are manifold: (1) They were both obtained with VVV data, therefore no photometric system conversion is needed. (2) The map of G12 provides very high resolution especially in the central regions (2 × 2 arcmin for −3°̣5 < b < 5°̣0, 4 × 4 arcmin for −7°̣0 < b < −3°̣5 and 6 × 6 arcmin for −10°̣0 < b < −7°̣0), while the “pixels” of the S14 map have a size of 6 × 6 arcmin everywhere. However, the map by S14 has the key advantage to be three-dimensional. It provides a grid of E(J − K_s) for 21 bins of distance, from 0 to 10.5 kpc. This is crucial for studying a complex structure such as that of the Galactic bulge and possibly intervening thick-disk populations. In the following, we indicate as E(J − K_s)_G12 and E(J − K_s)_S14 the two reddening values obtained from the G12 and S14 maps, respectively. As explained in Sect. 3, we adopt the 3D reddening map of S14 for our final estimates, but G12 serves as a comparison.

Proper motions. We retrieved relative proper motions for 894 of 924 targets from the publicly available VIRAC catalog (Smith et al. 2018), obtained with VVV data. According to their recommendations, we discarded all targets for which the flag reliable, based on the validation of the photometric solutions, is equal to zero. We point out that although the cross-match was made by unique VVV ID and not by coordinate, we found multiple (either double, triple, or quadruple) records for 120 targets in VIRAC. We checked the multiple identifications one by one on the basis of right ascension (α), declination (δ), and K_s-band magnitude. About two-thirds of the time, all of the records of a multiple identification had reliable = 0 and were thus all discarded. Of the remaining fraction, an a posteriori check revealed that the majority had reliable = 1 for the correct match and reliable = 0 for the incorrect matches. This validates our selection of the best match.

We also retrieved relative proper motions for 416 targets from PSF photometry, obtained from VVV data by the method explained in Contreras Ramos et al. (2017). The match was performed using a searching radius of 2″̣. This catalog of proper motions does not cover the entire VVV area, but only the low latitudes (−3°̣0 ≲ b ≲ 3°̣0), that is, the most crowded region, where PSF photometry has several advantages over aperture photometry. A comparison between the two sets of proper motions from the VVV is performed in Sect. 5.

Finally, we searched our targets within the recent Gaia DR2 (Gaia Collaboration 2016, 2018), using a searching radius of 4″. We retrieved matches for 920 targets. Because of the density of Gaia DR2, a search radius of 4″ means multiple records for almost all our targets. We note that for Gaia, we needed a larger searching radius because α and δ are at the epoch J2015.5. We selected the best matches on the basis of the separation in (α,δ) and the V − G and I − G color indexes, and retrieved 914 targets (for the remaining 6 targets, it was not possible to select a best match based on the adopted criteria). Of these, 868 have a five-parameter Gaia solution (coordinates, parallaxes, and absolute proper motions).

3. Light curves and properties of T2Cs

We have phased the VVV light curves using the periods provided by OGLE IV (Soszyński et al. 2017). We visually inspected all the VVV light curves and separated promising from noisy or poorly sampled light curves. We estimated the uncertainty on the mean magnitude (eK_s) as the sum in quadrature of the standard error of mean of the phase points around the fit (see below) plus the median photometric error of the phase points. We note that the standard error on the mean and median photometric error represents a statistical and systematic measure of uncertainties, respectively,

$$\mathit{e}{\mathit{K}}_{\mathit{s}}=\sqrt{\frac{1}{\mathit{n}}\frac{{\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{n}}{\left({\text{mag}}_{\mathit{i}}-{\text{mag}}_{\mathit{i}\left(\text{fit}\right)}\right)}^2}}{\mathit{n}-1}+\text{madian}{\left(\text{err}\right)}^2.}$$

We adopted eK_s as a quantitative criterion to select good-quality light curves, with a threshold at eK_s < 0.05. We performed visual inspection to validate this threshold and further selected promising light curves within the range 0.05 < eK_s < 0.10 mag. However, some of the brighter stars (〈K_s〉 ≲ 11.5 mag) were classified as poor-quality light curves in this range (see below), and all stars with eK_s > 0.10 mag were also included in the poor-quality sample. Examples of good- and poor-quality light curves are displayed in Fig. 1.

Fig. 1. Left: example of a good-quality light curve. Right: example of a poor-quality light curve. The periods are labeled at the top of the panels, while the names, as provided in the OGLE catalog, are labeled at the bottom. The PLOESS fit of the good-quality light curve is displayed as a red line.

For the 161 targets with poor-quality light curves, we adopted for the mean magnitude the median of the magnitudes, converted into flux, of the individual phase points. The uncertainty on the mean magnitude is the standard deviation of the median. We do not provide a light-amplitude estimate for these targets since the uncertainty would be on the same order as the amplitude itself. We point out that most of these targets are very bright (135 of 161 have 〈K_s〉 ≲ 11.5 mag), meaning that they are either saturated or in the nonlinear regime of the camera. However, we also point out that not all stars in this magnitude range have poor-quality light curves. The cause of poor photometry is a synergy between saturation and crowding by other nearby bright stars. The magnitude for saturated sources is derived from the most external apertures (those with the largest radii). When another bright source contaminates the external aperture, the photometric solution is worse. In these cases, PSF photometry would not help to improve the photometry because it would be even more affected by saturation.

To compute the mean magnitudes, amplitudes, and relative uncertainties of the good-quality targets, we adopted the fits of the light curves obtained with the PLOESS (Braga et al. 2018) fitting method, which is a variant of the GLOESS method (Persson et al. 2004). The mean magnitude is the integral of the fitting curve converted into flux, while the amplitude is the difference between the brightest and faintest points of the fit. The uncertainty on the mean magnitude was defined before as eK_s. The uncertainty on the amplitude was derived as the sum in quadrature of the median photometric errors of the phase points around the maximum and minimum, plus the standard deviation of these phase points around the fit of the light curve. The final value was weighted with the number of phase points around the maximum and the minimum. We also derived mean J- and H-band magnitudes by applying the light-curve templates of Bhardwaj et al. (2017a) to the single J and H phase points and using the time of maximum from OGLE (Soszyński et al. 2017). The photometric properties of our final sample are listed in Table 1.

Next, we corrected the mean magnitudes for extinction. As stated in Sect. 2, we derived two extinction values from two different reddening maps. We obtained E(J − K_s)_G12 and A_Ks(G12) for each target from the online tool BEAM². However, estimating E(J − K_s)_S14 is not straightforward because the distance of the target needs to be known in order to obtain E(J − K_s)_S14. Distance and reddening in S14 are degenerate. Therefore, we adopted the following method. First, we located the four pixels of the map that are closest to the position of our target. We weighted the E(J − K_s)_S14 of each pixel by their inverse angular distance of the pixel center from our target. We repeated this for each of the 21 bins and obtained 21 possible values for E(J − K_s)_S14 for each target. Using an iterative method (Bhardwaj et al. 2017b), we simultaneously found the most plausible values of E(J − K_s)_S14 and A_Ks(S14) and an estimate of the distance.

We adopted a new reddening law with both reddening maps to derive the extinction A_Ks from E(J − K_s). This law was derived by Alonso-García et al. (2017) using VVV data for the innermost regions of the bulge between |l| < 2°̣7 and |b| < 1°̣55. Specifically, we adopted as the ratio of total-to-selective extinction ${\mathit{R}}_{\mathit{J}\mathit{K}}=\frac{{\mathit{A}}_{\mathit{K}}{}_{{}_{\mathit{s}}}}{\mathit{E}\left(\mathit{J}-{\mathit{K}}_{\mathit{s}}\right)}$ the values shown in their Table 2, according to the quadrant, within their surveyed sky area. For targets outside the quoted area, we adopted R_JK = 0.428, which is the average suggested by Alonso-García et al. (2017).

We checked the differences ΔA_Ks = A_Ks(G12) − A_Ks(S14) for all our targets, and found that ΔA_Ks values follow an almost Gaussian distribution with a mean of −0.007 mag and σ = 0.032 mag. A tail of targets with 0.1< ΔA_Ks < 0.25 mag and one target with ΔA_Ks as high as 0.63 mag were found. Almost all of these targets are located at distance moduli (DM) smaller than 14.5 mag (∼8 kpc). This means that they are closer than the Galactic center. On the other hand, targets with ΔA_Ks < −0.1 mag are mostly located at DMs larger than 14.5 mag. This is expected when comparing a reddening map fixed at ∼8 kpc (G12) with one that takes distance into account (S14). This also suggests that we can adopt the S14 map to derive our final results.

However, since reddening and distance were derived simultaneously, A_Ks(S14) may depend on the calibration selected for the PL relation. We checked that for different calibrations (see Sect. 4.2), the differences between A_Ks(S14) from different calibrations are within 0.05 mag and are smaller than 0.01 mag for ∼80% of the targets.

4. PL relations and distances

4.1. Empirical PL relation

Figure 2 shows the targets in the log P−K_s0 plane. The BLHs and WVs were dereddened adopting the S14 map, but A_Ks(S14) values and distances were estimated simultaneously. Since we did not estimate the distances of pWVs and RVTs, as explained below, these targets were dereddened with A_Ks(G12), which is independent of distance.

Fig. 2. Top: bulge T2Cs in the period-luminosity plane. Dark blue, light blue, and cyan circles mark BLHs, WVs, and RVTs, respectively. Red crosses mark pWVs. Bottom: close-up of the empirical PLKs0 relation obtained by fitting BLHs and WVs in our sample (black solid line). Larger symbols display BLHs and WVs that were kept after an iterative 3σ clipping procedure. The purple, red, and green dashed lines display the literature PLKs0 shifted by 14.6 mag (∼8.3 kpc) of Matsunaga et al. (2006), Groenewegen et al. (2008), and Bhardwaj et al. (2017a), respectively.

In stellar system like the bulge, where stars are not all at the same distance (as in the case of globular clusters, external galaxies, etc.), extinction and distance are degenerate and a simple empirical PL does not provide precise insight into the structure, especially if the reddening is not constant. Moreover, there is a debate on whether the old population in the bulge is indeed spheroidal (Dékány et al. 2013; Kunder et al. 2016) or if it is ellipsoidal and tilted, similarly to the Galactic bar (Pietrukowicz et al. 2015). Therefore, as done before by Groenewegen et al. (2008), we fit the PL relation by adding the dependence on the Galactic longitude (l) and latitude (b). We selected all the BLHs and WVs from our sample, and after an iterative rejection of outliers at more than 3σ, we found K_s0 = (10.66 ± 0.02) − (2.21 ± 0.03) #x00b7; (log P − 1.2) − (0.020 ± 0.003) · l + (0.050 ± 0.008) · |b| mag, with a standard deviation of 0.07 mag. The positive coefficient in |b| means that fainter stars are located at higher distances from the Galactic plane, where the reddening is lower. This is an evidence that the T2C sample is biased by reddening. On the other hand, the nonzero dependence of K_s0 on l indicates that the T2C ellipsoid is tilted. If we ignore the l and |b| terms, the simple PL relation is K_s0 = (10.76 ± 0.02) − (2.23 ± 0.03) · (log P − 1.2) mag, and the standard deviation increases to 0.28 mag. Figure 2 shows our newly derived empirical, coordinate-independent PL as a black solid line. We set the zero of the independent variable at log P = 1.2 to facilitate comparison with Groenewegen et al. (2008). The coefficients of both the coordinate-dependent and the simple PL relations agree remarkably well with those obtained with an identical approach by Groenewegen et al. (2008): all of them agree within 1σ. However, while the error on the coefficient of the l coordinate in Groenewegen et al. (2008) was larger than the value itself (−0.028 ± 0.031 mag°, thus including the zero value within 1σ), our coefficient is more precise (−0.019 ± 0.003 mag°) because the set of variables is much larger, and it clearly rules out a PL relation that is independent of l. We also compared our empirical PL with that of Bhardwaj et al. (2017b), K_s0 = (10.749 ± 0.0056) − (2.189 ± 0.032) · (log P − 1.2) mag). We adjusted their zero-point since they had adopted log P = 1.0 as the zero of their independent variable. The slope and zero-point both agree with ours within 1σ.

The bottom panel of Fig. 2 also shows PLs from the literature, Galactic globular clusters (Matsunaga et al. 2006), bulge (Groenewegen et al. 2008), Large Magellanic Cloud (LMC; Bhardwaj et al. 2017a) as purple, red, and green dashed lines, respectively. The PL relations of Matsunaga et al. (2006) and Bhardwaj et al. (2017a) were placed at a DM of 14.6 mag (∼8.3 kpc), which is the best recommended value of R₀ in the literature in general (Bland-Hawthorn & Gerhard 2016; de Grijs & Bono 2016) and a very common value obtained from methods based on the PL of RRLs and T2Cs (Dékány et al. 2013; Pietrukowicz et al. 2015; Bhardwaj et al. 2017b; Majaess et al. 2018). We note that the slope of our relation is the same within the uncertainties as the slope found by (Bhardwaj et al. 2017a, b_B17 = −2.21 ± 0.02 mag).

Finally, we note that RVTs are systematically fainter than predicted by the PLs, which would be K_s0 < 10 mag. This is an effect of saturation, which in the VVV is severe at magnitudes brighter than K_s ∼ 10.0 mag (Minniti et al. 2010; Mauro et al. 2013). This is also supported by the fact that as we discuss in Sect. 4.2, RVTs are expected to be either brighter than predicted by the PL (Matsunaga et al. 2009, 2011; Ripepi et al. 2015) or to follow the PL (Matsunaga et al. 2006; Bhardwaj et al. 2017a), but not to be fainter. We also checked for possible selection biases that would favor the selection of the closest RVTs, but found none. This further supports the saturation scenario for these variables.

4.2. Individual distances

Before deriving individual distances, we discuss a few key points below.

(1) Metallicity dependence. There is a general consensus about the independence, or a very mild dependence, less than 0.1 mag dex⁻¹ (Di Criscienzo et al. 2007), of the PL relation of T2Cs on metal abundance, based on both empirical (Matsunaga et al. 2006, 2011) and theoretical (Bono et al. 1997; Di Criscienzo et al. 2007) arguments. The advantage of adopting a metal-independent PL relation to calibrate distances is straightforward because the metallicity of our targets is unknown. Since T2Cs belong to the same stellar population as RRLs, their metallicity dispersion should be similar, and it should be fairly small (0.25 dex, Pietrukowicz et al. 2012).

(2) RV Tauri and peculiar W Virginis. T2Cs are separated into BLHs (1 day < P < 5 days), WVs (5 days < P < 20 days), and RVTs (P > 20 days). The thresholds are those by Soszyński et al. (2011). They are based on the period distribution of OGLE III T2Cs, and the authors kept them unchanged in the OGLE IV catalog (Soszyński et al. 2017). However, there is a debate as to whether RVTs obey the same PL as shorter-period T2Cs: Matsunaga et al. (2009) using data from the InfraRed Survey Facility (IRSF) and Ripepi et al. (2015) using VMC data, showed that in the LMC, RVTs are overluminous with respect to the extrapolation at long periods of the PL of BLHs and WVs. In contrast, Bhardwaj et al. (2017a) used NIR data in the central bar of the LMC (Macri et al. 2015) and found that RVTs also fall on the linear PL fit to short-period T2Cs in the LMC. However, their photometry for short-period variables (that is, at fainter magnitudes) was more prone to be affected by crowding in the central regions than longer-period variables. On the other hand, RVTs in Globular clusters also follow the same PL relation as shorter period T2Cs (Matsunaga et al. 2006). We do not discuss this in detail since this is not the aim of this work, but it is crucial to remember that in addition to these empirical findings, there are two different evolutionary channels from which RVTs are generated: either from low-mass (∼0.50 M_⊙), very old (>10 Gyr) AGB stars (Wallerstein 2002), or from more massive (≳1 M_⊙) and younger objects (Willson & Templeton 2009). Recent findings about RVTs also involve binarity, to distinguish their evolutionary channel (Manick et al. 2018). This means that RVTs are not reliable as distance indicators, and we did not take them into account for our PL relations. Finally, following the classification in Soszyński et al. (2017), we also discarded pWVs, which are a subclass of WVs that likely belong to binary systems (Soszyński et al. 2008). They are overluminous (0.3– 0.5 mag in the K_s band Ripepi et al. 2015) when compared to WVs with similar periods. This means that they do not follow the PL relation of T2Cs and cannot be used as distance indicators. We note that pWVs in Fig. 2 are not more luminous than WVs at the same period. Nonetheless, we did not take them into account for the quoted reasons.

(3) No semi-empirical calibration. In principle, it is not possible, in our case, to adopt a semi-empirical calibration (empirical slope from our own sample and zero-point from literature) for the PLs because our targets are not all located at the same distance. This is true in principle, even though our empirical slope is identical within 0.01 mag to that by Bhardwaj et al. (2017a), as stated in Sect. 4.1.

(4) Selection of the calibration. Keeping these points in mind, we have searched the literature, where several calibrations of the PL relation of T2Cs in the NIR are available (Matsunaga et al. 2006, 2009; Di Criscienzo et al. 2007; Ripepi et al. 2015; Bhardwaj et al. 2017a; Gaia Collaboration 2017, henceforth PL_M06, PL_D07, PL_M09, PL_R15, PL_B17 and PL_C17). We rule out the preliminary calibration PL_C17, based on Gaia DR1 parallaxes, because it yields DM_LMC >18.9 mag, which does not agree with literature values nor with the Gaia DR1 calibration of RRLs and Classical Cepheids derived in the same work. We note that PL_M09, PL_R15 and PL_B17, which are all based on T2Cs in the LMC, do not provide absolute calibrations of their PLs. Therefore, we adopted the distance of the LMC based on eclipsing binaries (Pietrzyński et al. 2013, DM_LMC = 18.493 ± 0.008 ± 0.047 mag) to set the zero-point. In principle, we could calibrate the zero-point using Baade-Wesselink parallaxes for field T2Cs obtained by Feast et al. (2008), but this would be based only on two objects, with a strong difference among the different values of the parallaxes in the literature. After checking each of these calibrating PLs, we decided to adopt PL_B17. The choice was guided by the following reasons. First, the calibrating PL must be based on both BLHs and WVs. This rules out PL_D07 since they only used BLHs pulsation models. Second, of the three PLs from the LMC (PL_M09, PL_R15, and PL_B17), the latter is based on light curves with an average of 50 epochs, which is more than twice that of the other two together. Moreover, for targets that are outside their surveyed sky area, they include data from the previous works. Third, we excluded the M06 calibration because it is based on all the three subclasses of T2Cs, including RVTs, which we did not include in our distance analysis. In the end, we adopted PL_B17, which is based only on BLHs and WVs, which are the subclasses of variables for which we did estimate the distance.

As stated in Sect. 2, we simultaneously estimated distances and the 3D extinction E(J − K_s)_S14. This was made iteratively, following the method of Bhardwaj et al. (2017b). Normally, after the second or third iteration, the values of distance and reddening converge. We derived individual distances for 710 variables, which are listed in Table 2. Taking all sources of uncertainty into account (uncertainty on the mean magnitude, on the extinction, on the coefficients of the calibrating PL, and the intrinsic width of the PL), the relative uncertainties are in the range 8%–9% for 676 objects, 9%–20% for 38, and only one, at a distance of 2.87 kpc, has a relative uncertainty of ∼27%.

4.3. Spatial distribution of Type II Cepheids

Knowing the coordinates (l,b) and the distances d, we derived the coordinates x_GAL, y_GAL and z_GAL. We adopted a reference frame with the Sun at the origin (x_GAL = y_GAL = z_GAL = 0), x_GAL which increases toward the Galactic center, y_GAL which is on the Galactic plane and positive for l > 0, and z_GAL which is perpendicular to the Galactic plane and positive for b > 0.

Figure 3 clearly shows that several objects are located either in front of or beyond the bulge. The individual distances range from 2.0 to 111.7 kpc, with 11 objects more distant than 30 kpc, most likely belonging to the outer halo, which dominates at Galactocentric distances greater than ∼20 kpc (Carollo et al. 2007, 2018).

Fig. 3. Top: projection, on the Galactic plane, of T2Cs within 14 kpc. The Sun is at the vertex of the plot. Black circles with error bars represent the means and standard deviations of the distances projected onto the Galactic plane, in 2° wide angular areas (4° wide for the most peripheral one at positive l). Thirty-four T2Cs were not plotted because they are more distant than 14 kpc. Bottom: same as the top panel, but projected in latitude. The angular areas are 2° wide, except the central one, which is 3° wide.

The two panels of Fig. 3 also show as black circles with error bars the means and standard deviations of the distances projected onto the Galactic plane (top panel) or onto the longitude (l = 0) plane (bottom panel), in angular areas of 2°. We point out that in the top panel, the most peripheral area at positive l is 4° wide, and in the bottom panel, the most central area is 3° wide, to take into account the lower density of objects. These mean and standard deviations were derived by fitting a rescaled histogram with a Gaussian. A detailed explanation is provided below.

Galactic (b = 0) plane. The central angular areas share basically identical averages. However, the most peripheral areas (5°̣0 < |l| < 9°̣0, that is, between 0.7 and 1.25 kpc from the Galactic center) show a slight deviation similar to that of the barred distribution of Red Clump stars (Wegg & Gerhard 2013): closer at positive l and more distant at negative l. This agrees with our finding that the PL relation does depend on the l coordinate, and gives further evidence of the ellipsoidal symmetry of the old stellar component of the bulge, as has been outlined by Pietrukowicz et al. (2015), who found an ellipsoidal structure for bulge RRLs, with an inclination of the major axis with respect to the Sun of 20° ± 3°, similar to the orientation of the Galactic bar (∼30°Gonzalez et al. 2011; Wegg & Gerhard 2013).

Longitude (l=0) plane. A distinct yet fictitious trend of the average distance with b is displayed in the bottom panel of Fig. 3. The variables that are closer to the Galactic plane appear to have smaller distances. This is due to a selection effect caused by extinction in the original OGLE catalog, which more easily detects stars in the closest part of the bulge than stars in the farthest part, which are more heavily reddened. The extinction ratio between the I band of OGLE and the K_s band of VVV ranges within a factor of four to ten, and the lower limit of A_Ks(S14) for our targets is ∼0.3 mag. We have considered two other possible explanations for this trend: either an overestimate of reddening at low b or the effect of crowding on aperture photometry, but none can explain the quoted behavior. We discard the possibility that reddening is overestimated because other experiments that adopted either the G12 or the S14 map, which are consistent between themselves, provided estimates of R₀, all at about 8.3 kpc (Gonzalez et al. 2012; Bhardwaj et al. 2017b; Majaess et al. 2018). We also rule out the possibility that crowding affects the magnitudes from aperture photometry, making the targets brighter. A direct comparison of PSF versus aperture photometry does not reveal any clear trend with b, and the average difference of mean magnitudes is ΔK_{s(Aperture−PSF)} = 0.03 ± 0.013 mag.

With individual distances and coordinates for our targets, we can estimate the distance of the Galactic center R₀. However, the calculation is not straightforward, and cuts and resampling are needed to take the biases into account.

First, we selected only stars at ${\mathit{R}}_{\mathit{G}}=\sqrt{{\mathit{x}}_{\text{GAL}}^2+{\mathit{y}}_{\text{GAL}}^2,}$ which is the distance of the star, projected onto the Galactic plane, between 6 and 11 kpc, to avoid non-bulge stars within the sample. The choice is justified by the results of (Pietrukowicz et al. 2015, see their Fig. 5), who showed that the density of bulge stars is very low (lower than ∼10% of the peak) outside this distance range. Second, we only selected stars at b > 3°̣0 and b < −3°̣0 to avoid the OGLE selection bias. A similar cut was applied by Pietrukowicz et al. (2015), who only used RRLs at b < −2°̣7 to estimate R₀.

These selection criteria left us with 172 stars. However, their average longitude (〈 l 〉) is 0°̣275, which means that the sample is biased toward shorter distances. To overcome this bias, we performed a resampling of the data by randomly selecting 75 stars at negative l and 75 at positive l. If 〈l〉 of the 150 random targets is lower in absolute values than 0°̣1, we kept the sample, otherwise, we repeated the random target selection. We point out that this process is not a proper bootstrap method because we did not allow sampling the same element more than once.

With this set of 150 targets, we plot the distribution of R_G in bins of 0.25 kpc as shown in Fig. 4 (black histogram). However, this distribution is biased and shifted to greater distances. At fixed coordinates (l,b), the volume within a given sky area (Δl,_Δb) and a given depth range (Δd) increases with distance. This means that the number of stars within the volume (and therefore the probability of detecting a target in the volume) increases quadratically with distance. This causes a bias that shifts the distribution toward the more probable larger distances. To take this geometric effect into account, we scaled the distribution by d⁻². We fit the scaled distribution (red histogram in Fig. 4) with a Gaussian. We estimated the abscissa of the peak (x₀) and adopted it as our estimate of R_0(i) on the ith resampled set. Starting from the random extraction of 150 targets, this process was repeated 5000 times to avoid any selection bias.

Fig. 4. Top: overall histogram of the distances of 150 targets resampled 5000 times (750 000 in total, shown in black). d−2 scaled histogram (red). A Gaussian centered at R0 is displayed. The estimate of R0 is labeled with its uncertainty as derived by using percentiles, as described in the text. Bottom: d−2 scaled histogram of 75 targets at l < 0° resampled 5000 times (375 000 in total, shown in blue). Light blue: same as the blue, but for targets at l > 0°.

We adopted an overall average of the 5000 estimates of R_0(i) as our final estimate of R₀. Based on a sample with 〈l〉 = 0°̣016 ± 0°̣050, we obtain a final R₀ estimate of 8.46 kpc. The statistical uncertainty of both R₀ was derived as half of the range between the 15.8% and 84.1% percentiles of the distribution of R_0(i) (see Fig. 5). These thresholds were chosen to enclose 68.3% of the estimates provided by the simulations, like a ± 1σ range in a Gaussian distribution. We derived a range of 8.43–8.49 kpc for R₀, which means a statistical uncertainty of 0.03 kpc.

Fig. 5. Distribution of the 5000 estimates of R0 from the resampled sets of targets. The dashed lines display the percentiles at 15.8% and at 84.1% that we used to derive the uncertainty on R0.

We calculated the systematic uncertainty as the squared sum of the average uncertainty on the mean magnitude (0.020 mag), the average uncertainty on the extinction (0.021 mag), and the average propagation of the uncertainties of the calibrating PL coefficients (almost vanishing, 0.001 mag). These together are 0.028 mag, which is 0.11 kpc at 8.30 kpc. Our estimate of R₀ does not agree very well with estimates from similar works, either using T2Cs (8.34 ± 0.03[stat.] ± 0.41[syst.] kpc, Bhardwaj et al. 2017b) or RRLs (8.33 ± 0.05[stat.] ± 0.14[syst.] kpc, Dékány et al. 2013; 8.27 ± 0.01[stat.] ± 0.40[syst.] kpc, Pietrukowicz et al. 2015). However, all the quoted papers adopted the reddening law by Nishiyama et al. (2009), which provides a higher R_JK = 0.528 and, in turn, smaller distances. Had we adopted the Nishiyama et al. (2009) reddening law with our data, it would have provided R₀ = 8.30 ± 0.03(stat.) ± 0.11(syst.), which would agree perfectly well with the quoted papers. This is evidence of how crucial a correct understanding of the reddening law is. Our estimate of R₀ agrees within 1σ with the best overall recommended value from a recent review (∼8.3 ± 0.2[stat.] ± 0.4[syst.] kpc, de Grijs & Bono 2016).

Finally, a more detailed analysis of the data allows us to show further evidence of the asymmetrical distribution of T2Cs around the Galactic center. By resampling 5000 times only T2Cs at positive l, we obtain an average peak of the distribution of 8.29 ± 0.09 kpc. The same process on T2Cs at negative l provides a value of 8.68 ± 0.05 kpc. Together with the distribution of average distances in the Galactic plane (top panel of Fig. 3) and the dependence of the PL relation on l, this is strong evidence that T2Cs trace an old, ellipsoidal stellar population.

5. Kinematics

5.1. Proper motion of the center of mass

The proper motion of Sgr A^✻, the supermassive black hole at the center of the Milky Way, based on VLBA measures, is μ_{l∗(Sgr A^✻)} = −6.379 ± 0.026 mas yr⁻¹; μ_{b(Sgr A^✻)} = −0.202 ± 0.019 mas yr⁻¹ (Reid & Brunthaler 2004). Assuming that the center of mass of the old population traced by T2Cs overlaps with Sgr A^✻ and has the same proper motion, we adopt the proper motions of T2Cs to obtain an indirect estimate of the proper motion of the center of mass.

As discussed in Sect. 2, we collected proper motions from three different catalogs: VIRAC, PSF, and Gaia. For the analysis in this section, we rejected, from all three catalogs proper motions with a combined statistical error $\text{CSE}=\sqrt{\text{err}{\mathit{\mu }}_{\mathit{\alpha }✻}^2+\text{err}{\mathit{\mu }}_{\mathit{\delta }}^2}$ or $\text{CSE}=\sqrt{\text{err}{\mathit{\mu }}_{\mathit{l}✻}^2+\text{err}{\mathit{\mu }}_{\mathit{b}}^2}$ larger than 2 mas yr⁻¹, leaving 553, 343, and 837 targets from VIRAC, PSF, and Gaia, respectively. We point out that the error propagation for the Gaia proper motions, when converting from (μ_α∗,μ_δ) into (μ_l,μ_b), was performed taking into account the covariance terms as suggested by Luri et al. (2018). Gaia coordinates were precessed from their native J2015.5 epoch to J2000, the same reference epoch as for PSF, to perform the quoted conversion.

It is crucial to remember that while VIRAC and PSF provide relative proper motions in the frame of reference of the Galaxy, Gaia provides absolute proper motions in the practically inertial frame of reference defined by quasars. This allows an interesting comparison among the catalogs. As a first step, we left out VIRAC proper motions. As displayed in Figure 6, the distribution of Δμ_l∗ = μ_l∗(Gaia) − μ_l∗(PSF) is centered at Δμ_l∗(peak) = −6.41 ± 0.02 mas yr⁻¹. For the b component, we find Δμ_b∗(peak) = 0.12 ± 0.03 mas yr⁻¹. These numbers were derived using 251 T2Cs for which we have both Gaia and PSF proper motions and that are located within 2 kpc from the center of the Galaxy, as derived in Sect. 4.3. The latter criterion was adopted as a compromise to leave out possible thick-disk stars and to retain a large sample of targets.

Fig. 6. Left: distribution of Δμl∗ for 251 T2Cs located less than 2 kpc away from the center of the Galaxy. A Gaussian fit to the distribution is shown in red. The abscissa of the peak and the σ of the fit are labeled. Right: same as left panel, but for Δμb.

When we assume that the velocities of our targets are randomly distributed around the center of mass, which is reasonable because this is the behavior of old, low-mass bulge stars (Spaenhauer et al. 1992; Minniti 1996; Babusiaux et al. 2010; Kunder et al. 2016), then (Δμ_l∗(peak),Δμ_b(peak)) is an indirect estimate of the proper motion of the center of mass of the Galaxy. This is supported by the fact that Δμ_l∗(peak) is identical well within 1σ to the longitudinal component of the proper motion of Sgr A^✻. The latitudinal component Δμ_b(peak) does not agree with that of Sgr A^✻, but the mean uncertainty on the proper motions (0.82 mas yr⁻¹ for PSF and 0.47 for Gaia) is larger than the offset. We also checked whether there is a trend of μ_l∗ and μ_b with distance, and found nothing significant in either the Gaia or the PSF sample. This implies that there is no evidence of net rotation.

We have performed the same analysis using proper motions from VIRAC and Gaia and found that in this case, Δμ_l∗(peak) = −5.43 ± 0.03 mas yr⁻¹ and Δμ_l∗(peak) = −0.04 ± 0.02 mas yr⁻¹. While the b component is similar to that of Sgr A^✻, the component along l is different from μ_{l∗(Sgr A∗)} by as much as ∼1 mas yr⁻¹. Moreover, we found that for 265 targets in common between the VIRAC and PSF datasets, the difference between the medians of μ_l∗ and μ_b, are −0.97 ± 2.62 and +0.12 ± 2.62 mas yr⁻¹, respectively. No dependence of μ_l∗ on l or on b was found in either of the two catalogs. We conclude that the μ_PSF are more reliable than μ_VIR not only a priori, as discussed in Sect. 2, but also a posteriori. Because these are both relative measures of the proper motion, the assumption of randomly distributed motions around the Galactic center would imply a zero median μ_l∗. The comparisons with Gaia and that with the PSF sources both indicate that VIRAC proper motions are systematically shifted by ∼1 mas yr⁻¹ to the east on the Galactic plane.

5.2. Velocity dispersion

Tangential velocities can be used to test the kinematic properties of a stellar population. If bulge T2Cs trace a spheroidal population that is kinematically hot, then their distribution should be dominated by a velocity dispersion with negligible rotation such as the Bulge RRL population (Minniti 1996; Kunder et al. 2016; Marconi & Minniti 2018; Contreras Ramos et al. 2018; and references therein). In this case, the velocity ellipsoid should be fairly symmetric in the Galactic longitude (v_l∗) and latitude (v_b) components. To test this hypothesis, we derived the tangential velocity (v_t) and its two components v_l∗ and v_b, using the classical relation v_t = 4.74 · d · μ, where d is in kpc and and μ in mas yr⁻¹. We used d as derived in Sect. 4.2 and μ from Gaia because it is the most complete sample with the smallest uncertainties. While pWVs and RVTs might have valid values of μ,we cannot derive their v_t, since d is not available. Table 3 displays the proper motions and velocities.

Finally, we derived the standard deviation of v_l∗ and v_b and obtained 106.2 ± 3.0kms⁻¹ and 97.3 ± 5.5kms⁻¹, respectively. These values were corrected by subtracting, in quadrature, the average uncertainties on v_l∗ (20.5 km s⁻¹) and v_b (9.9 km s⁻¹). After this correction, we find σv_l∗ = 104.2 ± 3.0kms⁻¹, and σv_b = 96.8 ± 5.5kms⁻¹. The agreement between the two is better than 1.2σ, thus providing further evidence that the T2Cs belong to a kinematically hot population. We cannot rule out the possibility that despite our cuts, we still included some thick-disk objects. Stricter selections concerning the target distance from the center are hampered by the sample size.

6 Conclusions

We have retrieved K_s-band light curves from VVV aperture photometry for 894 of 924 T2Cs in the OGLE IV catalog (Soszyński et al. 2017). We calculated mean magnitudes and amplitudes based on PLOESS fits (Persson et al. 2004; Braga et al. 2018) to the light curves. For BLHs and WVs, we simultaneously estimated individual extinctions and distance moduli, based on a 3D reddening map (Schultheis et al. 2014) and on a PL relation. The calibration of the PL relation was based on the slope and zero-point of T2Cs in the LMC (Bhardwaj et al. 2017a), anchored with a late-type eclipsing binaries distance to the LMC (Pietrzyński et al. 2013). We found distances ranging from 2.0 to 111.7 kpc, which means that our objects are located in the bulge, in the inner and outer halo, and possibly in the thick disk. The mean individual relative uncertainty is 8.6%, independent of distance and with a small standard deviation of 1.2%.

The distribution of the individual distances, taking various geometric and selection biases into account, provides an estimate of the distance of the Galactic center R₀ of 8.46 ± 0.03(stat.) ± 0.11(syst.), which agrees with the recommended value of 8.3 ± 0.2(stat.) ± 0.4(syst.) kpc (de Grijs & Bono 2016). Our estimate of R₀ does not agree with other estimates with similar methods (R₀ ≈ 8.30 kpc Dékány et al. 2013; Pietrukowicz et al. 2015; Bhardwaj et al. 2017b), but the difference is consistent with the different reddening law that was adopted (Alonso-García et al. 2017 instead of Nishiyama et al. 2009).

Fig. 7. Left: distribution of the l component of vt for 164 T2Cs located at less than 2 kpc from the center of the Galaxy, as in Sect. 4.3. The simple and true standard deviations of the sample are labeled. Right: same as left panel, but for the b component of vt.

We provided solid evidence that the old stellar population in the bulge is ellipsoidal. First, we found a non-negligible dependence of the PL relation on the l coordinate. This has been described before by Groenewegen et al. (2008), but their limited sample hampered the precision of the coefficient (−0.028 ± 0.031 mag°), while ours is more precise (−0.019 ± 0.003 mag°). Second, we found that at l ≲ −5°, the average distance is larger, while at l ≳ 5°the average distance of T2Cs is smaller, on a map projected onto the Galactic plane. Third, which is a similar but more quantitative approach as the second point, we found that the distribution of T2Cs at positive l is centered at 8.29 ± 0.09 kpc, while that of T2Cs at negative l is centered at 8.68 ± 0.05 kpc.

We also adopted proper motions from both Gaia and VVV itself to constrain the kinematic properties of T2Cs in the bulge. The analysis was restricted to only the sources with a combined statistical error smaller than 2 mas yr⁻¹. The power of the synergy between Gaia and VVV astrometric data is clear when comparing the absolute proper motions from Gaia with relative proper motions from the VVV. The mean difference (−6.41 ± 0.02 and 0.12 ± 0.03 mas yr⁻¹ in the longitude and latitude direction, respectively) for T2Cs within 2 kpc from the Galactic center is similar within the uncertainties (0.82 mas yr⁻¹ for PSF and 0.47 mas yr⁻¹ for Gaia) to the VLBA estimate of the relative proper motion of Sgr A^✻ (−6.379 ± 0.026 and −0.202 ± 0.019 mas yr⁻¹). This is reasonable if we assume that the T2Cs of the bulge belong to the kinematically hot, old stellar population (Minniti 1996; Kunder et al. 2016). Another piece of evidence supporting this assumption is that the velocity dispersion in both the longitude and latitude directions agree within almost 1σ (σv_l∗ = 104.2 ± 3.0kms⁻¹, σv_b = 96.8 ± 5.5kms⁻¹). The difference may be due to contamination by thick-disk stars in the 2 kpc sphere around the Galactic center.

It is important to note that while the distribution and kinematics of metal-rich populations in the bulge, tracing the X-shaped structure, have been studied widely, the distributions of more metal-poor populations based on different tracers remain to be investigated in detail. This work on T2Cs provides results that are consistent with RRLs. The spectroscopic follow-up of these objects in the near future will allow us to confirm the differences in their spatial distributions and kinematics to those of metal-rich populations in the Galactic bulge.

Note added in proof. Recent investigations based on the position of T2Cs in the Magellanic Clouds, showed that WVs might not be exclusively old stars, but rather a mix of intermediate-age and old stars Iwanek et al. 2018. We have checked that the tilted distribution found in Sect. 4.3 does not substantially change if traced only by BLHs (purely old) and by WVs, although the tilt is more clear when traced by the latter.

Acknowledgments

We are grateful to the anonymous referee for the useful suggestions that helped us to improve our paper. We thank Prof. L. Smith, who helped us to point out important properties of VIRAC. 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.

References

All Tables

All Figures

	Fig. 1. Left: example of a good-quality light curve. Right: example of a poor-quality light curve. The periods are labeled at the top of the panels, while the names, as provided in the OGLE catalog, are labeled at the bottom. The PLOESS fit of the good-quality light curve is displayed as a red line.
In the text

Fig. 2. Top: bulge T2Cs in the period-luminosity plane. Dark blue, light blue, and cyan circles mark BLHs, WVs, and RVTs, respectively. Red crosses mark pWVs. Bottom: close-up of the empirical PLKs0 relation obtained by fitting BLHs and WVs in our sample (black solid line). Larger symbols display BLHs and WVs that were kept after an iterative 3σ clipping procedure. The purple, red, and green dashed lines display the literature PLKs0 shifted by 14.6 mag (∼8.3 kpc) of Matsunaga et al. (2006), Groenewegen et al. (2008), and Bhardwaj et al. (2017a), respectively.

In the text

Fig. 3. Top: projection, on the Galactic plane, of T2Cs within 14 kpc. The Sun is at the vertex of the plot. Black circles with error bars represent the means and standard deviations of the distances projected onto the Galactic plane, in 2° wide angular areas (4° wide for the most peripheral one at positive l). Thirty-four T2Cs were not plotted because they are more distant than 14 kpc. Bottom: same as the top panel, but projected in latitude. The angular areas are 2° wide, except the central one, which is 3° wide.

In the text

Fig. 4. Top: overall histogram of the distances of 150 targets resampled 5000 times (750 000 in total, shown in black). d−2 scaled histogram (red). A Gaussian centered at R0 is displayed. The estimate of R0 is labeled with its uncertainty as derived by using percentiles, as described in the text. Bottom: d−2 scaled histogram of 75 targets at l < 0° resampled 5000 times (375 000 in total, shown in blue). Light blue: same as the blue, but for targets at l > 0°.

In the text

	Fig. 5. Distribution of the 5000 estimates of R0 from the resampled sets of targets. The dashed lines display the percentiles at 15.8% and at 84.1% that we used to derive the uncertainty on R0.
In the text

	Fig. 6. Left: distribution of Δμl∗ for 251 T2Cs located less than 2 kpc away from the center of the Galaxy. A Gaussian fit to the distribution is shown in red. The abscissa of the peak and the σ of the fit are labeled. Right: same as left panel, but for Δμb.
In the text

	Fig. 7. Left: distribution of the l component of vt for 164 T2Cs located at less than 2 kpc from the center of the Galaxy, as in Sect. 4.3. The simple and true standard deviations of the sample are labeled. Right: same as left panel, but for the b component of vt.
In the text