Bayesian inflationary reconstructions from Planck 2018 data

TL;DR

This paper uses Bayesian methods to reconstruct the primordial power spectrum from Planck 2018 data, revealing potential oscillatory features and confirming the consistency of the standard power-law model across most scales.

Contribution

It introduces three non-parametric Bayesian reconstruction techniques for primordial spectra, highlighting potential oscillatory features and robustness of cosmological parameters.

Findings

Hints of oscillatory features in primordial spectrum at certain scales

No significant deviations from power-law spectrum across most scales

Reconstruction methods are consistent with standard cosmological parameters

Abstract

We present three non-parametric Bayesian primordial reconstructions using Planck 2018 polarization data: linear spline primordial power spectrum reconstructions, cubic spline inflationary potential reconstructions and sharp-featured primordial power spectrum reconstructions. All three methods conditionally show hints of an oscillatory feature in the primordial power spectrum in the multipole range $\ell\sim20$ to $\ell\sim50$, which is to some extent preserved upon marginalization. We find no evidence for deviations from a pure power law across a broad observable window ($50\lesssim\ell\lesssim2000$), but find that parameterizations are preferred which are able to account for lack of resolution at large angular scales due to cosmic variance, and at small angular scales due to Planck instrument noise. Furthermore, the late-time cosmological parameters are unperturbed by these extensions to the primordial power spectrum. This work is intended to provide a background and give more details of the Bayesian primordial reconstruction work found in the Planck 2018 papers.

Tables

Table 1. Table 1: The prior distributions on late-time cosmological parameters and Planck nuisance parameters for all analyses. The parameters of each Gaussian distribution are defined as [ μ , σ ] 𝜇 𝜎 [\mu,\sigma] , and the above distributions combine to make a truncated Gaussian distribution on the nuisance parameters A k ​ S ​ Z superscript 𝐴 𝑘 𝑆 𝑍 A^{kSZ} and A 143 t ​ S ​ Z subscript superscript 𝐴 𝑡 𝑆 𝑍 143 A^{tSZ}_{143} . The nuisance priors are the default ones in CosmoMC , whilst the cosmological priors are narrowed to speed up sampling, but remain sufficiently wide to effectively include the entire posterior mass. Also indicated is each parameter’s speed with respect to the CosmoMC fast-slow hierarchy.

Parameters	Prior type	Prior parameters	Speed
${\Omega }_b{h}^2$	uniform	$[0.019,0.025]$	slow
${\Omega }_c{h}^2$	uniform	$[0.095,0.145]$	slow
$100{\theta }_{MC}$	uniform	$[1.03,1.05]$	slow
$\tau$	uniform	$[0.01,0.4]$	slow
$y_{\mathrm{cal}}$	Gaussian	$1\pm 0.0025$	semi-slow
$A_{217}^{CIB}$	uniform	$[0,200]$	fast
${\xi }^{tSZ-CIB}$	uniform	$[0,1]$	fast
$A_{143}^{tSZ}$	uniform	$[0,10]$	fast
$A_{100}^{PS}$	uniform	$[0,400]$	fast
$A_{143}^{PS}$	uniform	$[0,400]$	fast
$A_{143\times 217}^{PS}$	uniform	$[0,400]$	fast
$A_{217}^{PS}$	uniform	$[0,400]$	fast
$A^{kSZ}$	uniform	$[0,10]$	fast
$A_{100}^{\mathrm{dust}TT}$	Gaussian	$8.6\pm 2$	fast
$A_{143}^{\mathrm{dust}TT}$	Gaussian	$10.6\pm 2$	fast
$A_{143\times 217}^{\mathrm{dust}TT}$	Gaussian	$23.5\pm 8.5$	fast
$A_{217}^{\mathrm{dust}TT}$	Gaussian	$91.9\pm 20$	fast
$A_{100}^{\mathrm{dust}TE}$	Gaussian	$0.13\pm 0.042$	fast
$A_{100\times 143}^{\mathrm{dust}TE}$	Gaussian	$0.13\pm 0.036$	fast
$A_{100\times 217}^{\mathrm{dust}TE}$	Gaussian	$0.46\pm 0.09$	fast
$A_{143}^{\mathrm{dust}TE}$	Gaussian	$0.207\pm 0.072$	fast
$A_{143\times 217}^{\mathrm{dust}TE}$	Gaussian	$0.69\pm 0.09$	fast
$A_{217}^{\mathrm{dust}TE}$	Gaussian	$1.938\pm 0.54$	fast
$c_{100}$	Gaussian	$1.0002\pm 0.0007$	fast
$c_{217}$	Gaussian	$0.99805\pm 0.00065$	fast
$A^{kSZ}+1.6A_{143}^{tSZ}$	Gaussian	$9.5\pm 3$	fast

Table 2. Table 2: The prior distributions on early-time cosmological parameters for the primordial power spectrum reconstruction.

Parameters	Prior type	Prior range
$N$	discrete uniform	$[1,9]$
${𝒫}_1,\cdots ,{𝒫}_N$	$\log$-uniform	${10}^{-10}[e^2,e^4]$
	sorted $\log$-uniform	$[{10}^{-4},{10}^{-0.3}]$

Table 3. Table 3: The prior distributions on early-time cosmological parameters for the inflationary potential reconstruction. ϕ ~ min subscript ~ italic-ϕ \tilde{\phi}_{\min{}} and ϕ ~ max subscript ~ italic-ϕ \tilde{\phi}_{\max{}} are defined by the observable window of the unperturbed potential. There is a further prior constraint in that we require that the inflaton should evolve in an inflating phase throughout the observable window and that the inflaton should be rolling downhill from negative to positive ϕ italic-ϕ \phi throughout.

Parameters	Prior type	Prior range
$N$	Discrete Uniform	$[0,8]$
$\ln V_{\ast }$	Uniform	$[-25,-15]$
$\frac{d\ln V_{\ast }}{d\varphi }$	Log-Uniform	$[{10}^{-3},{10}^{-0.3}]$
$\frac{d^2\ln V_1}{d{\varphi }^2},\dots ,\frac{d^2\ln V_N}{d{\varphi }^2}$	Uniform	$[-0.5,0.5]$
${\varphi }_1,\dots ,{\varphi }_N$	Sorted uniform	$[{\tilde{\varphi }}_{\min },{\tilde{\varphi }}_{\max }]$
$\ln {10}^{10}{𝒫}_{ℛ}\left(k\right)$	Indirect constraint	$[2,4]$

Table 4. Table 4: The prior distributions on early-time cosmological parameters for the sharp feature reconstructions.

Parameters	Prior type	Prior range
$N$	discrete uniform	$[0,8]$
$A_s$	uniform	${10}^{-10}[e^2,e^4]$
$n_s$	uniform	$[0.8,1.2]$
$h_1,\cdots ,h_N$	uniform	$[-1,1]$
	sorted $\log$-uniform	$[{10}^{-4},{10}^{-0.3}]$
${\Delta }_1,\cdots ,{\Delta }_N$	uniform	$[0,1]$
$\ln {10}^{10}{𝒫}_{ℛ}\left(k\right)$	Indirect constraint	$[2,4]$

Equations

$$\ddot{\varphi }+3H\dot{\varphi }+\frac{dV}{d\varphi }=0,$$

$$H^2=\frac{1}{3M_{\mathrm{p}}^2}\left(\frac{1}{2}{\dot{\varphi }}^2+V(\varphi )\right),$$

$${\mathcal{R}}_k^{\prime \prime }+2\frac{z^{\prime }}{z}{\mathcal{R}}_k^{\prime }+k^2{\mathcal{R}}_k$$

$${\mathcal{T}}_k^{\prime \prime }+2\frac{a^{\prime }}{a}{\mathcal{T}}_k^{\prime }+k^2{\mathcal{T}}_k$$

$$z=\frac{a\dot{\varphi }}{H},\eta =\int \frac{dt}{a}$$

$${\mathcal{P}}_{\mathcal{X}}(k)=\underset{aH\gg k}{\lim }\frac{k^3}{2{\pi }^2}{\left|{\mathcal{X}}_k\right|}^2,\mathcal{X}\in \{\mathcal{R},\mathcal{T}\}.$$

$$H=\sqrt{\frac{V({\varphi }_{\mathrm{init}})}{3M_{\mathrm{p}}^2}},\dot{\varphi }=-\frac{V^{\prime }({\varphi }_{\mathrm{init}})}{3H},$$

$${\mathcal{R}}_k=\frac{1}{z\sqrt{2k}}{e}^{-ik\eta },{\mathcal{T}}_k=\frac{1}{a\sqrt{2k}}{e}^{-ik\eta }.$$

$$P(\Theta |D,M)$$

$$P(D|M)$$

$$P(M_i|D)$$

$$P(D)$$

$$D_{KL}(P|Q)=\int \ln \left[\frac{P(x)}{Q(x)}\right]P(x)dx,$$

$$\mathrm{Lin}(x;{\Theta }_f)$$

$${\Theta }_f$$

$${\Theta }_{\mathrm{c}}=({\Omega }_b{h}^2,{\Omega }_c{h}^2,100{\theta }_{MC},\tau ).$$

$${\Theta }_{\mathrm{n}}=($$

$$A_{217}^{PS},A^{kSZ},A_{100}^{\mathrm{dust}TT},A_{143}^{\mathrm{dust}TT},A_{143\times 217}^{\mathrm{dust}TT},A_{217}^{\mathrm{dust}TT},$$

$$A_{100}^{\mathrm{dust}TE},A_{100\times 143}^{\mathrm{dust}TE},A_{100\times 217}^{\mathrm{dust}TE},A_{143}^{\mathrm{dust}TE},A_{143\times 217}^{\mathrm{dust}TE},$$

$$A_{217}^{\mathrm{dust}TE},c_{100},c_{217}).$$

$$\ln {\mathcal{P}}_{\mathcal{R}}(k)=\ln A_{\mathrm{s}}+(n_{\mathrm{s}}-1)\ln \left(\frac{k}{k_{\ast }}\right),$$

$$\ln 10^{10}{\mathcal{P}}_{\mathcal{R}}(k)=$$

$${\Theta }_{\mathcal{P}}=$$

$$\ln 10^{10}{\mathcal{P}}_1,\cdots ,\ln 10^{10}{\mathcal{P}}_N).$$

$${\theta }_n={\theta }_{\min }+({\theta }_{n+1}-{\theta }_{\min })x_n^{1/n},{\theta }_{N+1}={\theta }_{\max },$$

$$\ln V=$$

$$+{\int }_{{\varphi }_{\ast }}^{\varphi }d{\varphi }^{\prime }{\int }_{{\varphi }_{\ast }}^{{\varphi }^{\prime }}d{\varphi }^{\prime \prime }\mathrm{Lin}({\varphi }^{\prime \prime };{\Theta }_V)$$

$${\Theta }_V=$$

$${\epsilon }_V=\frac{1}{2}{\left(\frac{1}{V}\frac{d{V}_{\ast }}{d\varphi }\right)}^2,{\eta }_V=\frac{1}{V}\frac{d^2{V}_i}{d{\varphi }^2}\Rightarrow \frac{d^2\ln V}{d{\varphi }^2}={\eta }_V-2{\epsilon }_V.$$

$$\ln {\mathcal{P}}_{\mathcal{R}}(k)=$$

$$+\sum _{i=1}^N{h}_i\left[|{\log }_{10}k-{\log }_{10}{k}_i|<\frac{{\Delta }_i}{2}\right],$$

$${\Theta }_{\mathcal{P}}=$$

Full text

I Introduction

The final release of Planck satellite data Planck Collaboration (2018a, 2019, b, c) provides an unprecedented window onto the cosmic microwave background (CMB). These high-resolution CMB anisotropy data give constraints on the state of the Universe in its earliest observable stage. Assuming a theory of inflation, the primordial power spectrum of curvature perturbations provides an indirect probe of ultra-high energy physics.

This paper focuses on non-parametric reconstructions of primordial physics. The aim of such analyses is to provide information on quantities and functions of interest that are arguably model-independent. While unambiguous scientific detections will only ever result from a consideration of specific, physically motivated models, results from reconstructions such as these can be used to inform and guide observational and theoretical cosmology, providing insight and evidence for interesting features not clearly visible in the data when using standard modeling assumptions.

Throughout we adopt a fully Bayesian framework, treating our non-parametric reconstruction functions using priors, posteriors and evidences to marginalize out factors that are irrelevant to physical quantities of interest. We reconstruct both the inflationary potential and the primordial power spectrum directly using spline and feature-based reconstructions, in a manner related but not identical to the existing literature (Vázquez et al., 2012a; Aslanyan et al., 2014; CORE Collaboration, 2018; Planck Collaboration, 2014, 2016; Hee et al., 2017; Vázquez et al., 2012b; Millea and Bouchet, 2018; Olamaie et al., 2018; Gauthier and Bucher, 2012; Dvorkin and Hu, 2010a, b, 2011; Miranda et al., 2015; Obied et al., 2017, 2018).

In Sec. II we review the relevant background theory in primordial cosmology, Bayesian inference, non-parametric reconstruction and CMB data. Section III reconstructs the primordial power spectrum directly using a linear interpolating spline. Section IV takes the analysis one step back and reconstructs the inflationary potential using a cubic spline, treating the primordial power spectrum as a derived quantity. Section V works with a parameterization that is more suited for reconstructing sharp features in the primordial power spectrum as a complementary approach to that of Sec. III. Section VII draws conclusions from all three analyses.

II Background

II.1 Primordial cosmology

We begin by summarizing the background theory and establish notation. For a more detailed discussion of inflationary cosmology and perturbation theory, we recommend Mukhanov et al. (1992) or Baumann (2009, 2018).

The evolution equations for a spatially homogeneous, isotropic and flat universe filled with a scalar field $\phi$ with arbitrary potential $V(\phi)$ are

[TABLE]

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter, $a$ is the scale factor of the Universe, dots denote derivatives with respect to cosmic time $\dot{f}\equiv\frac{\operatorname{d}\!{f}}{\operatorname{d}\!{t}}$ and $M_{\mathrm{p}}$ is the reduced Planck mass. For most potentials $V(\phi)$, solutions to Eqs. 1 and 2 rapidly converge on the attractor slow roll state, satisfying $\dot{\phi}^{2}\ll V(\phi)$.

The evolution equations for the Fourier $k$-components of the gauge-invariant comoving curvature $\mathcal{R}$ and tensor $\mathcal{T}$ perturbations are

[TABLE]

where primes denote derivatives with respect to conformal time $f^{\prime}\equiv\frac{\operatorname{d}\!{f}}{\operatorname{d}\!{\eta}}$. These equations have the property that in-horizon solutions ($k\gg aH$) oscillate with time-varying amplitude and frequency, whilst out-of-horizon solutions ($k\ll aH$) freeze out. The dimensionless primordial power spectra of these perturbations are defined as

[TABLE]

Initial conditions for the background Eqs. 1 and 2 may be set using the slow roll approximation:

[TABLE]

where $V^{\prime}$ denotes the derivative of $V$ with respect to $\phi$. Whilst solutions set with these initial conditions do not lie precisely on the attractor state, they rapidly converge on it. Providing that $\phi_{\mathrm{init}}$ is chosen self-consistently (Mortonson et al., 2011; Easther and Peiris, 2012; Noreña et al., 2012) with enough additional evolution so that any transient effects are lost, these initial conditions are equivalent to choosing the background solution to be the attractor.

For the perturbation Eqs. 3 and 4, Bunch-Davies initial conditions are chosen such that the Mukhanov variables match onto the de-Sitter vacuum solutions

[TABLE]

Providing that the $k$-mode lies well within the horizon ($k\gg aH$), this is the canonical choice for initializing the perturbation spectrum, although other vacua are available Armendáriz-Picón and Lim (2003); Handley et al. (2016); Handley (2019a); Agocs et al. (2018).

Whilst Eqs. 1, 2, 3 and 4 take their simplest form using cosmic and conformal time as variables, for numerical stability it is more prudent to choose a time-like parameter which does not saturate during inflation, such as cosmic time $t$, or the number of $e$-folds $N=\log(a)$. We choose the logarithmic comoving horizon $\log(aH)$ as the independent variable for our analyses. In this form Eqs. 1, 2, 3 and 4 become complicated, so to avoid typographical errors we generate Fortran source code using the Maple Toronto: Maplesoft, a division of Waterloo Maple Inc (2018) computer algebra package. The numerical integration of all differential equations was performed using the NAG library The Numerical Algorithms Group , Oxford, United Kingdom() (NAG).

It should also be noted that Eqs. 3 and 4 are usually phrased in terms of the Mukhanov variables $v=z\mathcal{R}$ and $h=a\mathcal{T}$, but $\mathcal{R}$ and $\mathcal{T}$ prove to be more numerically stable as they have the attractive property that they explicitly freeze out.

II.2 Bayesian statistics

Once the primordial power spectra $\mathcal{P}_{\mathcal{R},\mathcal{T}}(k)$ have been determined, these form the initial conditions for Boltzmann codes Lewis et al. (2000); Blas et al. (2011). For a universe described by a cosmological model $M$ with corresponding late-time parameters $\Theta_{\mathrm{c}}$ and primordial power spectra $\mathcal{P}$, a Boltzmann code computes CMB power spectra $C_{\mathrm{\ell}}$ in both temperature and polarization. These CMB power spectra may then be fed into cosmological likelihood codes Planck Collaboration (2019), which typically depend on additional nuisance parameters $\Theta_{\mathrm{n}}$ associated with the experiment. The end result is a likelihood $P(D|\Theta,M)$ of the parameters $\Theta=(\mathcal{P}(k),\Theta_{\mathrm{c}},\Theta_{\mathrm{n}})$ given CMB data $D$, and a cosmological model $M$.

We may formally invert the conditioning on $\theta$ in the likelihood using Bayes theorem

[TABLE]

where the first expression above should be read as “posterior is likelihood times prior over evidence”, and the second expression indicates that the evidence is the normalizing constant of Eq. 9, and is a multidimensional marginalization of the likelihood over the prior. The evidence may also be used in a Bayesian model comparison, to assess the relative merits of a set of competing models $\{M_{i}\}$

[TABLE]

In the event of uniform priors over models, the evidence establishes the relative probability weighting to give to models describing the same data $D$.

Throughout this work, we use a modified version of CAMB Lewis et al. (2000) to compute $C_{\ell}$ power spectra and CosmoChord Handley (2019b) (a modified version of CosmoMC Lewis and Bridle (2002); Lewis (2013)) to interface the likelihoods. To sample the posterior and compute evidences we make use of the nested sampling Skilling (2006) algorithm PolyChord Handley et al. (2015a, b). The default Metropolis–Hastings sampler in CosmoMC is insufficient both due to the complexity of the posteriors that must be navigated, and the requirement of evidence computation. Furthermore, PolyChord is required in place of the previous nested sampling algorithm MultiNest Feroz and Hobson (2008); Feroz et al. (2009, 2013) due to the high dimensionality of the full Planck likelihood with nuisance parameters. As an added bonus, PolyChord has the ability to exploit the fast-slow cosmological hierarchy Lewis (2013), which greatly speeds up the sampling. Most importantly all parameters associated with the primordial power spectrum are “semi-slow”, given that one does not need to recompute transfer functions upon changing the primordial power spectrum.

II.3 Functional inference

For our reconstructions, the quantities of interest are functions $f(k;\Theta_{f})$ of wavenumber $k$, parameterized by a set of parameters $\Theta_{f}$, which presents a challenge in both plotting and quantifying our results.

We utilize two related techniques to plot the posterior of a function $f(k;\Theta_{f})$. First, we can generate equally-weighted samples of $\Theta_{f}$, and therefore of the function $f$, and plot each sample as a curve on the $(k,f)$ plane. In general, we simultaneously plot prior samples in red, and posterior samples in black. An example of such a plot can be found in the upper-left panel of Sec. III. For the second type of plot, we first compute the marginalized posterior distribution $P(f|k)$ of the dependent variable $f$ conditioned on the independent variable $k$ using Gaussian kernel density estimation. The iso-probability credibility intervals are then plotted in the $(k,f)$ plane, with their mass converted to $\sigma$-values via an inverse error function transformation. An example of this kind of plot can be seen in the upper-right panel of Sec. III. The code for producing such plots is published in Ref. Handley (2018).

To quantify the constraining power of a given reconstruction, we use the conditional Kullback-Leibler (KL) divergence Kullback and Leibler (1951) as exemplified by Hee et al. (2017). For two distributions $P(x)$ and $Q(x)$, the KL divergence is defined as

[TABLE]

and may be interpreted as the information gain in moving from a prior $Q$ to a posterior $P$ (Hosoya et al., 2004; Verde et al., 2013; Grandis et al., 2016; Raveri et al., 2016). For our reconstructions, we compute the KL divergence from prior to posterior for each distribution $P(f|k)$ conditioned on $k$. An example of such a plot can be found in the lower-right panel of Sec. III.

Throughout this work, plots use an approximate correspondence between wavenumber $k$ and multipole moment $\ell$ via the Limber approximation $\ell\approx k/D_{\mathrm{A}}$, where $D_{\mathrm{A}}=r_{*}/\theta_{*}$ is the Planck 2018 best-fit comoving angular distance to recombination at $r_{*}$.

II.4 Non-parametric reconstructions

Throughout this work, we explore various non-parametric functional forms for either the primordial power spectrum or the inflationary potential. “Non-parametric” is a slightly misleading terminology, as in general such reconstructions choose a function with a very large number of additional parameters. We prefer the terminology free-form Planck Collaboration (2018c); Olamaie et al. (2018), flexible Millea and Bouchet (2018) or adaptive Hee et al. (2017, 2016). The principle behind this is that the parameterization should have enough freedom to reconstruct any reasonable underlying function, independent of any underlying physical model.

For example, in this paper we work with variations on the linear spline, defined by parameters $\Theta_{f}$, producing a mapping from the independent variable $x$ to the dependent variable $y$ thus

[TABLE]

Here we have used a compact notation for denoting piecewise functions espoused by Graham et al. (1994) whereby $[R]$ is a logical truth function, yielding $1$ if the relation $R$ is true, and [math] if false. For consistency, we interpret the case $N=1$ as having a constant value of $y_{1}$ for all $x$.

In a Bayesian approach, one treats the additional degrees of freedom $\Theta_{f}$ of the non-parametric function as parameters in a posterior distribution, which one marginalizes out in order to obtain model-independent reconstructions. Typically there is a degree of choice as to how many parameters $N$ to use, and a penalty is applied for larger $N$ to avoid over-parameterization and noise fitting. In this work we treat $N$ in a Bayesian sense as well. Each reconstruction with a given number of parameters $N$ is treated as an independent model. We can then marginalize over the number of models using the Bayesian evidence.

II.5 Planck data and cosmology

In this paper in almost all cases we focus our efforts on using the pure Planck 2018 polarization data baseline, referred to in Planck Collaboration (2018a, 2019, b, c) as TT,TE,EE+lowE+lensing. Throughout we use a flat cold-dark-matter with dark energy ($\Lambda$CDM) late-time cosmology; therefore there are four associated cosmological parameters which in the default CosmoMC basis take the form

[TABLE]

Additionally there are $21$ nuisance parameters associated with Galactic foregrounds and the Planck instrumentation:

[TABLE]

We use the Planck PolyChord CosmoMC defaults as prior widths for all of these (indicated in Tab. 1), but as they are common to all models considered and sufficiently wide to encompass the entire posterior bulk, any prior effects from these parameters have been shown theoretically Trotta (2007) and in practice Handley and Lemos (2019) to cancel out.

II.6 Sampling strategy

Throughout this paper, we sample over the full parameter space $\Theta=(\Theta_{\mathrm{f}},\Theta_{\mathrm{c}},\Theta_{\mathrm{n}})$ of reconstruction, cosmological and nuisance parameters. The resulting posteriors are in general multimodal with complicated degeneracies between many parameters, particularly when there are a large number of reconstruction parameters $\Theta_{\mathrm{f}}$. In all cases, plots in this paper have the unmentioned parameters implicitly marginalized out. Marginalizing the likelihood over the prior in order to compute evidences is in general even more challenging. Nested sampling is ideally suited to performing such tasks, with PolyChord providing the cutting-edge of such technology in a cosmological context Handley et al. (2015a, b), proving to be essential for sampling over these complicated parameter spaces with up to $\sim\mathcal{O}(50)$ dimensions with fast-slow parameter hierarchies.

III Primordial power spectrum reconstruction

Traditionally in a $\Lambda$CDM cosmology, the primordial power spectrum $\mathcal{P}_{\mathcal{R}}(k)$ is modeled by a two-parameter function with an amplitude $A_{\mathrm{s}}$ and spectral index $n_{\mathrm{s}}-1$

[TABLE]

i.e.  a straight line in the $(\ln k,\ln\mathcal{P})$ plane. The tensor spectrum $\mathcal{P}_{\mathcal{T}}(k)$ may be parameterized by its own independent amplitude $A_{\mathrm{t}}$ and index $n_{\mathrm{t}}$, or via the tensor to scalar ratio $r=A_{\mathrm{t}}/A_{\mathrm{s}}$ and slow roll inflation consistency condition $n_{\mathrm{t}}=-r/8$ Baumann (2009). Extensions to parameterization (17) can be made by adding quadratic (running) and cubic (running of running) terms, but no evidence is found that these are required to describe the primordial power spectrum in the $k$-window which Planck probes.

Extending Eq. 17 with runnings of the spectral index creates a stiff parameterization, with no ability to account for sharper features, or large deviations at low or high-$k$. For our first primordial power spectrum reconstruction, we therefore parameterize as a logarithmic spline

[TABLE]

This represents a spline (Eq. 14) that is linear in the $(\ln k,\ln\mathcal{P})$ plane as shown in Fig. 1. For the tensor power spectrum we have analyzed the cases where $r$ is allowed to vary as a parameter, and also when $\mathcal{P}_{\mathcal{T}}(k)$ is given its own independent linear spline. Unsurprisingly, given that Planck measured an $r$ consistent with [math], the addition of a tensor power spectrum makes no difference to the scalar reconstructions. For simplicity we assume $r=0$ for the remainder of this section.

This technique has a history of being successfully applied to the primordial power spectrum (Vázquez et al., 2012a; Aslanyan et al., 2014; CORE Collaboration, 2018; Planck Collaboration, 2014, 2016), but has also been applied to dark energy equation of state by Hee et al. (2017) and Vázquez et al. (2012b), to the cosmic reionisation history by Millea and Bouchet (2018) and to galaxy cluster profiles by Olamaie et al. (2018). Our work differs from previous primordial power spectrum reconstructions in both the data we use, the styling of the priors, and in the application of more modern inference tools such as functional posterior plotting Handley (2018), conditional Kullback-Leibler divergences Hee et al. (2017) and PolyChord. The technique was applied to primordial power spectrum reconstruction from CORE simulated data in Section 6 of CORE Collaboration (2018), where it was shown that this approach accurately reconstructs complicated injected features (or lack thereof).

Priors

For priors on the vertical spline location parameters, we choose them to be independently uniform in $2<\log 10^{10}\mathcal{P}<4$. This spans an almost maximally wide range, increasing their width further has little effect due to CosmoMC discarding unphysically normalized spectra.

For priors on the horizontal spline location parameters, we choose the outermost knots to be fixed at $10^{-4}$ and $10^{-0.3}$. This corresponds roughly to a $C_{\ell}$ multipole range $1\lesssim\ell\lesssim 7000$, which fully encompasses the CMB window that Planck observes. For the remaining horizontal knots, we choose a prior which distributes the parameters logarithmically within this range, such that $k_{2}<\cdots<k_{N-1}$. This sorting procedure breaks the $(N-2)!$ implicit switching degeneracy, and is also termed a forced identifiability prior Handley et al. (2015b); Hall et al. (2018); Buscicchio et al. (2019).

To implement this sorted prior in the context of nested sampling, we need to define the transformation from the unit hypercube to the physical space. Coordinates in the unit hypercube $x_{1},\cdots,x_{N}$, can be transformed to coordinates in the physical space $\theta_{1},\cdots,\theta_{N}$, such that they are distributed uniformly in $[\theta_{\min{}},\theta_{\max{}}]$ and sorted so that $\theta_{1}<\cdots<\theta_{N}$ via the following reversed recurrence relation

[TABLE]

which is equivalent to saying that $\theta_{n}$ is marginally distributed as the largest of $n$ uniformly distributed variables within $[\theta_{\min{}},\theta_{n+1}]$. Another method for breaking the switching degeneracy is to exclude the region of the parameter space which does not satisfy the sorting criterion. This becomes exponentially small as more knots are added, which makes the initial sampling from the prior more challenging. It is also more in keeping with the nested sampling methodology to explicitly transform the full hypercube onto the space of interest. In our case, given that $k_{1}$ and $k_{N}$ are fixed, for $N\geq 4$, we sort the $N-2$ inner logarithmic coordinates $\log_{10}k_{2}<\cdots<\log_{10}k_{N-1}$.

We perform the reconstruction for $N=1,\cdots,9$ and then marginalize using Bayesian evidences with an implicit equal weighting for each $N$. This is equivalent to sampling from a full joint posterior with a uniform prior on $N$, and could alternatively be accomplished using the method described in Hee et al. (2016). Our priors on the reconstruction parameters are summarized in Tab. 2.

Results

We show results for our primordial power spectrum reconstruction using Planck 2018 TT,TE,EE+lowE+lensing data in Secs. III, III and III.

Sections III and III show the prior and posterior conditioned on each value of $N$. The case $N=1$ corresponds to a scale-invariant spectrum, whilst $N=2$ is equivalent (up to a small difference in prior) to the standard $\Lambda$CDM parameterization. As further knots are added, the reconstruction accounts for cosmic variance at low-$k$, and instrument noise at high-$k$. Furthermore, for large $N$, in a fraction of the samples there is a visibly clear oscillation characterized by a rise in power at $\ell\sim 50$, and a dip in power at $20<\ell<30$, as well as an overall suppression of power at low-$k$. For lower values of $\ell$, cosmic variance sets in, and few conclusions can be drawn from sampling differences between runs of different $N$ at these values.

To determine the statistical significance of these features, one should consider the Bayesian evidence, as indicated in the lower-left panel of Sec. III.

The first observation from Sec. III is that the $N=1$ scale-invariant power spectrum is completely ruled out, with a logarithmic difference of $\ln\mathcal{B}_{N=1}^{N=2}\sim\mathcal{O}(33)$. This represents overwhelming evidence for a tilted power spectrum, one of the key predictions of the theory inflation. A gambler could get odds of a quintillion to one against scale-invariance vs. $\Lambda$CDM.

The second observation is that the evidence for $N=3$ is greater than $N=2$, namely a model that is able to account for cosmic variance at low-$k$ and instrument noise at high-$k$ is, in a Bayesian sense, preferred to the simpler $\Lambda$CDM parameterization. Up until Planck 2018, the data had not been quite powerful enough for us to define the window that we observe in the primordial power spectrum in this Bayesian sense.

The third observation is that whilst $N=3$ is maximal in evidence, in fact $N=4,5,6$ are competitive, and $N=7,8,9$ are far from ruled out. With this lack of knowledge, the correct Bayesian approach is to marginalize over all models, using the Bayesian evidence as the relative weighting. Doing so, we can compute the marginalized spectrum and KL divergence as shown in Sec. III. We find that the observable window is now clearly defined, and hints of the low-$k$ features survive this marginalization.

Historical context

Section III shows reconstructions using the same methodology 111For the historical data prior to WMAP, we needed to significantly widen the priors on cosmological late-time parameters but now on data from a historical sequence of CMB experiments

COBE Bennett et al. (1996), 2. 2.

“pre-WMAP” (COBE Bennett et al. (1996), BOOMERANG Netterfield et al. (2002), MAXIMA Hanany et al. (2000), DASI Halverson et al. (2002), VSA Scott et al. (2003) and CBI Pearson et al. (2003)), 3. 3.

WMAP Bennett et al. (2013); Hinshaw et al. (2013), 4. 4.

Planck 2013 (TT+lowlike+lensing) Planck Collaboration et al. (2014) , 5. 5.

Planck 2015 (TT+lowTEB+lensing) Planck Collaboration et al. (2016) , 6. 6.

Planck 2018 (TT,TE,EE+lowE+lensing) Planck Collaboration (2018a).

The hint of an $20<\ell<30$ feature becomes visible after WMAP, is strengthened in switching to Planck, and remains stable as the Planck data are updated.

Examining the KL divergences in Sec. III, the vertical axis shows a greater overall constraint on the primordial power spectrum as improved cosmological constraints are obtained, and the horizontal axis shows that the $k$-window increases as the angular resolution of the experiments increases. The only alteration to this trend is the change from Planck 2013 to Planck 2015. In this case, the constraint on the primordial power spectrum is actually lowered, whilst the $k$-window increases. This is due to the fact that the $\tau$ constraint widened from 2013 to 2015, as can be seen it the top right panels of the Sec. III.

Planck 2018 provides the best constraints on the primordial power spectrum, both via its high-accuracy measurement of $\tau$, and in its small-scale angular resolution.

IV Inflationary Potential reconstruction

In contrast to the analysis from the previous section, instead of parameterising the primordial power spectrum directly, here we take the pipeline one stage backward and perform a non-parametric reconstruction of the inflaton potential $V(\phi)$. The scalar and tensor primordial power spectra $\mathcal{P}_{\mathcal{R},\mathcal{T}}$ are then derived from $V(\phi)$ via the procedure indicated in Section II.1.

To reconstruct the inflationary potential $V(\phi)$, it is more appropriate to work first with $\ln V$, as in general $V(\phi)$ can a-priori span a many scales, and it is typically $\frac{\operatorname{d}\!{}}{\operatorname{d}\!{\phi}}\ln V$ that drives much of the evolution of the inflaton during inflation.

More importantly, one cannot parameterize $V(\phi)$ via a linear interpolating spline as was the case in Sec. III. The equations of motion (Eqs. 1, 2, 3 and 4) in general depend on first (and sometimes second) derivatives of $V(\phi)$. Parameterizing the potential using a linear spline will typically yield primordial power spectra with (arguably) unphysical ringing effects.

One should therefore use a spline with continuous first derivatives, and it is natural to choose a cubic spline as the conceptually simplest smooth interpolator. It is tempting to try to do this directly by taking the locations of the knots of the spline as free parameters. Cubic splines, however, are stiff, yielding complicated posteriors that are very difficult to navigate, interpret, and set priors on Vázquez et al. (2012a).

Cubic splines have the property that their derivative is a smooth piecewise quadratic, and their second derivative is a piecewise linear spline. This suggests that the cleanest way to reconstruct the potential is to parameterize the second derivative as a linear spline, and then integrate this function twice to get the log-potential. Two additional parameters are created by this double integration, a gradient term $\frac{\operatorname{d}\!{\ln V_{*}}}{\operatorname{d}\!{\phi}}$ and an overall offset $\ln V_{*}$. These two free parameters function as an alternative constraint choice in comparison with natural or clamped splines. Our reconstruction function is therefore

[TABLE]

which can be viewed as working through Fig. 6 in reverse.

Priors

The priors in this analysis proved to be critically important for recovering sensible results. To harmonize with the analysis of the primordial power spectrum, our first requirement is that any primordial power spectrum generated from a potential $V(\phi)$ resides in the range ${2<\ln 10^{10}\mathcal{P}_{\mathcal{R}}(k)<4}$.

Consider the slow roll parameters Liddle et al. (1994), and their relation to the second derivative of the log potential

[TABLE]

We therefore take the priors on the second derivatives of the log potential $\frac{\operatorname{d}^{2}\!{\ln V_{i}}}{\operatorname{d}\!{\phi}^{2}}$ to be uniformly distributed, and the gradient $\frac{\operatorname{d}\!{\ln V_{*}}}{\operatorname{d}\!{\phi}}$ is taken to be negatively log-uniform. Negativity forces the inflaton to roll downhill from negative to positive $\phi$, breaking a symmetric degeneracy. We take the potential offset to vary across a wide range $\ln V_{*}$. Widening any of these priors detailed in Tab. 3 further has no effect, as any primordial power spectrum generated outside these bounds lies outside the range $[2,4]$.

Particular care must be taken with the horizontal locations of the knots. Any reconstruction of the potential will be sensitive only to the observable window of inflation $[\phi_{\min{}},\phi_{\max{}}]$, defined as when the largest and smallest observable scales $k_{\min{}}$ and $k_{\max{}}$ exit the horizon. As in Sec. III, we take $(k_{\min{}},k_{\max{}})=(10^{-4},10^{-0.3})\mathrm{Mpc}^{-1}$.

Unfortunately, the bounds of the window $[\phi_{\min{}},\phi_{\max{}}]$ are strongly dependent on the other primordial parameters. One cannot therefore take an arbitrarily wide range in $\phi$ for the horizontal locations, as the reconstruction is then dominated by the prior effect of unconstrained knot parameters. The locations $\phi_{1},\ldots,\phi_{N}$ of the reconstruction knots should instead be distributed throughout the observable window. Whilst the locations $\phi_{1},\ldots\phi_{N}$ and heights $\frac{\operatorname{d}^{2}\!{\ln V_{1}}}{\operatorname{d}\!{\phi}^{2}},\ldots,\frac{\operatorname{d}^{2}\!{\ln V_{N}}}{\operatorname{d}\!{\phi}^{2}}$ themselves influence the size of the observable window, a reasonable approach is to first estimate it using the unperturbed potential (i.e. setting $N=0$), giving an alternative window $[\tilde{\phi}_{\min{}},\tilde{\phi}_{\max{}}]$. In a similar manner as the horizontal knots in Sec. III, we take the $N$ horizontal $\phi$-knot locations to be sorted and uniform throughout this window.

Finally, we require that the inflaton should evolve in an inflating phase throughout the observable window, and that it should be rolling (not necessarily in slow roll) downhill from negative to positive $\phi$ throughout. These priors are summarized in Tab. 3.

Alternative methodologies for direct reconstruction of the potential exist in the literature. One approach is to to expand the potential $V(\phi)$ as a Taylor series (Lesgourgues and Valkenburg, 2007; Mortonson et al., 2011). Another is to expand $H(\phi)$ as a Taylor series (Lesgourgues et al., 2008; Noreña et al., 2012), and then derive the potential analytically via ${V(\phi)=3M_{\mathrm{p}}^{2}H^{2}-2M_{\mathrm{p}}^{4}{H^{\prime}}^{2}}$. Both of these approaches have been successfully applied in the Planck inflation papers Planck Collaboration (2014, 2016, 2018c).

Results

Due to the strong dependency of the $\phi$-window on the potential itself, it is not particularly illuminating to plot $V(\phi)$ directly. Instead, in the spirit of the other two sections we start by plotting the functional posterior of the primordial power spectrum $\mathcal{P}_{\mathcal{R}}(k)$, shown in Secs. IV and IV. Viewed in this manner, one can think of these primordial power spectrum reconstructions as having an alternative prior complementary to Sec. III, motivated by the assumption that the primordial power spectrum is derived from a smooth underlying potential.

In the same manner as Sec. III, Sec. IV and is consequently a form of exponential potential. Regardless, it recovers a primordial power spectrum with an appropriate amplitude and tilt and minimal running, almost identical to the traditional $A_{\mathrm{s}},n_{\mathrm{s}}$ parameterization. $N=1$ adds a constant second derivative term to the Taylor expansion, and produces a similar primordial power spectrum. As more knots are added, the potential has greater freedom, and the corresponding primordial power spectrum begins to gain similar features to the results in Sec. III; a loss of constraint at low and high-$k$. Intriguingly, there is also the same preference for an oscillation with a peak at $\ell\sim 50$ and trough at $20<\ell<30$. In the line plots of Sec. IV the oscillation is now smooth, on account of the physical potential-based prior created by this reconstruction. It should be noted that while such oscillations are characteristic of this integrated inflationary potential parameterization, these were also partially recovered a priori in the free-form approach from Sec. III.

Examining the evidences in Sec. IV, we can see that despite the similarities in primordial power spectra $N=1$ is preferred over $N=0$. The reason for this is that the restrictive form of potential for $N=0$ forces $r\approx 0.2$, which is now ruled out by Planck. Allowing a second derivative for the $N=1$ relaxes the $r$ constraint, resulting in a Bayesian preference for the $N=1$ case, consistent with the results of the Planck Collaboration (2018c). Adding further knots causes the evidence to drop, indicating that from a Bayesian standpoint, no further complexity is required by the data. The marginalized plots in Sec. IV show similar attributes to those of the primordial power spectrum in Sec. III, but in this case the stiffness of the primordial power spectrum reconstruction results in a slightly poorer recovery of the relative lack of power spectrum constraint at low and high-$k$.

Our second set of plots detail results for the inflationary slow roll parameter $\eta_{V}(k)$, shown in Secs. IV and IV. Instead of using $\phi$ as the independent variable (which suffers from the dependency of the window widths on the underlying potential), we use the effective wavenumber $k(\phi)$, which is monotonically related to $\phi$ in our reconstruction, defined to be the size of the comoving Hubble radius at that moment in the field’s evolution.

Section IV reveals that the oscillations in the primordial power spectrum at low $k$ are created by a partial breakdown in the slow roll conditions. For $N\geq 5$, $\eta_{V}\sim 0.5$, which are the same values of $N$ at which oscillations become apparent in Sec. IV. This will be of particular interest for just enough inflation models (Ramirez and Schwarz, 2012; Handley et al., 2014; Hergt et al., 2019a, b), models with singularities and discontinuities Starobinskij (1992); Adams et al. (2001); Joy et al. (2008), multi-field phase-transitions Silk and Turner (1987); Holman et al. (1991); Polarski and Starobinsky (1992); Adams et al. (1997); Hunt and Sarkar (2004), M-theory Lesgourgues (2000); Burgess et al. (2005) or supergravity Ashoorioon and Krause (2006) models, to name a few examples. There is a long history of confronting such models with data Peiris et al. (2003); Covi et al. (2006); Hamann et al. (2007); Hunt and Sarkar (2007); Joy et al. (2009); Martin and Ringeval (2004); Kawasaki et al. (2005); Jain et al. (2009). Section IV details the marginalized results.

Planck provides only a weak upper bound on the other slow roll parameter $\varepsilon_{V}\approx\frac{r}{16}$, meaning $\varepsilon_{V}$ is nearly indistinguishable from its logarithmic prior. Given the slow roll relations, we find that for our reconstructions $\eta_{V}\approx\frac{\operatorname{d}^{2}\!{\log V}}{\operatorname{d}\!{\phi}^{2}}$, so plots of $\eta_{V}$ are nearly identical to plots of the equivalent underlying linear second-derivative reconstruction parameters.

V Sharp feature reconstruction

In this section, we return to a direct analysis of the primordial power spectrum. Inspired by the oscillatory features present in both the linear spline primordial power spectrum reconstruction (Sec. III) and in the functional posterior of the primordial power spectrum from the inflationary potential reconstruction (Sec. IV), we now consider a parameterization of the primordial power spectrum which favors sharp features. Attempts at explaining these features has a long history in the literature, initially investigated in Refs. Hinshaw et al. (2003); Spergel et al. (2003); Peiris et al. (2003); Mortonson et al. (2009).

We introduce sharp features into the parameterization of the spectrum by placing a variable number $N$ of top-hat functions with varying widths $\Delta$, heights $h$, and locations $\log_{10}k$ on top of the traditional $A_{\mathrm{s}},n_{\mathrm{s}}$ parameterization

[TABLE]

where the square brackets (as in Section II.4) denote a logical truth function (Graham et al., 1994). This parameterisation is indicated schematically in Fig. 11, and Tab. 4 provides a summary of the priors that we use. Readers are referred to the priors sections of Secs. III and IV for further details.

Results

The results for the sharp feature reconstructions can be found in Secs. V, V and V. Section V shows that in the marginalized plots the oscillations and lack of power spectrum constraints at low- and high-$k$ are once again recovered. Visually, the features are even more striking in these reconstructions, on account of the ability for this parameterization to localize in $k$ more precisely. There are also hints of features at high-$k$ in this reconstruction, which were smoothed out by the parameterizations of the two previous sections.

Marginalization in Sec. V shows that there is little Bayesian evidence to support the introduction of more than two features, but the low-$k$ oscillation still comes through clearly in the fully marginalized plot.

Finally, we examine the effects of these reconstructions on the $C_{\ell}$ spectra by considering the functional posterior in Sec. V. By comparing the $\Lambda$CDM case $N=0$ with the case $N=8$ we see that the features at low-$k$ in the PPS correspond to both a suppression of power at low-$\ell$ in the TT spectrum, as well as a more specific localized reduction of power in the $20<\ell<30$ region. There seems to be no obvious correspondence with possible features seen in the polarization TE or EE spectra.

VI Cosmological parameter stability

Finally in Sec. V we show that the underlying cosmological parameter constraints remain effectively unchanged for all three of the analyses in Secs. III, IV and V, in spite of the additional degrees of freedom we have given to the primordial power spectrum.

The only exception is the PPS $N=1$ reconstruction. This model is highly-disfavored (Sec. III), since Planck rules out a Harrison-Zeldovich scale-invariant ($n_{s}=1$) spectrum Planck Collaboration (2014); Planck Collaboration et al. (2016); Planck Collaboration (2018c). Requiring $n_{s}=1$ gives a poorly-fit model which forces the cosmological parameters into locations discordant with $\Lambda$CDM.

This overall parameter stability for models consistent with the data demonstrates that one can explain features in CMB power spectra via modifications to the primordial cosmology, without the need to alter late-time cosmological parameters.

VII Conclusions

In this work, we have reconstructed the primordial Universe three ways. In Sec. III, we reconstructed the primordial power spectrum using a linear spline. In Sec. IV, we reconstructed the inflationary potential using a cubic spline. In Sec. V, we probed sharp features in the primordial power spectrum by superimposing top-hat functions on top of the traditional $\Lambda$CDM power spectrum.

We showed that the Bayesian odds of a scale-invariant power spectrum are around a quintillion to one against in comparison to the $\Lambda$CDM cosmology. This agrees with the Planck Collaboration’s conclusions Planck Collaboration (2018c) that there is decisive evidence for $n_{\mathrm{s}}\neq 1$ — one of the key predictions of the theory of inflation.

All methods reconstruct a featureless tilted power law consistent with a simple $(A_{\mathrm{s}},n_{\mathrm{s}})$ parameterization across a broad observable window $(50\lesssim\ell\lesssim 2000)$. In addition, all reconstructions demonstrate that in a Bayesian sense it is preferable to have models which are able to recover the lack of power spectrum constraints at low-$k$ due to cosmic variance, and at high-$k$ due to Planck instrument noise, reflected in the evidences and marginalized plots (Secs. III, IV and V).

All large $N$ conditional reconstructions partially recover oscillatory features in the primordial power spectrum, with a peak at $\ell\sim 50$ and a trough at $20<\ell<30$, which manifest themselves in the functional posteriors of the $C_{\ell}$ spectra (Sec. V). The inflationary potential reconstruction (Sec. IV) shows that this oscillation could be due to a breakdown in slow roll near the start of the inflationary window (Sec. IV), which is relevant for a wide variety of inflationary models (Ramirez and Schwarz, 2012; Handley et al., 2014; Hergt et al., 2019a, b; Starobinskij, 1992; Adams et al., 2001; Joy et al., 2008; Silk and Turner, 1987; Holman et al., 1991; Polarski and Starobinsky, 1992; Adams et al., 1997; Hunt and Sarkar, 2004; Lesgourgues, 2000; Burgess et al., 2005; Ashoorioon and Krause, 2006; Peiris et al., 2003; Covi et al., 2006; Hamann et al., 2007; Hunt and Sarkar, 2007; Joy et al., 2009; Martin and Ringeval, 2004; Kawasaki et al., 2005; Jain et al., 2009). However, the oscillations do not survive marginalization over $N$, indicating that the Bayesian evidence is not strong enough from Planck data to indicate a significant detection of such a feature.

The renewed upper bound on $r$ from Planck 2018 now has enough discriminatory power to begin reconstructing potentials, as shown by the preference for the $N=1$ case in the inflationary potential reconstructions.

As shown in Sec. V, in all cases, the distributions on the late-time cosmological parameters remain unperturbed by the additional degrees of freedom on the primordial cosmology provided by these reconstructions, indicating that any conclusions using late-time parameters are unlikely to be affected by modifying the primordial cosmology.

There is scope for inflationary models which a-priori predict these low-$k$ features to be preferred over the $\Lambda$CDM cosmology, particularly if such models are capable of producing sharper features in the $C_{\ell}$ spectra at $20<\ell<30$. Additionally, in light of further CMB data CORE Collaboration (2018), or failing that, strong $\tau$ characterization, it is likely that these hints of features will sharpen and provide further discriminatory power in constructing better models of the primordial Universe.

Acknowledgements.

WJH was supported via STFC PhD studentship RG68795, the European Research Council (ERC) under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 306478-CosmicDawn, and a Gonville & Caius College Research Fellowship. WJH and ANL thank Fabio Finelli, Martin Bucher and Julien Lesgourgues for numerous useful comments and suggestions for improvements over the course of this work. HVP acknowledges useful conversations with George Efstathiou, Matt Johnson and Keir Rogers. HVP was supported by the European Research Council (ERC) under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 306478-CosmicDawn. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation Grant PHY-1607611. This work was also partially supported by a grant from the Simons Foundation. This work was performed using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service, provided by Dell Inc. using Strategic Research Infrastructure Funding from the Higher Education Funding Council for England and funding from the Science and Technology Facilities Council, as well as resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service, provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council.