Accounting for Correlations When Fitting Extra Cosmological Parameters

TL;DR

This paper emphasizes the importance of considering correlations between additional cosmological parameters when fitting models, demonstrating how neglecting these correlations can lead to incomplete or misleading conclusions about model consistency.

Contribution

It introduces a method to quantify parameter correlations using simulations and MCMC, highlighting their significance in cosmological data analysis.

Findings

Correlations between parameters like $A_L$, $n_{run}$, and $Y_P$ can be substantial.

Accounting for correlations aligns Planck data with $ m extLambda$CDM expectations.

Neglecting correlations may obscure or misinterpret cosmological tensions.

Abstract

Current cosmological tensions motivate investigating extensions to the standard $\Lambda$CDM model. Additional model parameters are typically varied one or two at a time, in a series of separate tests. The purpose of this paper is to highlight that information is lost by not also examining the correlations between these additional parameters, which arise when their effects on model predictions are similar, even if the parameters are not varied simultaneously. We show how these correlations can be quantified with simulations and Markov Chain Monte Carlo (MCMC) methods. As an example, we assume that $\Lambda$CDM is the true underlying model, and calculate the correlations expected between the phenomenological lensing amplitude parameter, $A_L$, the running of the spectral index, $n_{\rm run}$, and the primordial helium mass fraction, $Y_P$, when these parameters are varied one at a time along with the $\Lambda$CDM parameters in fits to the $\textit{Planck}$ 2015 temperature power spectrum. These correlations are not small, ranging from 0.31 ($A_L-n_{\rm run}$) to $-0.93$ ($n_{\rm run}-Y_P$). We find that the values of these three parameters from the $\textit{Planck}$ data are consistent with $\Lambda$CDM expectations within $0.9\sigma$ when the correlations are accounted for. This does not explain the 1.8-2.7$\sigma$ $\textit{Planck}$ preference for $A_L>1$, but provides an additional $\Lambda$CDM consistency test. For example, if $A_L>1$ was a symptom of an underlying systematic error or some real but unknown physical effect that also produced spurious correlations with $n_{\rm run}$ or $Y_P$ our test might have revealed this. We recommend that future cosmological analyses examine correlations between additional model parameters in addition to investigating them separately, one a time.

Figures

Tables

Table 1. Table 1 The significance of the difference between the experimental ML values and the fiducial of A L subscript 𝐴 𝐿 A_{L} , n run subscript 𝑛 run n_{\rm{run}} and Y P subscript 𝑌 𝑃 Y_{P} , in terms of the χ 2 superscript 𝜒 2 \chi^{2} difference (with 3 degree of freedom) of difference, its PTE and Δ ​ σ Δ 𝜎 \Delta\sigma , the level of consistency in terms of the number of σ 𝜎 \sigma .

	MCMC	simulation
${\chi }^2$	4.94	4.78
PTE	0.18	0.19
$\Delta \sigma$	0.9	0.9

Table 2. Table 2 Parameter uncertainties and correlations for the extension parameters, estimated from the MCMC method and the simulations.

Parameter	MCMC	Simulations
parameter uncertainty
$A_L$	0.084	0.081
$n_{\mathrm{run}}$	0.0093	0.0093
$Y_P$	0.021	0.020
parameter correlation
$A_L$ v.s. $n_{\mathrm{run}}$	0.31	0.32
$A_L$ v.s. $Y_P$	-0.46	-0.47
$n_{\mathrm{run}}$ v.s. $Y_P$	-0.93	-0.93

Equations

$$\{{\Omega }_b{h}^2,{\Omega }_c{h}^2,100{\theta }_{MC},\tau ,\log A_s,n_s{\}}^{\mathrm{fid}}=$$

$$\{0.02216,0.1211,1.0407,0.07,3.0777,0.9601\}.$$

$${\chi }^2=({\mathbf{\theta }}_{\mathrm{dat}}^{\mathrm{ML}}-{\mathbf{\theta }}^{\mathrm{fid}}{)}^T{\mathbf{\Sigma }}_{\mathrm{ex}}^{-1}({\mathbf{\theta }}_{\mathrm{dat}}^{\mathrm{ML}}-{\mathbf{\theta }}^{\mathrm{fid}})$$

$$Cov(A,B{)}_{\mathrm{fitted}\mathrm{separately}}$$

$$\times (Var(A{)}_{\Lambda \mathrm{CDM}+A}$$

$$\log \mathbb{P}({\mathbf{d}}^{\mathrm{sim}}|{\mathbf{\theta }}^{\mathrm{fid}})=-\frac{1}{2}({\mathbf{d}}^{\mathrm{sim}}-\mathbf{\mu }({\mathbf{\theta }}^{\mathrm{fid}}){)}^T{\mathcal{C}}^{-1}({\mathbf{d}}^{\mathrm{sim}}-\mathbf{\mu }({\mathbf{\theta }}^{\mathrm{fid}})).$$

$$\log L(\mathbf{\theta })\equiv \log \mathbb{P}({\mathbf{d}}^{\mathrm{sim}}|\mathbf{\theta }).$$

$$\frac{\partial\log\mathbb{P}(\bm{d}^{\rm{sim}}|\bm{\theta})}{\partial\theta_{i}}\bigg{|}_{\bm{\theta}=\bm{\theta}^{\rm ML}}=0$$

$$\begin{split}\log&\mathbb{P}(\bm{d}^{\rm{sim}}|\bm{\theta})=\log\mathbb{P}(\bm{d}^{\rm{sim}}|\bm{\theta})\big{|}_{\bm{\theta}=\bm{\theta}^{\rm ML}}\\ &+\frac{1}{2}\sum_{ij}(\theta_{i}-\theta^{\rm ML}_{i})\left[\frac{\partial^{2}\log\mathbb{P}(\bm{d}^{\rm{sim}}|\bm{\theta})}{\partial\theta_{i}\partial\theta_{j}}\right]_{\bm{\theta}=\bm{\theta}^{\rm ML}}(\theta_{j}-\theta^{\rm ML}_{j})\\ &+\mathcal{O}(\bm{\theta}-\bm{\theta}^{\rm ML})^{3}.\end{split}$$

$$\mathbf{\Sigma }({\mathbf{\theta }}^{\mathrm{ML}}({\mathbf{d}}^{\mathrm{sim}}){)}_{ij}$$

$$\log \mathbb{P}({\mathbf{\theta }}^{\mathrm{ML}}({\mathbf{d}}^{\mathrm{sim}})|{\mathbf{\theta }}^{\mathrm{fid}})=$$

$$-\frac{1}{2}({\mathbf{\theta }}^{\mathrm{ML}}({\mathbf{d}}^{\mathrm{sim}})-{\mathbf{\theta }}^{\mathrm{fid}}{)}^T{\mathbf{\Sigma }}^{-1}({\mathbf{\theta }}^{\mathrm{ML}}({\mathbf{d}}^{\mathrm{sim}})-{\mathbf{\theta }}^{\mathrm{fid}})$$

$$\mathbb{P}(\mathbf{\theta }|\mathbf{d})=\frac{\mathbb{P}(\mathbf{d}|\mathbf{\theta })\mathbb{P}(\mathbf{\theta })}{\mathbb{P}(\mathbf{d})}\propto \mathbb{P}(\mathbf{d}|\mathbf{\theta })\mathbb{P}(\mathbf{\theta }),$$

$$\log \mathbb{P}(\mathbf{\theta }|\mathbf{d})=-\frac{1}{2}(\mathbf{\theta }-\bar{\mathbf{\theta }}{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{\theta }-\bar{\mathbf{\theta }}),$$

$$\mathbb{P}(\mathbf{\theta }|\mathbf{d}=\mathbf{\mu }({\mathbf{\theta }}^{\mathrm{fid}}))$$

$$\mathfrak{L}=-\frac{1}{2}(\mathbf{d}-\mathbf{\mu }(\mathbf{\theta }){)}^T{\mathcal{C}}^{-1}(\mathbf{d}-\mathbf{\mu }(\mathbf{\theta }))$$

$$\displaystyle\frac{\partial\mathfrak{L}}{\partial{\theta}_{i}}\bigg{|}_{\mathrm{ML}}$$

$${\mathcal{F}}_{ij}$$

$$\mathbf{\mu }(\mathbf{\theta })\approx {\mathbf{\mu }}^{\mathrm{fid}}+\sum _{i=1}^n{\mathbf{\mu }}_{,i}^{\mathrm{fid}}({\theta }_i-{\theta }_i^{\mathrm{fid}}),$$

$${\mathbf{\Sigma }}_{ij}=({\mathcal{F}}^{-1}{)}_{ij}.$$

$$({\mathcal{M}}^{-1}{)}_{ij}=(-1{)}^{i+j}\frac{|{\mathcal{M}}_{\setminus j\setminus i}|}{|\mathcal{M}|}$$

$$|\mathcal{M}|=\sum _{i=1}^n(-1{)}^{i+j}{\mathcal{M}}_{ij}|{\mathcal{M}}_{\backslash i\backslash j}|$$

$${\mathbf{\Sigma }}_{ij}=(-1{)}^{i+j}|{\mathcal{F}}_{\backslash j\backslash i}|/|\mathcal{F}|.$$

$$corr(A,B)=\frac{{\mathbf{\Sigma }}_{AB}}{({\mathbf{\Sigma }}_{AA}{\mathbf{\Sigma }}_{BB}{)}^{\frac{1}{2}}}.$$

$$corr(A,B)=\frac{-|{\mathcal{F}}_{\backslash A\backslash B}|}{(|{\mathcal{F}}_{\backslash A\backslash A}||{\mathcal{F}}_{\backslash B\backslash B}|{)}^{\frac{1}{2}}}.$$

$${\mathbf{\theta }}^{\mathrm{ML}}(\mathbf{d})-{\mathbf{\theta }}^{\mathrm{fid}}={\mathcal{F}}^{-1}{\mathbf{\mu }}_{,i}^T{\mathcal{C}}^{-1}(\mathbf{d}-\mathbf{\mu }).$$

$${\mathbf{\theta }}^{\mathrm{ML}}(\mathbf{d})={\mathcal{F}}^{-1}{\mathbf{\mu }}_{,i}^T{\mathcal{C}}^{-1}\mathbf{d}+\hat{\mathbf{\theta }},$$

$$\hat{\mathbf{\theta }}\equiv -{\mathcal{F}}^{-1}{\mathbf{\mu }}_{,i}^T{\mathcal{C}}^{-1}{\mathbf{\mu }}^{\mathrm{fid}}+{\mathbf{\theta }}^{\mathrm{fid}}.$$

$$\left(\begin{array}[]{ccc}\bm{\theta}^{\rm{ML}}_{1}\\ \bm{\theta}^{\rm{ML}}_{2}\\ \vdots\\ \bm{\theta}^{\rm{ML}}_{n}\\ A^{\rm{ML}}\end{array}\right)=(\mathcal{F}_{\backslash B\backslash B})^{-1}\left(\begin{array}[]{ccc}\bm{\mu}^{T}_{,1}\bm{\mathcal{C}}^{-1}\bm{d}\\ \bm{\mu}^{T}_{,2}\bm{\mathcal{C}}^{-1}\bm{d}\\ \vdots\\ \bm{\mu}^{T}_{,n}\bm{\mathcal{C}}^{-1}\bm{d}\\ \bm{\mu}^{T}_{,A}\bm{\mathcal{C}}^{-1}\bm{d}\end{array}\right)+\hat{\bm{\theta}}.$$

$$A^{\mathrm{ML}}(\mathbf{d})=$$

$$\displaystyle\bigg{[}\sum^{n}_{i=1}(\mathcal{F}_{\backslash B\backslash B})^{-1}_{Ai}\bm{\mu}^{T}_{,i}\bm{\mathcal{C}}^{-1}+(\mathcal{F}_{\backslash B\backslash B})^{-1}_{AA}\bm{\mu}^{T}_{,A}\bm{\mathcal{C}}^{-1}\bigg{]}\bm{d}+\hat{{\theta}}_{A}.$$

$$B^{\mathrm{ML}}(\mathbf{d})=$$

Full text

\pacs

95.36.+x, 98.80.Es, 98.80.-k

1. Introduction

Over the last decade, much progress has been made on putting precise constraints on cosmological parameters within the $\Lambda$ cold dark matter ($\Lambda\mathrm{CDM}$) model, particularly from Cosmic Microwave Background (CMB) experiments (e.g., Bennett et al., 2013; Planck Collaboration XIII, 2016; Sievers et al., 2013; Story et al., 2013). In the most recent release of $\it{Planck}$ results (Planck Collaboration VI, 2018), determinations of the standard $\Lambda\mathrm{CDM}$ parameters such as baryon density, Hubble constant and matter density have reached the percent level or below.

Although currently there is no convincing evidence for deviations from the standard $\Lambda\mathrm{CDM}$ model from any single experiment, tensions exist between the values of some parameters inferred from different datasets. The most severe one is the 4.4$\sigma$ disagreement between the Hubble constant measurements from the anchor-Cepheid-supernova distance ladder by the SH0ES collaboration (Riess et al., 2019) and from the $\it{Planck}$ CMB data (Planck Collaboration VI, 2018). The measurement of $H_{0}$ via strong lensing time delays (Bonvin et al., 2017; Birrer et al., 2019) is consistent with the distance ladder and in 2.5$\sigma$ tension with $\it{Planck}$. Addison et al. (2018) showed that the tension between early and late time universe measurements persists even without the inclusion of $\it{Planck}$ data, using baryon acoustic oscillation (BAO) scale measurements.

Extensions or alternatives to the $\Lambda\mathrm{CDM}$ model have been explored in an attempt to resolve the Hubble tension. For example, the effects of varying the effective number of neutrino species (e.g., Riess et al., 2016) and the equation of state parameter of dark energy (e.g., Joudaki et al., 2017) have been studied, though these extensions have not been able to effectively relieve the tension without including multiple turning points in the evolution of the dark energy equation of state (Zhao et al., 2017). Recently, new ideas have been proposed as more promising solutions, for example, with the introduction of early dark energy (Poulin et al., 2018) or self-interaction massive neutrinos (Kreisch et al., 2019).

If there is new physics, consistency tests within the $\Lambda\mathrm{CDM}$ model will eventually fail with sufficiently sensitive new data.

An example of a parameter used to test consistency is the lensing amplitude $A_{L}$. It was first introduced by Calabrese et al. (2008) as a phenomenological way to quantify the effect of weak gravitational lensing in the CMB power spectrum. By definition, $A_{L}=1$ is the physical value. However, the $\it{Planck}$ temperature power spectrum (TT) data have shown a persistent preference for $A_{L}>1$ at 1.8-2.7$\sigma$ depending on which datasets are included (Planck Collaboration XVI, 2014; Planck Collaboration XIII, 2016; Planck Collaboration VI, 2018). For discussion of the $A_{L}$ tension, see also Addison et al. (2016); Motloch & Hu (2018, 2019a). The cause of the deviation of $A_{L}$ from its physical value is unclear, however, varying $A_{L}$ is an example of the sort of test that might ultimately shed light on the origin of the distance ladder tension.

Typically results from fitting additional model parameters are presented one or two at a time, a series of separate tests (for recent examples, see Heavens et al., 2017; Joudaki et al., 2017; Planck Collaboration VI, 2018). Different parameters may have similar effects on the theory prediction so the constraints on extension parameters, even when fitted separately, from an experiment or particular combination of experiments may be correlated. As an example, Figure 1 shows the expected 2-D confidence contours for two parameters, $A$ and $B$, obtained from \Lambda\mathrm{CDM}$$+A and \Lambda\mathrm{CDM}$$+B fits. The green contours show the expected distribution (estimated from e.g. simulations) if $\Lambda\mathrm{CDM}$, with $A=0$ and $B=0$, is the true model. Only looking at their 1-D marginalized distributions would lead to an incorrect conclusion of consistency, as the values of $A$ and $B$ are each only 1 standard deviation away from their fiducial $\Lambda\mathrm{CDM}$ values, while in the 2-D space the shift from the fiducial point is orthogonal to the degeneracy direction and the values are actually outside the 99.7% contour. This result in the 2-D space would not mean that there is strong evidence for the $\Lambda\mathrm{CDM}$+$A$ or the $\Lambda\mathrm{CDM}$+$B$ model, since neither $A$ or $B$ disagrees with their fiducial values significantly in 1-D. Rather, it suggests that there is something else wrong with the assumed model, or there may be experimental systematics we have not accounted for.

Therefore, to more carefully assess whether the standard $\Lambda\mathrm{CDM}$ model can consistently describe data, the covariance of the set of extension parameters should be quantified.

In this paper, our goal is to answer the following questions:

How do we calculate the expected correlation between different extension parameters to the standard $\Lambda\mathrm{CDM}$ model when constrained by the same data? 2. 2.

How do we incorporate these correlations into a more stringent test of the $\Lambda\mathrm{CDM}$ model?

As an example, we use the $\it{Planck}$ 2015 temperature power spectrum likelihood code (Planck Collaboration XI, 2016). The outline of this paper is as follows. In Section 2 we present the theoretical basis of our work and two different methods to achieve our goals. We show results in Section 3, followed by a discussion in Section 4 of general recipes for similar analysis in the future and conclusions in Section 5.

2. Methodology

2.1. Estimating correlation between extension parameters using simulations

In general it is not straightforward to infer the correlation between parameters $A$ and $B$ directly from MCMC fits to the \Lambda\mathrm{CDM}$$+A and \Lambda\mathrm{CDM}$$+B models that have already been performed by experimental collaborations like $\it{Planck}$. To make progress, we assume that the data likelihood is Gaussian, with a covariance matrix $\bm{\mathcal{C}}$ that does not depend on the cosmological parameters, and that the posterior distribution for the cosmological parameters is also Gaussian.

The assumption of Gaussianity of the likelihood is made widely across different cosmological measurements, including CMB power spectra (e.g., Planck Collaboration XI, 2016; Louis et al., 2017; Henning et al., 2018), weak lensing shear (e.g., Krause et al., 2017; Hikage et al., 2019; Wright et al., 2018), galaxy clustering (e.g., Percival et al., 2014; Alam et al., 2017), supernovae distance moduli (e.g., Scolnic et al., 2018; DES Collaboration, 2018), and others. The likelihood can almost always be made more Gaussian through compression of the data, often with negligible loss of cosmological information (e.g., combining CMB power spectra over a range of multipoles into bins). Neglecting the cosmological parameter dependence in the covariance (e.g., assuming a fixed fiducial model and set of parameters for computing the cosmic variance contribution to the errors) has also been demonstrated to be a suitable approximation for $\it{Planck}$, and other current experiments (e.g., Hamimeche & Lewis, 2008; Krause et al., 2017). The fiducial model as well as the covariance are often estimated from some iterative process of fitting the actual data.

For the $\it{Planck}$ data, cosmological parameter posterior distributions can be well approximated as Gaussian for $\Lambda\mathrm{CDM}$, as well as for many one- and two-parameter extensions, although there are also cases (e.g., involving curvature or varying the dark energy equation of state), where $\it{Planck}$ alone does not provide Gaussian posteriors. Other experiments only constrain a portion of the $\Lambda\mathrm{CDM}$ parameter space sufficiently well to produce Gaussian constraints in one or two parameters. Many parameters are also physically required to be positive, truncating the available parameter space and causing departure from Gaussianity when the data do not constrain the parameter to significantly differ from zero. This is the case for current cosmological constraints on the neutrino mass, for example (see Figure 30 of Planck Collaboration XIII, 2016). We return to the handling of non-Gaussian cases in Section 4.3.

In certain cases, the Bayesian posterior parameter distribution and the distribution of the maximum likelihood (ML) parameter values from many realizations of the data are approximately equal. Specifically, $\mathbb{P}(\bm{\theta}|\bm{d})$, the Bayesian posterior distribution of parameters $\bm{\theta}$ sampled by MCMC given the experimental data, can well approximate $\mathbb{P}(\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})|\bm{\theta}^{\rm{fid}})$. The latter is the distribution of frequentist maximum-likelihood (ML) parameter estimation based on realizations of $\bm{d}^{\rm{sim}}$ from a fiducial model $\bm{\theta}^{\rm{fid}}$. The choice of $\bm{\theta}^{\rm{fid}}$ is usually physically motivated and is based on fits from actual data. For a detailed discussion, see e.g., Chapter 4 and Appendix B of Gelman et al. (2013).

We show in Appendix A that this correspondence is mathematically exact, if we: (1) assume the Gaussianity of both the data likelihood and the posterior distribution of parameters; (2) impose parameter priors that are far less constraining than the likelihood; (3) in the MCMC computation, replace the experimental data $\bm{d}$ with $\bm{\mu}(\bm{\theta}^{\rm{fid}})$, the theory prediction of the fiducial model. In other words, $\mathbb{P}(\bm{\theta}|\bm{d}=\bm{\mu}(\bm{\theta}^{\rm{fid}}))$ from MCMC and $\mathbb{P}(\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})|\bm{\theta}^{\rm{fid}})$ from simulations are the same mathematically. In this work, the fiducial model is from fitting the $\it{Planck}$ TT data. We will show in Section 3.2 that the exact choice of $\bm{\theta}^{\rm{fid}}$ is not very important.

This correspondence allows us to estimate the correlation between extension parameters $A$ and $B$ using simulated data (i.e., frequency sampling of the likelihood) in the following steps:

Generate many simulated data sets (in our case simulated $\it{Planck}$-like CMB power spectra) drawn from the likelihood in the form of a Gaussian distribution $\mathcal{N}(\bm{\mu}(\bm{\theta}^{\rm{fid}}),\bm{\mathcal{C}})$ for some choice of fiducial $\Lambda\mathrm{CDM}$ model parameters, $\bm{\theta}^{\rm fid}$, and using the covariance matrix $\bm{\mathcal{C}}$ provided by the experiment collaboration. 2. 2.

Calculate the maximum-likelihood parameters $\bm{\theta}^{\rm{ML}}$ for the \Lambda\mathrm{CDM}$$+A and \Lambda\mathrm{CDM}$$+B models for each simulated data set, and 3. 3.

Estimate the covariance between $A$ in the \Lambda\mathrm{CDM}$$+A fit and $B$ in the \Lambda\mathrm{CDM}$$+B fit using the sample covariance from the simulations.

2.2. Estimating correlation between extension parameters using MCMC chains

Alternatively, in the special case of Gaussian parameter posteriors, we can run MCMC chains to estimate the correlation between extension parameters directly. In Appendix B we show that all the information on the correlation between extension parameters $A$ and $B$ can be estimated from three sets of chains: \Lambda\mathrm{CDM}$$+A, \Lambda\mathrm{CDM}$$+B, and \Lambda\mathrm{CDM}$$+A+B. Specifically, we found that the correlation between $A$ and $B$ when they are fitted separately is equal to minus one times the correlation between $A$ and $B$ when they are fitted together. We refer to this property as Correlation Equivalence.

A complication of estimating correlation between extension parameters by performing MCMC on the real data in this way is that we clearly do not have control over the ‘underlying’ cosmological model in the way we do when generating simulations. It is also possible that the real data have imperfections or systematic biases that are not correctly accounted for in the likelihood. We therefore instead performed MCMC computations replacing the real data vector with a theory prediction computed using the same $\bm{\theta}^{\rm fid}$ as for the simulations.

2.3. An example using $\it{Planck}$ CMB spectra

Here we provide a worked example using $\it{Planck}$ data for both the simulation and MCMC approaches discussed above. At the time of writing, the $\it{Planck}$ 2018 likelihood code is not yet available. Therefore we perform a simple three-parameter test using the $\it{Planck}$ $\ell\geq 30$ 2015 TT data from the plik_lite likelihood111Can be downloaded from http://pla.esac.esa.int/pla/#cosmology (as described in Planck Collaboration XI, 2016). This simplified likelihood includes only CMB information, marginalizing over foreground template amplitudes and other nuisance parameters.

2.4. Fiducial model and extension parameters

We assume the fiducial model to be the best-fit $\Lambda\mathrm{CDM}$ model of $\ell\geq 30$ $\it{Planck}$ TT plik_lite data with the optical depth fixed at $\tau=0.07$, the $\it{Planck}$ calibration parameter calPlanck equal to 1 and other nuisance parameters marginalized. The $\ell<30$ likelihood is pixel based rather than power spectrum based. For simplicity we only include the power spectrum based likelihood of the $\ell\geq 30$ data. Moreover, because the $\ell\geq 30$ TT spectrum only well constrains the parameter combination $A_{s}e^{-2\tau}$, $\tau$ is fixed to break the strong degeneracy with $A_{s}$. As for the calibration parameter calPlanck, a Gaussian prior with mean equal to 1 and standard deviation 0.0025 was originally imposed on it. Its variation has minimal impact on other parameters. For simplicity, calPlanck is fixed at 1.

For the standard $\Lambda\mathrm{CDM}$ parameters, we use

[TABLE]

In addition, the three extension parameters of interest and their fiducial values for this example are $\{A_{L},n_{\rm run},Y_{P}\}^{\rm{fid}}=\{1,0,0.2453\}$. As mentioned in Section 1, $A_{L}$ is a phenomenological parameter that artificially scales the lensing power spectrum. It is worth investigating as $\it{Planck}$ has a curious preference for $A_{L}>1$. Running of the spectral index, $n_{\rm{run}}\equiv dn_{s}/d\ln k$, and the primordial Helium mass fraction $Y_{P}$ are an interesting pair of extension parameters to test because of their high correlation ($\sim 0.9$), that is, they produce similar changes in the power spectrum damping tail. Even a small deviation from their degeneracy direction in their expected 2-D distribution would be noticeable. In $\Lambda\mathrm{CDM}$, $Y_{P}$ is calculated from $\Omega_{b}h^{2}$ and the CMB temperature through big-bang nucleosynthesis (BBN) predictions (Planck Collaboration XIII, 2016). In \Lambda\mathrm{CDM}$$+Y_{P}, $Y_{P}$ is independent of BBN and decouples from $\Omega_{b}h^{2}$.

The fiducial parameters from fitting the $\it{Planck}$ data have uncertainties due to cosmic variance and experimental noise. A slightly different set of fiducial parameters, which are still consistent with the measured power spectrum, may produce different results in the covariance of the extension parameters as well as the significance of their deviations from the fiducial $\Lambda\mathrm{CDM}$ values. We found this effect to be small for the $\it{Planck}$ 2015 data. See Section 3.2 for details.

2.5. Simulations

We run simulations to sample the distribution of the extension parameters. We draw 1000 binned power spectrum samples with mean equal to the binned power spectrum of the fiducial model and covariance equal to the plik_lite CMB band power covariance. We use the same binning scheme as $\it{Planck}$ 2015 to convert power spectra computed by CAMB222https://camb.info/ to band powers. We replace the data spectrum in the plik_lite likelihood with our samples, thus forming the simulated likelihoods. Next, using the ML finding algorithm in CosmoMC (Lewis & Bridle, 2002), setting $\tau=0.07$ and the $\it{Planck}$ calibration parameter calPlanck = 1, we maximize the simulated likelihoods to obtain the best-fit parameters of a specific model for each realization. The models we explore are: the standard $\Lambda\mathrm{CDM}$, $\Lambda\mathrm{CDM}$+$A_{L}$, $\Lambda\mathrm{CDM}$+$n_{\rm run}$ and $\Lambda\mathrm{CDM}$+$Y_{P}$.

To quantify the overall shift between the values of the extension parameters estimated from the actual data and their fiducial $\Lambda\mathrm{CDM}$ values, we calculate the $\chi^{2}$ and its probability to exceed (PTE) for a $\chi^{2}$ distribution with degrees of freedom equal to the number of extension parameters:

[TABLE]

where $\bm{\Sigma}_{\rm{ex}}$ is the covariance matrix for the extension parameters only, $\bm{\theta}^{\rm{fid}}$ is an array of 3 elements equal to $\{A_{L},n_{\rm run},Y_{P}\}^{\rm{fid}}$ and $\bm{\theta}^{\rm{ML}}_{\rm{dat}}=(1.1245,0.00773,0.2306)$ is an array whose values are obtained from fitting models $\Lambda\mathrm{CDM}$+$A_{L}$, $\Lambda\mathrm{CDM}$+$n_{\rm run}$ and $\Lambda\mathrm{CDM}$+$Y_{P}$ respectively on the actual $\it{Planck}$ plik_lite TT data with $\tau$ fixed.

2.6. MCMC

As mentioned in Section 2.2, we can make use of Correlation Equivalence to estimate the expected covariance $\bm{\Sigma}_{\rm{ex}}$ between extension parameters fitted separately by running MCMC chains on the data likelihood, with the experimental power spectrum replaced by the fiducial one. Again we set $\tau=0.07$ and the $\it{Planck}$ calibration parameter calPlanck= 1.

The diagonal elements in $\bm{\Sigma}_{\rm{ex}}$ are estimated from the variances of the specific extension parameters from running the MCMC chains with the modified $\it{Planck}$ plik_lite TT likelihood, on the model $\Lambda\mathrm{CDM}$+$A_{L}$, $\Lambda\mathrm{CDM}$+$n_{\rm run}$ and $\Lambda\mathrm{CDM}$+$Y_{P}$.

For the off-diagonal elements, first we obtain the correlation between two extension parameters (again, denoting them generically as $A$ and $B$) from the MCMC runs that vary both of them in the \Lambda\mathrm{CDM}$$+A+B model with the same dataset. Then we calculate the covariance between $A$ and $B$, given their variances as described in the last paragraph:

[TABLE]

This way, we are able to estimate all of the elements in $\bm{\Sigma}_{\rm{ex}}$ without running simulations at all, and calculate the $\chi^{2}$ defined in Equation (1) and its PTE to quantify the shifts of the extension parameters from their fiducial values.

3. Results

3.1. Quantifying significance of deviations of extension parameters from their fiducial values

In Table 1, we show the $\chi^{2}$ of the difference, as defined in Equation (1), their corresponding PTE values and the level of consistency with the assumed model in terms of the number of $\sigma$. Our results imply no significant discrepancy (0.9$\sigma$) between the values of the three additional parameters estimated from the actual data in one-parameter extensions and their fiducial $\Lambda\mathrm{CDM}$ values. This is consistent with expectations if the standard $\Lambda\mathrm{CDM}$ is the correct model.

Taking a closer look at elements in the estimated covariance, we show in Table 2 the estimated uncertainties of one-parameter extensions to the base $\Lambda\mathrm{CDM}$ model, along with the estimated parameter correlations from the two methods. These results are consistent up to at most 5%, in the $Y_{P}$ uncertainty. The differences between using the two methods are likely due to numerical noise and have only a minimal impact on our main conclusion.

We also visualize the shifts of the ML extension parameters from their fiducial values compared to their expected distributions in Figure 2. Notice in the 2-D contour plots, how the ML points lie along the correlation/degeneracy direction of the estimated distributions, as expected if the base $\Lambda\mathrm{CDM}$ is the true model. We do not attempt to explain the 1.8$\sigma$ (from the MCMC method) or 1.9$\sigma$ (from simulations) difference between the value of $A_{L}$ inferred by the plik_lite TT data and its physical value $A_{L}=1$. Rather we note that there is no extra sign of discrepancy when we take into account its non-negligible correlation with $n_{\rm{run}}$ or $Y_{P}$: 0.31 and -0.46 respectively. Had we found inconsistency between the ML parameters from the data and the expected $\Lambda\mathrm{CDM}$ distribution, it could indicate the presence of systematic error or unknown physical effects producing spurious parameter correlations.

In addition, numerical differences between the simulated contours and those of MCMC are small in this figure, again providing confidence to both methods.

3.2. Testing stability of results against uncertainties in fiducial model

As mentioned in Section 2.4, the fiducial parameters are from fitting the actual plik_lite data and therefore have uncertainties. So the standard $\Lambda\mathrm{CDM}$ model we define is not just one point in the parameter space, but an ensemble of points whose shape is described by the parameter posterior distribution inferred from the data.

To assess the impact of different fiducial values on our results, we randomly draw 1000 sets of parameters from the MCMC chains fitted to the plik_lite data, with the base $\Lambda\mathrm{CDM}$ model, $\tau=0.07$ and calPlanck = 1. For each set of parameters, we approximate the parameter covariance matrix by computing the $C_{\ell}$ derivatives and Fisher matrices (see Equation (B3) in Appendix B). Then we use Correlation Equivalence to calculate the extension parameter covariance $\bm{\Sigma}_{\rm{ex}}$. We find that varying parameters causes less than 3% scatter in the matrix elements.

We also plug $\bm{\Sigma}_{\rm{ex}}$ into Equation (1) and find that the resulting $\chi^{2}$ ranges from 4.2 to 5.1 and the PTE from 0.17 to 0.24, corresponding to consistency within 0.7-1$\sigma$. These results show that the uncertainties in the fiducial model have very minimal impact on results from our consistency test. The stability of the results reflects the constraining power of the $\it{Planck}$ data on the model parameters in $\Lambda\mathrm{CDM}$.

4. Discussion

In Sections 2 and 3 we have shown an example of quantifying the level of consistency between extension parameters and their fiducial values given a specific dataset. In this section, we outline and discuss steps to perform this approach more comprehensively for the $\it{Planck}$ 2018 and other cosmology data sets.

4.1. MCMC method

In the ideal case where all extension parameters of interest are Gaussian, we can estimate their expected distribution from MCMC chains. We denote the individual extension parameter as ${\theta}_{\rm{ex},i}$ with $i$ runs over 1 to $n_{\rm{ex}}$, the total number of extension parameters. The suggested recipe is as follows:

Define a fiducial model (e.g. a fit from existing data) and calculate its theory prediction. 2. 2.

In the data likelihood of interest, substitute the experimental observables, e.g. power spectra in the CMB case, by the fiducial prediction. 3. 3.

Explore parameter space around the fiducial point by running MCMC chains on the modified data likelihood(s), fitting models \Lambda\mathrm{CDM}$$+{\theta}_{\rm{ex},i} and \Lambda\mathrm{CDM}$$+{\theta}_{\rm{ex},i}+{\theta}_{\rm{ex},j\neq i}. 4. 4.

Calculate the variance of ${\theta}_{\rm{ex},i}$ from the one-parameter extension, and the correlation between ${\theta}_{\rm{ex},i}$ and ${\theta}_{\rm{ex},j\neq i}$ from the two-parameter extension. 5. 5.

Using results from previous step, construct a covariance matrix $\bm{\Sigma}_{\rm{ex}}$ for $\bm{\theta}_{\rm{ex}}$, with the signs of all the correlations flipped. 6. 6.

Calculate $\chi^{2}$ defined in Equation (1) and its PTE to quantify deviation of experimental values to fiducial ones. Additionally, one can plot confidence ellipses such as in Figure 2 to visualize the deviations. 7. 7.

Check validity of the input fiducial model, for example by using a Fisher Forecast (e.g. Heavens, 2009) to approximate parameter covariance and estimating the shifts in the $\chi^{2}$ when varying fiducial parameter values.

As an example, on a computer cluster with 12-core 2.5 GHz processors333Our computation was conducted on the computer cluster of the Maryland Advanced Research Computing Center. See https://www.marcc.jhu.edu/cyberinfrastructure/hardware/ for descriptions of its system architecture., one MCMC run with 8 parallel chains usually take $\lesssim$ 24 hours to converge for the 2015 plik_lite likelihood. The CPU time for running one MCMC job is at the order of 100 hours and the total computing time is $\sim 5\times 10^{4}$ hours, for running all one-parameter and two-parameter extension fits for e.g. $n_{\rm{ex}}=11$, which is the number of extra parameters fitted in the publicly available 2015 $\it{Planck}$ chains.

4.2. Simulations

We can use simulations, as an alternative method to estimate the expected distribution of additional parameters around the fiducial point. As described in Section 2.1 and 2.5, one can perform the following procedure:

Generate simulated observables of the fiducial model using data covariance. 2. 2.

For each simulation, estimate the best-fit $\Lambda\mathrm{CDM}$+$\bm{\theta}^{ex}_{i}$ model for each extension parameter. 3. 3.

Calculate the covariance of extension parameter $\bm{\Sigma}_{\rm{ex}}$, the $\chi^{2}$ of difference defined in Equation (1) and its PTE to quantify the significance of difference. Confidence ellipses can also be plotted.

For a rough estimate of computing time for the simulation method, using CosmoMC and the same computer cluster with 12 cores per CPU, we find that the running time of the best-fit-finding algorithm is approximately 1 hour. For $n_{\mathrm{ex}}=11$, $n_{\mathrm{sim}}=2000$ and running 4 parallel best-fit-finding jobs for each simulation (to reduce numerical noise and avoid obtaining results from local minima), the total computing time is $9\times 10^{4}$ hours.

4.3. In the case of non-Gaussianity

To discuss the case of non-Gaussianity, first we need to clarify what exactly non-Gaussianity arises from. Recall that throughout this paper, we assume a Gaussian data likelihood and a prior that is approximately flat. If the data constrains the parameters well enough, changes in the model predictions can be treated as only linearly dependent on the parameters. This means the Taylor expansion of the log likelihood around the maximum point is significant up to quadratic terms of the parameters and therefore the parameters are Gaussian (see Appendix A).

In short, non-Gaussianity of parameters is a result of the data not being constraining enough for fitting parameters. When this is the case, there may not be a mathematical correspondence between the ML parameter distribution from the simulations and the posterior distribution of the MCMC chains, since our argument in Appendix A depends on the assumption of Gaussianity. Besides, since our proof of Correlation Equivalence in Appendix B also rests upon approximations of the log likelihood to only the second order derivatives, the property now becomes questionable. This means that we cannot simply estimate the correlation between parameters fit separately from the MCMC chains where they are fitted together. What is worse is that the means and the covariance matrix no longer carry all information of the parameter distributions and one might need to evaluate high order tensors (Sellentin et al., 2014). Additionally, the $\chi^{2}$ test is inappropriate for non-Gaussian parameters and we need new ways (such as one proposed in Appendix C of Motloch & Hu, 2019b) to quantify the significance of the overall difference between the experimental parameters and their fiducial values.

An example of non-Gaussian parameters are those with priors that are more informative than the data, e.g. the neutrino mass $\sum m_{\nu}$ and the tensor to scalar ratio $r$, physically defined as non-negative. Current CMB data are not sufficiently constraining to pull these parameters away from lower bounds (BICEP2 Collaboration et al., 2018; Planck Collaboration VI, 2018). Another non-Gaussian parameter given just the TT data, is the curvature density $\Omega_{k}$. Curvarture can only be weakly constrained as allowing it to be free worsens the existing degeneracy between the physical matter density $\Omega_{m}$ and the Hubble constant $H_{0}$ (Zaldarriaga et al., 1997; Percival et al., 2002; Kable et al., 2019). In the MCMC method, a two-parameter extension model with both $A_{L}$ and $\Omega_{k}$ results in very wide and non-Gaussian distributions, since they are highly correlated — a positive curvature has similar effect as an increased lensing signal (Planck Collaboration XIII, 2016).

To reduce non-Gaussianity, one can always include more datasets if available to have greater constraining power on the parameter, e.g. include the BAO data (Alam et al., 2017) in parameter fitting along with $\it{Planck}$. However, when extensions are used as a means to test the internal consistency of one dataset, adding extra data is not an option. Fortunately, there are existing methods of transforming non-Gaussian parameters into Gaussian ones. One such method is the Box-Cox transformation (Box & Cox, 1964). Joachimi & Taylor (2011) and Schuhmann et al. (2016) applied this bijective transformation to non-Gaussian cosmological parameters. Theoretically, Correlation Equivalence still holds for the transformed Gaussian parameters. Thus we may still apply the procedure outlined in Section 4.1 for transformed parameters, obtain the expected distribution for transformed ones, and then transform them back to the model parameters. One thing to note is that the Box-Cox transformation does not guarantee Gaussianity. So one should check for the Gaussianity of transformed parameters, e.g. calculate the skewness and kurtosis of the resulting distribution, and if needed, apply a second transformation. For consistency, it is also a good idea to compare the resulting distribution of model parameters to simulations, or compare the covariance of the transformed parameters with predictions from Fisher Forecast, keeping in mind that Fisher Forecast assumes Gaussianity.

Further understanding of non-Gaussian scenarios is left for future efforts.

5. conclusions

In this paper, we have presented a method to help identify potential new physics and/or systematic errors by calculating the correlations between additional parameters and performing a $\Lambda\mathrm{CDM}$ consistency test accounting for them.

Usually extension parameters are added to the model separately, one a time. However, different parameters may affect the theory prediction in a similar way, which means their values from a given data set can be correlated, even when they are fit separately. Examining the consistency of extension parameters with $\Lambda\mathrm{CDM}$ expectations, accounting for these correlations, provides an additional test of the model beyond looking at the results from individual one-parameter extensions.

Under the assumption of Gaussianity of both the likelihood and the posterior distribution of the parameters, one can fit a series of one-parameter extension models to simulated data and obtain the multivariate distribution of the extension parameters. With the base parameters marginalized over, the $\chi^{2}$ of difference and its PTE can be calculated to quantify the significance of deviations from $\Lambda\mathrm{CDM}$.

A more computationally economic alternative is to run MCMC, fitting the same series of one-parameter extension models and additionally two-parameter extensions across all the possible pairs in the set of parameters of interest. Using Correlation Equivalence as proven in Appendix B, the covariance matrix of the extension parameters fitted separately can be estimated from the results of these MCMC runs and so the expected multivariate Gaussian distribution can be obtained.

In an attempt to narrow down causes of the $A_{L}$ anomaly in the $\it{Planck}$ data and possibly shed light on tensions between cosmological measurements, we looked at the example parameter combination $\{A_{L},n_{\rm run},Y_{P}\}$, fitted with the $\it{Planck}$ 2015 plik_lite $\ell\geq 30$ TT data, as we do not yet have access to the $\it{Planck}$ 2018 likelihood. Results from MCMC and simulations show that the deviations of the three additional parameters from their fiducial $\Lambda\mathrm{CDM}$ values are consistent with statistical fluctuations within 0.9$\sigma$ when correlations are accounted for.

Although the cause of the reported 1.8-2.7$\sigma$ preference (depending on the the specific combination of datasets) for $A_{L}>1$ by the $\it{Planck}$ CMB data is yet to be understood, we find no further evidence for discrepancy when considering the correlations between $A_{L}$, $n_{\rm{run}}$ and $Y_{P}$. This is not a trivial test, as the correlations are significant: approximately 0.31 for $A_{L}$-$n_{\rm{run}}$, -0.46 for $A_{L}$-$Y_{P}$, and -0.93 for $n_{\rm{run}}$-$Y_{P}$. If the unphysical $A_{L}>1$ is a symptom of an underlying systematic error or some real but unknown physical effect that also produced spurious correlations with $n_{\rm run}$ or $Y_{P}$ our test could have revealed this.

We also tested the stability of results against the uncertainties in the parameters of the experimentally fitted fiducial model, and find that the change of the fiducial model has no impact on our conclusions, only shifting the PTE values from 0.17 to 0.24.

Furthermore, we discussed how our procedures depend on the assumption of Gaussianity of parameters. If the assumption is not valid, MCMC runs cannot simply be used to estimate correlations between parameter fitted separately, nor may there be a mathematical correspondence between parameter distributions from the simulations and the MCMC runs. Therefore, efforts might attempt to include Gaussianization of non-Gaussian parameters, such as using the Box-Cox transformation (Box & Cox, 1964).

The procedures developed in this paper can and should be applied to more extensive lists of extension parameters, with existing and future cosmological data, in order to provide a more stringent test and complete view of $\Lambda\mathrm{CDM}$ consistency.

Acknowledgments

This research was supported in part by NASA grants NNX16AF28G and NNX17AF34G. This work was partly based on observations obtained with $\it{Planck}$ (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. This research project was conducted using computational resources at the Maryland Advanced Research Computing Center (MARCC). The GetDist Python package was used to make Figures 1 and 2. We would like to thank Josh Kable for providing his code for Fisher matrix calculations, Haoyu Wang for his help with the derivation in Appendix B, and Janet Weiland for her comments on a draft of this paper.

Appendix A mathematical correspondence between frequentist maximum-likelihood parameters and bayesian parameter posterior

In this section, we will show that there is a mathematical correspondence between frequentist maximum likelihood (ML) parameters and Bayesian parameter posterior, under conditions described below. In other words, $\mathbb{P}(\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})|\bm{\theta}^{\rm{fid}})$, the distribution of ML parameters estimated from realizations of a fiducial model is mathematically the same as $\mathbb{P}(\bm{\theta}|\bm{d}=\bm{\mu}(\bm{\theta}^{\rm{fid}}))$, the Bayesian posterior distribution of parameters given a set of data that matches the theory prediction of the same fiducial model.

Given a set of fiducial parameters $\bm{\theta}^{\rm{fid}}$, we can calculate its theory prediction $\bm{\mu}(\bm{\theta}^{\rm{fid}})$. With $\bm{\mathcal{C}}$ as the data covariance estimated experimentally, we can then draw realizations of the data vector $\bm{d}$ from a multivariate Gaussian distribution $\mathcal{N}(\bm{\mu}(\bm{\theta}^{\rm{fid}}),\bm{\mathcal{C}})$

Up to some constants, the log probability density function of $\bm{d}$ is given by

[TABLE]

For a given set of simulated data $\bm{d}^{\rm{sim}}$, we can estimate the set of parameters that best describe it by maximizing the log-likelihood

[TABLE]

The ML parameter estimates, $\bm{\theta}^{\rm ML}$, are defined as the solutions to the simultaneous equations

[TABLE]

where the index $i$ runs over all parameters. Taylor expanding around $\bm{\theta}^{\rm ML}$ gives

[TABLE]

If higher-order terms can be neglected, Equation (A4) shows that given a set of $\bm{d}$, parameters around the ML point can be approximated by a Gaussian distribution, with the covariance $\bm{\Sigma}$ given by

[TABLE]

where the term inside the brackets is the Fisher matrix and $\bm{\theta}^{\rm ML}$ is written as $\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})$ to emphasize its dependence on $\bm{d}^{\rm{sim}}$.

Note that the negligibility of higher order terms in (A4) means that $\bm{\Sigma}$ is approximately constant for different values of $\bm{\theta}^{\rm ML}$. For simplicity we can set it to be $\bm{\Sigma}(\bm{\theta}^{\rm{fid}})$.

Furthermore, assuming that the zeroth order term \log\mathbb{P}(\bm{d}^{\rm{sim}}|\bm{\theta})\big{|}_{\bm{\theta}=\bm{\theta}^{\rm ML}} is also approximately a constant for all $\bm{\theta}^{\rm ML}$, we can interchange the positions of $\bm{\theta}$ and $\bm{\theta}^{\rm ML}$ in (A4) and set $\bm{\theta}=\bm{\theta}^{\rm{fid}}$, which gives $\mathbb{P}(\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})|\bm{\theta}^{\rm{fid}})$, the distribution of $\bm{\theta}^{\rm ML}$ from all simulated data $\bm{d}$. Up to some constants, the log of $\mathbb{P}(\bm{\theta}^{\rm ML}(\bm{d}^{\rm{sim}})|\bm{\theta}^{\rm{fid}})$ is:

[TABLE]

From the Bayesian viewpoint, given a single set of data $\bm{d}$, the posterior parameter distribution is given by Bayes’ theorem,

[TABLE]

where $\mathbb{P}(\bm{\theta})$ is the prior probability and $\mathbb{P}(\bm{d})$ is the evidence. If the posterior is Gaussian we can write, again up to some constants,

[TABLE]

where $\bm{\Sigma}$ is the parameter covariance matrix given by (A5). For flat priors $\bar{\bm{\theta}}=\bm{\theta}^{\rm ML}(\bm{d})$. If $\bm{\theta}^{\rm fid}$ mentioned above is close to $\bm{\theta}^{\rm ML}(\bm{d})$, then (A6) and (A8) is only different by a small offset. To make them exactly equal, we can choose to replace $\bm{d}$ with a theory prediction computed using $\bm{\theta}^{\rm{fid}}$ so that $\bm{\theta}^{\rm ML}$ becomes equal to $\bm{\theta}^{\rm{fid}}$ and the Bayesian posterior now describes the distribution of parameters around the fiducial value:

[TABLE]

which matches exactly the distribution of the frequentist ML parameter estimates given in (A6). This is asymptotically true even for non-flat priors when the data are sufficiently constraining, and we do not discuss the exact choice of priors further here. For more information we again direct the reader to, for example, Chapter 4 of Gelman et al. (2013).

Appendix B Correlation Equivalence

In this section, using the mathematical correspondence from Appendix A, we show that when using maximum likelihood estimation, the correlation between two parameters varied separately (e.g. the correlation between parameters $A$ and $B$ in model $\Lambda\mathrm{CDM}$+$A$ and $\Lambda\mathrm{CDM}$+$B$) is the same but with an opposite sign as the correlation between the same two parameters varying together (e.g. $\Lambda\mathrm{CDM}$+$A+B$).

B.1. Maximum likelihood estimation and parameter covariance

With $\mathfrak{L}\equiv\log L(\bm{\theta})$, the log likelihood of parameters given a specific sample of $\bm{d}$ can be written as

[TABLE]

up to some constants. $\bm{\mu}(\bm{\theta})$ is the theory prediction of the data given some set of parameters $\bm{\theta}$. For any $\bm{d}$, the goal is to find the $\bm{\theta}$ that maximize $\mathfrak{L}$.

As in Appendix A, at maximum likelihood, the partial derivative of the log likelihood with respect to (w.r.t.) ${\theta}_{i}$ is zero ($i$ runs over $n$, the number of parameters). That is

[TABLE]

Close to the ML point in parameter space, we can Taylor expand the log likelihood to the second order in $\Delta{\theta}_{i}\equiv{\theta}_{i}-{\theta}^{\rm{ML}}_{i}$ as shown in Equation (A4), where ${\theta}_{i}$ is well described by a multivariate Gaussian distribution. The parameter means are the values that maximize the likelihood.

The parameter covariance $\bm{\Sigma}$ can be calculated as the inverse of the Fisher matrix $\mathcal{F}^{-1}$, with the Fisher matrix, estimated from MCMC chains or calculated analytically, defined as:

[TABLE]

From now on we use the comma notation to denote partial derivatives w.r.t. $\bm{\theta}$, that is, $\bm{\mu}_{,i}\equiv\frac{\partial\bm{\mu}}{\partial\theta_{i}}$.

In the following proof, we continue to work under the assumption of Gaussianity for $\bm{\theta}$ and expand $\bm{\mu}$ only to the first order of $\bm{\theta}$, that is, assuming the Gaussian linear model (Raveri & Hu, 2018), around a chosen fiducial model $\bm{\theta}^{\mathrm{fid}}$:

[TABLE]

where we define $\bm{\mu}^{\mathrm{fid}}=\bm{\mu}(\bm{\theta}^{\mathrm{fid}})$. In the Gaussian linear model, the Jacobian $\bm{\mu}_{,i}$ may be taken as constant over changes in $\bm{\theta}$ and so $\bm{\mu}^{\mathrm{fid}}_{,i}$ is approximately equal to $\bm{\mu}^{\mathrm{ML}}_{,i}$. For simplicity, from here on, we just drop all the superscripts for $\bm{\mu}_{,i}$. In practice, the fiducial model is usually based on to the ML model given by the data and updated iteratively if necessary.

Next, we move on to calculate the elements in parameter covariance $\bm{\Sigma}$:

[TABLE]

Remember that for any invertible square matrix $\mathcal{M}$, its inverse can be calculated as

[TABLE]

where $|\mathcal{M}|$ is the determinant of $\mathcal{M}$. We use the notation “$\mathcal{M}_{\backslash j\backslash i}$”, to denote a smaller matrix, corresponding to $\mathcal{M}$ with the $j^{\rm{th}}$ row and the $i^{\rm{th}}$ column removed (the backslash symbol is borrowed from the notation for set difference in set theory). Then ${|\mathcal{M}_{\setminus j\setminus i}|}$ is a minor of $\mathcal{M}$. And the determinant for an $n\times n$ matrix $\mathcal{M}$ can be written as

[TABLE]

for any $1\leq j\leq n$.

Therefore

[TABLE]

If we choose to fix the $i^{\rm{th}}$ parameter, we can simply remove the $i^{\rm{th}}$ row and column from the Fisher matrix, and calculate a new parameter covariance matrix from the revised Fisher matrix.

B.2. Correlation between parameter $A$ and $B$, varying together

When two parameters $A$ and $B$ are both variables, we can calculate their covariance from the Fisher matrix that includes them along with the base $\Lambda\mathrm{CDM}$ parameters. So their correlation is just

[TABLE]

From here on, for clarity, we use $A$ and $B$ as subscripts for the rows and columns in the matrix for the two parameters of interest.

Following from equation (B8), (B9) becomes

[TABLE]

The minus sign is from setting the index of $B$ equal to that of $A$ plus one.

B.3. Correlation between parameter $A$ and $B$, varying separately

For a generic set of parameters, given the data array $\bm{d}$ (experimental or simulated) and the Fisher matrix $\mathcal{F}$, we can substitute (B5) into (B2) to obtain a relationship between the ML parameters $\bm{\theta}^{\rm{ML}}$ and the fiducial ones:

[TABLE]

So we can express $\bm{\theta}^{\rm{ML}}$ as:

[TABLE]

where we use the vector $\hat{\bm{\theta}}$ to represent all the terms in (B11) that are independent of data $\bm{d}$:

[TABLE]

As we will show below, $\hat{\bm{\theta}}$ does not contribute to the correlation between $A$ and $B$.

Again for clarity, we use $A$ and $B$ to index the rows and columns for parameter $A$ and $B$, respectively, while using the generic $i$ and $j$ to index $n$ base $\Lambda\mathrm{CDM}$ parameters. So $A$ is the $(n+1)^{\rm{th}}$ parameter and $B$ $(n+2)^{\rm{th}}$.

When fixing $B$ to its fiducial $\Lambda\mathrm{CDM}$ value $B^{\mathrm{fid}}$, the expression for the best estimate for other parameters is almost the same as equation (B12), except that $\mathcal{F}\rightarrow\mathcal{F}_{\backslash B\backslash B}$ and the elements associated with $B$ in $\bm{\mu}^{T}_{,i}$, $\bm{\theta}^{\rm{ML}}$ and $\hat{\bm{\theta}}$ are deleted. The expression for the rest of elements in $\bm{\theta}^{\rm{ML}}$ is:

[TABLE]

The last row of (B13) gives

[TABLE]

Similarly, fixing $A$ and letting $B$ vary leads to

[TABLE]

The variance of $A$ and $B$ can be obtained from the Fisher matrix (B3) and its minors (B7), as follows:

[TABLE]

And the covariance between $A^{\rm{ML}}$ and $B^{\rm{ML}}$ is:

[TABLE]

with the brackets here representing the averaging over all realizations of the data $\bm{d}$ .

Recall that the data covariance $\bm{\mathcal{C}}$ is assumed to be fixed and with our assumption of a Gaussian Linear Model, all partial derivatives w.r.t. the parameters are also constant. Then in (B14) and (B15), only $\bm{d}$ is a variable. Inserting (B14) and (B15) into (B18), we find that terms involving $\hat{{\theta}}_{A}$ and $\hat{{\theta}}_{B}$ cancel. In addition, all constant terms can be taken out of the brackets, leaving only $\langle(\bm{d}-\langle\bm{d}\rangle)^{T}(\bm{d}-\langle\bm{d}\rangle)\rangle$, which is just $\bm{\mathcal{C}}$. This simplifies to

[TABLE]

Notice that terms like $\bm{\mu}^{T}_{,i}\bm{\mathcal{C}}^{-1}\bm{\mu}_{,j}$ are just elements of the Fisher matrix and we can also use equation (B7) to express terms like $(\mathcal{F}_{\backslash B\backslash B})^{-1}_{Ai}$ in terms of determinants of minors of the Fisher matrix. Thus,

[TABLE]

To simplify (B20), first we combine term [1] and [2]:

[TABLE]

where we have used the property that the Fisher matrix is symmetric.

Compared to equation (B8), note that the terms inside the last bracket above sum up to the determinant of an $(n+1)\times(n+1)$ matrix, which is the same as $\mathcal{F}_{\backslash A\backslash A}$, except that its $(n+1)^{\rm{th}}$ column is replaced by its $i^{\rm{th}}$. So not all of the columns for this matrix are linearly independent, resulting in its determinant being zero. Thus $[1]+[2]=0$.

To simplify term [3] and [4], we have

[TABLE]

Insertng (B21) into (B20), we have

[TABLE]

Using Equation (B16) and (B17), the correlation between the best-fit values of parameter $A$ and $B$ estimated separately can be written as

[TABLE]

which is the same as equation (B10) up to a minus sign, thus completing our proof.