The Hubble Tension Revisited: Additional Local Distance Ladder Uncertainties

In a recent paper, we investigated possible systematic uncertainties related to the Cepheid color-luminosity calibration method and their influence on the tension between the Hubble constant as inferred from distances to Type Ia supernovae and the cosmic microwave background as measured with the Planck satellite. Here, we study the impact of other sources of uncertainty in the supernova distance ladder, including Cepheid temperature and metallicity variations, supernova magnitudes and GAIA parallax distances. Using Cepheid data in 19 Type Ia supernova host galaxies from Riess et al (2016), anchor data from Riess et al (2016, 2019, 2021) and a set of re-calibrated Milky Way Cepheid distances, we obtain $H_0 = 71.9\pm 2.2$ km/s/Mpc, $2.0\,\sigma$ from the Planck value. Excluding Cepheids with estimated color excesses $\hat{E}({\rm V} - {\rm I})=0.15$ mag to mitigate the impact of the Cepheid color-luminosity calibration, the inferred Hubble constant is $H_0 = 68.1\pm 2.6$ km/s/Mpc, removing the tension with the Planck value.


INTRODUCTION
Cepheid stars are crucial in building up the distance ladder to Type Ia supernovae (SNIa) when estimating the Hubble constant H 0 . To be used as standard candles, Cepheids need to be calibrated with respect to fact that • long period Cepheids are brighter, • red Cepheids are dimmer, • and Cepheids in high metallicity environments may be brighter.
The Cepheid color-luminosity (C-L) correlation may be understood as a combination of extrinsic dust extinction and intrinsic temperature variations Corresponding author: Edvard Mörtsell edvard@fysik.su.se (Pejcha & Kochanek 2012). Given the difficulty in separating these effects, Cepheids are usually calibrated using a phenomenological approach where a parameter, R, corrects for both dust and intrinsic variations (Madore 1982). The correction can be applied to observed colors, as by the SH0ES team in e.g. Riess et al. (2016), or estimated color excesses, as in Follin & Knox (2018), derived by subtracting a model for the mean intrinsic Cepheid color. Regardless of which method is employed, it will by necessity correct for both extrinsic and intrinsic color variations since the observed color and estimated color excess only differ by a linear term depending on the Cepheid period. When assuming a global value for R, the choice of calibration method is of minor importance since only the difference between Cepheid magnitudes in anchor and SNIa host galaxies matters for the H 0 inference. When allowing for the calibration parameter to vary between galaxies, motivated by the observed variation in dust properties (Mortsell et al. 2021), the choice of calibra-tion can potentially make a difference. As demonstrated in Mortsell et al. (2021), using Wesenheit H-band magnitudes calibrated with respect to the observed color R W (V − I) where a global R W is allowed to vary gives H 0 = 73.2 ± 1.3 (in units of km/s/Mpc implied from now on), whereas calibrating using estimated color excesses with individual galactic R EÊ (V − I) yields H 0 = 73.9 ± 1.8, a 3.4 σ tension with the value inferred from the cosmic microwave background (CMB) observations with the Planck satellite (Aghanim et al. 2020).
Arguably being closest to a physical Cepheid model as presented in Pejcha & Kochanek (2012), throughout the paper, we will use the color excess calibration as our default method for investigation 1 . We first point out that the Cepheid color excess distributions have substantial variations between galaxies, and that the inferred H 0 is sensitive to the introduction of color cuts in the data. Next, we investigate the impact of Cepheid temperature variations on the inferred H 0 uncertainty and discuss how systematic effects connected to the metallicity calibration may shift the inferred H 0 and increase the error budget. We express concerns about the reliability of the Milky Way (MW) Cepheid parallax distance calibration and show that the local Cepheid calibrated SNIa magnitude uncertainties may be underestimated. Finally, we combine the different systematic effects using Monte Carlo simulations and derive constraints on the Hubble constant.

Cepheid Calibration
We use the Hubble Space Telescope (HST) nearinfrared flux (H = F160W), color calibrated using optical (V = F555W and I = F814W) data, to derive Here,Ê (V − I) represents a proxy for the color excess the intrinsic Cepheid color. It is obtained by subtracting an estimate of the mean intrinsic colors, V − I 0 from the observed colors. We use mean intrinsic colors from Tammann et al. (2011), including the quoted uncertainties and a 0.075 dispersion in the mean intrinsic Cepheid color between galaxies, from the difference between Large Magellanic Cloud (LMC) and MW 1 The SH0ES team have presented an updated Hubble constant value of H 0 = 73.04 ± 1.04 (4.9 σ Planck tension) using an, yet not publicly available, expanded SNIa and Cepheid data set (Riess et al. 2022). As in Mortsell et al. (2021), we here make use of the data from Riess et al. (2016Riess et al. ( , 2019Riess et al. ( , 2021 and Breuval et al. (2020).
Cepheids, see Table 1. These colors are in good agree-  Pejcha & Kochanek (2012). We fit for the values of R E that minimize the scatter in m W H . The Wesenheit magnitude of the jth Cepheid in the ith galaxy, including the anchor galaxies MW, NGC 4258 and the Large Magellanic Cloud (LMC), is modeled as with [M/H] i,j the Cepheid metallicity, M W H the absolute Cepheid magnitude normalized to a period of P = 10 days and Solar metallicity and µ i the distance modulus to the ith galaxy. We use separate P-L relations for short and long period Cepheids where [P] s i,j = 0 for Cepheids with periods > 10 days and [P] l i,j = 0 for Cepheids with periods < 10 days.

Type Ia Supernovae
The apparent SNIa B-band peak magnitude in the ith host, corrected for the width-luminosity and C-L relations using the SALT2 model (Guy et al. 2007), is modeled by The errors on m B,i include the fitting uncertainty and a 0.1 magnitude contribution to take into account the intrinsic SNIa dispersion added in quadrature (Riess et al. 2016;Mortsell et al. 2021).

Data and Parameter Fitting
For Cepheids in M31 and beyond, we use data from Table 4 in Riess et al. (2016). For Cepheids in the LMC, we use data in Table 2 in Riess et al. (2019). Data for MW Cepheids, including GAIA parallax measurements are from Table 1 in Riess et al. (2021). 20 double eclipsing binaries (DEBs) observed using long-baseline near-infrared interferometry give a distance to the LMC of µ LMC = 18.477 ± 0.0263 mag (Paczynski 1996;Pietrzyński et al. 2019;Riess et al. 2019). Observations of mega-masers in Keplerian motion around its central super massive black hole give a distance to NGC 4258 of µ N4258 = 29.397 ± 0.032 mag (Reid et al. 2019).
Type Ia SN B-band magnitudes are from Table 5 in Riess et al. (2016), derived using version 2.4 of SALT II (Betoule et al. 2014).
Given the observed Cepheid magnitudes m H , color ex-cessesÊ (V − I), periods [P], metallicities [M/H], together with the SNIa magnitudes m B , the anchor distances µ k and the MW Cepheid parallaxes π, we fit simultaneously for R E , b W , Z W , the host galaxy distances µ i , the anchor distances µ k , the GAIA parallax offset zp, the Cepheid absolute magnitude M W H and the SNIa absolute magnitude M B . For linear parameters, the fit can be made analytically. The exception is the R E ; since the uncertainties in the observed Cepheid colors are non-negligible, the uncertainties in the derived Wesenheit magnitudes m W H will depend on the values of R E , requiring a non-linear treatment of these parameters using Markov chain Monte Carlo (MCMC) techniques, see Mortsell et al. (2021).

Default Results for Color excess Calibration
Using similar assumption as in Mortsell et al. (2021), allowing for R E to vary between galaxies, with weak flat prior constraints on their values R E = [0, 1], we obtain H 0 = 73.9 ± 1.5, in 4.0 σ tension with the Planck value. For a fixed global value of R E = 0.386, following Riess et al. (2016), we obtain H 0 = 73.0 ± 1.3 (4.0 σ). These are the default results to which we will relate the impact of other systematic effects, the latter having the advantage of only including linearly fitted parameters, substantially reducing the computational cost when comparing different scenarios. Note that the error bars do not yet include the intrinsic Cepheid color scatter, see Section 8.

COLOR EXCESS DISTRIBUTIONS AND CUTS
As discussed in Mortsell et al. (2021), the fact that Cepheid colors and periods are correlated (Tammann et al. 2011) may cause color selection effects related to the fact that longer period Cepheids are brighter. From Figure 1, showing the different distributions of the color excessÊ(V − I), large variations between the different galaxies employed in the Cepheid distance calibration are evident. This fact, possibly attributed to selection effects, may introduce systematic errors when comparing the derived Cepheid distances. MW are systematically redder compared to other anchors and SNIa host galaxies, and also have a strong M31 Host correlation of long period Cepheids having larger color excesses. For each individual MW Cepheid, we compare the estimated color excess with three dimensional dust maps based on GAIA, Pan-STARRS 1 and 2MASS data where applicable (Green et al. 2019;Green 2018). The generally good agreement as shown in Figure 2 confirms that the majority of the color excess is due to interstellar dust extinction, but with room for contributions from circumstellar dust and intrinsic temperature variations. Given that the statistical uncertainty contribution to H 0 is dominated by SNIa magnitude errors, one can afford to make rather severe cuts in the Cepheid data in order to minimize the impact of systematic effects, including color cuts to try to homogenize the different color excess distributions between galaxies.
Here, we investigate the impact on H 0 of mitigating the color calibration impact by removing the reddest Cepheids, most susceptible to dust extinction. In Figure 3, we show the fitted H 0 as a function of the cut in E(V − I) we apply, when calibrating using color excesses for individual R E in the range [0, 1]. We note that the size of the error bars do not significantly increase un-lessÊ(V − I) max 0.1 mag, and that the inferred H 0 drifts towards the Planck value when cutting out redder Cepheids.
The reason for H 0 to decrease when removing red Cepheids is that even after applying the color calibration of equation 1, in individual galaxies there are still left-over correlations between the m W H andÊ(V − I). This will be the case also when allowing for individual galactic R E since the fitted value will not only depend on the slope of the magnitude-color diagram, but also on its normalization. After color calibration, m W H andÊ(V − I) tend to be anti-correlated in SNIa host galaxies and correlated in anchor galaxies, see Table 2. Removing Cepheids with largeÊ(V − I) therefore increases the mean m W H in SNIa host galaxies and decreases the mean m W H in anchor galaxies. As discussed in Mortsell et al. (2021), both of these effects decrease the inferred H 0 . The exception is MW Cepheids that have m W H andÊ(V − I) with slightly stronger anticorrelation than SNIa hosts. However, the fact that the MW Cepheid color excess distribution is heavily biased toward high values means that when applying a cut in E(V − I) 0.3 mag, too few MW Cepheid will remain to contribute appreciably to the inferred H 0 .

TEMPERATURE VARIATIONS
Given our inability to separate dust and Cepheid temperature color effects, we will here treat the temperature magnitude and color variations as an additional source of uncertainty, following the physical Cepheid model from Pejcha & Kochanek (2012). Denoting the mean temperature deviation for a specific Cepheid by τ , the absolute X-band magnitude will shift according to where M X is the mean absolute magnitude at the specific period, color and metallicity. β X for a large range of filters are observationally constrained in Pejcha & Kochanek (2012), where also temperature variations are estimated to σ τ = 0.02 given the width of the Cepheid instability strip. The intrinsic color (ic) excess will accordingly be given bŷ For each Cepheid, we randomly assign a temperature from τ = 0 ± σ τ and adjust the H-magnitude and color excess according to assuming β H = 1.72, β I = 3.39 and β V = 5.14. The induced impact on the inferred H 0 is estimated generating random Monte Carlo samples, see Section 8. Note that it makes sense treating the temperature color excess as noise for the dust color excess since the former is subdominant, withÊ ic (V − I) 0.1 mag.

METALLICITIES
We correct for possible metallicity effects on the luminosity using δm  (2015), a value of H/O = 1500 ± 300 is given, corresponding to Z ⊙ = 8.824 ± 0.087 which we will use as our default value, with an 0.15 dispersion to take into account the uncertainty between Cepheid metallicities estimated from oxygen and iron abundances.

Impact on H 0
The derived H 0 decreases slightly with decreasing the LMC metallicity. For a fixed value R E = 0.386, shifting ∆ log(O/H) = −0.3 → −0.4079 decreases the inferred Hubble constant with δH 0 ∼ 0.26.
The inferred H 0 depends on k as δH 0 /δk ∼ −0.7, and changing from k = 0 to k = −0.5, we obtain H 0 = 73.3 ± 1.3 for R E = 0.386. A constant systematic shift in the Cepheid metallicities as inferred from iron and oxygen abundances can be parametrized by changing the assumed solar oxygen abundance. An increase δZ ⊙ will increase H 0 with δH 0 /δZ ⊙ ∼ 5. Shifting Z ⊙ = 8.824 → 8.674, for R E = 0.386 we get H 0 = 71.9 ± 1.4 (3.1 σ tension). The full impact on the inferred H 0 when allowing also for varying galactic R E is estimated using Monte Carlo techniques in Section 8.

MILKY WAY PARALLAX UNCERTAINTIES
Trigonometric parallaxes potentially provide the most direct calibration of the Cepheid absolute magnitude, M W H . We use data from Riess et al. (2021), with 68 MW Cepheids having estimated GAIA parallaxes. As described in Lindegren et al. (2021), there are systematic errors in the published GAIA parallax values, with the parallax bias depending on, for example, the magnitude, color and angular position of the source. Although the GAIA team provide tentative expressions for the parallax correction, given that these corrections are uncertain for sources as bright as the MW Cepheids used in this study, here, as well as in Riess et al. (2021), we will allow for the possibility of residual parallax biases to correct for. Possible systematic GAIA Cepheid parallax uncertainties were investigated using globular clusters in Vasiliev & Baumgardt (2021) and Maíz Apellániz et al. (2021). Based on MW Cepheid parallax uncertainties and metallicity effects, in Owens et al. (2022) and references therein, it was also concluded that the uncertainty in the derived Cepheid distance scale may be underestimated, suggesting a systematic error floor of ∼ 3 %.
For the jth Cepheid in the MW, where the distance modulii for each Cepheid is estimated using GAIA parallaxes, π, according to π j + zp = 10 −0.2(µj−10) , where zp is a residual parallax calibration offset that we fit for together with M W H , b W and Z W by writing µ j = 10 − 5 ln 10 ln π + ln 1 + zp π = 10 − 5 ln 10 Higher order terms, O(zp/π) 2 , are small and corrected for in an iterative manner. So far, with the exception that we fit for zp simultanously with all other parameters, this is similar to the approach in Riess et al. (2021). For the default case where we calibrate an individual extinction law for each galaxy in the interval R E = [0, 1], we obtain zp = −18.9 ± 6.1 µas. For a fixed global value of R E = 0.386, zp = −16.9 ± 5.1 µas.
After correcting for a constant residual parallax offset in the GAIA data, we may have a residual that correlates with magnitude and color, see  π → π + zp 1 + zp 2 · m H + zp 3 · (V − I). (15) Using the MW as the sole anchor galaxy with a fixed R E = 0.386, adding zp 2 increases the Hubble constant from H 0 = 73.8 ± 1.5 to H 0 = 77.0 ± 1.9. Adding also zp 3 gives H 0 = 72.5 ± 1.6. The fact that a change in the parallax offset parameterization can induce a shift of δH 0 = 4.5 introduces doubt about the uncertainty estimates.
An interesting possibility for calibrating MW Cepheid distances circumventing problems connected to their bright nature and variability is to employ Cepheids for which the distance can be estimated from the parallax of spatially resolved companions or their host open cluster. In Breuval et al. (2020), a sample of 36 MW Cepheids is constructed in this way, and assuming a fixed value for the GAIA calibration off-set, it was noted that the inferred H 0 is decreased comparing to the case of parallaxes measured from the Cepheids themselves. Here, we consider the impact of using the MW Cepheid sample from Breuval et al. (2020) when performing the full simultaneous parameter distance ladder fit, including the GAIA parallax off-set and assuming R E = 0.386. The inferred Hubble constant, H 0 = 72.5 ± 1.4 (3.5 σ Planck tension), is slightly decreased compared to using the Cepheids in Riess et al. (2021), with slightly larger uncertainties given the smaller MW Cepheid sample size. The value is very similar to the result obtained when skipping MW Cepheids altogether in which case H 0 = 72.7 ± 1.4 (3.5 σ), before taking other systematic uncertainties into account.
When allowing for individually fitted galactic R E , we obtain H 0 = 71.3 ± 1.6 (2.3 σ) using the Breuval et al. (2020) MW Cepheid sample, again neglecting additional systematic effects. of SNIae are in significant tension with the CMB result for H 0 , we also note that the scatter is slightly larger than expected given the individual error bars, possibly indicating underestimated SNIa magnitude uncertainties. Taking the full correlation between the inferred M B into account, the SNIa error bars should be increased by a factor of 1.22 in order for the χ 2 /dof = 1 assuming the individual M B have their origin in a common value. Increasing the σ(M B ) by this factor, for the default case we obtain H 0 = 73.0±1.5, decreasing the Planck tension slightly to 3.6 σ. Fitting galactic individual R E = [0, 1] increases the H 0 scatter and the SNIa errors bars need to be increased by by a factor of 1.37. Fitting for a global value Hubble constant for default parameter values, gives H 0 = 73.9 ± 1.8 (3.5 σ tension). For a set of more restricted priors, with flat priors R E = [0.15, 0.8] and Gaussian priors R E = 0.48 ± 0.1 (see Mortsell et al. 2021), the corresponding SNIa error factors are 1.17 and 1.23 respectively.

Impact of Type Ia Supernova Data Set
The fact that the statistical error on H 0 is dominated by SNIa data, and only a sub-set of the SNIa are in tension with the Planck H 0 may raise concerns about the reliability of the individual SNIa magnitude estimates. However, using the SNIa data set from Burns et al. (2018), gives very similar result, indicating that systematic effects related to the SNIa light curve fitting, stretch and color correction are subdominant, see figure 6 (assuming a fixed R E = 0.386). The possibility that Cepheid calibrated SNIae to a larger extent originate in star-forming environments, thus being dimmer compared to the average Hubble flow SNIae, may bias the Hubble constant measurements has been discussed in Rigault et al. (2015).  Figure 6. Comparing the distance modulii for SNIa host galaxies as derived using SNIa data from Burns et al. (2018) (B18) and Riess et al. (2016) Table 3. The Hubble constant obtained in this case is H 0 = 73.1 ± 1.8 (3.0 σ Planck tension). Table 3. Approximate individual contributions to the systematic error budget for a fixed RE = 0.386, with Int color being Cepheid intrinsic color uncertainties as given in Tammann et al. (2011), τ the intrinsic Cepheid temperature for which στ = 0.02 (Pejcha & Kochanek 2012), k and Z⊙ relating metallicities as inferred from iron and oxygen abundances as described in Section 5, and SNIa error being the effect of re-scaling SNIa magnitude uncertainties (see Section 7). Allowing for R E to vary between galaxies in the range [0, 1], we get H 0 = 73.9 ± 2.2 (2.9 σ) with a GAIA residual parallax calibration offset of zp = −19.3 ± 7.3 µas. Assuming flat priors R E = [0.15, 0.8] we obtain H 0 = 73.7±2.1 (3.0 σ) and for Gaussian priors R E = 0.48±0.1 H 0 = 73.5 ± 2.1 (2.9 σ).
We thus note that mitigating the color calibration impact by removing Cepheids for which dust extinction is expected to dominate the observed color excess, the H 0 inferred from supernovae agree with the Planck value re-gardless of Cepheid color calibration parameterization. Cutting out the blue Cepheid tail has to the opposite, although less pronounced, effect of increasing the inferred H 0 . RequiringÊ(V − I) > 0 mag, we obtain H 0 = 74.8 ± 2.3 for R E = [0, 1].

SUMMARY
In Mortsell et al. (2021), we investigated the sensitivity of the Hubble constant as inferred from SNIa distances to the choice of Cepheid color calibration method. Here, we complement this analysis by investigating the impact of Cepheid temperature variations, metallicity corrections, supernova magnitude uncertainties and Milky Way parallax systematics. The results are summarized in Figure 7 and Table 4. Only using MW Cepheids for which distances can be estimated from companions and host cluster parallaxes (Breuval et al. 2020) and flat priors R E = [0, 1], we obtain H 0 = 71.9 ± 2.2, decreasing the difference with the Planck value to 2.0 σ, and in good agreement with H 0 = 69.6 ± 1.6 obtained calibrating the absolute SNIa magnitude using the tip of the red giant branch observations (Freedman et al. 2019). Using more restricted prior values for R E yield similar results, see Table 4.
A possible caveat with ours, and all other H 0 measurements using Cepheid calibrated distances, comes from the fact that the observed Cepheid color excess distribution have large variations across galaxies, and that the inferred Hubble constant is sensitive to color excess cuts in the data. Also applying a color excess cut to Comparing results for H0. The solid black line is for individual galactic values of RE = [0, 1]. For the dotted black line, we have imposed an upper limit on the allowed estimated color excess, only including Cepheids for whichÊ (V − I) < 0.15 mag. For the dashed black line, we use the MW Cepheid sample from Breuval et al. (2020). The solid petrol line is fitted using the Wesenheit calibration with RW = 0.386 as in Riess et al. (2021) and the dashed brown region indicates the 1 σ region from Planck (Aghanim et al. 2020).
remove Cepheids for which dust extinction is expected to dominate the observed color excess,Ê(V − I) > 0.15 mag, we obtain H 0 = 68.1 ± 2.6.
We thus conclude that the current Hubble tension may be affected by systematic effects in the calibration of the local distance ladder, including MW Cepheid distances, as well as Cepheid selection biases.
We thank Vallery Stanishev for interesting discussions regarding stellar metallicities and the anonymous referee for useful comments on the manuscript. EM acknowledges support from the Swedish Research Council under Dnr VR 2020-03384. AG acknowledges support from the Swedish Research Council under Dnr VR 2020-03444, and the Swedish National Space Board, grant 110-18.