Acoustic Dark Energy: Potential Conversion of the Hubble Tension

TL;DR

This paper proposes a scalar field dark fluid model that converts potential to kinetic energy around matter-radiation equality, alleviating the Hubble tension by improving fit to CMB data and increasing the Hubble constant.

Contribution

It introduces a novel scalar field mechanism that balances acoustic oscillations in dark energy, providing a new way to address the Hubble tension within the mbda CDM framework.

Findings

Model improves fit to combined data by elta ext{ }.7 over mbda CDM.

Scalar field conversion mechanism is theoretically derived and constrained.

Future CMB polarization data can test this dark energy scenario.

Abstract

We discuss the ability of a dark fluid becoming relevant around the time of matter radiation equality to significantly relieve the tension between local measurements of the Hubble constant and CMB inference, within the $\Lambda$CDM model. We show that the gravitational impact of acoustic oscillations in the dark fluid balance the effects on the CMB and result in an improved fit to CMB measurements themselves while simultaneously raising the Hubble constant. The required balance favors a model where the fluid is a scalar field that converts its potential to kinetic energy around matter radiation equality which then quickly redshifts away. We derive the requirements on the potential for this conversion mechanism and find that a simple canonical scalar with two free parameters for its local slope and amplitude robustly improves the fit to the combined data by $\Delta\chi^2 \approx 12.7$ over $\Lambda$CDM. We uncover the CMB polarization signatures that can definitively test this scenario with future data.

Figures

Tables

Table 1. Table 1: Maximum likelihood (ML) parameters and constraints for the Λ Λ \Lambda CDM model, the canonical Acoustic Dark Energy (cADE) model, and the general ADE model. Δ ​ χ tot 2 = − 2 ​ Δ ​ ln ⁡ ℒ Δ subscript superscript 𝜒 2 tot 2 Δ ℒ \Delta\chi^{2}_{\rm tot}=-2\Delta\ln{\cal L} reflects the ratio between the maximum likelihood value and that of Λ Λ \Lambda CDM for the joint data.

Model	$\Lambda$CDM	cADE	ADE
$100{\theta }_{\mathrm{MC}}$	1.04115 (1.04110$\pm$0.00028)	1.04062 (1.04064$\pm$0.00031)	1.04072 (1.04081$\pm$0.00035)
${\Omega }_b{h}^2$	0.02246 (0.02241$\pm$0.00014)	0.02267 (0.02271$\pm$0.00022)	0.02270 (0.02263$\pm$0.00022)
${\Omega }_c{h}^2$	0.1170 (0.1174$\pm$0.0009)	0.1268 (0.1268$\pm$0.0032)	0.1274 (0.1242$\pm$0.0032)
$\tau$	0.082 (0.075$\pm$0.012)	0.064 (0.064$\pm$0.012)	0.064 (0.067$\pm$0.013)
$\ln \left({10}^{10}{A}_s\right)$	3.092 (3.079$\pm$0.022)	3.078 (3.078$\pm$0.023)	3.080 (3.081$\pm$0.023)
$n_s$	0.9726 (0.9701$\pm$0.0039)	0.9833 (0.9833$\pm$0.0065)	0.9873 (0.9832$\pm$0.0071)
$f_c$	-	0.082 (0.082$\pm$0.025)	0.086 (0.079$\pm$0.033)
${\log }_{10}{a}_c$	-	-3.45 (-3.46$\pm$0.06)	-3.52 (-3.50$\pm$0.15)
$w_f$	-	1 (fixed)	0.87 (1.89$\pm$0.86)
$c_s^2$	-	1 (fixed)	0.86 (1.07$\pm$0.25)
$H_0$	68.58 (68.35$\pm$0.42)	70.57 (70.60$\pm$0.85)	70.81 (70.20$\pm$0.88)
$\Delta {\chi }_{\mathrm{tot}}^2$	0	-12.7	-14.1

Table 2. Table 2: H 0 subscript 𝐻 0 H_{0} results for the ML cADE, ADE and EDE models and posterior constraints with the joint data set and with CMB polarization data removed (-POL). Δ ​ N Δ 𝑁 \Delta N is the number of additional parameters in addition to the Λ Λ \Lambda CDM ones. ∗ Note that ML ADE in the potential conversion case where c s 2 = w f superscript subscript 𝑐 𝑠 2 subscript 𝑤 f c_{s}^{2}=w_{\rm f} , is essentially the same as the general case but with Δ ​ N = 3 . Δ 𝑁 3 \Delta N=3. The total Δ ​ χ tot 2 Δ subscript superscript 𝜒 2 tot \Delta\chi^{2}_{\rm tot} relative to the Λ Λ \Lambda CDM model is broken down into contributions from the Planck CMB data sets and the local H 0 subscript 𝐻 0 H_{0} measurement.

Model (Data)	$\Delta N$	$H_0$	$\Delta {\chi }_{\mathrm{tot}}^2$	$\Delta {\chi }_{\mathrm{CMB}}^2$	$\Delta {\chi }_{\mathrm{H0}}^2$
cADE	2	70.57(70.60$\pm$0.85)	-12.7	-3.6	-8.8
ADE	3∗,4	70.81(70.20$\pm$0.88)	-14.1	-3.7	-9.6
EDE	4	71.92(71.40$\pm$1.09)	-16.6	-3.7	-12.5
cADE(-POL)	2	71.93(71.55$\pm$1.05)	-12.8	-0.4	-11.2
ADE(-POL)	3∗,4	72.27(71.30$\pm$1.03)	-15.1	-2.4	-11.8
EDE(-POL)	4	72.40(72.35$\pm$1.25)	-15.9	-2.9	-12.1

Equations

$$1+w_{\mathrm{ADE}}(a)=\frac{1+w_{\mathrm{f}}}{[1+(a_c/a{)}^{3(1+w_{\mathrm{f}})/p}{]}^p}.$$

$$f_{\mathrm{ADE}}(a)=\frac{{\rho }_{\mathrm{ADE}}(a)}{{\rho }_{\mathrm{tot}}(a)},$$

$${\sigma }_{\mathrm{CV}}=\{\begin{array}{ll}\sqrt{\frac{2}{2\ell +1}}{C}_{\ell }^{TT},& TT;\\ \sqrt{\frac{1}{2\ell +1}}\sqrt{C_{\ell }^{TT}{C}_{\ell }^{EE}+(C_{\ell }^{TE}{)}^2},& TE;\\ \sqrt{\frac{2}{2\ell +1}}{C}_{\ell }^{EE},& EE.\end{array}$$

$$P(X,\varphi )={\left(\frac{X}{A}\right)}^{\frac{1-c_s^2}{2c_s^2}}X-V(\varphi ),$$

$$\frac{d\varphi }{dN}\sim -\frac{V^{\prime }}{H^2},$$

$$V(\varphi )\approx V({\varphi }_c)+V^{\prime }({\varphi }_c)\Delta \varphi .$$

$${\left(\frac{V^{\prime }}{H}\right)}^2\gtrsim V,$$

$${ϵ}_V{f}_c\gtrsim 1,{ϵ}_V\equiv \frac{M_{\mathrm{Pl}}^2}{2}{\left(\frac{V^{\prime }}{V}\right)}^2.$$

$$V^{\prime \prime }({\varphi }_c)<-\frac{V^{\prime }({\varphi }_c)}{\Delta \varphi }$$

$$V^{\prime \prime }\lesssim H^2.$$

$${\eta }_V\lesssim 2{ϵ}_V,{\eta }_V\equiv M_{\mathrm{pl}}^2\frac{V^{\prime \prime }}{V},$$

$$V(\varphi )=\{\begin{array}{ll}A{\varphi }^m,& \varphi >0,\\ 0,& \varphi \le 0.\end{array}$$

$$2{ϵ}_V={\left(\frac{m}{\varphi }\right)}^2,{\eta }_V=\frac{m(m-1)}{{\varphi }^2},$$

$$V(\varphi )\propto [1-\cos (\varphi /f){]}^n,$$

$$w_{\mathrm{f}}=\frac{n-1}{n+1}.$$

$$c_s^2(a,k)=\{\begin{array}{ll}1,& a\le a_c,\\ \frac{2a^2(n-1){\varpi }^2+k^2}{2a^2(n+1){\varpi }^2+k^2},& a>a_c.\end{array}$$

$$\varpi (a)=\frac{\mathcal{G}}{\mathcal{P}}\sqrt{\frac{6\mathcal{P}+2}{8n}\frac{{\Theta }_i}{\sin {\Theta }_i}}H(a_c)a^{-3w_{\mathrm{f}}}.$$

$$\mathcal{P}(x=a_c/a_{\mathrm{eq}})=\frac{2}{3}\frac{x^2-x+2\sqrt{1+x}-2}{x^2},$$

$$\mathcal{G}(a_c,n)$$

$$V(\varphi )\propto {\varphi }^{2n},$$

Full text

I Introduction

The $\Lambda$CDM model of cosmology has been tested by an extensive number of independent probes, showing its general robustness and remarkable ability to explain a wide range of observational data with only six parameters. Despite its indubitable success, $\Lambda$CDM seems to fail in reconciling distance-redshift measurements when anchored at high redshift by cosmic microwave background (CMB) anisotropies to the same measurements anchored at low redshift by the local distance ladder.

This discrepancy is commonly quantified as tension between the inferences for the Hubble constant ($H_{0}$), and its statistical significance has been steadily increasing with increasing experimental precision. The most recent local estimate of the Hubble constant places its value at $H_{0}=74.03\pm 1.42$ km s*-1* Mpc*-1* Riess et al. (2019), showing a $4.4\sigma$ tension with the value inferred by the Planck 2018 CMB data, $H_{0}=67.4\pm 0.5$ km s*-1* Mpc*-1*, assuming the $\Lambda$CDM model Aghanim et al. (2018).

This tension is mainly a discrepancy between the anchors for the absolute distance scale rather than an indicator of missing physics between the anchors. Once anchored at one end, the same ladder of intermediate redshift measurements from baryon acoustic oscillations (BAOs) to supernovae Type IA (SN) predict the anchor at the other end, leaving little room for missing cosmological physics in between (for a recent assessment and discussion, see Refs. Bernal et al. (2016); Aylor et al. (2019) and references therein). On the high redshift side, the anchor is the CMB sound horizon $r_{s}$. Under $\Lambda$CDM, the shapes of the CMB acoustic peaks calibrate the sound speed and all of the energy densities of species relevant around recombination and thus determine the physical scale of $r_{s}$. Measurements of its angular scale in the CMB then fix the remaining parameter, the cosmological constant or equivalently $H_{0}$. Even beyond $\Lambda$CDM, this measurement determines the distance to recombination and anchors the inverse distance ladder from which $H_{0}$ can again be inferred so that even the most general dark energy or modified gravity model can only moderately alter its value Zhao et al. (2017); Lin et al. (2019); Raveri (2019); Benevento et al. (2019).

Altering the high redshift anchor requires modifying cosmological physics at recombination. Adding extra energy density raises the expansion rate before recombination and lowers the sound horizon $r_{s}$. For example, an extra sterile neutrino or other dark radiation has long been considered as a possible solution Wyman et al. (2014); Dvorkin et al. (2014); Leistedt et al. (2014); Ade et al. (2016a); Lesgourgues et al. (2016); Drewes et al. (2017); Di Valentino et al. (2016); Canac et al. (2016); Feng et al. (2017); Oldengott et al. (2017); Lancaster et al. (2017); Kreisch et al. (2019). However such a component would also change the driving of the acoustic oscillations and damping scale Hu and White (1996), which is now disfavored by increasingly precise CMB data, leaving little ability to raise $H_{0}$ (see Refs. Raveri et al. (2017); Aghanim et al. (2018) for recent assessments).

The problem with the damping scale arises because these additions affect the background expansion like radiation. As pointed out in Ref. Poulin et al. (2018a), this problem can be avoided by making the dark component only important transiently near the epoch of recombination. Specifically, Ref. Poulin et al. (2018a) introduces a component of so-called “early dark energy” (EDE) where a scalar field oscillates anharmonically around the minimum of a periodic potential and finds that the Hubble tension can be efficiently relieved. On the other hand, Ref. Agrawal et al. (2019) finds that for a monomial EDE model, which coincides with the periodic potential Poulin et al. (2018a) at the minimum, the Hubble tension is only partially relieved. Due to the behavior of its perturbations, these EDE scenarios change the amplitudes and phases of the CMB acoustic peaks in complex ways, leading to questions as to the robustness of this method for relieving the Hubble tension.

To build a more robust method to relieve the Hubble tension, we study the general phenomenology of perturbations in a dark fluid which similarly becomes transiently important, here around matter-radiation equality. Since its impact on the CMB comes through the gravitational effects of its own acoustic oscillations, we call this species acoustic dark energy (ADE) and in particular uncover the critical role its sound speed plays in relieving the Hubble tension. We find that the sound speed must vary with the equation of state in the background in a manner consistent with the conversion of potential to kinetic energy for a minimally coupled scalar with a general kinetic term Armendariz-Picon et al. (2001). Unlike the oscillatory EDE models, once released from Hubble drag, the scalar remains kinetic energy dominated until it redshifts away. Indeed for a simple canonical kinetic term, this allows for $H_{0}=70.60\pm 0.85$ with a better fit than $\Lambda$CDM even for the CMB alone and a better total $\chi^{2}$ by 12.7 for 2 extra parameters. This method is also robust and can be exactly realized in a wide class of potentials. We provide both the required conditions on the potential and explicit examples.

This paper is organized as follows. In §II, we introduce the ADE fluid model, its parameters, and the data sets that we use in the analysis. In §III, we discuss the phenomenological impacts of ADE, especially its sound speed, on acoustic driving and CMB polarization. In §IV, we show that ADE models that can relieve the Hubble tension correspond to scalars that convert potential to kinetic energy suddenly upon Hubble drag release and construct a canonical scalar model as proof of principle. In §V, we discuss the relation to the previous work and we conclude in §VI.

II Methodology

Acoustic dark energy is defined to be a perfect dark fluid and is specified by its background equation of state $w_{{\rm ADE}}=p_{\rm ADE}/\rho_{\rm ADE}$ and rest frame sound speed $c_{s}^{2}$ Hu (1998). The latter is only equivalent to the adiabatic sound speed $\dot{p}_{\rm ADE}/\dot{\rho}_{\rm ADE}$ for a barotropic fluid so that in the general case the acoustic phenomenology of linear ADE sound waves, which we shall see is crucial for relieving the Hubble tension, is defined independently of the background.

In order to have a transiently important ADE contribution, we model the ADE equation of state as

[TABLE]

The ADE component therefore changes its equation of state around a scale factor $a=a_{c}$ from $w_{{\rm ADE}}=-1$ to $w_{\rm f}$. Additionally, $p$ controls the rapidity of this transition, with small values corresponding to sharper transitions. Since this parameter does not qualitatively change our results, we use $p=1/2$ unless otherwise specified. We shall see in §IV that this corresponds to a simple quadratic potential for scalar field ADE. This is a generalization of the background of the EDE model Poulin et al. (2018a), discussed in §V, where $p=1$ and the fluid description is approximate.

The ADE background energy density is fully specified once its normalization is fixed, since $w_{\rm ADE}$ determines its evolution. Defining the ADE fractional energy density contribution

[TABLE]

we choose $f_{c}=f_{\rm ADE}(a_{c})$ as the normalization parameter.

The behavior of ADE perturbations is determined by their rest frame sound speed $c_{s}^{2}(a,k)$ which is, for an effective fluid, a function of both time and scale Hu (1998). In the context of a perfect fluid with a linear dispersion relation, it is scale independent. In particular this holds for K-essence scalar field models Armendariz-Picon et al. (2001), when treated exactly instead of in a time-averaged approximation. We shall return to this point in §V.

The equations of motion for ADE acoustic oscillations depend only on the value of $c_{s}^{2}$, not its time derivative. Since the impact of ADE on cosmological observables is extremely localized in time due to the parametrization of $w_{{\rm ADE}}$, we fix $c_{s}^{2}$ to be a constant, effectively its value at $a_{c}$. In §IV, we construct K-essence models where $c_{s}^{2}$ is strictly constant as a proof of principle, but our analysis holds more generally if we interpret the constant $c_{s}^{2}$ as matching a suitably averaged evolving one.

In our most general case, the ADE model is therefore characterized by four parameters $\{w_{\rm f},a_{c},f_{c},c_{s}^{2}\}$ once $p$ is fixed. When varying these parameters we impose flat, range bound priors: $-4.5\leq\log_{10}a_{c}\leq-3.0$, $0\leq f_{c}\leq 0.2$, $0\leq w_{\rm f}\leq 3.6$ and $0\leq c_{s}^{2}\leq 1.5$. We shall later see that the Hubble tension can be relieved by varying just two of these four parameters, fixing $c_{s}^{2}=w_{\rm f}=1$, corresponding to models where the ADE is a canonical scalar that converts its energy density from potential to kinetic around matter radiation equality (see §IV.2). We refer to this particular ADE model as cADE.

The full cosmological model also includes the six $\Lambda$CDM parameters: the cold dark matter density is characterized by $\Omega_{c}h^{2}$, baryon density by $\Omega_{b}h^{2}$, the angular size of the sound horizon by $\theta_{s}$, the optical depth to reionization by $\tau$, and the initial curvature spectrum by its normalization at $k=0.05$ (wavenumbers throughout are in units of Mpc*-1*, which we drop when no confusion should arise), $A_{s}$ and tilt $n_{s}$. These have the usual non-informative priors. We fix the sum of neutrino masses to the minimal value (e.g., Ref. Long et al. (2018)). We modify the CAMB Lewis et al. (2000) and CosmoMC Lewis and Bridle (2002) codes to include all the models that we discuss, following Ref. Hu (1998). We sample the posterior parameter distribution until the Gelman-Rubin convergence statistic Gelman and Rubin (1992) satisfies $R-1<0.02$ or better unless otherwise stated.

For the principal cosmological data sets, we use the publicly available Planck 2015 measurements of the CMB temperature and polarization power spectra at large and small angular scales and the CMB lensing potential power spectrum in the multipole range $40\leq\ell\leq 400$ (Ade et al., 2016b; Aghanim et al., 2016; Ade et al., 2016c). To expose the Hubble tension, we combine this with the latest measurement of the Hubble constant, $H_{0}=74.03\pm 1.42$ (in units of km s*-1* Mpc*-1* here and throughout) (Riess et al., 2019). To these data sets we add the Pantheon Supernovae sample Scolnic et al. (2018) and BAO measurements from BOSS DR12 Alam et al. (2017), SDSS Main Galaxy Sample Ross et al. (2015), and 6dFGS Beutler et al. (2011). These datasets prevent resolving the Hubble tension by modifying the dark sector only between recombination and the very low redshift universe Aylor et al. (2019).

Our baseline configuration thus contains: CMB temperature, polarization and lensing, BAO, SN and $H_{0}$ measurements. Unless otherwise specified all of our results will be for this combined data set. We include all the recommended parameters and priors describing systematic effects for these data sets.

As we shall see, the CMB polarization data provide an important limitation on the ability to raise $H_{0}$ and future polarization data can provide a definitive test of the ADE models that alleviate the Hubble tension. We therefore also consider the joint data set without CMB polarization data. We refer to this data set as -POL.

III ADE Phenomenology

In this section we discuss the phenomenology and observational implications of ADE and their dependence on its parameters.

At the background level, the addition of ADE increases the total energy density before recombination that changes the expansion history lowering the sound horizon $r_{s}$. This changes the calibration of distance measures not only for the CMB but also the whole inverse distance ladder through BAO to SN. Given the precise angular measurements of the sound horizon $\theta_{s}$, the inverse distance ladder scale is reduced and hence the inferred $H_{0}$ rises.

The prototypical example of this method for relieving the Hubble tension is an extra sterile neutrino that is at least mildly relativistic at recombination. Neutrinos, however, do not provide a good global solution (e.g., Ref. Raveri et al. (2017)) since they behave as free-streaming radiation before recombination and therefore change the phase of the CMB acoustic oscillations as well as the CMB damping scale, the distance a photon random walks through the ionized plasma before recombination, approximately as $\lambda_{D}\propto r_{s}^{1/2}$ Hu and White (1996). A more general dark fluid, on the other hand, can reduce the fraction of the dark component vs. normal radiation before matter radiation equality, allowing the two scales to change in a proportional way Poulin et al. (2018a).

Beyond these background effects, ADE and other dark sector candidates for relieving the Hubble tension, gravitationally drive photon-baryon acoustic oscillations changing the amplitudes and phases of the CMB peaks (e.g., Ref. Raveri et al. (2017); Lin et al. (2019)). ADE perturbations undergo their own acoustic oscillations under its sound horizon, leading to novel CMB driving phenomenology. As detailed in Table 1, this modified phenomenology leads to a maximum likelihood (ML) solution with $w_{\rm f}\approx c_{s}^{2}$ and $H_{0}=70.81$ which improves $\Delta\chi^{2}_{\rm tot}=-14.1$ over ML $\Lambda$CDM, and more generally a finite ADE fraction $f_{c}$ is detected at $2.4\sigma$. In the next section, we focus on the details of these physical effects. We then address the impact of Planck polarization data and the ability of future polarization data to test the ADE solutions to the Hubble tension.

Finally, note that although we do not consider measurements of the amplitude of local structure here, the ML and constraints for ADE are $\sigma_{8}\Omega_{m}^{1/2}$=0.4623 (0.4573$\pm$0.0073), whereas for $\Lambda$CDM with CMB TT only has nearly the same ML but larger errors 0.466 (0.466$\pm$0.013). For our combined data set, $\Lambda$CDM gives the value 0.4488 (0.4486$\pm$0.0056), where the ML is lower since raising $H_{0}$ in $\Lambda$CDM lowers $\sigma_{8}\Omega_{m}^{1/2}$ unlike in ADE. If the tension with weak lensing measurements of the amplitude increases in $\Lambda$CDM in the future, it will disfavor not only $\Lambda$CDM but these ADE models as well.

III.1 Acoustic Driving

Under the sound horizon or Jeans scale of the ADE, its density perturbations acoustically oscillate rather than grow, leading to changes in the decay of the Weyl potential $(\Psi+\Phi)/2$. This decay drives CMB acoustic oscillations, and the ADE impact is especially important for modes that enter the CMB sound horizon near $a_{c}$, roughly $k=0.04$ in the ML ADE model from Table 1. The excess decay is countered by raising the cold dark matter through $\Omega_{c}h^{2}$ since it remains gravitationally unstable on the relevant scales.

For ADE, at the parameter level, this effect is controlled by the sound speed $c_{s}$ in conjunction with the equation of state $w_{\rm ADE}(a)$ through $w_{\rm f}$. These two parameters are hence strongly correlated, as shown in Fig. 1, reflecting degenerate effects on the CMB when they are raised or lowered together. Near the ML solution, this requires $w_{\rm f}\approx c_{s}^{2}$.

We explore this degeneracy in Fig. 2 by showing the evolution of the Weyl potential for this mode in the ML model (red) from Table 1 relative to the same model with no ADE (ML, $f_{c}=0$) as a baseline (black). The Weyl potential is relatively suppressed at $a<a_{c}$, enhanced at $a\sim a_{c}$ and suppressed again at $a\gg a_{c}$ due to ADE. The enhancement and subsequent suppression correspond to the first acoustic compression extremum in the ADE density perturbation and the subsequent Jeans stabilization of the perturbations. The net impact is a reduction in the Weyl potential. This reduction is compensated by raising $\Omega_{c}h^{2}$. For comparison we also show the difference in reverting the value of $\Omega_{c}h^{2}$ in the baseline model to the ML $\Lambda$CDM value (cyan dashed). Since $a_{c}\sim a_{\rm eq}$, ADE becomes important around the same epoch when radiation driving has the maximal impact on the shape of the CMB acoustic peaks. Along with other adjustments in $\Lambda$CDM parameters, in particular $n_{s}$ and $\Omega_{b}h^{2}$, these effects compensate for each other.

This compensation leaves the CMB acoustic peaks nearly unchanged despite raising $H_{0}$ from 68.58 to 70.81. In Fig. 3, we show the data and model residuals relative to the $\Lambda$CDM ML in Table 1. As we can see, this effect does not exacerbate the oscillatory residuals in the data, where the acoustic peaks (vertical lines) are suppressed relative to troughs, which would occur if $H_{0}$ were raised in $\Lambda$CDM. Note that the residuals are scaled to the cosmic variance per $\ell$-mode for the ML $\Lambda$CDM model as:

[TABLE]

We can better understand the origin of the ADE effects and their impact on the CMB by varying $w_{\rm f}$ and $c_{s}^{2}$ independently. Fig 2 (upper) also shows a $+0.4$ variation in each with other parameters held fixed. Increasing $c_{s}^{2}$ makes the ADE acoustic oscillations and Jeans stability occur earlier and as a consequence also cuts into the enhancement. Raising $w_{\rm f}$ has two effects. Before the first compression peak and above the CMB sound horizon, the comoving ADE density perturbation grows adiabatically so that its amplitude grows relative to the photons approximately as $\delta_{\rm ADE}\propto(1+w_{\rm ADE})\delta_{\gamma}$. For $w_{\rm f}>1/3$ this first causes a dip at $a\lesssim a_{c}$ and then an enhancement in the Weyl potential for $a\gtrsim a_{c}$, especially approaching the first compression. A larger $w_{\rm f}$ then suppresses $f(a>a_{c})$ which also causes a relative enhancement at $a>a_{c}$. Combined, these effects imply that for a fixed amount of driving of acoustic oscillations through the decay of the Weyl potential, raising $c_{s}^{2}$ should be compensated by raising $w_{\rm f}$. This is the leading order degeneracy that we see in Fig. 1.

In terms of the residuals, shown in Fig. 3, a positive variation in $c_{s}^{2}$, when not compensated by $w_{\rm f}$, leads to a sharp TT feature around $\ell\sim 500$ near the second TT peak whereas along the degeneracy line the ML model TT residuals remain small. The modes responsible for higher multipoles are not sensitive to ADE perturbation parameters since the Weyl decay that drives them occurred before the ADE became important $a\ll a_{c}$.

The degeneracy is truncated at low $w_{\rm f}$ in Fig. 1. If $w_{\rm f}<1/3$, the ADE redshifts slower than the radiation and thus has a large impact on the driving of CMB acoustic oscillations between $a_{c}$ and recombination which cannot be balanced by the same variations in $c_{s}^{2}$.

III.2 CMB Polarization

Even for the ML ADE model, the compensation between $\Omega_{c}h^{2},c_{s}^{2}$, and $w_{\rm f}$ is imperfect for modes that enter the CMB sound horizon between $a_{c}$ and recombination. These modifications leave distinct imprints on the polarization spectra that already limit the ability of ADE to raise $H_{0}$ using the Planck data and, in the future, can definitively test this scenario. Polarization provides the cleanest signatures of driving on these scales given that it isolates the acoustic oscillations at recombination from the early integrated Sachs-Wolfe that smooths out its signatures in the TT spectrum.

In Fig. 2 (lower), we show the Weyl potential evolution for such a mode, $k=0.01$ in the ML model. While the qualitative behavior is similar to the higher $k$ mode, the balance between $w_{\rm f}$, $c_{s}^{2}$ and $\Omega_{c}h^{2}$ changes. First, the impact of the change in $\Omega_{c}h^{2}$ is relatively higher since the ADE redshifts faster than radiation. Second, the impact of $w_{\rm f}$ is also somewhat higher relative to $c_{s}^{2}$. These changes lead to uncompensated sharp features in the polarization residuals.

For the Planck data, where measurements of the EE spectrum are still noisy, this makes the TE spectrum the most informative for these features (see Fig. 3). In the $\Lambda$CDM model, this sensitivity provides an important constraint on $\Omega_{c}h^{2}$ and hence supporting evidence for a low $H_{0}$ from multipoles $\ell<1000$. The region around $\ell\sim 165$ (between the first EE and TT peaks) is particularly important due to the $2\sigma$ low point compared to the best fit $\Lambda$CDM model shown in Fig. 3 Obied et al. (2017). Note that in $\Lambda$CDM raising $H_{0}$ requires lowering $\Omega_{c}h^{2}$ which raises TE there, making the fit even worse.

In the ADE ML model, the impact of raising $\Omega_{c}h^{2}$ lowers TE in this region providing a better fit to the data. Even without direct $H_{0}$ data, the CMB data favor the ADE model (see Table 2). However raising $H_{0}$ further than the ADE ML would make TE too low at $\ell\lesssim 500$ for the data. Indeed for the -POL dataset, the ML model allows for a larger Hubble constant $H_{0}=72.27$ compared with $70.81$ for our joint dataset.

Raising $c_{s}^{2}$ has the impact of making the ADE more important for Weyl decay and counters the effect of raising $\Omega_{c}h^{2}$. As we can see from Fig. 3, this has the effect of raising TE in this region and degrading the fit. Thus the TE data are also important for disfavoring a canonical scalar field with $c_{s}^{2}=1$ if $w_{\rm f}$ is too low. We shall see in §V that this explains why previously considered models where a canonical scalar field oscillates in its potential must have its initial conditions set to avoid this region.

Finally, given the sharp features in the EE model residuals with up to $\sim 0.3-0.4$ amplitudes relative to cosmic variance per multipole, all of these cases where $H_{0}$ is raised by dark components that also change the driving of the acoustic peaks can be tested to high significance once EE measurements approach the precision of TT measurements today.

IV Potential-Kinetic Conversion

The ADE phenomenology favored by the Hubble tension can be concretely and exactly realized in the K-essence class of dark energy models, where the dark component is a perfect fluid represented by a minimally coupled scalar field $\phi$ with a general kinetic term Armendariz-Picon et al. (2001). More specifically, the class of constant sound speed $c_{s}$ models introduced in §II is given by the Lagrangian density Gordon and Hu (2004)

[TABLE]

where the kinetic term involves $X=-\nabla^{\mu}\phi\nabla_{\mu}\phi/2$ and $A$ is a constant density scale. For a scalar with a canonical kinetic term $c_{s}^{2}=1$, and more generally $w_{\rm ADE}\rightarrow c_{s}^{2}$ if the kinetic term dominates, whereas $w_{{\rm ADE}}\rightarrow-1$ if the potential $V(\phi)$ dominates. The fluid correspondence holds when $\nabla^{\mu}\phi$ remains always timelike; then $c_{s}^{2}=\delta p/\delta\rho$ in constant field gauge or rest frame, where the momentum density of the field vanishes and the potential energy is spatially constant.

The correlation shown in Fig. 1 implies that around the ML ADE model from Table 1, setting $w_{\rm f}=c_{s}^{2}$ provides a good fit to the combined data. Since $w_{\rm ADE}\rightarrow-1$ for $a\ll a_{c}$ and $w_{\rm ADE}\rightarrow w_{\rm f}$ for $a\gg a_{c}$, this suggests that the best fitting $P(X,\phi)$ models are those that suddenly convert nearly all of their energy density from potential to kinetic at $a_{c}$. If we focus on a model that has such a potential to kinetic energy conversion feature, we have $c_{s}^{2}=w_{\rm f}$ and the number of parameters of ADE reduces to three. The ML of this model gives $c_{s}^{2}=0.84$ , $H_{0}=70.84$ and $\Delta\chi^{2}_{\rm tot}=-14.1$ which is nearly identical to ML ADE with one fewer free parameter.

Indeed the preferred region includes $w_{\rm f}=c_{s}^{2}=1$ which corresponds to the case of a canonical field, cADE. From Table 1, we see that the ML cADE model has only a marginally smaller $H_{0}=70.57$ for a $\Delta\chi^{2}_{\rm tot}=-12.7$ and more generally a higher significance to the detection of a finite ADE fraction $f_{c}$ of $3.3\sigma$ given the smaller set of parameters. We shall now consider how to construct a corresponding potential $V(\phi)$.

IV.1 Canonical Conditions

A canonical scalar which converts its potential to kinetic energy around $a_{c}$ provides a simple, concrete example of ADE that alleviates the Hubble tension. To explicitly construct such a model that matches requirements on the two remaining quantities $a_{c}$ and $f_{c}$ we can determine the equivalent requirements for the potential $V(\phi)$.

At $a\ll a_{c}$, the $\phi$ field is stuck on its potential due to Hubble friction and rolls according to

[TABLE]

where $N=\ln a$ denotes e-folds. For the purposes of this qualitative discussion we drop factors of order unity. After $a_{c}$, we want the field to be released from Hubble drag and convert its potential energy to kinetic energy on the e-fold timescale $\Delta N\sim 1$. Defining $\phi_{c}=\phi(a_{c})$ and linearizing the change

[TABLE]

Therefore, around $\phi_{c}$, we want

[TABLE]

or in terms of $f_{c}\sim V/\rho_{\rm tot}$,

[TABLE]

This is the main condition for the potential to kinetic conversion.

For the linearization in Eq. (6) to be valid in the sense of the second order term $\frac{1}{2}(\Delta\phi)^{2}V^{\prime\prime}$ not preventing the conversion, we also want

[TABLE]

so

[TABLE]

Putting these two criteria together,

[TABLE]

where we have restored a factor of 2 so as to match the well-known condition for no tracking solution to exist Steinhardt et al. (1999). Tracking potentials do not work since the scalar field follows an equation of state that is determined by the dominant component of the total energy density rather than the kinetic energy dominated limit. A similar derivation applies to the $c_{s}^{2}\neq 1$ case with a modification to the Hubble drag evolution (5) Gordon and Hu (2004).

Finally we want the field to maintain kinetic energy domination until its energy density has largely redshifted away. This excludes models where the field oscillates around a minimum and so is different from those in Refs.  Poulin et al. (2018a); Agrawal et al. (2019) as we shall discuss in the next section. Furthermore, the fluid description is exact for our models whereas it is only approximate for oscillatory models.

Thus our requirements on the potential are fairly generic and correspond to setting the amplitude and slope of the potential at the desired point of Hubble drag release, along with the condition that the field remains kinetic energy dominated until most of the energy density has redshifted away. A wide class of potentials can satisfy these requirements and we shall give concrete examples next.

IV.2 Canonical Solution

To make these considerations concrete, consider the class of potentials:

[TABLE]

Then for $\phi>0$

[TABLE]

and any $m>0$ satisfies $\eta_{V}<2\epsilon_{V}$. The flat potential at $\phi\leq 0$ prevents the kinetic energy from converting back to potential energy. We choose $A$ and $\phi_{\rm initial}$ to give the desired $f_{c}$ and $a_{c}$.

In Fig. 4 (upper), we show a worked example of this matching. We fix cosmological parameters to the ML ADE model in Table 1 and take a quadratic potential with $m=2$. We find a good match to the form of Eq. (1) with $p=1/2$. This motivates our fixed fiducial choice in §II.

To showcase the robustness of the potential to kinetic conversion mechanism for relieving the Hubble tension, we also consider a quartic potential $m=4$ (Fig. 4, lower). The change in $f_{\rm ADE}(a)$ is itself small and, once a shift in $a_{c}$ is absorbed, corresponds to a slight broadening of the transition. Even this small change can be matched to the general ADE form of Eq. (1) by adopting $p=1$. For the linear $m=1$ case, $p\approx 0.1$ corresponding to a sharper transition. We have also tested that various values of $p$ in this range provide comparable ML solutions to our fiducial $p=1/2$ case. Finally we have tested that the correspondence between $p$ and $m$ holds for non-canonical values of $c_{s}^{2}$ with the Lagrangian (4).

These simple canonical or cADE models still provide good fits to the data as illustrated in Fig. 5 for the ML cADE model of Table 1. The main difference compared with ML ADE is the slight lowering of $H_{0}$ from 70.81 to 70.57. The total improvement over $\Lambda$CDM for 2 extra parameters is $\Delta\chi^{2}_{\rm tot}=-12.7$ with $-3.6$ actually coming from the improved fit to the CMB as shown in Table 2. More concretely, ML cADE makes CMB lensing a little bit worse by $\Delta\chi^{2}_{\rm lens}=+1.1$ but fits the TT and polarization spectrum better by $\Delta(\chi^{2}_{\rm plik}+\chi^{2}_{\rm lowTEB})=-4.7$. If compared to ML $\Lambda$CDM fit to CMB only, the ML cADE fits CMB lensing as well and fits the TT and polarization spectrum better by $\Delta(\chi^{2}_{\rm plik}+\chi^{2}_{\rm lowTEB})=-1.6$. 111 We have also explicitly checked that a direct solution for the scalar field Klein Gordon equation is nearly indistinguishable from a cADE model with the best matching parameters, e.g. $\Delta\chi^{2}=1.4$ for a cosmic variance limited TT, TE, EE measurement to $\ell\leq 2000$. This holds for these gravitationally sourced, or adiabatic, field perturbations whereas modeling isocurvature fluctuations from initial field perturbations from inflation would require matching the radiation dominated evolution in the equation of state $(1+w)\propto a^{4}$ implied by Hubble friction through Eq. (5) Gordon and Hu (2004), which Eq. (1) does not do.

Fig. 6 shows the parameter covariances and posteriors in the cADE model. The centered values for the ML parameters is indicative of the nearly Gaussian posteriors and reflects the fact that the parameters are constrained mainly by the data rather than the priors. The one exception is $a_{c}$ since if $a_{c}\rightarrow 0$ any $f_{c}$ is equivalent to $\Lambda$CDM so that the prior volume begins to matter. Even in this case, $\Lambda$CDM is sufficiently disfavored so that constraints on $a_{c}$ are data not prior driven. Correspondingly the ADE fraction is significantly detected with $f_{c}=0.082\pm 0.025$.

This model also illustrates the main compensation between raising the ADE fraction $f_{c}$ and raising the CDM density $\Omega_{c}h^{2}$ as well as adjusting $\Omega_{b}h^{2}$ and $n_{s}$ slightly higher to minimize the data residuals (see Fig. 6). The change in $\theta_{*}$ to lower values is also notable. The modifications in driving make a small change in the phasing of the CMB acoustic peaks relative to its sound horizon. Note that in $\Lambda$CDM, $\theta_{*}$ drifts lower once the high multipoles $\ell>800$ are included Aghanim et al. (2017).

Finally under the -POL data set, the ADE canonical model allows a higher $H_{0}=71.55\pm 1.05$ and ML value of 71.93. This is because of the limitations the TE spectrum around $\ell\lesssim 500$ places on these solutions as discussed in §III.2.

V Relation to Prior Work

In Ref. (Poulin et al., 2018a), a canonical scalar field component, referred to as early dark energy, with a potential

[TABLE]

plays a similar role as our ADE. Unlike ADE, EDE oscillates after being released from Hubble drag and Ref. Poulin et al. (2018b) finds that the time averaged background equation of state can be modeled by Eq. (1) with $p=1$ and

[TABLE]

As we have discussed in the previous sections, $w_{\rm f}$ is a relevant parameter for the resolution of Hubble tension, so its adjustment should be considered a parameter variation in the EDE model in spite of $n$ taking discrete integer values.

The time-averaged behavior of perturbations is described by a fluid approximation with a rest frame sound speed Poulin et al. (2018b):

[TABLE]

Unlike in our case, this fluid description is approximate, especially at $a_{c}$. The time dependence of $\varpi$ is fixed by the parameters $a_{c},w_{\rm f}$ and the initial field position $\Theta_{i}=\phi_{i}/f$ Poulin et al. (2018b):

[TABLE]

Here ${\mathcal{P}}=H(a_{c})t$ and we approximate it as:

[TABLE]

while

[TABLE]

Note that for $w_{\rm f}>0$, $\varpi$ decreases with $a$, and so the sound speed evolves from $1$ at $a\leq a_{c}$ back to $1$ at late times, with higher $k$ exhibiting smaller amplitude deviations, and a $k$-dependent minimum $c_{s}^{2}(a_{c},k)$.

The EDE model therefore has four parameters $a_{c},f_{c},\Theta_{i}$ and $n$. Following Ref. Poulin et al. (2018b), we choose the best value $n=3$, which corresponds to $w_{\rm f}=1/2$, and conduct an MCMC likelihood analysis on the remaining parameters. We treat $a_{c}$ and $f_{c}$ as in the ADE model and impose a flat prior on $0\leq\Theta_{i}/\pi\leq 1$. Because of the large parameter volume of degenerate models around $\Lambda$CDM, we only sample the posterior until $R-1<0.05$ which should give an adequate, but not perfect, estimate of parameter constraints out to $95\%$ C.L. The results are compared with our ADE model in Table 2. The EDE ML model allows a slightly higher $H_{0}=71.92$ and hence a better fit to the data $\Delta\chi^{2}_{\rm tot}=-2.5$ for one extra parameter over ML ADE with $c_{s}^{2}=w_{\rm f}$ or $\Delta\chi^{2}_{\rm tot}=-3.9$ for two extra parameters compared with the cADE ML. Note that with the -POL dataset, the EDE and ADE models have comparable performances.

The main phenomenological difference between the ADE and EDE models is the parametrization of the sound speed. The sound speed was indeed also varied in Ref. Poulin et al. (2018b) but its impact was not discussed. As we have seen in our ADE model a low $w_{\rm f}$ generally requires a low $c_{s}^{2}$. In the EDE model, this translates into specific requirements for the initial phase $\Theta_{i}$. In Fig. 7, we show the relationship between constraints on $\Theta_{i}$, $H_{0}$ and the minimum sound speed at $k=0.01,0.04$. As we can see, achieving a higher value of $H_{0}$ requires a large initial phase and its ML value is $\Theta_{i}/\pi=0.90$. The 68% confidence region is $0.72\leq\Theta_{i}/\pi\leq 0.94$ and $\Theta_{i}/\pi<1/2$ is excluded at 93% C.L. Note that the upper range exceeds the value required for the validity of the fluid mapping approximation, $\Theta/\pi\sim 0.96$ Poulin et al. (2018b).

The reason for this preference is that $\Theta_{i}$ controls the minimum sound speed. Following the degeneracy line in Fig. 1 to $w_{n}=1/2$, we would expect that for a constant sound speed $c_{s}^{2}\approx 0.77$. Given the effective EDE sound speed of Eq. (16), this represents an average over the relevant timescales and wavemodes. For the ML EDE model $c_{s}^{2}(a_{c},0.01)=0.51$ and $c_{s}^{2}(a_{c},0.04)=0.63$ as the minimum value for each $k$-mode.

The scale dependence of the sound speed also explains the slightly better fit to CMB data, specifically the TE data. In Fig. 5, we compare the EDE and cADE residuals for their respective ML models. Notice that the TT residuals are very similar. However, in TE, by allowing the sound speed to decrease in the $k$ range associated with $\ell<500$, the EDE model exposes more of the driving reduction at $\ell\sim 200$ from raising $\Omega_{c}h^{2}$ as discussed above but now without adverse consequences elsewhere. This in turn better fits the low TE residuals and allows $H_{0}$ to increase further relative to the ADE model. We also show the impact of reducing $\Theta_{i}$ in the ML EDE model, making the sound speed closer to one at all times. The most significant effect is localized to $\ell<500$ and in particular destroys the pattern of lower TE at $\ell\sim 200$ vs $\ell\sim 400$ compared with cADE.

We conclude that the small improvement of the EDE over ADE fit requires a specific range in the initial phase that lowers the sound speed in a scale-dependent way. Comparing to the canonical ADE mode, this improvement gives $\Delta\chi^{2}_{\rm tot}=-3.9$ for two extra parameters $w_{n},\Theta_{i}$ and is therefore marginal. In the future, polarization measurements that approach the cosmic variance limit can distinguish between the EDE and ADE classes. For example, we forecast that with cosmic variance TT,TE,EE measurements to $\ell\leq 2000$, the current best fit EDE model differs from the closest ADE model by $\Delta\chi^{2}=22.4$. Furthermore the ADE model provides a general class of exact solutions where the potential energy is converted quickly to kinetic, whereas the EDE model requires a specific set of initial conditions to achieve a similar phenomenology with an approximation to an oscillating field.

Relatedly, Ref. Agrawal et al. (2019) considers a model where the scalar field oscillates in a monomial potential

[TABLE]

with parameters adjusted to reproduce the EDE phenomenology. This coincides with Eq. (14) only near the bottom of the potential, $\Theta_{i}\ll\pi/2$ where the potential is convex rather than concave. As pointed out in Ref. Agrawal et al. (2019) the model has significantly worse performances than the EDE model. We identify here that this is related to the initial field being in the concave rather than the convex part of the potential which raises the sound speed.

Finally, while this work was nearing completion, Ref. Alexander and McDonough (2019) proposed that a fast-roll or kinetic energy dominated period in a two-field model might relieve the Hubble tension.

VI Discussion

Acoustic dark energy, appearing around the epoch of matter radiation equality, can substantially relieve the tension between CMB inference of $H_{0}$ and local measurements, exhibited in the $\Lambda$CDM model. The presence of extra energy density lowers the CMB sound horizon that anchors the inverse distance ladder for BAO and SN, while its disappearance before and after equality allows for a good fit to CMB data in the damping tail. Furthermore by introducing ADE at equality, the gravitational effects of raising the cold dark matter density can be balanced by the acoustic oscillations in the ADE itself.

Our main findings regarding the Hubble tension are summarized in Table 2 and Fig. 8. In all cases relieving the Hubble tension requires ADE to be a $\sim 8\%$ contribution to the total energy density around matter radiation equality, leading to at least a two parameter extension to $\Lambda$CDM.

In the general ADE class of models, the acoustic phenomenology is controlled by two additional parameters, the asymptotic equation of state $w_{\rm f}$ at late times, and the sound speed $c_{s}^{2}$. The sound speed plays a crucial role in the gravitational driving of CMB acoustic oscillations through its impact on the Weyl potential, leading to a strong correlation between the two, consistent with $w_{\rm f}=c_{s}^{2}\approx 1$ around the maximum likelihood model.

Fixing these two parameters to the canonical model, cADE, leads to only a minor degradation in the ability of ADE to relax the Hubble tension in the ML model, $H_{0}=70.57$ vs. $H_{0}=70.81$. In fact, the two parameter model has the advantage of providing a more significant detection of the ADE fraction $f_{c}=0.082\pm 0.025$ and thus allows less parameter volume around the $\Lambda$CDM limit, producing a posterior centered around the ML: $H_{0}=70.60\pm 0.85$ in cADE vs. the lower $H_{0}=70.20\pm 0.88$ shift in ADE.

Note that the fit to the CMB data themselves improves and combined with the higher $H_{0}$ provided for an improvement of $\Delta\chi^{2}_{\rm tot}=-14.1$ for 3 parameters ($w_{\rm f}=c_{s}^{2}$) and $\Delta\chi^{2}_{\rm tot}=-12.8$ for 2 parameters ($w_{\rm f}=c_{s}^{2}=1$) respectively.

This class of $w_{\rm f}=c_{s}^{2}\sim 1$ ADE models corresponds to scalar fields which convert their potential to kinetic energy efficiently around their release from Hubble drag. By setting this epoch to be around matter radiation equality, we obtain a robust mechanism for relieving the Hubble tension in a wide class of potentials.

For canonical scalar fields, we explicitly determine the requirements on the potential: that its slope allows for Hubble drag release around equality where its amplitude is set to the $\sim 8\%$ fraction required by the data. Any potential that obeys this property and efficiently converts potential to kinetic energy until the latter redshifts away will satisfy these requirements. As a proof of principle, we explicitly construct an example where the potential is locally quadratic around its release. In this model, the timing of the release to equality is not explained but the identification of this coincidence may lead to more sophisticated models where it is.

The robustness and generality of this potential-kinetic conversion mechanism for relieving the Hubble tension separates it from similar models in the literature. In Refs. Poulin et al. (2018a); Agrawal et al. (2019), the EDE scalar field oscillates after Hubble drag release leading to an effective fluid described by time-averaged values of $w_{\rm f}$ and $c_{s}^{2}$. Converting our requirements on the relationship between the two, we find that in Ref. Poulin et al. (2018a), where the potential is periodic, the initial field must be on the concave part of the potential, and near the maximum to best relieve the Hubble tension. This also explains the poorer fits in Ref. Agrawal et al. (2019), where the potential is convex and matches the periodic potential only near the minimum.

The periodic EDE model Poulin et al. (2018a); Agrawal et al. (2019) allows for a slightly higher $H_{0}=71.40\pm 1.09$ and better $\Delta\chi^{2}_{\rm tot}=-16.6$ for four parameters as compared with ADE. Most of this improvement comes from the fit to the Planck polarization data. As shown in Fig. 8, without these data, EDE and ADE perform similarly at ML with $H_{0}=72.27$ vs $72.40$ respectively. The reason is that between equality and recombination, changes in acoustic driving between raising the CDM density and adding the dark component no longer cancel. This is in fact a beneficial feature of both models since the Planck TE data show low residuals with respect to the ML $\Lambda$CDM model around $\ell\sim 200$. By allowing the effective sound speed to depend on wavenumber, the EDE model fits this region better without violating constraints elsewhere. However, when compared with our cADE model this extra improvement of $\Delta\chi^{2}_{\rm tot}\approx 3.9$ comes at the cost of two extra parameters and a less robust mechanism for relieving Hubble tension.

Finally, in all cases, the predicted deviations in the EE power spectrum, while not at a level testable by Planck data, are highly significant compared with cosmic variance. Future polarization data can provide key tests for these and other dark component explanations of the Hubble tension.

Acknowledgements.

We thank Nicola Bartolo, Daniel Grin, Tanvi Karwal, Macarena Lagos, Michele Liguori, Samuel Passaglia, Vivian Poulin for useful discussions. MXL, WH and MR are supported by U.S. Dept. of Energy contract DE-FG02-13ER41958 and the Simons Foundation. Computing resources were provided by the University of Chicago Research Computing Center through the Kavli Institute for Cosmological Physics at the University of Chicago. GB acknowledges financial support from Fondazione Ing. Aldo Gini.