Validating the Fisher approach for stage IV spectroscopic surveys

A&A 649, A52 (2021) https://doi.org/10.1051/0004-6361/201937312 c S. Yahia-Cherif et al. 2021 Astronomy & Astrophysics Validating the Fisher approach for stage IV spectroscopic surveys S. Yahia-Cherif1, A. Blanchard1, S. Camera2,3,4, S. Casas5, S. Ili´ c6,7,1, K. Markoviˇ c8, A. Pourtsidou9, Z. Sakr1,10, D. Sapone11, and I. Tutusaus12,13,1 1Université de Toulouse, UPS-OMP, IRAP and CNRS, IRAP, 14, avenue Edouard Belin, 31400 Toulouse, France e-mail: [email protected] 2Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy 3INFN – Istututo Nazionale di Fisica Nucleare, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy 4INAF – Istituto Nazionale di Astrofisica, Osservatorio Astrofisico di Torino, Strada Osservatorio 20, 10025 Pino Torinese, Italy 5AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France 6Université PSL, Observatoire de Paris, Sorbonne Université, CNRS, LERMA, 75014 Paris, France 7CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2, Praha 8, Czech Republic 8Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr, Pasadena, CA 91109, USA 9School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK 10 Université St Joseph; UR EGFEM, Faculty of Sciences, Beirut, Lebanon 11 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile 12 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain 13 Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain Received 13 December 2019 /Accepted 28 July 2020 ABSTRACT In recent years, forecasting activities have become an important tool in designing and optimising large-scale structure surveys. To predict the performance of such surveys, the Fisher matrix formalism is frequently used as a fast and easy way to compute constraints on cosmological parameters. Among them lies the study of the properties of dark energy which is one of the main goals in modern cosmology. As so, a metric for the power of a survey to constrain dark energy is provided by the figure of merit (FoM). This is defined as the inverse of the surface contour given by the joint variance of the dark energy equation of state parameters {w0,wa}in the Chevallier-Polarski-Linder parameterization, which can be evaluated from the covariance matrix of the parameters. This covariance matrix is obtained as the inverse of the Fisher matrix. The inversion of an ill-conditioned matrix can result in large errors on the covariance coefficients if the elements of the Fisher matrix are estimated with insufficient precision. The conditioning number is a metric providing a mathematical lower limit to the required precision for a reliable inversion, but it is often too stringent in practice for Fisher matrices with sizes greater than 2×2. In this paper, we propose a general numerical method to guarantee a certain precision on the inferred constraints, such as the FoM. It consists of randomly vibrating (perturbing) the Fisher matrix elements with Gaussian perturbations of a given amplitude and then evaluating the maximum amplitude that keeps the FoM within the chosen precision. The steps used in the numerical derivatives and integrals involved in the calculation of the Fisher matrix elements can then be chosen accordingly in order to keep the precision of the Fisher matrix elements below this maximum amplitude. We illustrate our approach by forecasting stage IV spectroscopic surveys cosmological constraints from the galaxy power spectrum. We infer the range of steps for which the Fisher matrix approach is numerically reliable. We explicitly check that using steps that are larger by a factor of two produce an inaccurate estimation of the constraints. We further validate our approach by comparing the Fisher matrix contours to those obtained with a Monte Carlo Markov chain (MCMC) approach – in the case where the MCMC posterior distribution is close to a Gaussian – and finding excellent agreement between the two approaches. Key words. dark energy – cosmological parameters – large-scale structure of Universe – galaxies: statistics 1. Introduction Since the discovery of the acceleration of the expansion of the Universe in the late 1990s (Riess et al. 1998;Perlmutter et al. 1999), the Lambda Cold Dark Mattery (ΛCDM) model remains the most successful model in cosmology. It provides a simple and accurate description of the properties of the Universe, with a very limited number of parameters that are now constrained at the level of a few percent by using measurements from the Planck CMB mission (Planck Collaboration VI 2020) and other experiments (see e.g., Percival et al. 2001;Blake et al. 2011; Alam et al. 2017;Abbott et al. 2018). The accelerated expansion of the Universe can be explained by postulating a new form of matter or the mechanism dubbed “dark energy”, which would be present alongside dark matter. Dark energy dominates the energy budget of the Universe today, making up for ∼70% of its total energy density. Overthepast twodecades,arangeofdarkenergyexperiments have been proposed in order to study the observed cosmic acceleration of the Universe through various probes: optical spectroscopic and photometric galaxy clustering, and weak gravitational lensing, supernovae, with galaxy clusters being the most common. In order to quantify the performance of a survey and, in particular, its ability to constrain the properties of dark energy, the Dark Energy Task Force (DETF) defined a metric known as the figure of merit (FoM, Albrecht et al. 2006): a figure inversely proportional to the surface bounded by the confidence contours for the w0and waparameters from the Chevalier-Polarski-Linder Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. A52, page 1 of 19 A&A 649, A52 (2021) (CPL) parameterization (Chevallier & Polarski 2001;Linder 2003).TheapplicationoftheFoMisnowstandardfor quantifying the performance of a particular experiment that aims to constrain the dark energy parameters. We are currently entering the era of high-precision cosmology with a suite of stage IV experiments that are about to be deployed. Currently, ongoing and forthcoming ground-based experiments as well as space missions include: DESI1(DESI Collaboration 2016), Euclid2(Laureijs et al. 2011;Amendola et al. 2018;Euclid Collaboration 2020), LSST3(LSST Science Collaboration 2009;Ivezi´ c et al. 2019;Chisari et al. 2019), WFIRST4(Akeson et al. 2019;Dore et al. 2019), and the SKA5 (Square Kilometre Array Cosmology Science Working Group 2020). Forecasting activities are key for predicting future constraints on cosmological parameters and the corresponding dark energy FoM. Indeed, forecasts not only allow us to predict the performance of a future survey but also aid in its design and optimization. There are a number of different ways to compute forecasts. For instance, it’s possible to use a simulated realization of the data and a Monte Carlo Markov chain (MCMC) method to sample the likelihood of the parameters of interest. This approach is known to provide reliable estimations of the constraints on the cosmological parameters (see e.g., Dunkley et al. 2005). However, it is often quite time-consuming, especially if there is a wide variety of cases involved. The Fisher matrix formalism, as adopted by the DETF, is a computationally fast approach that allows for constraints to be obtained within a few seconds or minutes. However, this method is dubious from the technical point of view as the covariance matrix is obtained by inversion of the Fisher matrix, whose coefficient calculations often require the computation of derivatives and integrals, which can lead to inaccurate results in cases where, for example, the step sizes (as we show in subsequent sections of this paper) are not chosen properly to ensure the appropriate precision. Moreover, the Fisher formalism can miss contour shapes and, thus, miss potential degeneracies (see e.g., Wolz et al. 2012) as it assumes a Gaussian likelihood in a frequentist approach or a Gaussian posterior of the parameters in a Bayesian approach; that might not be true if, for example, the parameters are not tightly constrained by the experiment(s). Several methods have been proposed to handle these complex situations (Joachimi & Taylor 2011;Sellentin et al. 2014;Sellentin 2015;Amendola & Sellentin 2016). In the present study, we propose a numerical approach that is independent of the specific problem treated, with the purpose of testing the numerical reliability of constraints obtained within the Fisher matrix formalism. As a working example, we consider the spectroscopic galaxy clustering probe for the following stage IV surveys: DESI, Euclid, and WFIRST-2.4. The paper is organized as follows. In Sect. 2, we give a brief description of the aforementioned surveys. In Sect. 3, we describe the different tests performed in order to validate the Fisher approach: namely, we determine the required precision to compute a valid Fisher matrix, carry out several stability and convergence tests, and perform a comprehensive comparison with the MCMC approach. In Sect. 4, we present the results and address possible ways of numerically validating the Fisher matrix approach. Finally, we report our conclusions in Sect. 5. 1www.desi.lbl.gov 2www.euclid-ec.org 3www.lsst.org 4https://wfirst.gsfc.nasa.gov/ 5www.skatelescope.org/ 2. Surveys 2.1. Euclid Euclid is a medium-class ESA mission that is aimed at gaining an improved understanding of the expansion of the Universe, as well as dark energy and dark matter. Euclid will measure the shape of more than one billion galaxies and provide tens of millions of spectroscopic redshift. It will use two main cosmological probes, namely galaxy clustering and weak gravitational lensing. These will allow us to probe the dark sector of the Universe with unprecedented precision. Euclid (Laureijs et al. 2011) targets a FoM of 400, while also aiming to distinguish between different theories of gravity, that is, testing general relativity against alternative models by measuring the exponent of the growth factor, γ, with a precision higher than 0.02 at 1σ.Euclid will also measure the neutrino masses with a 1σprecision better than 0.03 eV. Combined with Planck,Euclid will also probe some inflation models through the measurement of the non-Gaussianity of initial conditions. So far, two surveys are planned: a wide one of 15 000 deg2and a deep one of 40deg2. The 1.2m telescope contains two instruments that will allow to observe around 2 billion galaxies and produce 50 million precise redshift measurements. Euclid’s large format visible imager (VIS), will probe the photometric galaxy clustering in the optical wavelength range. The near-infrared imaging and spectroscopy (NISP) instrument will perform photometry imaging and produce spectra in the near infrared range. 2.2. Dark Energy Spectroscopic Instrument (DESI) Known as the successor of the stage III BOSS and eBOSS surveys (Dawson et al. 2013,2016), DESI is a Stage IV ground-based dark energy experiment whose primary objectives are measuring the baryonic acoustic oscillations (BAO) and constraining the growth of structure through redshift space distortions measurements. It will allow for the estimation of the deviations from Gaussianity through the fNL parameter. DESI will also measure the neutrinos masses with a precision of 0.02 eV at 1σ, for a maximum scale of kmax <0.2hMpc−1. It will scan a 14 000 deg2sky area and measure more than 30 million spectra from four galaxy populations: the bright galaxies (BGS) at low redshift (until z≤0.6), the luminous red galaxies (LRGs) composed of highly biased objects at intermediate redshift z<1, the bright emission lines galaxies (ELGs) up to z=1.7 (which can only be detected at high resolution since the OII doublet must be well-resolved), and the quasi-stellar objects (QSOs). The latter group allows us to trace both the matter distribution at high redshift and the neutral hydrogen by the Ly −α absorption in the spectra. DESI is now underway and will conduct observations over five years (2019–2023) using the Mayall four-meter telescope at Kitt Peak. 2.3. Wide Field Infrared Survey Telescope (WFIRST) WFIRST is a NASA near-infrared imaging and low-resolution spectroscopy observatory that, similarly to Euclid, aims to address fundamental questions about the accelerated expansion of the Universe. WFIRST will determine the expansion history of the Universe and structure growth in order to constrain dark energy and modified gravity models. To estimate the dark energy equation of state, WFIRST will perform weak lensing, along with measurements of supernovae distance and baryonic A52, page 2 of 19 S. Yahia-Cherif et al.: Validating Fisher approach Table 1. Redshift range, binning choices, and density of galaxies for the surveys used in this work. Surveys zmin zmean zmax dN(zmean)/dΩdz[deg−2] DESI(BGS) 0 0.125 0.25 3029 0.25 0.375 0.5 528 DESI(ELGs) 0.5 0.6 0.7 309 0.7 0.8 0.9 2124 Euclid 0.9 1.05 1.2 2287 1.2 1.35 1.5 1574 1.5 1.65 1.8 764 WFIRST-2.4 1.8 1.95 2.1 6718 2.1 2.25 2.4 1368 2.4 2.55 2.7 781 Notes. For more details on survey specifications, see Sect. 3. acoustic oscillations. In this work, we use the spectroscopic probe of WFIRST-2.4, a 2000deg2survey that will observe more than 20 million galaxies in the redshift range of 1 ≤z≤3. 3. Methodology 3.1. The surveys combination specification In general, using a single probe of a stage IV survey does not constrain cosmological parameters well enough to get a Gaussian posterior distribution, in which case testing simply the Fisher formalism is meaningless. However, when the constraints are very tight, the posterior is likely to be close to gaussian. For this reason, we chose to investigate the Fisher formalism in a case where we can anticipate a Gaussian posterior (an assumption that is checked by the MCMC forecast used for the full validation of the method). In order to substantially improve the constraints, we therefore combine Euclid, DESI, and WFIRST-2.4. However, we note that DESI and WFIRST-2.4 cover Euclid’s entire redshift range. Hence, if we want to fully and consistently combine these three surveys, computing the cross-correlation power spectra between those surveys is necessary. For the purposes of our study, we do not wish to consider the cross-spectra, therefore, we have to remove all the DESI and WFIRST-2.4 redshift bins that overlap with Euclid. In practice, we perform the full spectroscopic Euclid forecasts in the redshift range of 0.9≤z≤1.8, the DESI forecasts at low redshift 0<z<0.9,and the WFIRST-2.4 forecasts at high redshift 1.8<z<2.7. Moreover, between z=0.6 and z=0.9 the ELGs, LRGs, and QSOs populations from DESI also overlap. Each of these populations has a different bias and their cross-power spectra would have to be computed as well. For simplicity, we only consider the ELGs population at z={0.6,0.9}because it has the highest density of galaxies and potentially offers the best constraints on the parameters. At z<0.6,we consider the BGS population. We consider four redshift bins for the DESI survey with the following binning: 2 redshift bins for the BGS population centered at z={0.125,0.375}with a redshift size ∆z=0.25 and 2 redshift bins for the ELGs population centered at z={0.6,0.8} with ∆z=0.2. Both Euclid and WFIRST-2.4 utilise 3 redshift bins with ∆z=0.3, centered at z={1.05,1.35,1.65}and z={1.95,2.25,2.55}, respectively. The full binning procedure is summarized in Table 1. Table 2. Fiducial values for the cosmological parameters (M1 approach). Cosmological parameter Fiducial value Ωb0.05 h0.67 Ωm0.32 ns0.96 ΩDE 0.68 w0−1 wa0 σ80.83 3.2. Two approaches towards Fisher matrix cosmological constraints In this work, we follow two approaches for the Fisher matrix forecast of cosmological constraints using the galaxy clustering probe. We describe these two approaches, dubbed M1 and M2, below. 3.2.1. M1 approach In the first approach, M1, the cosmological model considered, is an extension of the ΛCDM model in which we assume a variable dark energy equation of state parameter described by the CPL parameterization. Furthermore, we allow for non-flatness by letting ΩDE vary. We also assume massless neutrinos. The full set of cosmological parameters is: θcosmo :{Ωb,h,Ωm,ns,ΩDE,w0,wa, σ8}. The cosmological parameters are the quantities of direct interest here. Ωb,h,Ωm,nsare parameters that modulate the shape of the power spectrum. Ωbis the baryon density parameter. The his the reduced Hubble constant defined as h=H0/(100 km s−1Mpc−1), where H0is the present day Hubble parameter. Ωmis the matter density parameter, and nsthe spectral scalar index. The ΩDE is the dark energy content directly linked to the curvature of the universe and (w0,wa) represent, respectively, the dark energy equation of state at z=0 and its derivative with respect to scale factor at z=0, in the so-called CPL form (Chevallier & Polarski 2001;Linder 2003). The equation of state is then given by: w(z)=w0+wa z 1+z.(1) The σ8parameter gives the amplitude of matter fluctuations in the linear regime on a scale of 8h−1Mpc. Furthermore, at each redshift bin, we consider two nuisance parameters: the logarithm of the product between the bias and σ8(z), namely, ln(bσ8), and the shot noise, Pshot(z),due to the finite self-pair counts in two-point statistics. The shot noise affects the power spectrum by adding a white systematic component increasing thus the statistical uncertainties. We consider 10 redshift bins, zi, thus, the total number of nuisance parameters amounts to 20: θnuisance ={ln(bσ8(zi)),Pshot(zi)}. The cosmological parameters fiducial values are summarized in Table 2and the nuisance parameters are given in Table 3. A52, page 3 of 19 A&A 649, A52 (2021) Table 3. Fiducial values for the nuisance parameters. Nuisance parameter Fiducial value ln(bσ8(0.125)) −0.157865 ln(bσ8(0.375)) 0.052047 ln(bσ8(0.6)) −0.830210 ln(bσ8(0.8)) −0.947904 ln(bσ8(1.05)) −0.305297 ln(bσ8(1.35)) −0.293715 ln(bσ8(1.65)) −0.302905 ln(bσ8(1.95)) −0.321325 ln(bσ8(2.25)) −0.342024 ln(bσ8(2.55)) −0.364400 Pshot(z) 0.0 Notes. They are strictly identical for both approaches: M1 and M2. The bias values are taken from DESI Collaboration (2016), Euclid Collaboration (2020), Green et al. (2011) and Spergel et al. (2013). The shot noise fiducial value is the same at each redshift bins. 3.2.2. Approach M2 In the second approach, M2, adopted in the Euclid Collaboration (2020), four quantities which determine the shape of the power spectrum are treated as independent of the cosmological parameters. We call these four quantities the shape parameters. The other quantities, the angular diameter distance, Da(z), the Hubble rate, H(z), and the growth rate of cosmic structure, f(z) are called the redshift dependent (RD) parameters. With the addition of three redshift-dependent parameters and 2 redshift-dependent nuisance parameters (the same as in M1) per redshift bin, for ten redshift bins, this amounts to 30 redshift dependent parameters and 20 nuisance parameters. The full set consists of 54 parameters (see Table 4): θshape :{ωb,h, ωm,ns} θrd+nuisance :{ln(Da(zi)),ln(H(zi)),ln(fσ8(zi)),ln(bσ8(zi)),Pshot(zi)}. The nuisance parameters are identical in the two approaches. As in the Euclid Collaboration (2020), we estimate the constraints of the physical baryon density (ωb= Ωbh2), and the physical matter density (ωm= Ωmh2) parameters. In the M2 approach, the constraints on the final cosmological parameters (θcosmo) are obtained from the 54 parameters stated above by projection. We also note that taking the logarithm of the parameters helps to stabilize matrix inversion by reducing the dynamic range between minimum and maximum eigenvalues, as described in Euclid Collaboration (2020). 3.3. The Fisher matrix formalism The Fisher matrix formalism (see Coe 2009 for a pedagogical introduction) is commonly used in cosmology to forecast the Gaussian uncertainties of a set of model parameters under some constraints. The Fisher approach relies on the assumption of Gaussian posterior distribution. The Fisher matrix corresponding to the information on cosmological parameters provided by a galaxy clustering survey is given by (Tegmark 1997;Tegmark et al. 1998): Fαβ =1 8π2Z1 −1 dµZkmax kmin k2dk ×∂ln Pobs(z,k,µ) ∂α ∂ln Pobs(z,k,µ) ∂β Veff(z,k,µ),(2) Table 4. Fiducial values of the shape and redshift dependent parameters for the approach M2. Parameter Fiducial value ωb0.022445 h0.67 ωm0.143648 ns0.96 ln(Da(0.125)) 6.177874 ln(H(0.125)) 4.268284 ln(fσ8(0.125)) −0.757651 ln(Da(0.375)) 7.00907 ln(H(0.375)) 4.411368 ln(fσ8(0.375)) −0.714737 ln(Da(0.6)) 7.264910 ln(H(0.6)) 4.548941 ln(fσ8(0.6)) −0.730402 ln(Da(0.8)) 7.378986 ln(H(0.8)) 4.672001 ln(fσ8(0.8)) −0.767895 ln(Da(1.05)) 7.452405 ln(H(1.05)) 4.821969 ln(fσ8(1.05)) −0.830894 ln(Da(1.35)) 7.488258 ln(H(1.35)) 4.992418 ln(fσ8(1.35)) −0.916615 ln(Da(1.65)) 7.494160 ln(H(1.65)) 5.150878 ln(fσ8(1.65)) −1.004727 ln(Da(1.95)) 7.483737 ln(H(1.95)) 5.297446 ln(fσ8(1.95)) −1.090957 ln(Da(2.25)) 7.463977 ln(H(2.25)) 5.432989 ln(fσ8(2.25)) −1.173509 ln(Da(2.55)) 7.438762 ln(H(2.55)) 5.558599 ln(fσ8(2.55)) −1.251760 where αand βrun over the parameters we vary. kis the total wave vector magnitude in Mpc−1,µthe cosine of the angle to the line of sight, and Veff(z,k,µ) the effective volume of the survey: Veff(z,k,µ)=Vs(z)"n(z)Pobs(z,k,µ) n(z)Pobs(z,k,µ)+1#2 ,(3) with n(z) the number density of galaxies in each redshift bin and Vs(z) the volume of each redshift bin. The observed power spectrum Pobs(z,k,µ), taking into account nonlinear effects is given by: Pobs(z,k,µ)=1 q2 ⊥qk [bσ8(z)+fσ8(z)µ2]2 1+[f(z)kµσp(z)]2! ×Pdw(z,k,µ) σ2 8(z)Fz(z,k,µ)+Pshot(z).(4) The first term represents the Alcock-Paczynski (AP) volume dilation effect that reflects spurious anisotropies in the power spectrum which arises when assuming a cosmology (that might be different from the true cosmology) to convert redshifts into distances. The coefficient qk=Href(z)/H(z) is given by the ratio A52, page 4 of 19 S. Yahia-Cherif et al.: Validating Fisher approach between the Hubble parameter in the true cosmology and the one in the reference cosmology. The coefficient q⊥=Da(z)/Da,ref(z) is given by the ratio between the angular distance in the reference cosmology and the one in the true cosmology. The term, between parentheses, contains the linear redshift space distortions due to the growth of structure and finger-of-God effects that aim to describe the nonlinear damping due to the velocity dispersion of satellite galaxies inside host halos. Then, σp= 1/(6π2)RPm(k,z)dkrepresents the linear theory velocity dispersion that can be computed from the linear matter spectrum Pm(k,z). The third term, Pdw, is the dewiggled power spectrum (Eisenstein & Hu 1998,1999) that takes into account the nonlinearities in the matter power spectrum. The fourth term, Fz, represents the uncertainties coming from the spectroscopic precision of the survey under consideration and Pshot(z) the residual shot noise. For a comprehensive description of the observed power spectrum model introduced above, we recommend the recent Euclid paper on validated forecasts (Euclid Collaboration 2020), where the specific form of the function Fzas well as the rescaling of kand µdue to the AP effect are given in Eqs. (74)–(79). In terms of nonlinearities, we should note that the model used accounts for the finger-of-God effect as well as for the damping of the BAO feature. More details can be found in Sect. 3.2.2 of Euclid Collaboration (2020). We move on to describing the estimation of the uncertainties using the Fisher matrix formalism. The inverse of the Fisher matrix gives the covariance matrix that represents the likelihood curvature evaluated at the fiducial values of αand β. According to the Cramér-Rao inequality, the square root of each covariance matrix diagonal elements gives the 1σlower bound constraint for each parameter (marginalized over all other parameters) where the posterior distribution is Gaussian: σα=p(F−1)αα.(5) In order to compute the Fisher matrix we use SpecSAF (Spectroscopic Super Accurate Forecast), a modified and improved version of a code first used in Tutusaus et al. (2016). SpecSAF computes high point stencil derivatives with a very high level of precision. This code is linked to CAMB (Lewis et al. 2000) and CLASS (Lesgourgues 2011) to compute the matter power spectra and allows the user to directly compute the FoM and plot the contours. We note that SpecSAF is one of the validated codes used in Euclid Collaboration (2020). 3.4. Precision of the Fisher matrix formalism In order to obtain reliable final constraints from the covariance matrix, we have to estimate the precision needed for the Fisher elements themselves to achieve a given precision (for instance, a relative precision lower than 10% for the square root of elements of the diagonal of the covariance matrix). This problem is related to the conditioning of the Fisher matrix that needs to be inverted. The condition number ,CN,of a matrix, A,is known to provide an upper limit to the precision, δx,of the determination of the solution, x,of the equation, Ax =b, where bis a vector. That is, the elements of band Ahave to be known with enough precision such that δbi/bi≤(1/CN)δx/x. That is, if the condition number is large, even a very small error in bwill result in a large error in x. The condition number (Belsley et al. 2005) is given by the ratio between the largest eigenvalue to the smallest eigenvalue. This is however an upper limit that can be far too stringent in practice: as an example, in our case, the condition number in approach M1 is 1.3×105, while in the M2 approach, the condition number is 4.4×1011. This formally calls for a precision of ∼10−6and ∼10−13 to ensure an accuracy of 10% in the final constraints. However, a diagonal matrix (provided that no term in the diagonal vanishes) is well-conditioned in practice, even if the condition number is very large. Therefore, for a matrix with size significantly greater that 2 ×2, the condition number does not provide useful information on the conditioning of the matrix. This observation has important consequences in practice, as the elements of the Fisher matrix are generally computed from numerical derivatives whose required precision has to be determined in order to achieve reliable results on the covariance. To better quantify the precision needed for reliable constraints we introduce a simple method that consists of “vibrating” each element of the Fisher matrix as follows: Fvibrated αβ =Fαβ(1 +N(0,1)),(6) with αand βas the indices which run over all unique pairs of parameters. N(0,1) is a normal distribution centered at 0 with variance 1 and the amplitude of the perturbations. In order to keep the Fisher matrix symmetric, we only perturb the lower triangle part of the Fisher and compute the symmetric part across the diagonal. The coefficients in the Fisher matrix are computed from Boltzmann codes and are therefore subjects to correlated noise. In our approach, this correlation is lost during the Fisher matrix vibration. In the approach M1, we apply the perturbations for 2000 amplitudes, , regularly logspaced in the range 10−12– 10−1. To ensure enough statistical data and reduce the Poisson noise, we produce 10 000 vibrated Fisher matrix per epsilon values. In the M2 approach, this procedure is more expensive in terms of computational time. We therefore reduced the number of values to 600 in the range 10−6–10−1and we produced 2000 vibrated Fisher matrices per value. The constraints obtained with each vibrated Fisher matrix are then compared with the original Fisher matrix by computing the relative error. For each parameter, we compute the number of draws that gives an agreement level on the constraints at the chosen precision (10% in our case) or better. The tolerated precision is the largest value of for which 68% of the constraints on each parameter remain lower than the chosen precision (10% in our case). The results are illustrated in Fig. 1for the cosmological parameters in the M1 approach. In the M2 approach, the size of the matrix is too large, so in Fig. 2 we only show the shape parameters and the redshift-dependent parameters for a single redshift bin (z=1.35). 3.5. Stability and convergence tests The stability and convergence of numerical derivatives can be an important issue and testing it is an essential validation step. Fortunately, testing the stability of the derivative of the observed power spectrum over any parameter is straightforward and determines the appropriate step sizes to use and which ones to avoid. However, the stability itself does not prove that a numerical derivative is correctly computed since it only shows to what extent the derivative changes as a function of step size. Ideally, a convergence test, which consists in computing the relative error between the analytical derivative and the numerical one, also has to be performed. Unfortunately, directly exploiting the convergence with a relative error test is not possible when the derivatives cannot be computed analytically. This is the case for most of the parameters we consider in this study. For the cosmological parameters, we need to use a Boltzmann code and the semi-analytical form of the derivatives over Da(z) and H(z) A52, page 5 of 19 A&A 649, A52 (2021) Fig. 1. Vibration matrix: approach M1. The values are expressed in absolute values. Each cell containing 1 has to be taken as a lower bound value of the limit because we take 1 as the threshold perturbation amplitude. contain derivatives of the matter power spectrum over the scale k (dP/dk). Here, we follow an alternative approach to test derivative convergence by considering three different numerical methods to compute the derivatives: the 3, 5,and 7 point-centered schemes. With the 3 (respectively 5, 7) scheme being the method in which we use 3 (respectively 5, 7) perturbed points of the function f(x) to determine the numerical derivative at x, we assume that the 7 points derivative is very near the true value when an appropriate step size is chosen. We compare the 7-point with the 3and 5-point stencil for different step sizes and look for the best agreement between the two derivatives; in other words, we carry out convergence tests between the 3and the 5-point stencil towards the 7 points stencil by comparing the relative error between the 7-point stencil and the 3and 5-point stencils. More generally, the convergence is always faster with the high-stencil derivatives method than with the low-stencil derivatives, which means that the optimum for the 5-point stencil often occurs at larger step values than the 3-point stencil. However, this does not hold in all cases. For instance, applying this method to a highly noisy or oscillatory function can lead to erroneous results and finding stable behaviors would not be possible. It is therefore very useful to plot the functions under consideration in order to identify any such features. 3.6. The MCMC sampling In order to consolidate the reliability of our Fisher matrix constraints, we compare them with MCMC constraints. As is well known, the MCMC sampling is much more time-consuming than the Fisher matrix computation due to the computational time needed to obtain the power spectrum at each iteration. Moreover, by construction, a Markov chain cannot be parallelized, as each draw depends on the previous displacement; however, it is possible to launch several chains at the same time to sample the parameter space more effectively. Considering a set of data and a model equipped with a set of parameters Θ, the MCMC allows us to use Bayesian inference and model comparison: P(Θ|data) =P(data|Θ)P(Θ) P(data) ,(7) with P(Θ|data) the posterior distribution, which provides the constraints and the joint covariance between all the parameters considered, P(data|Θ) the Likelihood, P(Θ) the prior that we take to be flat in the present study, and P(data) the Bayesian evidence that can be safely ignored as it is a constant of no interest in our A52, page 6 of 19 S. Yahia-Cherif et al.: Validating Fisher approach Fig. 2. Vibration matrix: approach M2. We note that some values are identical because we decreased the resolution of the perturbation amplitude range. situation. Thus, the previous equation can be reduced to: P(Θ|data) ∝P(data|Θ).(8) We compute the posterior distribution with the MCMC module of SpecSAF using the Metropolis-Hastings algorithm (Metropolis & Ulam 1949;Robert 2015). This is described in more detail in Appendix A. We test the chains convergence using the Gelman-Rubin diagnostic (Gelman & Rubin 1992;Brooks & Gelman 1998) described in Appendix B. Our data are modeled by the a realization of our fiducial model given by Eq. (4). 3.7. MCMC Comparison The full validation of the Fisher matrix approach is performed with a direct comparison with the results from the MCMC chains. We sample the two parameter spaces of the two approaches and compare the relative error given by the covariance matrix of the sampling and the inverse of the Fisher matrix. We also visualize the Fisher contours versus the MCMC sampling. As we use exactly the same specifications for both methods, we expect to get the same constraints with a relative error lower than 10% as long as the posterior distribution remains gaussian and doesn’t present any degeneracies. 4. Results 4.1. Precision requirements The vibration matrices for M1 and M2 are presented in Figs. 1 and 2, respectively. We can clearly see that the sensitivity of the constraints depends on the element vibrated. According to Fig. 1, the diagonal elements are the most sensitive, more specifically, Ωm,ΩDE and wa. For our working requested precision (10%), these three parameters require a precision between 0.06% and 0.09%; for the other parameters, the diagonal elements required precision spreads between 0.1% and 0.4%; the nondiagonal elements are less cumbersome: some of them can even change by a factor of 2 without affecting the constraints (for instance Ωb vs [ΩDE,w0,wa,σ8]). Most of the nondiagonal elements can be computed with a precision worse than 1%. The Fisher matrix produced following the M2 approach (Fig. 2) shows a greater sensitivity to vibration than the M1 approach, which is in line with their different condition numbers. The diagonal element, h,seems to be the most sensitive parameter with a required precision around 0.0012%, which is, however, much less stringent than what the condition number would suggest (∼10−11%). The background quantities and ωmrequire a precision around 0.015%, while ωband nsare one order of magnitude more tolerant, with a target precision of 0.2%. The A52, page 7 of 19 A&A 649, A52 (2021) Fig. 3. Stability and convergence tests towards the 7-point stencil for waand ln(Da(1.35)) at a given kand µat z=1.35. Upper left panel: stability of wa;upper right panel: stability of ln(Da(1.35)); bottom left panel: convergence of ωm;bottom right panel: convergence of ln(Da(1.35)). Blue: 3-point stencil, red: 5-point stencil, green: 7-point stencil. off-diagonal elements show a higher tolerance for perturbations, however, it is only the background elements versus [Ωb,ns] that demonstrate a vibration tolerance higher than one percent. To summarize, compared to the M1 approach, the Fisher matrix produced in the M2 approach is generally at least two orders of magnitude more sensitive for the most sensitive parameters and one order of magnitude more sensitive for the less sensitive ones. These final figures are poorly reflected by the condition number of each approach. The above vibration matrix approach allows us to quantify the precision requested for the elements in the Fisher matrix to achieve the requested precision on the constraints of interest, namely, the FoM in our case. This effective precision is hard to anticipate otherwise: the condition number provides a mathematical lower limit to it. However, for a simple 2×2 diagonal matrix, it is clear that the precision needed for each element is just the one sought after with regard to the constraints, while the condition number (being the product of the eigenvalues) can be arbitrary large. The precision requested for the elements in the Fisher matrix is therefore intimately related to the internal structure of the Fisher matrix. In the next section, we examine the step sizes necessary to achieve the requested accuracy on the Fisher matrix elements. 4.2. Stability and convergence In order to compute the elements of the Fisher matrix (Eq. (2)), we have to compute several integral quantities involving the derivatives. The precision on these elements is therefore directly related to the precision that can be achieved on these computations, determined by the step sizes used and the numerical schemes used in these derivations. In the following, all the parameters step sizes shown are relative to their corresponding fiducial value, except for waand Pshot and that is because their fiducial values is zero. Moreover, we chose to illustrate one shape parameter and one redshift dependent parameter in the following plots because the derivatives behavior between the shapes parameters is similar as well as the derivative behavior between the redshift-dependent parameters. We have therefore examined the precision obtained on each of these elements by first examining the precision reached when the step size varies. In order to illustrate this, we provide illustration for derivatives against step size, which are representative of general behavior. The upper panels of Fig. 3illustrate the stability of the square of the derivative of ln P(k,µ) with respect to waat (k,µ)=(0.0121,0.5) and lnDaat (k,µ)=(0.098,−0.1) (blue: 3-point stencil, red: 5-point stencil, green: 7-point stencil), at z=1.35. The upper left panel demonstrates the cosmological parameters stability behavior. The stability of the derivative is reached when the slope of the derivative value over the step size is close to zero (when the curve shows a horizontal behavior) while decreasing the step size. These derivatives are usually relatively stable in the step size range [10−4,10−1] for each derivative. The truncation errors remain small in all of the range for most of the derivatives. We note that for the low step values, instabilities occur due to rounding errors. The upper-right panel is representative of the typical stability behavior of ln(H(z)) and ln(Da(z)). For A52, page 8 of 19 S. Yahia-Cherif et al.: Validating Fisher approach 10 310 210 1 k 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 b : 3pts stencil 10 7 10 6 10 5 10 4 10 3 10 2 10 1 Optimal step size 10 310 210 1 k 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 b : 5pts stencil 10 7 10 6 10 5 10 4 10 3 10 2 10 1 Optimal step size 10 310 210 1 k 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 b : 3pts stencil 10 7 10 6 10 5 10 4 10 3 10 2 Critical step size 10 310 210 1 k 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 b : 5pts stencil 10 8 10 7 10 6 10 5 10 4 10 3 10 2 Critical step size Fig. 4. Optimal/Limit steps histograms for Ωb(M1). Upper left panel: optimal steps for the 3-point derivatives; upper right panel: optimal steps for the 5-point derivatives; bottom left panel: critical steps for the 3-point derivatives; bottom right panel: critical steps for the 3-point derivatives. step sizes higher than 3×10−3, instabilities begin to arise for the 3-point stencil. The truncation errors dominate and increase the total error budget. The instability range for the 5and 7-point derivatives is smaller: typically around 1 ×10−2. The lower panels (Fig. 3) represent the corresponding derivative convergence tests with respect to the 7-point derivatives; namely, the relative error between the 7-point stencil and the (3, 5) point stencil in (blue, red). Generally, convergence tests are carried out to compare the relative error between the analytical solution of a derivative and solutions obtained numerically when we vary the step size. The minimum relative error between the analytical form and the numerical solutions corresponds to the optimal step 6that ensures the most accurate derivative. Because we can6See Appendix Cfor a summary of the main definitions relevant for this work. not compute the analytical form of most of the derivatives presented here, we carry out convergence tests of the 5and the 3-point stencil using the 7-point stencil results as a reference since the latter is theoretically more accurate. In the bottom-left panel, the 3-point derivative achieves its convergence when the step is equal to 10−2. The optimal step for the 5-point stencil is around 1 ×10−1at most. Concerning the angular diameter distance (bottom right panel), the 3-point convergence level is achieved down to a step of 3 ×10−7, whereas for the 5-point derivatives, it is around 5×10−5. Figure 3shows that the convergence is always reached after the stability if we decrease the step size. For instance, let’s take a closer look on Fig. 3right panel (ln Da(1.35)). If we observe the stability and the corresponding convergence plots at a step of 10−1, we see that the convergence and the stability aren’t reached yet. Now if we start to decrease A52, page 9 of 19 A&A 649, A52 (2021) 0.000 0.500 1.000 / max -48.832 0.000 48.832 Ps (1.35) 0.650 0.670 0.690 h 0.130 0.144 0.157 m 0.932 0.960 0.988 ns 7.418 7.488 7.559 ln ( Da (1.35)) 4.919 4.992 5.066 ln ( H (1.35)) -1.011 -0.917 -0.822 ln ( f 8(1.35)) -0.388 -0.294 -0.200 ln ( b 8(1.35)) 0.019 0.022 0.026 b -48.832 0.000 48.832 Ps (1.35) 0.650 0.670 0.690 h 0.130 0.144 0.157 m 0.932 0.960 0.988 ns 7.418 7.488 7.559 ln ( Da (1.35)) 4.919 4.992 5.066 ln ( H (1.35)) -1.011 -0.917 -0.822 ln ( f 8(1.35)) -0.388 -0.294 -0.200 ln ( b 8(1.35)) Optimal Fisher contours: DESI + Euclid + WFIRST-2.4 Fig. 10. MCMC vs Fisher optimal (red) contours: approach M2. resulting in a much larger matrix; the final constraints are projected on the same parameter space for both approaches. We also found that the larger Fisher matrix is globally more sensitive to perturbations (as reflected by its high condition number). For instance, the Hubble parameter, h,is two orders of magnitude more sensitive in the second approach. We studied the stability of the numerical derivatives over the cosmological parameters and found that the truncation error is particularly small within the step ranges considered. Cosmological parameters remain stable in the step range from 1 ×10−1 to 3 ×10−4. Lower step values increase the rounding errors and result in higher total errors. Taking “large” steps is the safest way to perform the derivatives with respect to the cosmological parameters. The background quantities behave differently. The rounding errors remain small in all the step ranges considered. However, steps larger than 1×10−2for a 5-point stencil derivative give unstable results. Taking “large” steps for the background quantities is not safe because the differentiated function can oscillate on a scale smaller than the step size; this is exactly what must be avoided when computing the numerical derivatives. The convergence tests with respect to the 7-point stencil derivatives are performing better for the 5-point stencil. The convergence is reached between one and two orders of magnitude before the 3point stencil case, and similarly for the convergence level. The safe steps window range of the 5-point stencil is then wider. A comparison with the MCMC sampling shows the efficiency of the Fisher formalism even with non-adaptive derivatives, at least in the case where the posterior distribution is close to Gaussian. Although the optimal steps change with kand µ, choosing a unique step size for all the parameters does not produce appreciably larger errors that would significantly alter the forecasts. Furthermore, we considered a 1000 ×1000 (k,µ) grid which, multiplied by the number of redshift bins, amounts to computing 107derivatives per parameter. One of the reasons for the stability comes from the fact that the contribution of the large scales on the derivatives is very small compared to the small scales contribution. Thus, taking optimal steps from high kvalues remains the best option. We have also seen that the transition between stable and unstable contours is quite sharp. Changing a step size by a factor of 2 in the wrong direction can give unstable contours. A52, page 16 of 19 S. Yahia-Cherif et al.: Validating Fisher approach Identifying and explaining the derivatives behavior and the effects on the constraints is difficult without the tests presented in this work. Our strategy allows for a comprehensive and robust validation of the constraints derived from the Fisher matrix formalism through a numerical method. Acknowledgements. S. C. is funded by Miur through Rita Levi Montalcini project ‘Prometheus – Probing and Relating Observables with Multiwavelength Experiments To Help Enlightening the Universe’s Structure’. S. C. acknowledges support from the ‘Departments of Excellence 2018-2022’ Grant (L. 232/2016) awarded by the Italian Ministry of Education, University and Research (Miur). K. M.’s part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). I. T. acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838-C3-1-R, and the H2020 programme of the European Commission through Grant 776247. A. P. is a UK Research and Innovation Future Leaders Fellow, Grant MR/S016066/1. D. S. acknowledges financial support from the Fondecyt Regular project number 1200171. References Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526 Akeson, R., Armus, L., Bachelet, E., et al. 2019, ArXiv e-prints [arXiv:1902.05569] Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617 Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, ArXiv e-prints [arXiv:astro-ph/0609591] Amendola, L., & Sellentin, E. 2016, MNRAS, 457, 1490 Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Liv. Rev. Rel., 21, 2 Belsley, D. A., Kuh, E., & Welsch, R. E. 2005, Regression Diagnostics: Identifying Influential Data and Sources of Collinearity (John Wiley & Sons) Blake, C., Kazin, E. A., Beutler, F., et al. 2011, MNRAS, 418, 1707 Brooks, S. P., & Gelman, A. 1998, J. Comput. Graph. Stat., 7, 434 Chevallier, M., & Polarski, D. 2001, Int. J. Mod. Phys. D, 10, 213 Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2 Coe, D. 2009, ArXiv e-prints [arXiv:0906.4123] Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10 Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44 DESI Collaboration (Aghamousa, A., et al.) 2016, ArXiv e-prints [arXiv:1611.00036] Dore, O., Hirata, C., Wang, Y., et al. 2019, BAAS, 51, 341 Dunkley, J., Bucher, M., Ferreira, P. G., Moodley, K., & Skordis, C. 2005, MNRAS, 356, 925 Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 Eisenstein, D. J., & Hu, W. 1999, ApJ, 511, 5 Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191 Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 Green, J., Schechter, P., Baltay, C., et al. 2011, ArXiv e-prints [arXiv:1108.1374] Ivezi´ c, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111 Joachimi, B., & Taylor, A. N. 2011, MNRAS, 416, 1010 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193] Lesgourgues, J. 2011, ArXiv e-prints [arXiv:1104.2932] Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301 LSST Science Collaboration (Abell, P.A., et al.) 2009, ArXiv e-prints [arXiv:0912.0201] Metropolis, N., & Ulam, S. 1949, J. Am. Stat. Assoc., 44, 335 Percival, W. J., Baugh, C. M., Bland-Hawthorn, J., et al. 2001, MNRAS, 327, 1297 Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 Planck Collaboration VI. 2020, A&A, 641, A6 Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 Robert, C. P. 2015, ArXiv e-prints [arXiv:1504.01896] Sellentin, E. 2015, MNRAS, 453, 893 Sellentin, E., Quartin, M., & Amendola, L. 2014, MNRAS, 441, 1831 Seo, H.-J., & Eisenstein, D. J. 2003, ApJ, 598, 720 Spergel, D., Gehrels, N., Breckinridge, J., et al. 2013, ArXiv e-prints [arXiv:1305.5422] Square Kilometre Array Cosmology Science Working Group (Bacon, D. J., et al.) 2020, PASA, 37, e007 Tegmark, M. 1997, Phys. Rev. Lett., 79, 3806 Tegmark, M., Hamilton, A. J. S., Strauss, M. A., Vogeley, M. S., & Szalay, A. S. 1998, ApJ, 499, 555 Tutusaus, I., Lamine, B., Blanchard, A., et al. 2016, Phys. Rev. D, 94, 123515 Wang, Y. 2006, ApJ, 647, 1 Wang, Y., Percival, W., Cimatti, A., et al. 2010, MNRAS, 409, 737 Wolz, L., Kilbinger, M., Weller, J., & Giannantonio, T. 2012, JCAP, 2012, 009 A52, page 17 of 19 A&A 649, A52 (2021) Appendix A: The Metropolis-Hastings sampler Here we describe the steps to sample the parameter space from the SpecSAF MCMC module: 1. We initialize each parameter at its fiducial value and compute the corresponding observed power spectrum (Eq. (4)). 2. We draw a new parameter set from a multivariate normal distribution (mean =current parameter set) and compute the new observed power spectrum. The multivariate normal distribution is the one given by the inverse of the Fisher matrix for the same specifications (i.e., the covariance matrix). 3. We perform a χ2test in order to get the sum of the square errors: χ2 new =X z,k,µ (ln Pobs,new(z,k,µ)−ln Pobs,previous(z,k,µ))2 σp(z,k,µ)2,(A.1) with σp(k,µ,z) the errors of the logarithm of the observed power spectrum given by Seo & Eisenstein (2003): σp=2πs2 Veff(z,k,µ)k2∆k∆µ 1+n(z)Pobs,fid n(z)Pobs,fid !,(A.2) with n(z) the mean density of galaxies at a given redshift bin. 4. We finally compute the likelihood ratio: LR=exp(−0.5(χ2 new −χ2 previous)).(A.3) 5. We draw a random number RNin the interval [0,1] from a uniform distribution. 6. We compare LRand RN: – if LR≥RNthe parameter set is accepted and the new parameter set is updated. – else the parameter set is not accepted and the new parameter set is not updated. 7. We loop over 2–6. Finally we stop the chains using the Gelman-Rubin diagnostic, described in Appendix B. Appendix B: The Gelman-Rubin diagnostic The Gelman-Rubin diagnostic aims to monitor MCMC outputs for parallelized chains. The convergence test computes a constant, (R−1), proportional to the ratio of the parameters variance of one chain to the mean parameters variances of the independent chains. The mean of the empirical variance within each chain W is defined as: W=1 M(N−1) N X i=1 P X α=1 (vα i−ˆvα)2,(B.1) with Mthe number of chains, Pthe number of parameters and N the number of accepted points for each chain; ˆvαrepresents the mean value drawn of a parameter αfor one chain and vα ithe i-th αvalue drawn. The between-chain variance Bis given by: Table B.1. (R−1) values for approach M1. Parameter (R−1) value Ωb4.518e−4 h1.900e−3 Ωm3.657e−4 ns6.984e−4 ΩDE 1.807e−3 w04.007e−5 wa3.765e−4 σ81.528e−4 ln(bσ8(0.125) 1.316e−3 Pshot(0.125) 4.903e−5 ln(bσ8(0.375) 1.091e−3 Pshot(0.375) 7.877e−4 ln(bσ8(0.6) 5.615e−4 Pshot(0.6) 8.707e−5 ln(bσ8(0.8) 7.964e−4 Pshot(0.8) 1.920e−4 ln(bσ8(1.05) 7.551e−4 Pshot(1.05) 1.213e−3 ln(bσ8(1.35) 8.478e−3 Pshot(1.35) 1.458e−3 ln(bσ8(1.65) 5.064e−4 Pshot(1.65) 1.003e−3 ln(bσ8(1.95) 4.007e−4 Pshot(1.95) 1.635e−3 ln(bσ8(2.25) 2.594e−4 Pshot(2.25) 7.626e−4 ln(bσ8(2.55) 1.823e−2 Pshot(2.55) 9.062e−4 B N=1 (M−1) P X α=1 (vα i−ˆv)2.(B.2) The weighted sum of Wand Bis written as: V=N−1 NW+B N 1+1 N!.(B.3) Finally, the Rnumber is defined as: R=sˆ d+3 ˆ d+1 V W≥1,(B.4) with ˆ dthe degrees of freedom estimate of a given distribution. The convergence is reached once (R−1) <0.03 “(i.e., R must be sufficiently close to 1)”. We provide the (R−1) values obtained for the approaches M1 and M2 in Tables B.1 and B.2, respectively. A52, page 18 of 19 S. Yahia-Cherif et al.: Validating Fisher approach Table B.2. (R−1) values for approach M2. Parameter (R−1) value ωb1.811e−2 h1.820e−2 ωm1.820e−2 ns1.827e−2 ln(Da(0.125)) 1.812e−2 ln(H(0.125)) 1.827e−2 ln(fσ8(0.125)) 1.819e−2 ln(bσ8(0.125)) 1.819e−2 Pshot(0.125) 1.817e−2 ln(Da(0.375)) 1.808e−2 ln(H(0.375)) 1.817e−2 ln(fσ8(0.375)) 1.802e−2 ln(bσ8(0.375)) 1.819e−2 Pshot(0.375) 1.808e−2 ln(Da(0.6)) 1.821e−2 ln(H(0.6)) 1.825e−2 ln(fσ8(0.6)) 1.811e−2 ln(bσ8(0.6)) 1.828e−2 Pshot(0.6) 1.809e−2 ln(Da(0.8)) 1.807e−2 ln(H(0.8)) 1.810e−2 ln(fσ8(0.8)) 1.816e−2 ln(bσ8(0.8)) 1.813e−2 Pshot(0.8) 1.815e−2 ln(Da(1.05)) 1.818e−2 ln(H(1.05)) 1.814e−2 ln(fσ8(1.05)) 1.825e−2 ln(bσ8(1.05)) 1.818e−2 Pshot(1.05) 1.811e−2 ln(Da(1.35)) 1.814e−2 ln(H(1.35)) 1.819e−2 ln(fσ8(1.35)) 1.811e−2 ln(bσ8(1.35)) 1.821e−2 Pshot(1.35) 1.822e−2 ln(Da(1.65)) 1.822e−2 ln(H(1.65)) 1.815e−2 ln(fσ8(1.65)) 1.825e−2 ln(bσ8(1.65)) 1.821e−2 Pshot(1.65) 1.804e−2 ln(Da(1.95)) 1.820e−2 ln(H(1.95)) 1.820e−2 ln(fσ8(1.95)) 1.821e−2 ln(bσ8(1.95)) 1.823e−2 Pshot(1.95) 1.812e−2 ln(Da(2.25)) 1.820e−2 ln(H(2.25)) 1.820e−2 ln(fσ8(2.25)) 1.826e−2 ln(bσ8(2.25)) 1.820e−2 Pshot(2.25) 1.827e−2 ln(Da(2.25)) 1.820e−2 ln(H(2.55)) 1.816e−2 ln(fσ8(2.55)) 1.811e−2 ln(bσ8(2.55)) 1.823e−2 Pshot(2.55) 1.830e−2 Appendix C: General definitions In this last part of the appendix, we summarize the main definitions relevant for this work. Optimal step. Step size providing the minimum relative error between the analytical form and the numerical solutions. It ensures the most accurate derivative. Critical step. Step size that can potentially lead to an error of 10% on the final constraints. Hypercritical step. Step size that can potentially lead to catastrophic errors (several dozen percent) on the final constraints. Truncation error. Error introduced by considering too large step sizes. The approximation of computing the derivative with a finite amount of terms is in this case no longer valid. Rounding error. Error introduced by considering too small step sizes. The finite precision floating point numbers used on computers cannot distinguish between numbers that are too close. Vibration matrix. Matrix providing the largest value of the vibration for which 68% of the constraints on each parameter remain smaller than the chosen precision. A52, page 19 of 19