Sensitivity of the Hubble Constant Determination to Cepheid Calibration

Motivated by the large observed diversity in the properties of extra-galactic extinction by dust, we re-analyse the Cepheid calibration used to infer the Hubble constant, $H_0$, from Type Ia supernovae, using Cepheid data in 19 Type Ia supernova host galaxies from Riess et al (2016) and anchor data from Riess et al (2016, 2019, 2021). Unlike the SH0ES team, we do not enforce a fixed universal color-luminosity relation to correct the Cepheid magnitudes. Instead, we focus on a data driven method, where the optical colors and near infrared magnitudes of the Cepheids are used to derive individual color-luminosity relations for each Type Ia supernova host and anchor galaxy. We present two different analyses, one based on Wesenheit magnitudes resulting in $H_0=73.2\pm 1.3$ km/s/Mpc, a $4.2\,\sigma$ tension with the value inferred from the cosmic microwave background. In the second approach, we calibrate an individual extinction law for each galaxy with non-informative priors using color excesses, yielding $H_0=73.9\pm 1.8$ km/s/Mpc, in $3.4\,\sigma$ tension with the Planck value. Although the two methods yield similar results, in the latter approach the Hubble constant inferred from the individual Cepheid absolute distance calibrator galaxies range from $H_0=68.1\pm 3.5$ km/s/Mpc to $H_0=76.7\pm 2.0$ km/s/Mpc. Taking the correlated nature of $H_0$ inferred from individual anchors into account and allowing for individual extinction laws, the Milky Way anchor is in $2.1\,\sigma - 3.1\,\sigma$ tension with the NGC 4258 and the Large Magellanic Cloud anchors, depending on prior assumptions regarding the color-luminosity relations and the method used for quantifying the tension.


INTRODUCTION
As is well known, there is a tension between the value of the Hubble constant as inferred from small and large distance measurements, most significantly between the values inferred from Type Ia supernova (SNIa) distances to redshifts z ∼ 0.1 calibrated by Cepheid observations, as measured by the SH0ES team, and the distance to the cosmic microwave background (CMB) decoupling surface at z ∼ 1090, as measured by the Planck satellite. The former yields H 0 = 73.2±1.3 (in units of km/s/Mpc used from now on) (Riess et al. 2021a) and the latter H 0 = 67.4 ± 0.5 (Aghanim et al. 2020); a 4.1 σ tension 1 The tension between other measurements is not as significant: Calibrating the absolute SNIa magnitude using the tip of the red giant branch (TRGB) observations gives H 0 = 69.6 ± 1.6 (Freedman et al. 2019), right between the Cepheid calibrated SNIa and the CMB inferred values.
Another local estimate of the expansion rate of the Universe is derived from the amplitude of the gravitational wave signal GW170817, the merger of a binary neutron-star system located to the galaxy NGC 4993 at z = 0.010 through the electromagnetic counterpart, AT2017gfo, yielding a Hubble constant with relatively large uncertainties of H 0 = 70.0 +12.0 −8.0 (Abbott et al. 2017).
A different route involves gravitational lensing. Time delays in the TDCOSMO sample of seven lensed quasars yield H 0 = 74.5 +5.6 −6.1 (Birrer et al. 2020). Constraining the galaxy lens mass profiles, using kinematics observations of an independent set of gravitational lenses in the Sloan Lens ACS sample (SLACS), lowers the value to H 0 = 67.4 +4.1 −3.2 , assuming that the TDCOSMO and SLACS galaxies are drawn from the same parent population. These results are illustrated in Figure 1, from which is evident that only the Cepheid calibrated SNIa distance scale from SH0ES is in definite tension with the CMB inferred distance. . The Hubble constant as inferred from distances measured to the CMB, strongly lensed quasars, SNIa calibrated with Cepheids (SH0ES) and the TRGB, and the gravitational wave signal GW170817. The redshift corresponds to the mean redshift of the distances employed in the method, slightly shifted where needed to avoid overlap. The major tension is between the values inferred from Cepheid calibrated SNIa by the SH0ES team and CMB observations.
The value inferred from the CMB depends on the entire expansion history of the Universe, whereas the SNIa measurement only depends on the local expansion rate it sets out to measure. On the other hand, the inference from SNe Ia depends on a combination of a larger number of astrophysical probes. Therefore, attempts to modify the CMB inferred H 0 usually rely on modifications of the cosmological model, whereas the SNIa value is usually studied with emphasis on possible systematic effects concerning the local distance measurements.
At least in principle, the CMB value can be increased by various departures from the concordance cosmological constant and cold dark matter, ΛCDM model, see e.g. Mörtsell & Dhawan (2018); Knox & Millea (2020). Options include decreasing the physical size of the sound horizon used to measure the distance to z = 1090. This can be accomplished by adding sources of energy present before CMB photon decoupling, e.g., new thermal relativistic species or early dark energy, or by reducing the sound speed. However, such modifications are severely constrained when taking the full CMB power spectrum into account. Attempts to shift the CMB value also involve changing the expansion history at redshifts z < 1090, with modest success since the expansion rate is tightly constrained by SNIa and baryonic acoustic oscillation (BAO) observations. Given that the proposals mentioned above require substantial modifications of the current concordance cosmological model, and still fail in relieving the full tension, we investigate the Cepheid-SNIa value and its uncertainties, see also Follin & Knox (2018);Efstathiou (2020). In particular, we concentrate on dust extinction, affecting all astronomical observations in the optical and near infrared (NIR) regime. We focus on a a very specific assumption made by the SH0ES team throughout their series of publications, namely that there is a fixed universal reddening law in all galaxies.

REVISITING DUST EXTINCTION CORRECTIONS
In spite of the critical importance for precision cosmology, the current understanding of light attenuation in the interstellar medium (ISM) of galaxies remains very limited. In comparison, the Milky Way (MW) ISM has been studied in great detail, including the properties of dust grains responsible for dimming of light(see Draine 2003, for a review). In particular, several MW reddening laws have been devised, among these Cardelli et al. (1989) (CCM), O'Donnell (1994 and Fitzpatrick (1999) (F99). They have in common the use of a single parameter, the total to selective extinction coefficient, R BV V , as a proxy for the grain composition and size distribution, where the attenuation in the optical V-band relates to the color excess , here referred to as the "CCMrelationship". Low values of R BV V indicate a steep wavelength dependence, while large R BV V correspond to gray extinction.
While an average R BV V = 3.1 for the MW is found in most studies, significant variations are found in individual lines of sight in the galaxy, ranging from R BV V ≈ 2 in some diffuse sight lines, to R BV V ≈ 6 in dense molecular clouds (Fitzpatrick 1999). The diversity in the MW has been confirmed by Nataf et al. (2016), who find significantly lower values of R BV V in the Galactic bulge. Moving the scope outside the MW, a study by Gordon et al. (2003) of the extinction in the Magellanic Clouds found that a small number of Large Magellanic Cloud (LMC) extinction curves are consistent with the CCM relationship, but the majority of the LMC and all the SMC curves are not. Fausnaugh et al. (2015) report a gray extinction law for NGC 4258 in the line of sight of the Cepheids, R BV V ∼ 4.9, although they caution that this could be the result of unresolved systematics. For more distant galaxies, observed SNIa colors highlight the observed diversity in extinction properties, ranging from R BV V ∼ 1 to values consistent with the MW average (see Krisciunas et al. 2006;Nobili & Goobar 2008;Goobar et al. 2014;Amanullah et al. 2014Amanullah et al. , 2015Burns et al. 2018, and references therein). For the SNe Ia in the Hubble flow, color corrections are based on the SALT2 lightcurve fitter (Guy et al. 2007), which again differ from the CCM parameterization, but are most consistent with values of R BV V ∼ 2.5 (see e.g., Biswas et al. 2021, and references therein), although dust extinction differences between host galaxy environments has been suggested as an explanation for a systematic "mass step" in the derived distances (Brout & Scolnic 2021;Johansson et al. 2021).
To minimize the impact from extinction correction uncertainties, the SH0ES team use flux measurements in the NIR H-band, centered at 1.6 µm, where extinction by dust, gauged using the observed color V − I, is significantly smaller. Adopting the CCM-like relationship from F99 and the extinction correction A H = R VI H E(V − I), the value corresponding to the MW average is R VI H ∼ 0.4. However, there is no theoretical, nor any empirical studies of extinction suggesting that a universal value of R VI H can be assumed. On the contrary, a recent study by Fitzpatrick et al. (2019) finds considerable variations between lines of sight for the extinction curves in the NIR in the MW. Based on a parameterization fitting extinction laws from ultraviolet to NIR for 72 well-measured stars, a very wide range in R VI H can be inferred, as shown in Figure 2. Hence, assuming a narrow range in R VI H for the anchor and Cepheid hosts in not warranted by current observations. In this paper, we investigate to what degree relaxing this assumption affects the inferred value of H 0 and its corresponding uncertainties.  (Gordon et al. 2003, LMC1 and LMC2) and in the water mega-maser anchor NGC 4258 (Fausnaugh et al. 2015, FA15).

DIFFERENCE IN METHODOLOGY
The local distance ladder uses the difference between Cepheid magnitudes in SNIa hosts and anchor galaxies. Therefore, it is relatively insensitive to changes in the global properties of the extinction law. In Riess et al. (2016), a global change of the parameter R W parameterizing a Cepheid color-luminosity (C-L) correction with respect to the observed color R W (V−I) of R W = 0.39 → 0.35 changed the Hubble constant by δH 0 ∼ 0.5. When multiplied with the observed color V − I, the parameter R W is argued to also partly correct for an intrinsic Cepheid C-L relation.
In Follin & Knox (2018), a slightly different C-L correction with respect to an estimated color excess of R EÊ (V − I) was employed, whereÊ (V − I) ≡ (V − I) − V − I 0 with V − I 0 an estimate of the mean intrinsic Cepheid color. Here, R E is interpreted as the dust total to selective extinction ratio, R VI H . Imposing a prior of R E = 0.39 ± 0.1, a value of H 0 = 73.3 ± 1.7 was derived when allowing R E to vary between galaxies, in good agreement with H 0 = 73.2 ± 1.3 from Riess et al. (2021a). At face value, this result seems to suggest that the method used for calibrating the C-L relation and/or the possibility of varying this calibration between galaxies have a small impact on the inferred Hubble constant.
Since dust extinction is uncertain at NIR wavelengths, we investigate the effect of allowing for R W and R E to be fitted by the Cepheid data. Since dust properties are also known to vary between galaxies, when calibrating with respect to color excesses, we allow R E to do the same. Given the lack of solid independent constraints on R E , we employ less restrictive priors than Follin & Knox (2018). Any systematic difference in dust properties or the intrinsic C-L between SNIa hosts and anchor galaxies will shift the inferred value of H 0 .

METHOD AND DATA
The apparent magnitude of a source at redshift z with absolute magnitude M is given by where D(z) is the luminosity distance in units of Mpc. By combining observed magnitudes of Cepheids in SNIa host and anchor galaxies, m host Ceph and m anch Ceph with SNIa magnitudes in host galaxies and in the Hubble flow, m host SN and m flow SN , we can derive where r(z) ≡ H 0 D(z) can be approximated by r(z) ≈ cz in the close Hubble flow and we have defined Apart from getting the SNIa redshifts and anchor distances right, we thus need to make sure there are no systematic offsets in the Cepheid and SNIa magnitudes between host, anchor and cosmic flow galaxies. Ignoring the weak cosmology dependence of r(z) (Dhawan et al. 2020), the inferred value of H 0 will decrease (increase) if we: 1. Increase (decrease) the independent anchor distances, D anch .
Here, we focus on option 3. With regards to option 2, ∆m SN will increase if SNIa in Cepheid host galaxies are systematically made brighter than in the Hubble flow, e.g., if there is additional dust extinction not accounted for in the host galaxies, or if the effect that SNIa in high mass hosts are systematically brighter than in low mass galaxies, such as Cepheid hosts, have been underestimated (see Rigault et al. 2020, and references therein). In terms of option 3, if there is additional dust extinction not accounted for in the anchor galaxies, or if we have over-corrected for dust extinction in the host galaxies, the inferred value of H 0 will decrease, and vice versa. The fractional shift in the Hubble constant is (5) and a lower limit to the precision in H 0 is set by the precision of the anchor distance measurements. Shifting δ(∆m Ceph ) = ±0.1 will shift the Hubble constant by δH 0 ∓ 4.6 %.

Cepheid calibration
We use the Hubble Space Telescope (HST) flux in the NIR filter (H = F160W) band, color calibrated using optical (V = F555W and I = F814W) data, to derive Wesenheit magnitudes In the last step, we see that m W H is corrected both for dust extinction, identifying the first R W with the total to selective extinction ratio R VI and for a possible intrinsic C-L relation, identifying the second R W with β VI H , as parameterized in e.g. Madore (1982) and Madore et al. (2017). The term β VI H (V − I) 0 corresponds to the intrinsic magnitude-color relation at a fixed Cepheid period, whereas the correlation between the intrinsic color with period is included in the periodluminosity (P-L) calibration parameterized by b W in equation 9 below. As an alternative approach, also employed in Follin & Knox (2018), we calibrate the C-L relation using Here,Ê (V − I) represents a proxy for the color excess obtained by subtracting an estimate of the mean intrinsic colors, V − I 0 from the observed colors 2 , The estimated color excessÊ (V − I) also represents a combination of dust extinction and intrinsic color, since the mean intrinsic colors, V − I 0 does not take into account individual variations in Cepheid temperature along the width of the Cepheid instability strip, see e.g. Madore & Freedman (1991); Sandage & Tammann (2006); Pejcha & Kochanek (2012). Using multi wavelength data, one can in principle attempt to distinguish the contribution from dust and intrinsic color variations, see e.g. Pejcha & Kochanek (2012); Madore et al. (2017) and calibrate them separately. In the following, we will follow standard practice and assume that R W and R E effectively corrects both for dust and intrinsic color variations, noting that the former will dominate since temperature variations are subdominant given the narrow width of the instability strip. Taking an empirical approach, for the calibration with respect to the observed color we will fit for the value of R W that minimizes the scatter in m W H . When calibrating with respect to color excesses, we will allow for individual galactic R E , representing a variation in the dust properties between galaxies.
We model the Wesenheit magnitude of the jth Cepheid in the ith SNIa host as where [M/H] i,j is a measure of the metallicity of the Cepheid, [P] i,j ≡ log P i,j − 1 where P i,j is the period measured in days, 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. In what follows, we will allow for separate P-L relations for short and long period Cepheids using where [P] s i,j = 0 for Cepheids with periods > 10 days and [P] l i,j = 0 for Cepheids with periods < 10 days, see Section A.
Similarly for the jth Cepheid in the kth anchor galaxy, here MW, NGC 4258 and the LMC,

Milky Way Cepheids
Trigonometric parallaxes potentially provide the most direct calibration of the Cepheid absolute magnitude, M W H . We use data from Riess et al. (2021a), with 75 MW Cepheids, out of which 68 have reliable GAIA parallaxes. For the jth Cepheid in the MW, where the distance modulii for each Cepheid is estimated using GAIA parallaxes, π, according to where zp is a residual parallax calibration offset that we fit for together with M W H , b W and Z W . In Riess et al. (2021a), M W H and zp are fit for using only MW data setting b W = −3.26 and Z W = −0.17 (as fitted to all Cepheids), finding zp = −14 ± 6 µas. Since we want to fit for zp simultaneously with all parameter, we write µ j = 10 − 5 ln 10 ln π + ln 1 + zp π = 10 − 5 ln 10 effectively transforming zp into a linear parameter, and Higher order terms, O(zp/π) 2 , are small and corrected for in an iterative manner.

Type Ia Supernovae
The calibrated SNIa B-band peak magnitude in the ith host is modelled by The SNIa peak apparent magnitudes need to be corrected for the width-luminosity and C-L relations.
There are several lightcurve fitting algorithms for deriving the SNIa peak magnitude, lightcurve shape and color, the most widely used for cosmology being the SALT2 model (Guy et al. 2007). The derived lightcurve widths and colors are used to correct the peak magnitude m B in equation 16. The errors on the corrected peak magnitude include the fitting error and a 0.1 mag term from the SNIa model added in quadrature.

Data
For the extra-galactic (M31 and beyond) Cepheids, including Cepheids in the anchor galaxy NGC 4258, we use the data set from Table 4 in Riess et al. (2016). This table is restricted to Cepheids passing a best-fit, global 2.7 σ outlier rejection, the impact of which is claimed to be small in Riess et al. (2016), but not possible to confirm independently by us.
For Cepheids in the LMC, we use data in Table 2 in . Data for MW Cepheids, including GAIA parallax measurements are from Table 1 in Riess et al. (2021a).
Double eclipsing binaries (DEBs) provide a means to measure distances by determining the physical sizes of the member stars via their radial velocities and light curves (Paczynski 1996). 20 DEBs observed using longbaseline near-infrared interferometry give a distance to the LMC of µ LMC = 18.477 ± 0.0263 (Pietrzyński et al. 2019;. We use the updated distance to NGC 4258 of µ N4258 = 29.397 ± 0.032 (Reid et al. 2019), using observations of mega-masers in Keplerian motion around its central super massive black hole.
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).

Parameter Fitting
Given the observed Cepheid magnitudes m H , colors V − I, periods [P], metallicities [M/H], together with the SNIa magnitudes m B , the anchor distances µ k and the MW Cepheid parallaxes π, we can fit simultaneously for R W or 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 as described in Section A. Although R W and R E appear as linear parameters in equation 6 and 7, since the uncertainty in the observed Cepheid colors are non-negligible, the uncertainty in the derived Wesenheit magnitudes m W H will depend on the values of R W or R E , requiring a non-linear treatment of these parameters, see Section A. Effectively, a proper color uncertainty treatment will assign a larger Wesenheit magnitude uncertainty for larger values of R W and R E , whereas not taking the color correction uncertainty into account will bias the result towards too low values of R W and R E . Whether this will bias the inferred Hubble constant low or high depends on the relative bias between R E in SNIa host and Cepheid distance calibrator (anchor) galaxies. For illustrative purposes, we will present an example of the magnitude of this bias in the Results sections.
Given M B , the Hubble constant is calculated as where a B is the intercept of the SNIa magnitude-redshift relation Using the same data with R W = 0.386, using double P-L relations, doing a full count-rate non-linearity correction following , identifying [O/H] = [Fe/H] using Z ⊙ = 8.824, and fitting for the residual GAIA offset zp simultaneously with all other parameters, we obtain H 0 = 73.1 ± 1.3. Here, we have added a scatter in the P-L relation of σ(m W H ) = 0.0682, to give a χ 2 /dof = 1. Despite slight differences in the analysis method, this value is in very good agreement with the value in Riess et al. (2021a), and in 4.1 σ tension with the Planck value of H 0 = 67.4 ± 0.5 (Aghanim et al. 2020).

Fitting for R W
The fact that dust extinction extrapolated to the Hband is very uncertain naturally opens up for the option of fitting for the value of R W . Assuming a global value common for all galaxies gives R W = 0.37 ± 0.02 with H 0 = 73.2 ± 1.3 (4.2 σ tension). Here, we have used a flat prior of R W = [0, 1], but given the small posterior uncertainty on R W , results are insensitive to the specific choice of this prior.

Individual P-L relations
So far, we have assumed that all Cepheids can be described by a global P-L relation, described by b s W and b l W . However, since there is evidence that the P-L relation can vary between galaxies (Tammann et al. 2011;Efstathiou 2020), in a similar spirit to our approach of allowing R W to vary between galaxies, we investigate to what extent relaxing this assumption will affect the inferred Hubble constant. We will allow for individual galactic values of b W,i to be fitted for simultaneously with all other parameters, in this case restricting to the same P-L relation for short and long period Cepheids. For a fixed global value of R W = 0.386, the resulting b W,i are shown in Figure 3. The fact that the fitted b W,i are systematically higher in SNIa host galaxies compared to hosts, the inferred Hubble constant is increased from H 0 = 73.1 ± 1.3 to H 0 = 76.5 ± 1.9. If we allow for both individual P-L and a globally fitted C-L relation, the Hubble constant is again H 0 = 76.5 ± 1.9, with Planck tension 4.7 σ. The inferred value of H 0 will shift if there is a systematic offset in R E between anchor and SNIa host galaxies. We derive values for H 0 when varying the value of R E in the anchor(s) and the SNIa hosts, see Figure 4. For a common R E (indicated by the dotted line), the inferred value of H 0 decreases when R E is increased. Also, H 0 is decreased when R E is larger in anchor than in host galaxies. Given the result in Figure 4, one could argue that a simple solution to the Hubble tension, would be a systematic shift of R E between anchor and host galaxies, or an increased overall value of R E bringing H 0 as inferred from SNIa down to the Planck value. This could be argued for, e.g., if dust properties in SNIa host galaxies have a systematically steeper extinction law than the anchor galaxies. In lack of solid independent evidence for such a systematic shift, or the value of R E and possible variations of it, we will next include R E as model parameters to be constrained by the available data. Fitting for a global value of R E (again using a flat prior R E = [0.1], common for all galaxies, gives R E = 0.37 ± 0.02 with H 0 = 73.2 ± 1.7 with Planck tension 3.3 σ. Allowing also for individual galactic P-L relations, we obtain H 0 = 76.0 ± 2.2 (3.9 σ tension).
We next allow for R E to vary between galaxies. Since R E primarily represent corrections to dust, with large uncertainties in the NIR, we use wide priors of R E = [0, 1] as a default to investigate the full impact of using the data at hand to constrain R E . We obtain H 0 = 73.9 ± 1.8, in 3.4 σ tension with the Planck value, see Figure 5. If not properly taking into account the dependence of the Wesenheit magnitude uncertainties on the R E , as discussed in Section 4.5, one obtains H 0 = 72.2±1.7 (2.7 σ tension), i.e., a bias of ∆H 0 = 1.7.
Imposing a set of more restrictive flat priors R E = [0.15, 0.8], the inferred Hubble constant is H 0 = 73.6 ± 1.8. Our most restrictive set of priors is derived by assuming that the (mean) reddening of extra-galactic Cepheids is well-represented by the reddening distribution of MW stars. Using the sample of R BV At face value, these results suggest that the inferred Hubble constant is quite insensitive to the Cepheid color calibration scheme. In terms of the results for individual anchors however, when allowing for individually fitted R E = [0, 1], the Hubble constants inferred from each individual anchor distance show a substantial spread with the MW anchor distance being dominant in pushing the combined Hubble constant to the high value in tension with Planck measurements, see Section B. Since the different H 0 derived for each individual anchor distance are correlated, the significance of the tension can not be immediately read out from the H 0 and their corresponding uncertainties.
Instead, we quantify the significance of differences between anchors in the following ways. First, using only MW Cepheids, we derive the Cepheid absolute magnitude, M W H . We next compare each of the independently measured µ LMC and µ N4258 to the values derived using the other two anchors. For example, using NGC 4258 and LMC as anchors but no MW Cepheids, we compare the derived value of M W H to that previously obtained using only MW Cepheids. Also, using the MW Cepheids value for M W H and the LMC distance as anchors, we compare the derived µ N4258 to the independently derived value. In order to account for the "look elsewhere" effect, we make use of a Monte Carlo approach where we simulate three independent data points from the same mean value, and arbitrary variances. Using a large number of realisations, we compare each of the three points to the weighted mean of the other two points. We compute the fraction of realizations that has a certain maximum deviation for the three data points to quantify how a given value for the observed significance corresponds to a slightly lower significance of the results. For the case of R W = 0.386, NGC 4258 shows the largest deviation of 2.0 σ, corresponding to a modest 1.6 σ significance accounting for the "look elsewhere" effect. For individually fitted R E = [0, 1], the MW anchor has a deviation of 3.0 σ, corresponding to a 2.6 σ significance.
A drawback of this method when quantifying the tension of the MW anchor is that using only the 67 MW Cepheids, the uncertainty of the inferred M W H is fairly large (we typically obtain M W H = −5.88 ± 0.04 for the Wesenheit calibration and M W H = −5.57 ± 0.05 when color calibrating with respect to estimated intrinsic colors). As an alternative approach, using only the MW Cepheids as distance anchors, we infer the distances to NGC 4258 and the LMC and compare to their independently measured values. Despite the slight differences between the two approaches, they yield broadly consistent results for the MW anchor tension. Results are summarized in Figure 11 and Table 1 for different priors on the color-luminosity relation. For the color excess calibration, the MW anchor is in 2.1 σ − 3.1 σ tension with the NGC 4258 and the LMC anchors, depending on the assumed prior distributions of R E . Note again that we are here using data from Riess et al. (2016Riess et al. ( , 2021a. In Riess et al. (2021b), an expanded data set is used to show good consistency between the individual anchor distances for a fixed value of R W = 0.386, attributed to a refined Cepheid metallicity dependence decreasing the inferred correction term from Z W = −0.13 ± 0.07 (Riess et al. 2016) to Z W = −0.22 ± 0.05.

SUMMARY AND DISCUSSION
We have investigated the sensitivity of the Hubble constant inferred from SNIa distance measurements to the choice of Cepheid calibration method. Specifically, we have compared results when color calibrating the Cepheid magnitudes with respect to observed colors and estimated color excesses. Guided by the lack of independent evidence for dust extinction properties at NIR wavelengths, we have allowed for the color calibration parameters to be determined by the Cepheid data; either a global value or in the case of the color excess calibration, individual galactic values.
For the color excess calibration, we derive H 0 = 73.9 ± 1.8 in 3.4 σ tension with the value inferred from Planck, using priors of R E = [0, 1]. Using more restrictive priors on R E only has small effects on the inferred Hubble constant. Calibrating with respect to observed 70 75 H 0 Figure 6. Results for H0 for the two main analyses employed in this paper. The solid black line is for a fitted global value of RW using Wesenheit magnitudes, overlapping with the dashed petrol line fitted using the Wesenheit calibration with RW = 0.386 as in Riess et al. (2021a). The dark petrol line is for individual RE using color excesses, assuming prior values RE = [0, 1]. The dashed brown region indicates the 1 σ region from Planck (Aghanim et al. 2020). colors yields H 0 = 73.2 ± 1.3, see Figure 6. Allowing also for individual galactic P-L relations, the corresponding H 0 values are increased to H 0 = 75.2 ± 2.4 and H 0 = 76.5 ± 1.9, respectively.
Results for different calibration choices are summarized in Table 2 in Section C. From the quality of the fits, there is no clear preference for any of the calibration methods with different information criteria showing preference for different amount of freedom in the model.
Regardless of calibration method, since Cepheid colors and periods are correlated (Tammann et al. 2011), so will the inferred P-L and C-L relations, i.e., the parameters b W and R W or R E . This may be of importance if there are color selection effects related to the fact that longer period Cepheids are brighter, see Figure 7. We have tested the possible impact of such an effect by imposing cuts on the observed color V − I. In Figure 8, we show the fitted H 0 as a function of the cut in (V − I) we apply, when calibrating using color excesses for a fixed R E = 0.386. With (V − I) max = [2, 1.5, 1.25, 1], [98 %, 85 %, 65 %, 29 %] of the original Cepheids remain. As evident from Figure 8, only for (V − I) max ≈ 1 does the cut significantly degrade the statistical uncertainty in H 0 . The trend of obtaining a lower H 0 when cutting out redder Cepheid is a common feature for all calibration methods and color-luminosity priors we have tested.
Finally, we note that when estimating distances using the I-band for which color corrections are larger, the sen-  With the limited information at hand regarding dust extinction for Cepheids at NIR wavelengths, the color calibration of Cepheid magnitudes could potentially introduce large uncertainties in the local distance ladder.
Allowing for a global R W or individually fitted values of R E in the anchor(s) and the SNIa hosts does not significantly change the inferred value of H 0 , although the H 0 as derived from the individual anchors will shift. In the case of individually fitted values of R E , the inferred values range from H 0 = 68.1 ± 3.5 for the NGC 4258 anchor to H 0 = 76.7 ± 2.0 for the Milky Way. Neither approach employed in this paper is in one-to-one correspondence with an underlying physical model and there is no clear evidence in the data for any of them.
We thank the anonymous referee for the insightful and thorough reviews of the manuscript, in particular for pointing out the importance of a proper treatment of color errors having substantial impact on the results when fitting for the values of R W and R E and the corresponding inferred values of H 0 . 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.  Riess et al. (2016), we collect all data points and their corresponding uncertainties and possible correlations in the matrices Y and C. This includes the Wesenheit magnitudes of all Cepheids, m W H,i,j , including the anchors. The exception is the MW Cepheids for which we use m π,j ≡ m W H,j − 10 + 5 ln 10 ln π. Next, we have the measured anchor distances, µ N4258 , µ LMC and possibly µ M31 . Finally, data points include the B-band SNIa magnitudes in the Cepheid hosts, m B,i : Collecting the model parameters in the matrix X we can relate data and parameters through Y = AX where in schematic form Note that the first N rows correspond to equation 15 and so forth, so that A is a n param × n data matrix. We can solve for the parameter matrix X and its covariance matrix Σ analytically )  Table 1. The tension between the three different anchors. For the first three columns, the independently measured distance to an anchor is compared to the value derived using the Cepheid and SNIa distance ladder calibrated using the other two anchors.
In the case of the MW, the distance is represented by the Cepheid absolute magnitude, M W H derived using MW Cepheids only. In the last column (Milky Way 2), we employ an alternative approach for the Milky anchor where we compare the distances to NGC 4258 and LMC derived using the MW anchor, calibrated using the full Cepheid sample, to their independently measured values. For the anchor with the largest deviation from the other two anchors, we also (in parenthesis) give the significance of deviation corrected for the "look elsewhere" effect. We use the emcee Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) (Foreman-Mackey et al. 2013) to obtain constraints on R W and R E . At each step, we recompute the Wesenheit magnitude uncertainties for the specific values of R W and R E , and use the formalism described above to solve for the linear parameters in the fit. The exceptions are for fixed values of R W , as well as cases where we investigate the bias induced from not correctly taking the color uncertainties into account, in which cases we can extend the linear formalism to include also R W and R E . We have checked that in these cases, the linear treatment and the MCMC formalism give consistent results.

B. RESULTS FOR INDIVIDUAL ANCHORS
In Figure 10, we show the inferred Hubble constant from each individual anchor. Here, we have used a globally fitted R W for the observed color calibration and individual R E for each galaxy for the color excess calibration. A global P-L relation is assumed in both cases. Dotted lines are for a global value R W = 0.386 using the Wesenheit calibration, closest resembling the case in Riess et al. (2021a). The solid lines are for the case of fitting for a global R W and individual R E .  Figure 11. Upper row: Individual anchor distances derived using the Cepheid and SNIa distance ladder using the other two anchors compared with their independently measured distances in Reid et al. (2019), Pietrzyński et al. (2019) and using MW Cepheids only to derive M W H (solid black lines, except M W H derived with the Wesenheit calibration using observed colors that is depicted with the solid brown line). Lower row: The distances to the NGC 4258 and the LMC anchors as inferred using MW data to calibrate the Cepheid magnitudes, compared to their independently measured distances (solid black lines).
Wesenheit magnitude errors are slightly decreased because of the slightly shifted values of R W and R E 4 . In terms of model selection, results are ambiguous, see Table 2. Here, the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) are defined by AIC ≡ 2n par + χ 2 , BIC ≡ n par ln n dat + χ 2 , with the latter more penalizing for models with extra parameters for n dat > 8. In terms of the AIC, fitting for individual values of R E with or without individual b W constitute the preferred models, whereas in terms of BIC, fixed values of R W = R E = 0.386 together with a global P-L relation yields the best result. A more careful analysis, beyond the scope of this paper, would entail calculating the Bayesian evidence by the integrating the prior times the likelihood over the entire parameter space of the model. In the last column of Table 2, we also account for the p-value of the fit. As it stands, it is unclear to what degree these results provide guidance on the choice of calibration model.