© ESO 2019

1. Introduction

The past few decades have witnessed the advent of an array of galaxy redshift surveys, with the state-of-the-art catalogues mapping millions of galaxies with precision positioning and accurate redshifts. The Sloan Digital Sky Survey (SDSS; York et al. 2000; Abazajian et al. 2009; Ahn et al. 2014; Alam et al. 2015) and the Six Degree Field Galaxy Redshift Survey (6dFGRS; Jones et al. 2009) are two notable examples. Future cutting-edge surveys from the Euclid (Laureijs et al. 2011; Racca et al. 2016; Amendola et al. 2018) and Large Synoptic Survey Telescope (LSST; Ivezic et al. 2008) missions, currently under construction, further highlight the wealth of galaxy redshift data sets which would be available within a five to ten year time frame. Sophisticated and optimal data analysis techniques, in particular large-scale structure analysis methods, are in increasing demand to cope with the present and upcoming avalanches of cosmological and astrophysical data, and therefore optimise the scientific returns of the missions.

With the metamorphosis of cosmology into a precision (and data-driven) science, the three-dimensional (3D) large-scale structures have emerged as an essential probe of the dynamics of structure formation and evolution to further our understanding of the Universe. The two-point statistics of the 3D matter distribution have developed into key tools to investigate various cosmological models and test different inflationary scenarios. Various techniques to measure the power spectrum and several reconstruction methods attempting to recover the underlying density field from galaxy observations are described in literature (e.g. Bertschinger & Dekel 1989; Bertschinger et al. 1991; Hoffman et al. 1994; Lahav et al. 1994; Fisher et al. 1995; Sheth 1995; Webster et al. 1997; Bistolas & Hoffman 1998; Schmoldt et al. 1999; Saunders et al. 2000; Zaroubi et al. 1999; Zaroubi 2002; Erdoǧdu et al. 2004, 2006), with the recent focus being on large-scale Bayesian inference methods (e.g. Kitaura & Enßlin 2008; Kitaura et al. 2009; Jasche & Kitaura 2010; Jasche et al. 2010a; Jasche & Wandelt 2012, 2013b; Jasche & Lavaux 2015, 2018). A formal and rigorous Bayesian framework provides the ideal setting to solve the ill-posed problem of inferring signals from noisy observations, while quantifying the corresponding statistical uncertainties.

The potential of such Bayesian algorithms to jointly infer cosmological constraints, nevertheless, has not yet been exploited. We present, for the first time, a non-linear Bayesian inference framework for cosmological parameter inference from galaxy redshift surveys via an implementation of the Alcock-Paczyński (AP, Alcock & Paczynski 1979) test. We extend the hierarchical Bayesian inference machinery of Bayesian Origin Reconstruction from Galaxies (BORG; Jasche & Wandelt 2013a), originally developed for the non-linear reconstruction of large-scale structures, to constrain cosmological parameters. BORG encodes a physical model for gravitational structure formation, yielding a highly non-trivial Bayesian inverse problem. This consequently allows us to reformulate the standard problem of present 3D density field reconstruction as an inference problem for initial conditions at an earlier epoch from current galaxy observations. BORG builds upon the implementation of the Hamiltonian Monte Carlo (HMC) method (Neal 1993), initially introduced in the HAmiltonian Density Estimation and Sampling (HADES) algorithm (Jasche & Kitaura 2010), for efficiently sampling the high dimensional and non-linear parameter space of possible initial conditions at an earlier epoch.

In this work, the conceptual framework is to constrain the comoving-redshift coordinate transformation and therefore infer the appropriate cosmology which would result in isotropic correlations of the galaxy density field. The key aspect of this application of the AP test consequently lies in its robustness to a misspecified model and the approximations therein, yielding a near-optimal exploitation of the model predictions, without relying on its accuracy in modelling the scale dependence of the correlations of the density field. Here, we employ Lagrangian perturbation theory (LPT) as a physical description for the non-linear dynamics and perform a joint inference of initial conditions, and consequently the corresponding non-linearly evolved density fields and associated velocity fields, and cosmological parameters, from incomplete observations. This augmented framework with cosmological applications is designated as ALTAIR (ALcock-Paczyński consTrAIned Reconstruction).

The paper is organised as follows. In Sect. 2, the underlying principles of the AP test are outlined, followed by a description of the forward modelling approach and data model implemented in Sect. 3. We then test the algorithm in Sect. 5 on an artificially generated galaxy survey, with the mock generation procedure described in the preceding Sect. 4, by investigating its performance via statistical efficiency and consistency tests. In Sect. 6, we summarise the main aspects of our work and discuss further possible extensions to our algorithm in order to fully exploit its potential in deriving cosmological constraints. In Appendix A, we describe the LPT-Poissonian posterior implemented in this work, followed by the computation of the Jacobian of the comoving-redshift transformation in Appendix B. We provide a brief overview of the Hamiltonian sampling approach in Appendix C, and follow up by deriving the required equations of motion in Appendix D, with the numerical implementation outlined in Appendix E. We subsequently describe how we increase the efficiency of our cosmological parameter sampler via a rotation of the parameter space in Appendix F. Finally, we outline the derivation of the adjoint gradient for a generic 3D interpolation scheme in Appendix G.

2. The Alcock-Paczyński test

The Alcock-Paczyński (AP) test (Alcock & Paczynski 1979) is a cosmological test of the expansion of the Universe and its geometry. The main advantage of this test is that it is independent of the evolution of galaxies but depends only on the geometry of the Universe. The assumption of incorrect cosmological parameters in data analysis produces distortions in the appearance of any spherical object or isotropic statistical distribution. The AP test provides a pathway to exploit this resulting spurious anisotropy to constrain the cosmological parameters. Here, we invoke the AP test to ensure that the underlying geometrical properties of isotropy of the Universe (Friedmann 1922, 1924; Lemaître 1927, 1931, 1933; Robertson 1935; 1936a,b; Walker 1937; Saadeh et al. 2016) are maintained. As such, the key underlying assumption adopted in this work relies purely on the geometrical symmetries of the cosmological principle. As a result, such a test does not employ the growth of structures to constrain cosmology, unlike cluster abundance (e.g. Wang & Steinhardt 1998).

The AP test, and various formulations thereof, have been studied extensively in the context of galaxy and quasar surveys (e.g. Phillipps 1994; Ryden 1995; Ballinger et al. 1996; Matsubara & Suto 1996; Popowski et al. 1998; de Laix & Starkman 1998; López-Corredoira 2014; Li et al. 2014, 2016). Variants of the AP test have also been successfully applied to cosmic voids (e.g. Sutter et al. 2012, 2014; Lavaux & Wandelt 2012; Hamaus et al. 2014, 2015, 2016) and also to other cosmological observables like supernovae (Blake et al. 2011), the Lyman-α forest (Hui et al. 1999) and 21 centimetre emission maps (Nusser 2005; Barkana 2006).

With baryon acoustic oscillations (BAOs) being a robust standard ruler, the AP test has been utilised for the simultaneous measurement of the Hubble parameter and angular diameter distance of distant galaxies (e.g. Seo & Eisenstein 2003; Blake & Glazebrook 2003; Glazebrook & Blake 2005; Padmanabhan & White 2008; Shoji et al. 2009). In Fig. 1, we depict the 1σ confidence region of the cosmological constraints inferred via our implementation of the AP test. As a comparison, we also indicate the corresponding confidence region obtained via BAO measurements from the SDSS-III (Data Release 12; Alam et al. 2017). These BAO constraints have not been combined with Planck measurements, which would significantly tighten the constraints. Nevertheless, this highlights the significant potential constraining power of our AP test compared to standard BAO analyses, while being at least as robust. While this improvement is extremely substantial for the mock SDSS-III survey considered here, we will investigate to what extent the above promise holds when applied to actual SDSS-III data in a follow-up work, as unknown systematics represent a potential caveat. We discuss Fig. 1 in more depth in Sect. 5.

Fig. 1. Comparison of cosmological constraints from BAO measurements and our implementation of AP test in ALTAIR. The grey and green lines denote the 1σ confidence region, centred on the fiducial cosmological parameters, obtained from our AP test and BAO constraints from SDSS-III (DR 12; Alam et al. 2017), respectively. The BAO constraints have not been combined with Planck CMB measurements. This demonstrates the potentially unprecedented constraining power of our AP test compared to standard BAO analyses, as discussed in Sect. 2, with the inset focusing on the ALTAIR constraints where the fiducial cosmology is depicted in dashed lines. This error forecast is validated on a simulated analysis (cf. Fig. 6).

The crucial aspect of our AP test is that it does not assume that the correlation function is correctly modelled. This robustness to a misspecified model is illustrated explicitly in Sect. 5, where we demonstrate that the shape of the prior power spectrum adopted in the inference framework does not impact on the inferred cosmological constraints (cf. Fig. 10). As a result of this robustness, our AP test has a definite edge over standard approaches. Moreover, it has been pointed out that other cosmological tests, such as the luminosity distance – redshift relation, can be considered as generalised formulations of the AP test (Mukherjee & Wandelt 2018), further underlining the strength of the approach presented in this work.

3. The forward modelling approach

The large-scale structure (LSS) posterior implemented in this work, based on the BORG framework (Jasche & Wandelt 2013a), is described in depth in Appendix A. A key component of the inference framework is the forward model ℳ_p which links the initial conditions ${\delta }_p^{\mathrm{ic},(r)}$ to the redshift space representation of the evolved density field ${\delta }_p^{\mathrm{f},(z)}$ as follows:

$$\begin{array}{cc}\hfill {\delta }_p^{\mathrm{f}(z)}& ={\mathcal{M}}_p(1+{\delta }_p^{\mathrm{ic}(r)})={\mathcal{M}}_p^{(1)}°{\mathcal{M}}_p^{(2)}(1+{\delta }_p^{\mathrm{ic}(r)})={\mathcal{M}}_p^{(1)}\left({\rho }_p^{\mathrm{f}(r)}\right)\hfill \\ \hfill & ={\mathcal{J}}_p\left({\displaystyle \sum _{i,j}{\mathcal{E}}_{\mathit{ij}}^{-1}{\rho }_j^{\mathrm{f}(r)}{x}_p^{\alpha (i)}{y}_p^{\beta (i)}{z}_p^{\gamma (i)}}\right)-1,\hfill \end{array}$$(1)

where ${\rho }_p^{\mathrm{f}(r)}\equiv 1+{\delta }_p^{\mathrm{f}(r)}$ is the final density field in comoving space. The forward model consists of two components, ${\mathcal{M}}_p={\mathcal{M}}_p^{(1)}°{\mathcal{M}}_p^{(2)}$. The first component, ${\mathcal{M}}_p^{(2)}\equiv {\mathcal{G}}_p(a,\{{\delta }_p^{\mathrm{ic}}\})$, contains a physical description of the non-linear dynamics, and consequently propagates the initial conditions forward in time using LPT, yielding a non-linearly evolved final density field in comoving space, ${\delta }_p^{\mathrm{f}(r)}$.

To encode the AP test, we incorporate another component in the forward model that takes care of the coordinate transformation from comoving (r) to redshift (z) space, encoded in ${\mathcal{M}}_p^{(1)}$ (cf. Fig. 2). Schematically, we construct a second grid in redshift space, which involves a triquintic interpolation (fifth order interpolation scheme in three dimensions) on the comoving grid. This interpolation scheme is described in Appendix G, with the notation employed in Eq. (1) clearly laid out. The corresponding Jacobian factor of this transformation, $|{\mathcal{J}}_{\mathcal{r}}^{\mathcal{z}}|$ (cf. Appendix B), entails cosmological dependence and is consequently included in the AP test as well as through the direct coordinate dependence ℰ_ij.

Fig. 2. Schematic representation of the reconstruction pipeline. The forward model consists of a chain of various components for the non-linear evolution from initial conditions and the subsequent transformation from comoving to redshift space for the application of the AP test. This consequently transforms the initial density field into a set of predicted observables, i.e. a galaxy distribution in redshift space, for comparison with data via a likelihood or posterior analysis.

The redshift space representation then allows for comparison with data via the likelihood or posterior. The essence of this AP test to constrain cosmological parameters can be summarised as follows: The Bayesian inference machinery explores the various cosmological expansion histories and selects the cosmology-dependent evolution pathways which result in isotropic correlations of the galaxy density field.

Figure 2 illustrates the reconstruction scheme implemented in ALTAIR. First, galaxies are projected from the survey onto a 3D grid, such that we have a distribution of galaxies in redshift space and this constitutes our observable. We then generate a 3D density field according to Gaussian initial conditions (homogeneous prior) with a reference power spectrum, typically ΛCDM cosmology. The forward model subsequently transforms the initial density field into a set of predicted observables which are then compared to data via a likelihood or posterior analysis. And conversely, given the position of galaxies, we can infer this density field.

While we implement LPT to approximately describe gravitational non-linear structure formation in this work, other more adequate physical descriptions such as 2LPT or the non-perturbative particle mesh (see recent upgrade of BORG in Jasche & Lavaux 2018) can be straightforwardly employed, within the flexible block sampling approach described in Appendix A.3, by upgrading the first component, ${\mathcal{M}}_p^{(2)}\equiv {\mathcal{G}}_p(a,\{{\delta }_p^{\mathrm{ic}}\})$, of our forward model (cf. Fig. 2). Nevertheless, our implementation of the AP test exploits essentially the isotropy of the correlation function, such that there is no explicit dependence on the accuracy of modelling the scale dependence of the correlations, rendering this method robust to a misspecified model and the approximations therein.

3.1. The galaxy data model

Galaxies can be considered as tracers of the matter fluctuations since they follow the gravitational potential of the underlying matter distribution, with the statistical uncertainty due to the discrete nature of the galaxy distribution usually modelled by a Poissonian distribution (e.g. Layzer 1956; Peebles 1980). Poissonian likelihoods have emerged as the standard for non-linear LSS inference (e.g. Jasche et al. 2010b; Jasche & Kitaura 2010; Kitaura et al. 2010). The Poissonian likelihood distribution implemented in ALTAIR, for multiple subcatalogues or galaxy observations labelled by the index g, can be expressed as follows:

$$\begin{array}{c}\hfill \mathcal{L}\left(\{N_p^{\mathit{g}}\}|{\lambda }_p^{\mathit{g}}\right)={\displaystyle \prod _p}\frac{{({\lambda }_p^{\mathit{g}})}^{N_p^{\mathit{g}}}{\mathrm{e}}^{-{\lambda }_p^{\mathit{g}}}}{N_p^{\mathit{g}}!},\end{array}$$(2)

where $N_p^{\mathit{g}}$ is the observed galaxy number counts in redshift space in the given voxel p. ${\lambda }_p^{\mathit{g}}$ is the expected number of galaxies at this given position and is related to the final density field ${\delta }_p^{\mathrm{f}}$, in redshift space, via

$$\begin{array}{c}\hfill {\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\})=R_p^{\mathit{g}}({\theta }_i)\mathcal{T}[1+{\delta }_p^{\mathrm{f}}],\end{array}$$(3)

where $R_p^{\mathit{g}}$ is the overall linear response operator of the survey that incorporates the survey geometry and selection effects, and θ_i corresponds to a set of cosmological parameters. 𝒯 is the galaxy biasing model which accounts for the fact that galaxies do not trace exactly the underlying matter distribution, and are therefore biased tracers with clustering properties that do not exactly mirror those of dark matter (Kaiser 1984). This is currently one of the most challenging and unresolved issues hindering the analysis of galaxy distributions in non-linear regimes (e.g. see the review by Desjacques et al. 2018; Schmidt et al. 2018).

In this work, we adopt the standard approach of a local, but non-linear bias function, in particular, the phenomenological model proposed by Neyrinck et al. (2014), such that the above Poisson intensity field can be expressed as

$$\begin{array}{c}\hfill {\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})=R_p^{\mathit{g}}{\overline{N}}^{\mathit{g}}{[1+{\delta }_p^{\mathrm{f}}]}^{\beta }{\mathrm{e}}^{-{\rho }^{\mathit{g}}{[1+{\delta }_p^{\mathrm{f}}]}^{-{ϵ}^{\mathit{g}}}},\end{array}$$(4)

where the bias function, described by four parameters, N̄^g, the mean density of tracers, and $\{b_i^{\mathit{g}}\}=\{\beta ,{\rho }^{\mathit{g}},{ϵ}^{\mathit{g}}\}$, is a truncated power law bias model with the additional exponential function suppressing galaxy clustering in under dense regions. This bias model, with a power law and an exponential at low densities, were found to be in good agreement with standard excursion set and local-growth-factor models (for more details, see Neyrinck et al. 2014). The main limitation of this bias model is that it is purely local. Nevertheless, it is more adequate than a simplistic linear bias model and mitigates in practice the deficiencies of our physical model (LPT) at the considered resolution. The expected number of galaxies can subsequently be related to the initial conditions ${\delta }_p^{\mathrm{ic}}$ via the forward model ℳ_p, as described above, due to the deterministic nature of structure formation, as implied by the Dirac delta function in Eq. (A.2).

The logarithm of the likelihood from Eq. (2) can therefore be expressed, in terms of the initial conditions, as

$$\begin{array}{cc}& ln\mathcal{L}[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\}),\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}]\hfill \\ \hfill & =-{\displaystyle \sum _p\{{\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})}\hfill \\ \hfill & -N_p^{\mathit{g}}ln\left[{\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})\right]+ln\left(N_p^{\mathit{g}}!\right)\}.\hfill \end{array}$$(5)

We therefore have a likelihood distribution that encodes the statistical process describing the generation of galaxy observations given a specific realisation of 3D initial conditions. This data model is inherently non-linear as a result of the galaxy biasing model employed and also due to the signal dependence of Poissonian noise, which does not behave as an additive nuisance.

3.2. The augmented joint posterior distribution

The augmented joint posterior distribution corresponds to the following:

$$\begin{array}{cc}& \mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\},\{{\theta }_i\}|\{N_p^{\mathit{g}}\},\mathbf{S})\hfill \\ \hfill & \propto \mathcal{L}[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\}),\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}]\hfill \\ \hfill & \times \mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})\mathrm{\Pi }(\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})\mathrm{\Pi }(\{{\theta }_i\}),\hfill \end{array}$$(6)

where the Π’s correspond to the respective priors for each parameter. Hence, given our forward model ℳ, which incorporates sequential components of structure formation and coordinate transformation, as described above, the complex task of modelling accurate priors for the statistical behaviour of present-day matter fluctuations can be recast into a Bayesian inference problem for the initial conditions. From this joint posterior distribution, we can construct the various conditional posterior distributions for each parameter of interest. The modular statistical programming approach adopted in outlined in Appendix A.3.

Since the non-linear LSS analysis has been reformulated as an initial conditions statistical inference problem, as described by the joint posterior distribution Eq. (6), this method depends solely on forward evaluations, and consequently has a definite edge over traditional approaches of initial conditions inference that require backward integration of the equations of motion or the inversion of the flow of time (e.g. Nusser & Dekel 1992). The latter methods are prone to erroneous fluctuations in the initial density and velocity fields, resulting from spurious growth of decaying modes. Moreover, such schemes are hindered by survey incompleteness which requires the knowledge of the complex and, as yet, unknown multi-variate probability distribution for the matter fluctuations, to render the backward integration of non-linear models physically meaningful via constrained realisations. In comparison, the forward modelling approach here conveniently accounts for survey masks and statistical uncertainties in the initial conditions, which amounts to modelling straightforward uncorrelated Gaussian processes, to generate data constrained realisations of the initial and evolved density fields.

3.3. The cosmological parameter posterior distribution

In this work, we sample the present-day values of matter density and dark energy equation of state parameters, {θ_i}={Ω_m, w₀}, via the following conditional posterior distribution,

$$\begin{array}{cc}& \mathcal{P}(\{{\theta }_i\}|\{N_p^{\mathit{g}}\},\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\},\mathbf{S})\hfill \\ \hfill & =\mathcal{P}(\{{\theta }_i\}|\{N_p^{\mathit{g}}\},\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}),\hfill \end{array}$$(7)

assuming conditional independence of the cosmological power spectrum, that is, Fourier transform of the covariance matrix S, once the density field is known. This assumption holds as we are only probing the cosmological expansion in this work, with the power spectrum anchored with a fiducial cosmology. We quantitatively demonstrate the validity of this assumption in Sect. 5 by comparing the entropy of prior information against that of posterior information (cf. Fig. 9). We defer power spectrum sampling to a future work. Applying Bayes’ identity, and using the joint posterior distribution from Eq. (6), we obtain

$$\begin{array}{cc}\hfill \mathcal{P}(\{{\theta }_i\}& |\{N_p^{\mathit{g}}\},\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})\hfill \\ \hfill & =\frac{\mathcal{P}(\{{\theta }_i\},\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}|\{N_p^{\mathit{g}}\})}{\mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}|\{N_p^{\mathit{g}}\})}\hfill \\ \hfill & =\frac{\mathcal{L}[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\}),\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}]}{\mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}|\{N_p^{\mathit{g}}\})}\hfill \\ \hfill & \times \mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})\mathrm{\Pi }(\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})\mathrm{\Pi }(\{{\theta }_i\})\hfill \\ \hfill & \propto \mathcal{L}[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\}),\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}]\times \mathrm{\Pi }(\{{\theta }_i\}),\hfill \end{array}$$(8)

after omitting the terms without any cosmological parameter dependence.

Since this work attempts to demonstrate the capabilities of ALTAIR to constrain the cosmological parameters, as proof of concept, we set uniform prior distributions on θ_i. To sample from the above marginal posterior distribution, we make use of a slice sampling procedure (e.g. Neal 2000, 2003). After obtaining a realisation of θ_i, we need to update the comoving-redshift coordinate transformation in the forward model.

We adopt a dynamical dark energy model, in particular, the standard Chevallier-Polarski-Linder (CPL) parameterisation (Chevallier & Polarski 2001; Linder 2003), where the evolution of the dark energy equation of state parameter w is a linear function of the scale factor a, as follows:

$$\begin{array}{c}\hfill \mathit{w}={\mathit{w}}_0+(1-a){\mathit{w}}_a.\end{array}$$(9)

In this work, we set w_a = 0 and infer the present-day value w₀. Moreover, we impose the assumption of flatness, that is, Ω_k = 0, such that the dark energy density is Ω_de = 1 − Ω_m.

Due to the correlation between Ω_m and w₀, we perform a rotation of the (Ω_m, w₀) parameter space, using orthonormal basis transformations derived from the covariance matrix, to improve the efficiency of the slice sampler. This procedure is outlined in Appendix F, with the significant gain in efficiency illustrated in Fig. F.1. The corresponding sampler for the bias parameters is outlined in Appendix A.4.

4. Generation of a mock galaxy catalogue

We describe the generation of an artificial galaxy survey using as template the Sloan Digital Sky Survey (SDSS-III), consisting of four galaxy subcatalogues, to validate the methodology described in the previous sections. The procedure implemented here for the mock generation is essentially based on the descriptions provided in Jasche & Kitaura (2010) and Jasche & Wandelt (2013a).

At first, we generate a realisation for the initial density contrast ${\delta }_p^{\mathrm{ic}}$ from a normal distribution with zero mean and covariance corresponding to an underlying cosmological power spectrum for the matter distribution. This power spectrum, including baryonic wiggles, is computed following the prescription provided in Eisenstein & Hu (1998, 1999), assuming a standard Λ cold dark matter (ΛCDM) cosmology with the set of cosmological parameters (Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.0486, h = 0.6774, σ₈ = 0.8159, n_s = 0.9667) from Planck (Planck Collaboration XIII 2016). This yields a 3D Gaussian initial density field in an equidistant Cartesian grid with N_side = 128, consisting of 128³ voxels, where each voxel corresponds to a discretised volume element with a physical voxel size of 15.6 h⁻¹ Mpc, and comoving box length of 4000 h⁻¹ Mpc. As a result, this implies a total of ∼2.1 × 10⁶ inference parameters, corresponding to the amplitude of the primordial density fluctuations at the respective grid nodes.

To ensure a sufficient resolution of inferred final density fields, we oversample the initial density field by a factor of eight, thereby evaluating the LPT model with 256³ particles in every sampling step. The following step is to scale this 3D distribution of initial conditions to a cosmological scale factor of a_init = 0.001 via multiplication with a cosmological growth factor D⁺(a_init). These initial conditions are subsequently evolved forward in time, in non-linear fashion, with LPT providing an approximate description of gravitational LSS formation. We then construct a final 3D non-linearly evolved density field ${\delta }_p^{\mathrm{f}}$ from the resulting particle distribution via the cloud-in-cell (CIC) method (e.g. Hockney & Eastwood 1988).

To generate the mock galaxy survey, we essentially need to simulate the inhomogeneous Poissonian process described by Eq. (2), by drawing random samples from the distribution, on top of the final density field ${\delta }_p^{\mathrm{f}}$. In this work, we generate a mock data set with realistic features emulating the highly structured survey geometry and selection effects of the SDSS-III survey, with the observed sky completeness depicted in the left and right panels of Fig. 3, respectively, for the CMASS and LOW-Z components of the SDSS-III survey. To account for the different selection effects in the northern and southern galactic planes, each component is further divided into two subcatalogues, corresponding to north/south galactic caps (NGC/SGC). The respective radial selection functions for these four subcatalogues are illustrated in Fig. 4. Here, the selection functions are numerical estimates obtained by binning the corresponding distribution of tracers N(d) in the CMASS and LOW-Z components (e.g. Ross et al. 2017), where d is the comoving distance from the observer.

Fig. 3. Observed sky completeness, used to generate and analyse the mock catalogue in this work, are illustrated in the left and right panels, corresponding to the CMASS and LOW-Z components of the SDSS-III survey, respectively.

Fig. 4. Radial selection functions for the CMASS sample, in solid and dashed lines for the north galactic cap (NGC) and south galactic cap (SGC), respectively. The corresponding radial selection functions for the LOW-Z sample are depicted in dash-dotted (NGC) and dotted lines (SGC). These selection functions are used to generate the mock data to emulate features of the actual SDSS-III BOSS data.

The projection of the completeness functions into the 3D volume produces the 3D observation mask $C_p^{\mathit{g}}$. The survey properties are described by the galaxy selection function and the 3D completeness function, with the product of these two functions yielding the survey response operator $R_p^{\mathit{g}}$. More specifically, it is the average of the product of the 2D survey geometry $C^{\mathit{g}}(\widehat{x})=C_{p(\widehat{x})}^{\mathit{g}}$ and the selection function f(x) at each volume element of the 3D grid:

$$\begin{array}{c}\hfill R_p^{\mathit{g}}=\frac{1}{|{\mathcal{V}}_p|}{\displaystyle {\int }_{{\mathcal{V}}_p}{\mathrm{d}}^3\mathit{x}{C}^{\mathit{g}}(\widehat{x})f^{\mathit{g}}(x)},\end{array}$$(10)

where the volume occupied by the pth voxel is indicated by 𝒱_p, with |𝒱| the volume of the set 𝒱.

Finally, we generate the four artificial galaxy subcatalogues, labelled by g, by Poisson sampling on the grid with ${\overline{N}}_p^{\mathit{g}}=\{110.42,122.94,71.43,205.48\}$, resulting in a total number of 997828 galaxies. The values of ${\overline{N}}_p^{\mathit{g}}$ are chosen such that the mock catalogue reflects the characteristics of the actual SDSS-III data which contains around one million tracers.

5. Results

We perform a series of tests to evaluate the performance of our algorithm in a realistic context by applying ALTAIR on the simulated SDSS-III-like galaxy catalogue described in the previous section. In particular, the focus is on the burn-in and convergence behaviour of our method, which are key indicators of the overall numerical feasibility and statistical efficiency for real data applications. To validate the conceptual framework for cosmological parameter inference and the robustness of our implementation of the AP test, the reproducibility of the input cosmology and the correlations with the other inferred quantities such as galaxy bias are also of interest.

The Markov chains for the cosmological parameters, displayed in Fig. 5, were initialised with an over-dispersed state. This figure consequently illustrates an initial burn-in phase, lasting ∼250 MCMC steps, where the Markov chains follow a persistent drift towards the high probability region of the parameter space. The rotation of the (Ω_m, w₀) parameter space before slice sampling, as described in Appendix F, reduces the burn-in period significantly, by roughly a factor of five, as shown in the right panel of Fig. F.1, resulting in improved sampling efficiency.

Fig. 5. MCMC chains for the cosmological parameters, for the first 1000 samples, with the reference cosmology employed in the mock generation indicated by the horizontal dashed lines. An initial burn-in phase lasting ∼250 Markov transitions is illustrated by the coherent drift of the Markov chain towards the preferred region in parameter space.

The corresponding marginal and joint posterior distributions for the cosmological parameters are displayed in Fig. 6, demonstrating the capability of ALTAIR to infer tight constraints from galaxy redshift surveys. This robust AP test fully exploits the high information content from the cosmic expansion as a result of probing a deep redshift range, where the distortion is more pronounced, yielding the following cosmological constraints: Ω_m = 0.3080 ± 0.0036 and w₀ = −0.998 ± 0.008. As a comparison, the SDSS-III (DR12, BAO + Planck) constraints are as follows: Ω_m = 0.310 ± 0.005 and w₀ = −1.01 ± 0.06 (Alam et al. 2017), further highlighting the significant constraining power of our AP test. We acknowledge the significant difference in the size of uncertainties. A back of the envelope computation of the information gain is as follows: Considering a sphere of 100 h⁻¹ Mpc for BAO against all voxels in 4000 h⁻¹ Mpc, N_BAO = 4000³/(4π × 100³/3) = 15 278, compared to N_vox = 128³, yields an improvement of $\sqrt{N_{}\text{vox}/N_{}\text{BAO}}=11.7$, which provides an order of magnitude of our uncertainties on the cosmological parameters. This is an approximate attempt to quantify the information gain from including smaller scales (in our work, ∼0.17 h⁻¹ Mpc) than the BAO scale by essentially counting the number of modes. However, the above argument does not imply that employing finer resolutions will result in an infinite gain of information. There is a saturation of information at a certain resolution due to the slow variation in the density fields across neighbouring voxels, such that further refinement beyond this limit will not yield any additional information.

Fig. 6. Marginal posteriors for Ωm (left panel) and the dark energy equation of state, w0 (middle panel), for ∼3000 MCMC realisations, ignoring the burn-in phase of ∼250 Markov steps. The corresponding mean and standard deviation for each parameter are indicated in the top right corner of each plot. The joint posterior (right panel) for Ωm and w0, depicting the high level of correlation between these two parameters. The highly informative distortion due to the cosmic expansion, as a result of probing a deep redshift range, yields extremely tight constraints on the above cosmological parameters. As a consistency test, this validates our implementation of the AP test to correctly recover the input cosmology.

To verify that the cosmological information stems purely from the geometric distortion due to the cosmic expansion, that is, the AP test, we perform the following experiment: We deactivate the cosmic expansion component in our forward model and sample the cosmological parameters using only the LPT evolved density field. We consequently recover the corresponding prior distributions for the marginal posteriors of the cosmological parameters. As a result, this implies that the information derives purely from the geometry and not from the clustering of the non-linearly evolved density field, at least for the test case with the physical voxel size considered in this work.

The mean initial and final density fields computed from ∼4000 realisations (ignoring the burn-in phase) are illustrated in the left panel of Fig. 7 for a particular slice of the 3D field. The maps reveal that on average we can recover highly detailed and well-defined structures in the observed regions. In particular, the filamentary nature of the non-linearly evolved density field can be clearly seen, while the Gaussian nature of the initial conditions can also be deduced visually. However, in poorly or not observed regions, the ensemble mean density field approaches the cosmic mean density, as expected, since regions lacking any observational information should on average reflect the cosmic mean. The uncertainty in these regions is accurately accounted for in the inference process, as demonstrated by the right panel, since each sampled density field is a constrained realisation, which means that these regions are augmented with statistically correct information.

Fig. 7. Mean and standard deviation maps for the initial (top panels) and final density fields (bottom panels), computed from the MCMC realisations, with the same slice through the 3D density fields being illustrated above. In unobserved or masked regions, the density fields are not constrained by data, and they average out to the cosmic mean density. However, in observed regions, the Gaussian nature of the initial conditions and the filamentary nature of the non-linearly evolved density field is manifest. The corresponding variance is therefore higher in regions devoid of data.

In the top panel of Fig. 8, we illustrate the power spectra reconstructed from all the realisations of 3D initial density field, obtained from the posterior via the HMC sampler. This is a self-consistency test to confirm that the sampled density fields are in agreement with the reference (prior) power spectrum adopted for the mock generation, as substantiated quantitatively in the bottom panel, where the ratio of the a posteriori power spectra to the prior distribution is shown. The measured power spectra therefore demonstrate that the individual realisations possess the correct power throughout the entire domain of Fourier modes considered in this work. Moreover, they do not exhibit any spurious power artefacts typically resulting from erroneous treatment of survey characteristics, such as survey geometry and selection effects, and galaxy biases, implying that such effects have been properly accounted for.

Fig. 8. Top panel: reconstructed power spectra from the inferred posterior initial density field realisations for all sampling steps of the Markov chain. The MCMC samples are colour-coded according to their sample number, as indicated by the colour bar. This is a self-consistency test to verify whether the sampled density fields are in accordance with the reference power spectrum, depicted in dashed lines, employed in the mock generation. Bottom panel: deviation from the true underlying solution, with the thick dashed lines representing the prior power spectrum. The thin grey and black dotted lines correspond to the Gaussian 1σ and 2σ limits, respectively.

In order to verify the impact of the prior power spectrum on the actual inference of cosmological parameters, we illustrate, in Fig. 9, the distributions of power spectra computed using the inferred cosmology and the prior analytic prescription (Eisenstein & Hu 1998, 1999) and the reconstructed power spectra from Fig. 8, via their respective summary statistics, normalised with respect to the fiducial power spectrum. The scatter in the latter a posteriori power spectra reconstructed from the sampled density field realisations is significantly higher than distribution of the former prior spectra. This implies that the entropy of prior information is much lower than that of posterior information, thereby demonstrating that the prior power spectrum does not influence the cosmological parameter inference via the AP test and justifying the assumption made in Sect. 3.3.

Fig. 9. Summary statistics of the reconstructed power spectra from the inferred posterior initial density field realisations depicted in Fig. 8, with the ensemble mean indicated by a solid green line. The solid blue line corresponds to the ensemble mean of the power spectra realisations generated using the inferred cosmological parameters. The shaded regions indicate their respective 1σ confidence region, i.e. 68% probability volume. This plot shows that the prior information entropy is inferior to the posterior information entropy, due to the narrower distribution of the former. The prior power spectrum adopted, as a result, does not impact significantly on the cosmological parameter inference via the AP test.

We further demonstrate this robustness of our AP test by employing a modified prior power spectrum in the inference procedure. By adopting a different cosmology (Ω_m = 0.40, w₀ = −0.85), we modify the shape of the power spectrum, and subsequently apply ALTAIR on the same mock catalogue from Sect. 4. As shown in Fig. 10, we recover the fiducial cosmological parameters employed in the mock generation, although with slightly larger uncertainties than for the original run by roughly 15%. This test case therefore explicitly highlights the robustness of our implementation of the AP test to a misspecified model since it does not optimise the information from the scale dependence of the correlations of the density field, but rather from the isotropy of the field.

Fig. 10. Same as Fig. 6, but employing a different prior power spectrum (Ωm = 0.40, w0 = −0.85), for ∼3000 MCMC realisations. By recovering the fiducial cosmological parameters employed in the mock generation, this test case explicitly highlights the robustness of our approach to the shape of the prior power spectrum adopted. The corresponding uncertainties are slightly larger than for the original run by around 15%.

From the correlation matrix of Ω_m, w₀ and the galaxy bias parameters for each of the four subcatalogues, illustrated in Fig. 11, we deduce the extremely weak correlation between the cosmological constraints and the bias. This is a key positive aspect of our method, as galaxy bias remains nevertheless a highly active and challenging field of research (see, for e.g., Desjacques et al. 2018), due to its complex non-linear behaviour on intermediate and small scales, which may potentially limit the effectiveness of traditional methods of cosmological parameter inference (Pollina et al. 2018).

Fig. 11. Correlation matrix of the galaxy bias and cosmological parameters, normalised using their respective variance. This illustrates the weak correlation between the inferred cosmological constraints and the galaxy bias. The lack of dependence on the currently unknown phenomenon of galaxy biasing is therefore a key highlight of our implementation of the AP test for cosmological parameter inference.

Moreover, this insensitivity to the galaxy bias implies that our method does not rely on absolute density fluctuation amplitudes, but on the actual location of matter. This entails that our AP test exploits the geometrical structure of the density field and not its absolute amplitude since the power spectrum does not influence the inferred cosmological constraints. To the best of our knowledge, this is a novel aspect of cosmological inference, with most of our current understanding of cosmology based on measurements of density contrast amplitudes. We present, therefore, one of the first methods which extracts a large fraction of information from statistics other than that of direct density contrast correlations, without relying on the power spectrum or bispectrum. Our method consequently yields complementary information to state-of-the-art methods.

6. Summary and conclusions

We presented the implementation of a robust AP test that performs a detailed fit of the cosmological expansion via a non-linear and hierarchical Bayesian LSS inference framework. This forward modelling approach employs LPT as a physical description for the non-linear dynamics and sequentially encodes the cosmic expansion effect for joint inference of the cosmological parameters and underlying 3D density fields, while also fitting the mean density of tracers and bias parameters. In essence, this inference machinery explores the various cosmological expansion histories and selects the cosmology-dependent evolution pathways which yield isotropic correlations of the galaxy density field, thereby constraining cosmology.

We demonstrated the application of our algorithm ALTAIR on an artificially generated galaxy catalogue, consisting of four subcatalogues, that emulates the highly structured survey geometry and selection effects of SDSS-III. We performed a series of statistical efficiency and consistency tests to validate the methodology adopted and showcased its potential to yield tight constraints on cosmological parameters from current and future galaxy redshift surveys. The main strength of our implementation of the AP test lies in its robustness to a misspecified model and its inherent approximations, thereby near-optimally exploiting the model predictions, without relying on its accuracy in modelling the scale dependence of the correlations of the field.

Moreover, another key aspect of our approach, resulting from the robustness to a misspecified model, is that the cosmological constraints show extremely weak dependence on galaxy bias. This yields two crucial advantages. First, this is especially interesting as the lack of a sufficient physical description of this bias remains a potential limiting factor for standard approaches of cosmological parameter inference from such redshift surveys. Furthermore, this lack of sensitivity to the bias also implies that our method does not depend on the absolute density fluctuation amplitudes. This is therefore among the first methods to extract a large amount of information from statistics other than that of direct density contrast correlations, without relying on the power spectrum or bispectrum, thereby providing complementary information to state-of-the-art techniques.

There is scope for further development of the ALTAIR framework, such as incorporating power spectrum inference, which is a highly non-trivial undertaking. We also intend to augment the current formalism to include the treatment of redshift space distortions, which is key for unbiased constraints on the cosmological parameters, and apply ALTAIR on state-of-the-art galaxy redshift catalogues for cosmological inference¹.

Acknowledgments

We express our appreciation to the anonymous reviewer for their constructive comments that helped us to improve the manuscript. We are grateful to Franz Elsner, Florent Leclercq and Natalia Porqueres for interesting discussions and/or for their comments on a draft version of the paper. We thank Tom Charnock for his help with coining the name ALTAIR in a coherent context. We also thank Jeremy Tinker for his assistance in interpreting the SDSS-III (DR12) science products. We acknowledge financial support from the ILP LABEX (under reference ANR-10-LABX-63) which is financed by French state funds managed by the ANR within the Investissements d’Avenir programme under reference ANR-11-IDEX-0004-02. This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. BDW is supported by the Simons Foundation. This work has made use of the Horizon/Beyond cluster at the Institut d’Astrophysique de Paris. This work is supported, through the maintenance of Horizon/Beyond cluster, under the name ORIGIN by the Domaine d’Intérêt Majeur (DIM) Astrophysique et Conditions d’Apparition de la Vie (ACAV), and received financial support from Région Île-de-France. We thank Stéphane Rouberol for running smoothly this cluster for us. This work is done within the Aquila Consortium.

References

Appendix A: The LPT-Poissonian posterior

In this section, we describe the large-scale structure (LSS) posterior distribution implemented in this work. We demonstrate how the complex problem of exploring the high dimensional joint posterior distribution can be reduced to a set of distinct steps via a multiple block sampling scheme.

A.1. The density posterior distribution

The primary objective is to fully characterise the 3D cosmic LSS in observations via a numerical representation of the corresponding LSS posterior using sophisticated Markov chain Monte Carlo (MCMC) techniques, in particular to provide data constrained realisations of a set of plausible 3D density contrast amplitudes underlying a given set of galaxy observations. The posterior distribution for the evolved density field fluctuation ${\delta }_p^{\mathrm{f}}$ can be formulated, in a general context, via Bayes’ identity as:

$$\begin{array}{c}\hfill \mathcal{P}\left(\{{\delta }_p^{\mathrm{f}}\}|\{N_p^{\mathit{g}}\}\right)=\frac{\mathcal{L}\left(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\}\right)\mathrm{\Pi }(\{{\delta }_p^{\mathrm{f}}\})}{\mathrm{\Pi }(\{N_p^{\mathit{g}}\})},\end{array}$$(A.1)

where $N_p^{\mathit{g}}$ is the observed number of galaxies in voxel p, at position x_p in the sky, in redshift space, for the gth galaxy sample, with the prior $\mathrm{\Pi }(\{{\delta }_p^{\mathrm{f}}\})$ incorporating our a priori knowledge of the present-day matter fluctuations in the Universe, the likelihood $\mathcal{L}(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\})$ and the normalizing factor given by the evidence $\mathrm{\Pi }(\{N_p^{\mathit{g}}\})$.

A major stumbling block consequently arises, as discussed extensively in Jasche & Wandelt (2013a), since the inference framework requires a suitable prior $\mathrm{\Pi }(\{{\delta }_p^{\mathrm{f}}\})$ which adequately describes the physical behaviour of the gravitationally evolved density field. Nevertheless, as elaborated in the following section, most attempts made in this direction so far have been based on heuristic approximations and the absence of a closed form description of the present day matter fluctuations encoded in a multi-variate probability density distribution still persists.

However, Jasche & Wandelt (2013a) proposed an elegant approach to circumvent this key impediment based on the following assertions: There is substantial evidence that primordial seed fluctuations at redshifts z ∼ 1000 can be modelled as a Gaussian random field to great accuracy (e.g. Linde 2008; Komatsu et al. 2011; Planck Collaboration XVII 2016), consistent with inflationary theories and CMB observations. Moreover, the evolution of the initial conditions relies solely on deterministic gravitational structure formation processes. Therefore, a conceptually reasonable alternative to modelling the complex statistical behaviour of the non-linear matter distribution is to formulate the inference problem at the level of the initial conditions adequately described by Gaussian statistics. This constitutes the conceptual foundation of the BORG framework (Jasche & Wandelt 2013a).

Given a forward model ℳ_p that connects the initial conditions ${\delta }_p^{\mathrm{ic},(r)}$, in comoving (r) space, to the redshift (z) space representation of the final density field ${\delta }_p^{\mathrm{f},(z)}$, we can therefore express the conditional posterior for the evolved density field as

$$\begin{array}{c}\hfill \mathcal{P}\left(\{{\delta }_p^{\mathrm{f}}\}|\{{\delta }_p^{\mathrm{ic}}\}\right)={\displaystyle \prod _p}{\delta }^{\mathrm{D}}[{\delta }_p^{\mathrm{f}}-{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})],\end{array}$$(A.2)

where δ^D(x) denoting the Dirac delta distribution encapsulates the assumption that the structure formation process is deterministic. Within this generic framework, the forward model may be generalised to a chain of arbitrary components linking the primordial density fluctuations to the present-day density contrast. Nevertheless, at its crux lies a cosmic structure formation model ${\mathcal{G}}_p(a,\{{\delta }_p^{\mathrm{ic}}\})$ for the non-linear evolution of initial conditions into a final density field at a given scale factor a, that is, ${\delta }_p^{\mathrm{f},(r)}={\mathcal{G}}_p(a,\{{\delta }_p^{\mathrm{ic},(r)}\})$. The forward model implemented in this work to encode the AP test is described in Sect. 3 above.

We can then obtain a prior distribution for ${\delta }_p^{\mathrm{f}}$ via a two-step sampling procedure: First, a realisation of ${\delta }_p^{\mathrm{ic}}$ is generated from the prior distribution $\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\})$ and subsequently evolved with a given forward model ${\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})$. This essentially implies generating samples from the joint prior distribution of ${\delta }_p^{\mathrm{ic}}$ and ${\delta }_p^{\mathrm{f}}$:

$$\begin{array}{cc}\hfill \mathcal{P}(\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\})& =\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\})\mathcal{P}\left(\{{\delta }_p^{\mathrm{f}}\}|\{{\delta }_p^{\mathrm{ic}}\}\right)\hfill \\ \hfill & =\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}){\displaystyle \prod _p}{\delta }^{\mathrm{D}}[{\delta }_p^{\mathrm{f}}-{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})],\hfill \end{array}$$(A.3)

after plugging in the conditional posterior distribution from Eq. (A.2).

Assuming a normal distribution with zero mean and covariance S corresponding to an underlying cosmological power spectrum for the initial conditions ${\delta }_p^{\mathrm{ic}}$, the joint prior distribution can be expressed as

$$\begin{array}{cc}\hfill \mathcal{P}(\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})& =\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S}){\displaystyle \prod _p}{\delta }^{\mathrm{D}}[{\delta }_p^{\mathrm{f}}-{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})]\hfill \\ \hfill & =\frac{{\mathrm{e}}^{-\frac{1}{2}{\sum }_{p,q}{\delta }_p^{\mathrm{ic}}{\mathbf{S}}^{-1}{\delta }_q^{\mathrm{ic}}}}{\det (2\pi \mathbf{S})}{\displaystyle \prod _p}{\delta }^{\mathrm{D}}[{\delta }_p^{\mathrm{f}}-{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})].\hfill \end{array}$$(A.4)

Reformulating the statistical inference problem in terms of the initial conditions ${\delta }_p^{\mathrm{ic}}$ results in the following joint posterior distribution of ${\delta }_p^{\mathrm{ic}}$ and ${\delta }_p^{\mathrm{f}}$, conditional on the observed galaxy number counts $N_p^{\mathit{g}}$,

$$\begin{array}{c}\hfill \mathcal{P}(\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\}|\{N_p^{\mathit{g}}\},\mathbf{S})=\frac{\mathcal{L}(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\})\mathcal{P}(\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})}{\mathrm{\Pi }(\{N_p^{\mathit{g}}\}|\mathbf{S})},\end{array}$$(A.5)

after making the dependence on the underlying power spectrum explicit. Assuming that the galaxy observations are conditionally independent of the initial conditions once the final density field is known, that is, $\mathcal{L}(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\})=\mathcal{L}(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\})$ and using the joint prior distribution from Eq. (A.4) leads to the LSS posterior distribution:

$$\begin{array}{cc}\hfill \mathcal{P}(\{{\delta }_p^{\mathrm{f}}\},\{{\delta }_p^{\mathrm{ic}}\}|\{N_p^{\mathit{g}}\},\mathbf{S})& =\mathcal{L}\left(\{N_p^{\mathit{g}}\}|\{{\delta }_p^{\mathrm{f}}\}\right)\frac{\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})}{\mathrm{\Pi }(\{N_p^{\mathit{g}}\}|\mathbf{S})}\hfill \\ \hfill & \times {\displaystyle \prod _p}{\delta }^{\mathrm{D}}[{\delta }_p^{\mathrm{f}}-{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})].\hfill \end{array}$$(A.6)

By marginalizing over the final density field ${\delta }_p^{\mathrm{f}}$, we finally obtain our posterior distribution as follows:

$$\begin{array}{c}\hfill \mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\}|\{N_p^{\mathit{g}}\},\mathbf{S})=\mathcal{L}\left[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})\right]\frac{\mathrm{\Pi }(\{{\delta }_p^{\mathrm{ic}}\}|\mathbf{S})}{\mathrm{\Pi }(\{N_p^{\mathit{g}}\}|\mathbf{S})},\end{array}$$(A.7)

thereby connecting present galaxy observations $N_p^{\mathit{g}}$ to their corresponding primordial density fields ${\delta }_p^{\mathrm{ic}}$ via a given forward model, ${\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\})$. We must therefore sample from the highly non-linear and non-Gaussian posterior above to obtain realisations of the 3D initial density fields conditioned on the galaxy observations via the sophisticated HMC method described below in Appendix C.

A.2. Choice of density prior

Standard Wiener filtering approaches employ isotropic Gaussian priors on the present-day density contrast which is justified for inference on the largest scales (e.g. Zaroubi et al. 1999; Zaroubi 2002; Erdoǧdu et al. 2004, 2006; Kitaura & Enßlin 2008; Kitaura et al. 2009; Jasche et al. 2010a; Jasche & Lavaux 2015), where the density field can be reasonably approximated as a Gaussian random field. Although the exact probability distribution for the density field in the mildly and strongly non-linear regimes is not known, the deviation from Gaussianity is well-established and we therefore require a non-Gaussian prior to effectively capture the details of the highly complex filamentary nature of the present day cosmic web.

Lognormal density priors were subsequently proposed (Coles & Jones 1991) and this proved to be an adequatephenomenological choice (e.g. Hubble 1934; Peebles 1980; Gaztanaga & Yokoyama 1993; Kayo et al. 2001), albeit withsome limitations, for modelling the evolved matter field in the mildly non-linear regime (Kitaura et al. 2009; Jasche et al. 2010; Jasche et al. 2010b). The use of Edgeworth expansions to construct prior distributions involving third order moments has also been proposed in the literature (Kitaura 2012).

In this work, we rely on Lagrangian perturbation theory (LPT) to model cosmic structure formation, which accounts for non-local effects of gravitational mass transport from initial to final positions. It has been extensively shown that LPT provides a sufficient description of the cosmic LSS on large scales, where it is capable of reproducing the exact one-, two- and three-point statistics of the evolved field, while still being a reasonable approximation to the higher order statistics (e.g. Moutarde et al. 1991; Buchert et al. 1994; Bouchet et al. 1995; Scoccimarro 2000; Scoccimarro & Sheth 2002).

The essence of the above approach is that it provides a pathway to recover the high-order statistics of the matterdistribution using only the 2-point statistics of the initial conditions, obviating the need for additional parameters to describe higher order statistics for the matter inference problem.Moreover, since our model encodes the dynamics, the reconstruction of the large-scale velocity field is automatically generated without incurring the expense of an increased parameter space.

A.3. Modular statistical programming

In this work, we exploit the modular statistical programming framework inherent in the BORG algorithm to encode the AP test as an additional component to the original block sampling machinery. This approach allows us to model any Bayesian hierarchical problem to take into account any observational systematics via data models of higher complexity. Here, we account for the unknown parameters $\{b_i^{\mathit{g}}\}$ of the galaxy bias model described in Sect. 3.1 and unknown normalisations {N̄^g} for distinct galaxy samples, as illustrated in Fig. A.1. These last normalisations, in practice, scale as the noise amplitudes as we are using a per-voxel Poissonian likelihood.

Conceptually, within such a framework, the overall aim is to characterise fully the augmented joint posterior $\mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\},\{{\theta }_i\}|\{N_p^{\mathit{g}}\},\mathbf{S})$ of different tracer populations labelled by the index g. Since it is not computationally feasible and advisable to sample directly from the high dimensional joint posterior distribution, we make use of an important theorem on Metropolis-Hastings block sampling which allows us to break the high dimensional sampling problem into a series of lower dimensional ones (Hastings 1970). We therefore dissect the exploration of the full joint parameter space into a sequence of conditional sampling procedures. The block sampling approach consists of drawing samples from the following conditional probability distributions:

$$\begin{array}{cc}& (1){\{{\delta }_p^{\mathrm{ic}}\}}^{s+1}\mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\}|{\{{\theta }_i\}}^s,{\{N_p^{\mathit{g}}\}}^s,{\{{\overline{N}}^{\mathit{g}}\}}^s,{\{b_i^{\mathit{g}}\}}^s,\mathbf{S})\hfill \\ \hfill & (2){\{{\overline{N}}^{\mathit{g}}\}}^{s+1}\mathcal{P}(\{{\overline{N}}^{\mathit{g}}\}|{\{{\theta }_i\}}^s,{\{N_p^{\mathit{g}}\}}^s,{\{{\delta }_p^{\mathrm{ic}}\}}^s,{\{b_i^{\mathit{g}}\}}^s,\mathbf{S})\hfill \\ \hfill & (3){\{b_i^{\mathit{g}}\}}^{s+1}\mathcal{P}(\{b_i^{\mathit{g}}\}|{\{{\theta }_i\}}^s,{\{{\overline{N}}^{\mathit{g}}\}}^s,{\{N_p^{\mathit{g}}\}}^s,{\{{\delta }_p^{\mathrm{ic}}\}}^s,\mathbf{S})\hfill \\ \hfill & (4){\{{\theta }_i\}}^{s+1}\mathcal{P}(\{{\theta }_i\}|{\{{\overline{N}}^{\mathit{g}}\}}^s,{\{N_p^{\mathit{g}}\}}^s,{\{{\delta }_p^{\mathrm{ic}}\}}^s,{\{b_i^{\mathit{g}}\}}^s,\mathbf{S}).\hfill \end{array}$$(A.8)

Fig. A.1. Schematic representation of the multi-step iterative block sampling procedure implemented in ALTAIR. Initially, a realisation of the 3D density contrast is obtained conditional on galaxy observations, followed by the other parameters being sampled conditional on the respective previous samples. Subsequent iterations of this procedure yield samples from the full joint posterior distribution within a modular statistical programming framework.

In the above expressions, s denotes the sampling step and the symbol ↶ indicates sampling from the distribution on the right hand side. A series of iterations of the individual sampling steps above will converge to the target distribution which corresponds the full joint posterior distribution (Hastings 1970). Hence, by simultaneously exploring the 3D initial density field $\{{\delta }_p^{\mathrm{ic}}\}$, the galaxy bias parameters $\{b_i^{\mathit{g}}\}$, the galaxy density contrast normalisations {N̄^g} and the cosmological parameters {θ_i} via an implementation of high dimensional MCMC methods in a multiple block sampling procedure, we can obtain samples from the desired joint posterior distribution.

A.4. The bias posterior distribution

The formalism for the data model is presented in a generic context in Sect. 3.1, such that it can be implemented for two or more different galaxy surveys, distinct subsamples of a given catalogue or even different cosmological probes of the LSS. The advantage of splitting a galaxy sample into various subsamples is that we can treat the respective systematic and statistical uncertainties of each subsample separately, thereby allowing us to account for the distinct clustering behaviour of galaxy populations in the LSS sample via their respective bias parameters.

These additional parameters can be trivially incorporated in the flexible block sampling approach adopted here, as described in Sect. A.3. The mean numbers of galaxies, N̄^g, for the various subsamples are essential for defining the density contrast of the galaxy distribution, with an erroneous value of N̄^g resulting in a non-zero value of the mean, yielding spurious large-scale power in the inferred density fields (Jasche & Wandelt 2013b). Due to a lack of a priori knowledge, we must perform a joint inference of the set of four N̄^g and $b_i^{\mathit{g}}$ bias parameters, together with initial and final density fields, to take into account possible cross-correlations and interdependencies. Unlike traditional approaches (e.g. Tegmark et al. 2004), here we infer the bias factors directly from the relation between the data and the density field.

The conditional posterior distribution for the bias parameters, given a realisation of the density field and galaxy number counts for the respective subcatalogues, can be expressed as:

$$\begin{array}{cc}\hfill \mathcal{P}(\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}|\{{\theta }_i\},\{N_p^{\mathit{g}}\},\{{\delta }_p^{\mathrm{ic}}\})& \propto \mathcal{L}[\{N_p^{\mathit{g}}\}|{\mathcal{M}}_p(\{{\delta }_p^{\mathrm{ic}}\}),\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}]\hfill \\ \hfill & \times \mathrm{\Pi }(\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\}),\hfill \end{array}$$(A.9)

following an analogous reasoning to that described in Sect. 3.3. Adopting a standard maximally agnostic philosophy, we set uniform prior distributions for the bias parameters:

$$\begin{array}{c}\hfill \mathrm{\Pi }(\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})=\mathrm{\Theta }(\{{\overline{N}}^{\mathit{g}}\})\mathrm{\Theta }(\{b_i^{\mathit{g}}\}),\end{array}$$(A.10)

where the Heaviside function Θ(x) ensures that the parameters are positive, as required by the bias model. The non-linear shape of the galaxy biasing function, as given in Eq. (4), poses a particular challenge as no straightforward sampling procedure can be derived. We therefore sample from the above conditional bias posterior distribution via sequential slice sampling steps (e.g. Neal 2000, 2003) to ensure unit acceptance rates.

Appendix B: Jacobian of comoving-redshift transformation

The derivation of the Jacobian, ${\mathcal{J}}_r^z$, of the coordinate transformation between comoving (r) and redshift (z) transformation is described below.

Since we desire the redshift space representation of the density field,

$$\begin{array}{c}\hfill {\rho }_z[\mathit{z}(\mathit{r})]={\rho }_r(\mathit{r}){\left|\frac{\partial z}{\partial r}\right|}^{-1}={\rho }_r(\mathit{r})|{\mathcal{J}}_z^r{|}^{-1}={\rho }_r(\mathit{r})|{\mathcal{J}}_r^z|,\end{array}$$(B.1)

or in terms of the density contrast,

$$\begin{array}{c}\hfill {\delta }_z[\mathit{z}(\mathit{r})]=[1+{\delta }_r(\mathit{r})]|{\mathcal{J}}_r^z|-1,\end{array}$$(B.2)

where ${\mathcal{J}}_z^r$ is the corresponding Jacobian for the converse (redshift to comoving) transformation.

The set of three indices {i, j, k} below are spatial coordinates {1, 2, 3}. With the functional dependence of z = f(r), we have

$$\begin{array}{c}\hfill z_i=\frac{r_i}{r}f(\mathit{r})=\overline{z}\left(r\right)\frac{r_i}{r},\end{array}$$(B.3)

where we defined $r=\sqrt{\sum _j{r}_j^2}$, and z̄ is the cosmological redshift. Taking the derivative yields

$$\begin{array}{cc}\hfill \mathrm{d}{z}_i& =\mathrm{d}\overline{z}\left(r\right)\frac{r_i}{r}+\overline{z}\left(r\right)\mathrm{d}\left(\frac{r_i}{r}\right)\hfill \\ \hfill & ={\displaystyle \sum _j\frac{\partial \overline{z}}{\partial r}\frac{\partial r}{\partial r_j}\mathrm{d}{r}_j\frac{r_i}{r}+\overline{z}\frac{\mathrm{d}{r}_i}{r}+\overline{z}{r}_i\mathrm{d}\left(\frac{1}{r}\right)},\hfill \end{array}$$(B.4)

after applying the chain rule. Differentiating with respect to the comoving coordinate, and using d(1/r) = −∑_j(1/r³)r_jdr_j,

$$\begin{array}{cc}\hfill \frac{\partial z_i}{\partial r_k}& =\frac{\partial \overline{z}}{\partial r}\frac{\partial r}{\partial r_k}\frac{r_i}{r}+\frac{\overline{z}}{r}{\delta }_{\mathit{ik}}^{\mathrm{k}}-\frac{z{r}_i{r}_k}{r^3}\hfill \\ \hfill & =(\frac{\partial \overline{z}}{\partial r}-\frac{\overline{z}}{r})\frac{r_i{r}_k}{r^2}+\frac{\overline{z}}{r}{\delta }_{\mathit{ik}}^{\mathrm{k}},\hfill \end{array}$$(B.5)

where we made use of rdr = ∑_kr_kdr_k which follows from the definition $r=\sqrt{\sum _k{r}_k^2}$. To compute ∂z̄/∂r, we start from the definition of comoving distance,

$$\begin{array}{c}\hfill r=\frac{c}{H_0}{\displaystyle {\int }_0^{\overline{z}}\frac{\mathrm{d}z}{E(z)}},\end{array}$$(B.6)

with the conventional definition for ℰ(z),

$$\begin{array}{c}\hfill E(z)=\sqrt{{\mathrm{\Omega }}_{\mathrm{de}}{(1+z)}^{3(1+w)}+{\mathrm{\Omega }}_{\mathrm{k}}{(1+z)}^2+{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{r}}{(1+z)}^4},\end{array}$$(B.7)

where c is the speed of light and H₀ is the Hubble constant today, such that

$$\begin{array}{c}\hfill \frac{\partial \overline{z}}{\partial r}=\left(\frac{H_0}{c}\right)\mathcal{E}(\overline{z}).\end{array}$$(B.8)

Finally, the Jacobian matrix can be expressed as follows:

$$\begin{array}{c}\hfill {\mathcal{J}}_{z\{ik\}}^r=\frac{\partial z_i}{\partial r_k}=(\frac{H_0}{c}E(\overline{z})-\frac{\overline{z}}{r})\frac{r_i{r}_k}{r^2}+\frac{\overline{z}}{r}{\delta }_{\mathit{ik}}^{\mathrm{k}},\end{array}$$(B.9)

and we obtain the desired $|{\mathcal{J}}_{r\{ik\}}^z|$ by taking the reciprocal of $|{\mathcal{J}}_{z\{ik\}}^r|$.

Appendix C: Hamiltonian sampling

The exploration of the LPT-Poissonian posterior, described by Eq. (6), requires highly non-linear reconstruction methods and is therefore numerically intensive. From the Bayesian viewpoint considered in this work, a single estimate of the density field is not of particular interest, with the desired scientific output being a sampled representation of the multi-dimensional LSS posterior. Extracting any relevant statistical summary, such as mean, mode or variance, given this posterior representation, is then straightforward. Furthermore, the propagation of uncertainties to the final inferred quantities is smooth and coherent.

In the absence of a known procedure to directly sample from a LPT-Poissonian distribution, a Metropolis-Hastings sampling scheme with an accept-reject step must be implemented. But such a mechanism suffers from the well-known drawback of possible high rejection rate, where a low acceptance rate of proposed samples will render the method numerically inefficient. This is especially significant since the inference of 3D density fields typically involves extremely high number of inference parameters. For instance, in this work, there are over 2 × 10⁶ free parameters, which correspond to primordial density fluctuation amplitudes ${\delta }_p^{\mathrm{ic}}$ at respective grid nodes. Exploring this high dimensional parameter space via a random walk, as in conventional MCMC methods, consequently results in a high rejection rate. This is an understatement, as in practice, the acceptance rate is virtually zero, with one chance over 10⁶ of going in the right direction, without taking into account the step size.

To limit the number of samples and alleviate the numerical scaling of the method, we implement a Hamiltonian Monte Carlo (HMC) scheme. This method guarantees an acceptance rate of unity in the absence of numerical errors. The numerical efficiency of the Hamiltonian sampling lies in deterministically proposing new samples to the Metropolis-Hastings algorithm, based on techniques to follow dynamical particle motion in potentials. The HMC algorithm therefore maintains a reasonable sampling efficiency in high dimensional spaces by suppressing the dominant random walk behaviour.

The HMC method has been successfully implemented for the inference of density fields in the non-linear regime, and has been found to be very efficient (e.g. Jasche & Kitaura 2010; Jasche et al. 2010b; Kitaura et al. 2012; Jasche & Wandelt 2013a,b). Moreover, this approach is less prone to spurious effects induced by numerical inaccuracies due to the final accept/reject step of the Metropolis-Hastings sampler ensuring correctness of the target density. Numerical inaccuracies would only be detrimental to the efficiency, without comprimising the correctness.

We briefly review the underlying framework and rationale of the HMC method. An excellent in-depth review of this approach is provided in Duane et al. (1987) and Neal (1993), with its application to LSS inference illustrated in Jasche & Kitaura (2010), Jasche et al. (2010b), Jasche & Wandelt (2013a,b) and Jasche & Lavaux (2018).

The essence of the HMC method is as follows: If we want to generate random variates according to a desired probability distribution 𝒫({x_i}), with {x_i} being a set of N elements x_i, then we may interpret the negative logarithm of the distribution as a potential ψ(x),

$$\begin{array}{c}\hfill \psi (x)=-ln(\mathcal{P}(x)).\end{array}$$(C.1)

We can subsequently formulate a corresponding Hamiltonian that describes the dynamics in the multi-dimensional phase space by adding a kinetic term to the above potential. To this end, we introduce a “momentum” variable p_i and a “mass matrix” M as nuisance parameters, as follows:

$$\begin{array}{c}\hfill H=\frac{1}{2}{\displaystyle \sum _{i,j}{p}_i{M}_{\mathit{ij}}^{-1}{p}_j+\psi (x)}.\end{array}$$(C.2)

This Hamiltonian, as in classical mechanics, describes the dynamics in a high dimensional parameter space. We then obtain the new target probability distribution by exponentiating the above Hamiltonian:

$$\begin{array}{c}\hfill e^{-H}=\mathcal{P}(\{x_i\})exp(\frac{1}{2}{\displaystyle \sum _{i,j}}{p}_i{M}_{\mathit{ij}}^{-1}{p}_j).\end{array}$$(C.3)

The form of the Hamiltonian ensures that the new joint distribution can be separated into a Gaussian distribution for the momenta {p_i} and the target distribution 𝒫({x_i}). This implies that the two sets of variables {x_i} and {p_i} are independent, and hence, we can obtain samples from the target distribution 𝒫({x_i}) simply by marginalizing over the auxiliary momenta.

The next step is to generate a random variate from the joint distribution, which is proportional to exp(−H). We must first draw a set of momenta from the distribution defined by the kinetic energy term, which is an N-dimensional Gaussian with covariance M. To obtain a new sample from the joint distribution, we allow the system to evolve deterministically, from a given initial point in the high dimensional parameter space ({x_i},{p_i}), according to Hamilton’s equations:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{x}_i}{\mathrm{d}t}& =\frac{\partial H}{\partial p_i}\hfill \end{array}$$(C.4)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_i}{\mathrm{d}t}& =\frac{\partial H}{\partial x_i}=-\frac{\partial \psi (x)}{\partial x_i}·\hfill \end{array}$$(C.5)

Integrating the above equations of motion for the time variable t up to a fixed pseudo-time τ results in the new position ($\left\{x_i^{\prime }\right\},\left\{p_i^{\prime }\right\}$) in phase space. This new point in phase space is accepted according to the standard Metropolis-Hastings acceptance rule:

$$\begin{array}{c}\hfill {\mathcal{P}}_A=\min \{1,exp[-(H(\{x_i^{\prime }\},\{p_i^{\prime }\})-H(\{x_i\},\{p_i\}))]\}.\end{array}$$(C.6)

Since the equations of motion provide a solution to a Hamiltonian system, the energy or the Hamiltonian described by Eq. (C.2) must be conserved. A direct consequence of this conservation of the Hamiltonian is that the HMC approach yields an acceptance rate of unity. In practice, however, the acceptance rate may be lower due to numerical inaccuracies in the integration scheme. After accepting a new sample, the mechanism proceeds by discarding the momentum counterpart and restarts the sampling procedure by randomly drawing a new set of momenta. The individual momenta {p_i} do not have to be stored in memory which is numerically convenient, and marginalisation simply implies discarding them.

The HMC sampling scheme can therefore be summarised in two steps. The first step is essentially a Gibbs sampling step which results in a new set of Gaussian distributed momenta. The second step involves computing the deterministic dynamical trajectory on the posterior surface. We obtain samples of the desired target distribution by marginalizing over the nuisance parameters, that is, by discarding the auxiliary momenta {p_i}.

Appendix D: Equations of motion for LSS inference

As described in Appendix C, the Hamiltonian Monte Carlo (HMC) approach allows the exploration of the non-linear LSS posterior via the study of the Hamiltonian dynamics in the high dimensional parameter space. To ensure a high acceptance rate, the HMC exploits the gradient of the logarithmic posterior distribution to optimally explore this parameter space, such that the algorithm also requires the derivatives of the posterior distribution, as outlined below.

From the LSS posterior illustrated in Eq. (6), we can derive the corresponding forces required to evaluate the Hamiltonian trajectories. The Hamiltonian potential ψ(δ^ic) can be written as

$$\begin{array}{cc}\hfill \psi ({\delta }^{\mathrm{ic}})& =-ln\left[\mathcal{P}(\{{\delta }_p^{\mathrm{ic}}\},\{{\overline{N}}^{\mathit{g}}\},\{b^{\mathit{g}}\},\{{\theta }_i\}|\{N_p^{\mathit{g}}\},\mathbf{S})\right]\hfill \\ \hfill & ={\psi }_{\mathrm{likelihood}}({\delta }^{\mathrm{ic}})+{\psi }_{\mathrm{prior}}({\delta }^{\mathrm{ic}}),\hfill \end{array}$$(D.1)

where the potential ψ_prior(δ^ic) corresponds to

$$\begin{array}{c}\hfill {\psi }_{\mathrm{prior}}({\delta }^{\mathrm{ic}})=\frac{1}{2}{\delta }^{\mathrm{ic},†}{\mathbf{S}}^{-1}{\delta }^{\mathrm{ic}}=\frac{1}{2}{\displaystyle \sum _{\mathit{ij}}{\delta }_i^{\mathrm{ic},\ast }{S}_{\mathit{ij}}^{-1}{\delta }_j^{\mathrm{ic}}},\end{array}$$(D.2)

while ψ_likelihood(δ^ic) is given by

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{likelihood}}({\delta }^{\mathrm{ic}})& ={\displaystyle \sum _p\{{\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})}\hfill \\ \hfill & -N_p^{\mathit{g}}ln\left[{\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})\right]\}.\hfill \end{array}$$(D.3)

Using the definition of the Hamiltonian potential ψ(δ^ic) from Eq. (D.1), we obtain the corresponding equations describing the Hamiltonian forces by differentiating with respect to δ^ic:

$$\begin{array}{c}\hfill \frac{\partial \psi ({\delta }^{\mathrm{ic}})}{\partial {\delta }_i^{\mathrm{ic}}}=\frac{\partial {\psi }_{\mathrm{prior}}({\delta }^{\mathrm{ic}})}{\partial {\delta }_i^{\mathrm{ic}}}+\frac{\partial {\psi }_{\mathrm{likelihood}}({\delta }^{\mathrm{ic}})}{\partial {\delta }_i^{\mathrm{ic}}}·\end{array}$$(D.4)

The prior term is trivially obtained as

$$\begin{array}{c}\hfill \frac{\partial {\psi }_{\mathrm{prior}}({\delta }^{\mathrm{ic}})}{\partial {\delta }_i^{\mathrm{ic}}}={\displaystyle \sum _j{\delta }_j^{\mathrm{ic},\ast }{S}_{\mathit{ij}}^{-1}}.\end{array}$$(D.5)

The corresponding likelihood term, however, cannot be derived in straightforward fashion. We compute this adjoint gradient as a sequential application of linear operations, as follows:

$$\begin{array}{c}\hfill \frac{\partial {\psi }_{\mathrm{likelihood}}({\delta }^{\mathrm{ic}})}{\partial {\delta }_s^{\mathrm{ic}}}=\sum _q\frac{\partial {\rho }_q^{(r)}}{\partial {\delta }_s^{\mathrm{ic}}}\sum _{\stackrel{\sim }{p}}\frac{\partial {\rho }_{\stackrel{\sim }{p}}^{(z)}}{\partial {\rho }_q^{(r)}}\sum _{p,\mathit{g}}\frac{\partial {\lambda }_p^{\mathit{g}}}{\partial {\rho }_{\stackrel{\sim }{p}}^{(z)}}\frac{\partial {\psi }_{\mathrm{likelihood}}}{\partial {\lambda }_p^{\mathit{g}}},\end{array}$$(D.6)

where we explicitly expressed the coordinate transformation from comoving (r) to redshift (z) space, and ${\rho }_p=1+{\delta }_p^{\mathrm{f}}$ is the density field. Equation (D.6) constitutes a sequence of matrix vector applications, as follows:

$$\begin{array}{cc}& {\stackrel{\sim }{\mathcal{G}}}_{\mathit{qs}}\equiv \frac{\partial {\rho }_q^{(r)}}{\partial {\delta }_s^{\mathrm{ic}}}\hfill \end{array}$$(D.7)

$$\begin{array}{cc}& {\mathcal{Q}}_{\stackrel{\sim }{p}q}\equiv \frac{\partial {\rho }_{\stackrel{\sim }{p}}^{(z)}}{\partial {\rho }_q^{(r)}}={\displaystyle \sum _i}{\mathcal{E}}_{i\stackrel{\sim }{p}}^{-1}{\displaystyle \sum _{\stackrel{\sim }{p}}}{\mathcal{J}}_{\stackrel{\sim }{p}}{x}_{\stackrel{\sim }{p}}^{\alpha (i)}{y}_{\stackrel{\sim }{p}}^{\beta (i)}{z}_{\stackrel{\sim }{p}}^{\gamma (i)}\hfill \end{array}$$(D.8)

$$\begin{array}{cc}& {\mathcal{K}}_{p\stackrel{\sim }{p}}\equiv {\displaystyle \sum _{\mathit{g}}}\frac{\partial {\lambda }_p^{\mathit{g}}}{\partial {\rho }_{\stackrel{\sim }{p}}^{(z)}}={\displaystyle \sum _{\mathit{g}}}{\delta }_{p,\stackrel{\sim }{p}}^{\mathrm{k}}\frac{{\lambda }_p^{\mathit{g}}}{{\rho }_{\stackrel{\sim }{p}}^{(r)}}[{\rho }^{\mathit{g}}{ϵ}^{\mathit{g}}{\rho }_{\stackrel{\sim }{p}}^{-{ϵ}_{\mathit{g}},(r)}+\beta ]\hfill \end{array}$$(D.9)

$$\begin{array}{cc}& {\upsilon }_p\equiv \frac{\partial {\psi }_{\mathrm{likelihood}}}{\partial {\lambda }_p^{\mathit{g}}}=1-\frac{N_p^{\mathit{g}}}{{\lambda }_p^{\mathit{g}}(\{{\delta }_p^{\mathrm{f}}\},\{{\theta }_i\},\{{\overline{N}}^{\mathit{g}}\},\{b_i^{\mathit{g}}\})},\hfill \end{array}$$(D.10)

where 𝒥 is the Jacobian of the comoving-redshift transformation, $\stackrel{\sim }{\mathcal{G}}$ indicates the derivative of LPT (cf. Appendix D in Jasche & Wandelt 2013a) and ${\delta }_{p,\stackrel{\sim }{p}}^{\mathrm{k}}$ is the Kronecker delta. The derivation of the adjoint gradient for the triquintic interpolation is illustrated in Appendix G.

Using the above results, the two equations of motion (C.4) and (C.5) for the Hamiltonian system can be expressed as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{\delta }}_i^{\mathrm{ic}}}{\mathrm{d}t}=\sum _j{M}_{\mathit{ij}}^{-1}{p}_j,\end{array}$$(D.11)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\widehat{p}}_i}{\mathrm{d}t}=-{\displaystyle \sum _j}{\delta }_j^{\mathrm{ic},\ast }{S}_{\mathit{ij}}^{-1}-{\displaystyle \sum _{q,\stackrel{\sim }{p},p}}{\stackrel{\sim }{\mathcal{G}}}_{\mathit{qi}}{\mathcal{Q}}_{\stackrel{\sim }{p}q}{\mathcal{K}}_{p\stackrel{\sim }{p}}{\upsilon }_p,\end{array}$$(D.12)

where the hats denote Fourier space representation. We can now obtain new samples from the LSS posterior by following the dynamical evolution of the Hamiltonian system, governed by Eqs. (D.11) and (D.12), in phase space.

Appendix E: Numerical implementation

Our numerical implementation of the augmented version of BORG, that incorporates cosmological parameter inference, while assuming a more realistic non-linear bias model, is designated as ALTAIR. It employs the FFTW3 library for fast Fourier transforms (Frigo & Johnson 2005), whose feature of parallel transforms allows for straightforward parallelisation of our code. For random number generation, we make use of the GNU scientific library (gsl; Galassi et al. 2009), in particular, the Mersenne Twister MT19937, with 32-bit word length, from the gsl_rng_mt19937 routine. The use of the Mersenne Twister algorithm as a pseudo-random number generator for Monte Carlo simulations has been proven to be efficient and dependable (Matsumoto & Nishimura 1998).

E.1. The leapfrog scheme

In order to obtain samples from the LSS posterior, we must integrate the equations of motion (D.11) and (D.12) numerically over a pseudo-time interval τ by means of a leapfrog scheme. Ideally, we want to solve the equations exactly for an optimal acceptance rate. As such, the choice of this integrator is motivated by several reasons. Primarily, the leapfrog scheme is convenient as it is a symplectic integrator and is therefore exactly reversible, thereby maintaining detailed balance (Duane et al. 1987) of the Markov chain. The high accuracy of such an integration scheme is conducive to the conservation of the total Hamiltonian along a given trajectory, within the limit of numerical errors, resulting in high acceptance rates. It is also numerically stable and allows for simple propagation of errors. To avoid resonant trajectories, τ is drawn from a uniform distribution. The numerical implementation of the “kick-drift-kick” leapfrog scheme here follows closely the descriptions provided in Jasche & Kitaura (2010) and Jasche & Wandelt (2013a), which we refer the reader to for more elaborate and complementary explanations.

E.2. Hamiltonian mass

The mass matrix M contains a large number of parameters that must be tuned to optimise the performance of the HMC sampler. Conceptually, this Hamiltonian mass characterises the inertia of individual parameters as they move through the parameter space. As a result, the choice of M is a compromise to limit slow exploration efficiency due to extremely large masses and avoid large numerical errors of the integration scheme due to extremely light masses, and is therefore a trade-off between efficiency and acceptance rate.

A general prescription for the Hamiltonian masses, when the density fluctuations ${\delta }_p^{\mathrm{f}}$ can be characterised as a Gaussian random field, is to set them inversely proportional to the variance of that specific ${\delta }_p^{\mathrm{f}}$ (Taylor et al. 2008). However, this prescription has been found to be effective even for the case of non-Gaussian distributions, such as the LPT-Poissonian posterior encountered in this work (Jasche & Wandelt 2013a; Jasche et al. 2015; Lavaux & Jasche 2016). Since the leapfrog scheme requires the inversion of M, a diagonal representation of M is numerically convenient, given the extremely high dimensionality of the problem. We therefore choose M to be inversely proportional to the cosmological power spectrum, M = S⁻¹, which is diagonal in its Fourier basis. This choice, moreover, improves the efficiency of the HMC sampler and results in faster convergence since it accounts for the correlation structure of the underlying density field (Jasche & Wandelt 2013a).

Appendix F: Rotation of cosmological parameter space

To increase the efficiency of the cosmological parameter sampler, we rotate the (Ω_m, w₀) parameter space, using orthogonal basis transformations, before slice sampling. The corresponding covariance matrix can be orthonormally decomposed as C = Q^† ΔQ, where Q^† Q = 1.

We perform a rotation of the θ vector of cosmological parameters, about the mean $\overline{\mathit{\theta }}$, using the transformation:

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathit{\theta }}={\mathbf{\Delta }}^{-\frac{1}{2}}\mathbf{Q}(\mathit{\theta }-\overline{\mathit{\theta }}),\end{array}$$(F.1)

where the tilde denotes the transformed variable. A key point is that we must rotate back to the original frame when computing the logarithm of the posterior via

$$\begin{array}{c}\hfill \mathit{\theta }={\mathbf{Q}}^{†}{\mathbf{\Delta }}^{\frac{1}{2}}\stackrel{\sim }{\mathit{\theta }}+\overline{\mathit{\theta }}.\end{array}$$(F.2)

This is implemented at each step of the Markov chain, such that the covariance used to compute the orthogonal basis operators must be learnt progressively as follows:

$$\begin{array}{c}\hfill {\mathbf{C}}_{n+1}=\frac{n}{n+1}{\mathbf{C}}_n+\frac{{({\mathit{\theta }}_n-\overline{\mathit{\theta }})}^{†}({\mathit{\theta }}_n-\overline{\mathit{\theta }})}{n},\end{array}$$(F.3)

at the given nth step. The above transformation helps to decorrelate the variables involved when sampling, thereby reducing the burn-in period significantly. This procedure is especially effective with highly correlated parameters, as is the case for (Ω_m, w₀) (cf. Fig. 6).

We investigate the statistical efficiency of our cosmological parameter sampler by computing the correlation length of the Markov chain. Subsequent samples in the chain are generally correlated and hence do not qualify as independent samples of the target posterior distribution. The correlation length characterises the statistical efficiency of generating independent samples of a given parameter θ, as follows:

$$\begin{array}{c}\hfill C{(\theta )}_n=\frac{1}{N-n}{\displaystyle \sum _{i=0}^{N-n}\left(\frac{{\theta }^i-⟨\theta ⟩}{\sqrt{\mathrm{Var}(\theta )}}\frac{{\theta }^{i+n}-⟨\theta ⟩}{\sqrt{\mathrm{Var}(\theta )}}\right)},\end{array}$$(F.4)

where n is the number of transition steps in the chain, ⟨θ⟩ and Var(θ) correspond to the mean and variance, respectively, ⟨θ⟩ = 1/N∑_i θⁱ and Var(θ) = 1/N∑_i(θⁱ − ⟨θ⟩)², while N is the total number of realisations in the chain. This correlation length therefore indicates the number of independent realisations that can be drawn from the MCMC chain with a given number of transition steps. For an unbiased comparison, we computed the correlation length using the same chain length for both scenarios, while ignoring their respective burn-in phases. From the left panel of Fig. F.1, we deduce that the correlation length for the Ω_m chain is of the order of 25 samples. Without the rotation of (Ω_m, w₀) parameter space, the correlation length is longer by nearly a factor of five, further demonstrating the gain in efficiency obtained with the more sophisticated sampler.