Cosmological inference from Bayesian forward modelling of deep galaxy redshift surveys

We present a large-scale Bayesian inference framework to constrain cosmological parameters using galaxy redshift surveys, via an application of the Alcock-Paczy\'nski (AP) test. Our physical model of the non-linearly evolved density field, as probed by galaxy surveys, employs Lagrangian perturbation theory (LPT) to connect Gaussian initial conditions to the final density field, followed by a coordinate transformation to obtain the redshift space representation for comparison with data. We generate realizations of primordial and present-day matter fluctuations given a set of observations. This hierarchical approach encodes a novel AP test, extracting several orders of magnitude more information from the cosmological expansion compared to classical approaches, to infer cosmological parameters and jointly reconstruct the underlying 3D dark matter density field. The novelty of this AP test lies in constraining the comoving-redshift transformation to infer the appropriate cosmology which yields isotropic correlations of the galaxy density field, with the underlying assumption relying purely on the cosmological principle. Such an AP test does not rely explicitly on modelling the full statistics of the field. We verify in depth via simulations that this renders our test robust to model misspecification. This leads to another crucial advantage, namely that the cosmological parameters exhibit extremely weak dependence on the currently unresolved phenomenon of galaxy bias, thereby circumventing a potentially key limitation. This is consequently among the first methods to extract a large fraction of information from statistics other than that of direct density contrast correlations, without being sensitive to the amplitude of density fluctuations. We perform several statistical efficiency and consistency tests on a mock galaxy catalogue, using the SDSS-III survey as template.


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. 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 , for efficiently sampling the high dimensional and nonlinear 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.

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  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).
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(Friedmann , 1924Lemaître 1927Lemaître , 1931Lemaître , 1933Robertson 1935Robertson , 1936aWalker 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).
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 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. 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.

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 M p which links the initial conditions δ ic,(r) p to the redshift space representation of the evolved density field δ f,(z) p as follows: where ρ f(r) is the final density field in comoving space. The forward model consists of two components, , 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, δ f(r) p . 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 M (1) p (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, |J ‡ ∇ | (cf. Appendix B), entails cosmological dependence and is consequently included in the AP test as well as through the direct coordinate dependence E i j .
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 cosmologydependent 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 nonperturbative 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, M (2) p ≡ G p (a, {δ ic p }), 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.

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;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: where N g p is the observed galaxy number counts in redshift space in the given voxel p. λ g p is the expected number of galaxies at this given position and is related to the final density field δ f p , in redshift space, via where R g p 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. T 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 where the bias function, described by four parameters,N g , the mean density of tracers, and {b g i } = {β, ρ g , 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 δ ic p via the forward model M 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 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.

The augmented joint posterior distribution
The augmented joint posterior distribution corresponds to the following: where the Π's correspond to the respective priors for each parameter. Hence, given our forward model M, 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 presentday 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.

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 0 }, via the following conditional posterior distribution, 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 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 2000Neal , 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: In this work, we set w a = 0 and infer the present-day value w 0 . 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 0 , we perform a rotation of the (Ω m , w 0 ) 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.

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  and Jasche & Wandelt (2013a).
At first, we generate a realisation for the initial density contrast δ ic p 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, assuming a standard Λ cold dark matter (ΛCDM) cosmology with the set of cosmological parameters (Ω m = 0.3089, Ω Λ = 0.6911, Collaboration XIII 2016). This yields a 3D Gaussian initial density field in an equidistant Cartesian grid with N side = 128, consisting of 128 3 voxels, where each voxel corresponds to a discretised volume element with a physical voxel size of 15.6 h −1 Mpc, and comoving box length of 4000 h −1 Mpc. As a result, this implies a total of ∼2.1 × 10 6 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 3 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 δ f p 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 δ f p . 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.
The projection of the completeness functions into the 3D volume produces the 3D observation mask C g p . 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 g p . More specifically, it is the average of the product of the 2D survey geometry C g (x) = C g p(x) and the selection function f (x) at each volume element of the 3D grid: where the volume occupied by the pth voxel is indicated by V p , with |V| the volume of the set V.

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 0 ) parameter space before slice sampling, as described in Appendix F, reduces the burnin period significantly, by roughly a factor of five, as shown in the right panel of Fig. F.1, resulting in improved sampling efficiency.
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 = −0.998 ± 0.008. As a comparison, the SDSS-III (DR12, BAO + Planck) constraints are as follows: Ω m = 0.310 ± 0.005 and w 0 = −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 −1 Mpc for BAO against all voxels in 4000 h −1 Mpc, N BAO = 4000 3 /(4π × 100 3 /3) = 15 278, compared to N vox = 128 3 , yields an improvement of √ N vox /N 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 −1 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.
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 for Ω m and w 0 , 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. 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.
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 selfconsistency 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.
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 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. 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.
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 = −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.
From the correlation matrix of Ω m , w 0 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).
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.

Summary and conclusions
We presented the implementation of a robust AP test that performs a detailed fit of the cosmological expansion via a nonlinear 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 Mean posterior power spectrum Mean prior power spectrum 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. 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 1 .  Fig. 10. Same as Fig. 6, but employing a different prior power spectrum (Ω m = 0.40, w 0 = −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%.  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.
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. . 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. 2004Erdoǧdu et al. , 2006Kitaura & 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 adequate phenomenological choice (e.g. Hubble 1934;Peebles 1980;Gaztanaga & Yokoyama 1993;Kayo et al. 2001), albeit with some limitations, for modelling the evolved matter field in the mildly non-linear regime (Kitaura et al. 2009;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 threepoint 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 matter distribution 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 g i } 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 pervoxel Poissonian likelihood.
Conceptually, within such a framework, the overall aim is to characterise fully the augmented joint posterior P({δ ic p }, {N g }, {b g i }, {θ i }|{N g p }, 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: 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 {δ ic p }, the galaxy bias parameters {b g i }, 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 ofN 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 fourN g and b g i 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: 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: 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 2000Neal , 2003 to ensure unit acceptance rates. 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  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 δ f p can be characterised as a Gaussian random field, is to set them inversely proportional to the variance of that specific δ f p (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;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 −1 , 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 0 ) 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θ, using the transformation: 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 θ = Q † ∆ 1 2θ +θ.

(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: 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 0 ) (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: where n is the number of transition steps in the chain, θ and Var(θ) correspond to the mean and variance, respectively, θ = 1/N i θ i and Var(θ) = 1/N i (θ i − θ ) 2 , 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 0 ) 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.

Appendix G: Adjoint gradient for generic 3D interpolation
The 3D interpolation can be expressed generally, for any arbitrary nth order of interpolation, as: where P is the interpolated surface and E is a (n + 1) 3 × (n + 1) 3 matrix required for computing the vector of interpolation coefficients a. In this work, we implement a triquintic interpolation scheme, which corresponds to n = 5, with the indices {i, j, k}, {α, β, γ} = {0, 1, 2, 3, 4, 5} constituting the particular indexing scheme employed to denote the voxels involved. E is therefore a 216 × 216 matrix for the triquintic scheme.
The above system of equations can be reformulated as a matrix for the linear system described by P = Ea, such that the vector of interpolation coefficients can be computed in straightforward fashion through matrix inversion, a = E −1 P. The advantage of this approach is that E −1 can computed only once and stored, and then used for interpolation at any location inside the cube (e.g. Lekien & Marsden 2005).
The derivative of the generic 3D interpolation is required as a component of the adjoint gradient at the core of the HMC sampler described in Appendix C. The 3D interpolated density field, in redshift space (cf. Eq. (B.1)), can be written explicitly as: where J is the Jacobian of the comoving-redshift transformation from Appendix B, evaluated at the location of the mesh element k. Here, we employed a compact notation, where {x, y, z} correspond to the fractional steps with respect to a reference node in the grid, that is, interpolation weights, and {α, β, γ} label the powers of the elements of E via: i = α(n + 1) 2 + β(n + 1) + γ, (G.3) with the index m ∈ {i, j, k} from Eq. (G.1). Using the chain rule, the resulting gradient w, after application of the incoming gradient u from the previous components in the forward model, is obtained as follows: represents the contribution to the adjoint gradient from a specific voxel labelled by the indexing scheme adopted. This is convenient as E −1 is already available from the execution of the 3D interpolation in the forward model.