© A. Kostić et al. 2022

1. Introduction

One of the most outstanding questions of astrophysical research asks where we might look next to find something new in the universe. This reflects the overall wish to confirm our current physical understanding of the world and find new evidence that could shift prevailing scientific paradigms. The latest advances in cosmological observations and data analysis now unlock the opportunity to perform machine-aided targeted searches for cosmological physics across the Universe. This has become feasible due to the availability of large-scale inferences that are informed by physics and causality of the cosmic large-scale structures and their evolution with time.

State-of-the-art cosmological surveys gather scientific information from tracers of the cosmic large-scale structures, such as galaxies, via large-area homogeneous scans of the sky. The historical origin of such search strategies is inherently related to the theoretical formulation of cosmology. Without any detailed knowledge about the actual spatial distribution of cosmic matter, theorists have resorted to predicting the mean value of statistical summaries, constituting ensemble averages over different realizations of the unknown cosmic matter distribution. Until today, most cosmological analyses have been performed on the basis of two-point statistics (Tegmark 2004; Percival et al. 2007; Porredon et al. 2021), although there are ongoing efforts to go beyond two- or three-point analyses. The initial need to use statistical summaries was twofold. First, the calculation of detailed matter field realizations by evaluating gravitational structure growth via simulations was costly and complex. Therefore, cosmological perturbation theory provided an efficient alternative to make predictions on the average statistical properties of the cosmic structures (Buchert et al. 1994; Bouchet et al. 1995). Secondly, little was known about the particular realization of the cosmic matter distribution of our Universe. The best way forward was to test statements that would be true, on average, using the quantitative diagnostics provided by statistical summaries.

The above situation has witnessed a dramatic improvement in the era of modern cosmology. The abundance of cosmological observations and the availability of computational resources have now contributed significantly to constructing a much more refined picture of the actual spatial configuration of the cosmic structures beyond characterizations made via statistical summaries. In particular, we have proposed the Bayesian physical forward-modeling framework (Jasche & Wandelt 2013; Lavaux & Jasche 2016; Jasche & Lavaux 2019) as a novel method enabling us to reconstruct the full 3D density field underlying observed galaxies in surveys with high fidelity. This machinery performs a causal inference that is informed by physics of the cosmic large-scale structures, their initial conditions and dynamical evolution. By characterizing the full 3D density field, our inference framework exploits the information on the phase distribution and higher-order statistics in the data, which is inaccessible when carrying out traditional approaches. While analyses relying on the use of the information content contained in the phase distribution of the large-scale structures have been proposed for a long time (see, for e.g., Chiang & Coles 2000; Byun et al. 2020), it was not propagated in the same manner as in the aforementioned physical forward-modeling machinery. Therefore, the inference carried out in previous studies by Jasche & Wandelt (2013), Lavaux & Jasche (2016), Jasche & Lavaux (2019) transcends the limitations of contemporary analyses based on statistical summaries in order to yield detailed characterizations of individual 3D structures. Presently, apart from our framework, several research groups are developing the technology to perform full 3D characterizations of the cosmic structures probed by galaxy surveys (Wang et al. 2014; Modi et al. 2018; Kitaura et al. 2021; Porqueres et al. 2020, 2021).

So far, most of the work in the literature is focused on constraining cosmological physics from existing data sets, such as observations of the cosmic microwave background (CMB) (see, for e.g., Planck Collaboration VI 2020, and references therein) and large-scale structure (e.g., Einasto et al. 2011; Byun et al. 2020; Porredon et al. 2021). Here, we want to go beyond the standard task of constraining models to answering the question of how to optimally acquire future data that will be most informative to update our cosmological knowledge and to study physics. In doing so, we use existing information on cosmic structures to identify regions in the sky that promise the highest discovery potential as quantified by the Fisher information. Given the advent of new data analysis technologies, we think that now it is especially timely to start a debate on whether the standard cosmological search strategy should be revised.

We note that some past studies aimed at understanding if certain regions of the sky are particularly informative with regard to cosmology (e.g., Mukherjee & Wandelt 2018) or if one could use such information to design better survey geometries in the future (Bassett 2005). The method proposed here goes beyond these studies by using the full forward model of the 3D matter field for quantifying the information content, weighted by the posterior of plausible realizations of the large-scale structures in the Universe.

In summary, we would like to address the question of whether it is true that a homogeneous scan of the sky is the optimal search strategy for testing fundamental physics with cosmological surveys. We also consider whether there are some particular regions in the sky that are more informative than others with respect to specific research questions and whether it is possible to answer certain research questions faster with small, cheap, and targeted searches. These questions are relevant since, in turn, they raise questions on the optimal use of scientific resources. For example, two upcoming large-scale cosmological surveys, namely, the Vera Rubin Observatory (Ivezić et al. 2019) and Euclid (Racca et al. 2016), constitute a total financial investment of over a billion dollars. As such, to optimize the scientific returns of these missions, we must ensure that the limited survey resources are adequately managed and optimally utilized.

With the ideas outlined in this work, we wish to trigger a discussion in the cosmological community on whether scientific progress in terms of information gain on fundamental physics can be sped up by using more refined cosmological searches, akin to active learning strategies. We also hope to initiate new technological developments for targeted searches of cosmological physics. To make advances in this direction, we address the question of what parts of the Universe ought to be observed to optimally retrieve information about an area of interest for a particular research focus. Specifically, in this work, we are interested in identifying the regions in the Universe that are most relevant to obtain new information about the cosmological parameters underlying our standard model of cosmology. Nonetheless, our approach is equally applicable to other research questions of interest.

We present, for the first time, maps detailing the expected information gain provided by cosmic structures of the nearby Universe on cosmological parameters, based on reconstructions of the 2M++ galaxy catalog (Lavaux & Hudson 2011). Our proposed methodology is not limited to the nearby Universe, however, and is indeed applicable to any existing cosmological data sets.

2. Results

In our quest for mapping out the cosmologically sensitive large-scale structures in the sky, we adopted the Fisher information (Fisher 1925) methodology to provide a quantitative measure of the information encoded by the cosmological parameters on the observed large-scale structures. Although computing the Fisher information in this context is a highly non-trivial task, this ambitious undertaking is rendered feasible by employing a physical forward-modeling machinery, such as our BORG algorithm (Jasche & Wandelt 2013; Jasche & Lavaux 2019), which constitutes a causal model of structure formation with a fully non-linear treatment of the dark matter clustering. It allows for making a connection between the cosmic initial conditions and the observed galaxy distribution. A schematic view of the distinct components of the BORG forward model is illustrated in Fig. 1. In our study, we marginalize over the ensemble of plausible realizations of the 3D primordial matter fluctuations of the very early Universe, as inferred by BORG via a hierarchical Bayesian statistical inference framework, from the 2M++ galaxy catalog (Lavaux & Hudson 2011) that traces the matter distribution of the nearby Universe. We refer the reader to the appendices for more detailed information pertaining to the BORG algorithm, 2M++ galaxy catalog, mathematical formalism underlying the Fisher information map, and numerical implementation thereof. We stress that our approach for deriving the Fisher information map, based on the constrained realizations of cosmic structures conditional on galaxy observations, bears a stark contrast to the standard Fisher analyses for obtaining forecasts on cosmological constraints from forthcoming galaxy surveys.

Fig. 1. Schematic of the Bayesian physical forward-modeling framework of BORG, which solves a large-scale Bayesian inverse problem by fitting a dynamical structure formation model to galaxy observations and subsequently inferring the primordial initial conditions (ICs) that lead to the formation of the presently observed cosmic structures via gravitational evolution. The BORG forward-modeling approach naturally marginalizes over unknown galaxy bias and accounts for all relevant physical effects, such as redshift space distortions (RSDs) resulting from the peculiar velocities of galaxies, as well as instrumental selection effects.

The availability of such inferred primordial matter fluctuations enables us to test the causal sensitivity of cosmic structures when forward-modeled with perturbed cosmological parameters. To illustrate this possibility, in Fig. 2, we illustrate the response of cosmic structures in the Universe to perturbative changes in six parameters of the so-called concordance ΛCDM cosmological model, namely, the matter density, Ω_m, baryon density, Ω_b, cosmic curvature, Ω_k, Hubble constant, h, amplitude of matter fluctuations, σ₈, and scalar spectral index, n_s, of the primordial power spectrum. This entails computing the gradient of the 3D galaxy field with respect to the cosmological parameters. This gradient quantifies the sensitivity of the galaxy distribution to changes in the cosmological parameters. As such, the gradient is described by a 6D vector for each volume element in the 3D grid. We computed the gradient using finite differencing by executing BORG forward model evaluations on the ensemble of data constrained realizations of the cosmic initial conditions from the BORG 2M++ analysis (Jasche & Lavaux 2019), while varying the cosmological parameters about their corresponding fiducial values, as given by the latest best-fit values from the Planck Collaboration (Planck Collaboration VI 2020).

Fig. 2. Components of the gradient of the matter density field with respect to the cosmological parameters for a spherical slice of thickness ∼2.65 h−1 Mpc at a comoving distance of 100 h−1 Mpc from the observer, with the corresponding galaxies from the 2M++ catalog denoted via black dots. Visually, we find that the baryon density Ωb has the largest influence on the spatial distribution of the cosmic structures, with the regions surrounding the filamentary galaxy distribution being particularly sensitive to changes in the baryon density. Conversely, the least significant response emanates from the cosmic curvature Ωk, with only the vicinity of the dense galaxy clusters reacting to changes in the geometry of the Universe.

In Fig. 2, we visualize the individual gradient components of the six cosmological parameters by representing them as “cosmological sensitivity maps”. To this end, we consider a spherical slice of thickness ∼2.65 h⁻¹ Mpc at a comoving distance of 100 h⁻¹ Mpc from the observer that is projected onto a HEALPIX map (Górski et al. 2005). The observed galaxies from the 2M++ catalog lying in the projected spherical slice are also indicated. From a visual comparison, we find that the density of baryonic matter, as characterized by the Ω_b parameter, induces the most significant response in the observed distribution of cosmic structures, with the Ω_m and σ₈ parameters also having a notable influence. The astrophysical regions surrounding the filamentary patterns traced by the galaxies respond effectively to changes in the latter cosmological parameters. In contrast, the cosmic matter distribution displays a relatively minimal sensitivity to the Ω_k parameter, with the most noticeable effect of the cosmic curvature on the large-scale structures being on the densest regions, such as the galaxy clusters. The underdense regions of the matter distribution, the so-called cosmic voids, are most impacted by changes in the Ω_b and σ₈ parameters. We clarify that our structure formation model does not account for baryonic physics, so that Ω_b only modifies the amplitude of growth of matter perturbations and the shape of the linear power spectrum. Nevertheless, it is straightforward to repeat our analysis with hydrodynamical simulations to properly account for the presence of baryons. Given the significantly higher computational cost of such simulations, one possible option would be to use physics emulators, which have recently been developed on a large scale (Villaescusa-Navarro et al. 2021).

In addition to displaying the relative strength of the causal response of the cosmic structures to changes in our underlying cosmological model, the array of sensitivity maps in Fig. 2 reveals some interesting features. The σ₈ map, for example, illustrates how the clustering of matter occurs at the expense of the neighbouring regions, while the h map shows that a modified Hubble flow yields different clustering features. Moreover, the Ω_m and Ω_b maps display distinct anti correlated signatures. This may be attributed to the fact that the sound horizon distance scale due to baryon acoustic oscillations (BAOs) increases with both Ω_m and Ω_b, with the effect induced by Ω_b being much stronger, such that a fixed BAO scale encoded in the data yields this anti correlation between these two maps. Similarly, the Ω_m and σ₈ maps depict a striking anti correlation. This can be naturally explained on the basis that the primary constraining power emanates from the combination of the growth rate, f, of cosmological perturbations and σ₈, where $f{\sigma }_8\approx {\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.55}$ for the ΛCDM model.

A crucial ingredient in the Fisher information formalism derived in our study, as described by Eq. (C.12), lies in computing the above gradient. The desired 3D Fisher information field, obtained using Eq. (C.12), can be represented as a “Fisher information map” in the same way as the cosmological sensitivity maps. This Fisher information map represents the combined information gain on all six cosmological parameters considered in this study. The Fisher information map for a spherical slice of thickness ∼2.65 h⁻¹ Mpc at a distance of 100 h⁻¹ Mpc from the observer is displayed in Fig. 3, along with the observed galaxy distribution in this particular slice. The Fisher information map indicates that the regions with the highest information gain are those in the vicinity of the filamentary distribution traced by galaxies. One plausible interpretation is that the regions of gravitational infall, where there is an accumulation of matter, are the most informative with regard to our cosmological model. This is particularly interesting as these infalling regions are still relatively obscure and must be properly understood to make further progress.

Fig. 3. Fisher information map for the same spherical slice as in Fig. 2. The observed galaxy distribution from the 2M++ catalog, lying in the corresponding spherical shell centered around the observer, is represented by the red dots. We find that the regions in the vicinity of the massive cosmic structures, as traced by the galaxy distribution, are the most informative according to the Fisher information map. These regions correspond to the regime of gravitational infall of the galaxy clusters.

The BORG analysis of the 2M++ galaxy catalog showcased the capacity of the physical forward-modeling machinery to resolve the key features of prominent cosmic structures in the present Universe (Jasche & Lavaux 2019). The inferred mass profile of the Coma cluster, in particular, was found to be in remarkable agreement with state-of-the-art weak lensing measurements. Therefore, it is also interesting to study the characteristic features underlying the source of cosmological information for this well-known cluster. Figure 4 displays the Fisher information map of the Coma cluster and its corresponding mass density. We note that the figure depicts the central slice of thickness ∼5 h⁻¹ Mpc through a 3D patch extending over 40 h⁻¹ Mpc. The features present in the Fisher map of the Coma cluster are in accordance with those from the sky projected Fisher map from Fig. 3, thereby supporting the interpretation that the information that can be gleaned from the regions surrounding the filaments, namely, the gravitationally infalling regions of the cluster, provides the most significant information gain. Recent studies relying on conceptually distinct methodologies have also indicated that accretion filaments and the surroundings of voids are highly sensitive to predictions of dark energy (Leclercq et al. 2016) and gravity (Lam et al. 2012) models. In contrast, the central core of the cluster provides a relatively lower information gain as quantified by Fisher information, with the regions devoid of matter displaying insubstantial information. A key point worth stressing here is that our relatively simple physics informed algorithm was capable of pinning down potential regions of cosmological interest in the sky.

Fig. 4. Fisher information map of the Coma cluster, with its corresponding mass density overlaid as a contour, for the central slice of thickness ∼5 h−1 Mpc through a 3D region centered on the cluster that extends over 40 h−1 Mpc. According to the Fisher information map, the central region containing the core of the cluster encodes fairly limited information gain, whilst the peripheral regions of the filaments and cluster core, where mass accretion occurs via gravitational infall, constitute the greatest proportion of information gain.

Now, it is possible to use the maps displayed in Figs. 3 and 4 as a guide for where to look in order to optimally collect data for testing our physical model of structure formation. For the particular case demonstrated here, the data consists of galaxy counts and the corresponding information gain predictions are geared at the galaxy clustering data. The idea behind the targeted search approach is to recursively search for galaxies in the high information gain regions, as quantified by the Fisher map, and then to use this newly acquired data to constrain the model parameters and repeat the procedure. This is schematically depicted in Fig. 5. The justification pertaining to why this search strategy is optimal is further elaborated in Appendix D. It should be understood that this approach can be similarly extended to other observables for constructing the corresponding information gain maps.

Fig. 5. Flowchart of the targeted search approach. The core idea is to use the existing knowledge as quantified by the posterior p(θ, ϕ|d) in order to calculate the Fisher map marginalized over this posterior, ⟨ℐ(θ|d0)⟩(ϕ|θ, d0), while keeping the observables of interest. This allows this Fisher map to provide us with regions of the sky with highest information gain potential for acquiring new data that are optimal for testing our model predictions, ⟨ℐ(θ|d0)⟩(ϕ|θ, d0): d0 → d. This data will in turn be used to update our knowledge about the model through updating the posterior and the procedure would repeat until the information content is fully depleted.

3. Discussion

The results presented in this work demonstrate the feasibility of machine-aided targeted searches for cosmological physics signals. These have become feasible through recent developments of physics-informed causal inference frameworks to study the 3D cosmic large-scale structures, their origin, and evolution over time. We used 3D initial conditions, inferred with the BORG algorithm, and a physics simulator of cosmological structure formation to chart the response of observed structures in the Universe and their corresponding information gain with respect to cosmological parameters. Our results are the first of their kind and elucidate the inhomogeneous distribution of cosmological information in the Universe. This study paves a new way forward to perform efficient targeted searches for the fundamental physics of the Universe, where search strategies are progressively refined with new cosmological data sets within an active learning framework.

We have further illuminated the response of the cosmic large-scale structures with respect to individual cosmological parameters. These results suggest that different features of the cosmic structures respond differently to perturbations in the physics determined by the cosmological parameters. For instance, we find that the vicinity of the filamentary cosmic structures, corresponding to gravitationally infalling regions, are highly sensitive to changes in the baryon density, with the cosmic curvature impacting only the surroundings of the dense galaxy clusters. Cosmic voids are primarily affected by the amplitude of matter fluctuations and the baryon density. Our results demonstrate the value of going beyond state-of-the-art analyses of the cosmic large-scale structures that are limited to summary statistics and ignore this richness of the 3D cosmic structures. Even though we consider one particular observable, namely galaxy counts, in this study, our proposed framework can be seamlessly applied to other tracers, such as Lyman-α forest (Porqueres et al. 2019a, 2020), which would bear complementary information.

Our findings further suggest that optimal targeted searches for research questions have become feasible. Given a specific research question, our approach enables us to propose targets for optimal information retrieval. This raises the question if traditional survey strategies should be revised and if significant scientific progress could also be driven by targeted searches with smaller, cheaper, and faster instrumentation. We hope that our contribution will trigger a discussion and new technological advances in the field.

While we studied the influence of cosmological parameters on the Universe in this work, our proof of concept reveals the potentially far-reaching and groundbreaking implications when considering the physical effects induced by modified gravity (Koyama 2016), dynamical dark energy (Zhao et al. 2017), and exotic dark matter models, such as self-interacting (Carlson et al. 1992) and fuzzy dark matter (Hu et al. 2000) or massive neutrinos (Villaescusa-Navarro et al. 2014), on the cosmic large-scale structures. Once we identify the astrophysical region(s) of interest, based on the Fisher information map, for a particular model, we may subsequently proceed by computing accurate theoretical predictions for the spatial distribution of matter and luminous tracers for the given model via cosmological N-body or hydrodynamical simulations with extremely high resolution, thereby providing highly detailed physical features. The final step would then entail observing the relevant region(s) and comparing the theoretical predictions with the galaxy observations in a likelihood or posterior analysis. In essence, we are proposing a novel way of doing science that optimizes scarce resources for maximal scientific returns via an efficient observational strategy.

Acknowledgments

We would like to express our appreciation to Fabian Schmidt, Benjamin Wandelt, Eleni Tsaprazi, Natalia Porqueres, Minh Nguyen, Radek Wojtak, Florent Leclercq and Harry Desmond for their remarks on our manuscript. A.K. acknowledges support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” of the European Research Council at MPA. J.J. acknowledges support by the Swedish Research Council (VR) under the project 2020-05143 – “Deciphering the Dynamics of Cosmic Structure”. DKR is a DARK fellow supported by a Semper Ardens grant from the Carlsberg Foundation (reference CF15-0384). This work was supported by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence Nationale de la Recherche. This work has made use of the Horizon Cluster hosted by Institut d’Astrophysique de Paris. This work has been done within the activities of the Domaine d’Intérêt Majeur (DIM) Astrophysique et Conditions d’Apparition de la Vie (ACAV), and received financial support from Région Ile-de-France. We thank Cepheid Studio (https://www.cepheidstudio.com) for providing us with the telescope illustration used in Fig. 1. This work is done within the Aquila Consortium (https://aquila-consortium.org). Author contributions. A.K.: led the project; methodology; software; obtained, validated, and interpreted results; writing – editing. J.J.: project conceptualization; methodology; validation and interpretation; supervision; writing - editing. D.K.R.: methodology; visualization; validation and interpretation, writing – editing. G.L.: methodology; validation and interpretation; supervision; provided resources; funding acquisition, editing.

References

Appendix A: 2M++ galaxy catalog

The 2M++ catalog (Lavaux & Hudson 2011) is a compilation galaxy redshifts derived from the Two-Micron-All-Sky-Survey (2MASS) Redshift Survey (2MRS) (Huchra et al. 2012), the Six-Degree-Field Galaxy Redshift Survey Data Release 3 (6dFGRS-DR3) (Jones et al. 2009), and the Sloan Digital Sky Survey Data Release 7 (SDSS-DR7) (Abazajian et al. 2009). The resulting catalog has a greater depth and a higher sampling rate relative to the Infrared Astronomical Satellite (IRAS) Point Source Catalog Redshift Survey (PSCZ) (Saunders et al. 2000). The photometry is based on the 2MASS Extended Source Catalog (2MASS-XSC) (Skrutskie et al. 2006), an all-sky survey in the J, H and K_S bands, with redshifts in the K_S band of the 2MRS complemented by those from the SDSS-DR7 and 6dFGRS-DR3.

Since the 2M++ catalog is a combination of several surveys, the galaxy magnitudes from all sources were first recomputed by measuring the apparent magnitude in the K_S band within a circular isophote at 20 mags arcsec⁻². The apparent K_S band magnitudes were subsequently corrected by taking into account Galactic extinction, cosmological surface brightness dimming, stellar evolution and k-corrections, while masking the Galactic Plane. To account for the incompleteness due to fibre collisions in 6dFGRS and SDSS, the redshifts of nearby galaxies within each survey region are cloned. The final 2M++ catalog contains 69190 galaxies in total, and is fairly sampling the galaxies of the cosmic volume up to a distance of 200 h⁻¹ Mpc for the area covered by the 6dFGRS and SDSS, and up to 125 h⁻¹ Mpc for the region mapped by 2MRS. For a more in-depth description of the construction of the 2M++ catalog, we refer the interested reader to the original compilation (Lavaux & Hudson 2011), with the computation of radial selection functions and target selection completeness as required for the BORG 2M++ analysis (Jasche & Lavaux 2019) detailed in Sect. 4.1 thereof.

Appendix B: BORG

The Bayesian Origin Reconstruction from Galaxies (BORG) algorithm (Jasche & Wandelt 2013; Jasche & Lavaux 2019) constitutes a hierarchical Bayesian statistical inference framework to infer the initial conditions of the Universe from the observed galaxy positions in the sky as provided by galaxy redshift surveys. By encoding a physical forward model at its core, as illustrated in Fig. 1, BORG allows for a detailed 3D reconstruction and characterization of observed cosmic large-scale structures.

The forward model entails the sequential application of several components that causally relate the cosmic initial conditions to the observed galaxy distribution, starting with a physical description of the non-linear dynamics involved in gravitational structure formation. While the BORG machinery incorporates several models of structure growth, such as the approximate first- and second-order Lagrangian perturbation theory and fully non-linear particle mesh models, we employ the latter option to evolve the primordial density fluctuations to their corresponding 3D dark matter distribution via cosmological N-body simulations. The predicted dark matter fields must subsequently be connected to a galaxy population via an adequate galaxy bias model. Although the galaxy biasing effect is a presently challenging and as yet unresolved issue in cosmology, the 2M++ analysis made use of a local but non-linear truncated power law bias model (Neyrinck et al. 2014) as motivated by numerical simulations. Specifically, BORG inherently accounts for the physical effects such as the geometric distortion due to the cosmic expansion and the redshift space distortions due to the peculiar velocities of galaxies. The final component in the forward model deals with all the relevant observational and instrumental effects, such as survey geometry, selection effects, and foreground contamination, to yield the desired galaxy distribution that can be compared to observations via a suitable likelihood, with a Poissonian distribution adopted in the 2M++ analysis (Jasche & Lavaux 2019).

The above physical forward model results in a highly non-trivial Bayesian inverse problem. To efficiently sample the high-dimensional and non-linear parameter space of plausible initial conditions at an earlier epoch, with typically 𝒪(10⁷) free parameters corresponding to the discretized volume elements of the observed domain, BORG relies on a Hamiltonian Monte Carlo (HMC) technique (Jasche & Kitaura 2010; Jasche & Wandelt 2013; Jasche & Lavaux 2019). In essence, BORG performs the joint inference of the respective posterior distributions of the initial conditions and the bias parameters given the observed galaxy distribution via a modular Markov Chain Monte Carlo (MCMC) sampling framework. Feeding the MCMC samples of the initial conditions to the forward model therefore yields physically plausible realizations of the 3D cosmic structures, namely, non-linearly evolved density fields and associated velocity fields underlying the spatial distribution of observed galaxies.

Recent extensions to the BORG framework have led to substantial improvements in cosmological parameter inference (Kodi Ramanah et al. 2019; Elsner et al. 2020; Schmidt et al. 2020), Lyman-α (Porqueres et al. 2019a, 2020) and cosmic shear (Porqueres et al. 2021) reconstructions, with the field-level treatment transcending the capabilities of conventional cosmological analyses. Novel sophisticated additions to the forward model include a robust likelihood to account for unknown foreground contamination and systematics (Porqueres et al. 2019b) and machine learning-based galaxy bias models (Charnock et al. 2020). The reconstructed density and velocity fields from BORG analyses have been employed in investigating extensions of the standard model of particle physics and cosmology (Desmond et al. 2018, 2019; Pardo et al. 2019; Desmond & Ferreira 2020; Bartlett et al. 2021) and for unbiased and accurate measurements of the Hubble constant from gravitational wave events (Mukherjee et al. 2021), as well as those of the kinematic Sunyaev-Zel’dovich effect from CMB observations (Nguyen et al. 2020).

Appendix C: Fisher information maps

The proposed targeted search approach requires identifying regions of the observational domain whose observations can thus provide optimal information gain to update existing knowledge. To quantify this information gain, we chose to employ the Fisher information (Fisher 1925) approach. This approach measures the amount of information that unseen observable data can carry about uncertain model parameters. More explicitly, Fisher information estimates the squared norm of the likelihood score for given model parameters averaged over all possible future data realizations permitted by the likelihood. It therefore measures the expected strength with which the likelihood will respond to changes in the model parameters once new data becomes available. For the specific case considered in this work, we assume that prior knowledge on the spatial cosmic matter configuration is available from previous observations. More specifically, we assume that the white noise realizations (the phases) of the initial conditions ϕ are provided by the BORG reconstruction of the 2M++ survey (Jasche & Lavaux 2019). Assuming this realization of initial conditions we can express the conditional Fisher information as:

$$\begin{array}{cc}\hfill \mathcal{I}(\mathit{\theta }|\mathit{\varphi })& ={\mathbb{E}}_{(\mathit{d}|\mathit{\theta },\mathit{\varphi })}\left[{\left(\frac{\partial ln(\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }))}{\partial \mathit{\theta }}\right)}^2\right]\hfill \\ \hfill & ={\displaystyle \int }\mathcal{D}\mathit{d}{\left(\frac{\partial ln(\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }))}{\partial \mathit{\theta }}\right)}^2\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }),\hfill \end{array}$$(C.1)

where θ corresponds to the set of cosmological parameters parametrizing our forward model, with the galaxy observations within the voxels denoted by d = [N_i]_{i = 1, …Nbox} and N_box corresponds to the number of the voxels inside the 3D volume considered in the 2M++ BORG reconstruction (see Jasche & Lavaux 2019, for further details). We note that we express the vector nature of these random variables using boldface symbols. For the sake of this work, we assume a Poisson likelihood to describe the galaxy clustering data:

$$\begin{array}{c}\hfill \mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi })={\displaystyle \prod _i}{\mathrm{e}}^{-{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{{\lambda }_i{(\mathit{\theta };\mathit{\varphi })}^{N_i}}{N_i!},\end{array}$$(C.2)

where λ_i(θ; ϕ) denotes the rate of the specific realization of the Poisson process as dictated by the initial phases ϕ and cosmological parameters θ, while N_i represents the galaxy counts within the i-th cell of the 3D volume. Now, taking the derivative with respect to the cosmological parameters θ gives:

$$\begin{array}{cc}\hfill \frac{\partial ln(\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }))}{\partial \mathit{\theta }}& ={\displaystyle \sum _i}(-\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}+\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }})\hfill \\ \hfill & ={\displaystyle \sum _i}(\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}-1)\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}.\hfill \end{array}$$(C.3)

The square of the above derivative is given by:

$$\begin{array}{cc}\hfill {\left(\frac{\partial ln(\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }))}{\partial \mathit{\theta }}\right)}^2& ={\displaystyle \sum _{\mathit{ij}}}\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\frac{\partial {\lambda }_j(\mathit{\theta },\mathit{\varphi })}{\partial \mathit{\theta }}\hfill \\ \hfill & \times (\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{N_j}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}-\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}-\frac{N_j}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}+1)\hfill \\ \hfill & ={\displaystyle \sum _i}{\left(\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2(\frac{N_i^2}{{\lambda }_i^2(\mathit{\theta },\mathit{\varphi })}-2\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}+1)\hfill \\ \hfill & +{\displaystyle \sum _{\mathit{ij}}}(1-{\delta }_{\mathit{ij}}^K)\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\frac{\partial {\lambda }_j(\mathit{\theta },\mathit{\varphi })}{\partial \mathit{\theta }}\hfill \\ \hfill & \times (\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{N_j}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}-\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}-\frac{N_j}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}+1),\hfill \end{array}$$(C.4)

where we split the sum in diagonal and off-diagonal terms. Again, the indices i and j denote the pixels in the HEALPIX projection. Since given a density field, the individual Poisson realizations are independent, the off-diagonal terms will vanish when computing the Fisher information as given by Eq. (C.1). In order to show this, we consider the following data average of Eq. (C.4):

$$\begin{array}{c}\hfill {\langle {\left(\frac{\partial ln(\mathcal{L}(\mathit{d}|\mathit{\theta },\mathit{\varphi }))}{\partial \mathit{\theta }}\right)}^2\rangle }_{(\mathit{d}|\mathit{\theta },\mathit{\varphi })}={\displaystyle \sum _{\mathit{ij}}}\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\frac{\partial {\lambda }_j(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\\ \hfill \times {\langle (\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{N_j}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}-\frac{N_i}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}-\frac{N_j}{{\lambda }_j(\mathit{\theta };\mathit{\varphi })}+1)\rangle }_{(\mathit{d}|\mathit{\theta },\mathit{\varphi })}\\ \hfill =\sum _i{\left(\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2(\frac{{\lambda }_i{(\mathit{\theta };\mathit{\varphi })}^2+{\lambda }_i(\mathit{\theta };\mathit{\varphi })}{{\lambda }_i^2(\mathit{\theta };\mathit{\varphi })}-2\frac{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}+1)\\ \hfill +\sum _{\mathit{ij}}(1-{\delta }_{\mathit{ij}}^K)\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\frac{\partial {\lambda }_j(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\\ \hfill \underset{=0}{\underbrace{(\frac{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\frac{{\lambda }_j(\mathit{\theta };\mathit{\varphi })}{{\lambda }_j(\mathit{\theta },\mathit{\varphi })}-\frac{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}-\frac{{\lambda }_j(\mathit{\theta };\mathit{\varphi })}{{\lambda }_j(\mathit{\theta };\mathit{\varphi })}+1)}},\end{array}$$(C.5)

where we used the following identities:

$$\begin{array}{cc}\hfill {⟨N_i⟩}_{(\mathit{d}|\mathit{\theta },\mathit{\varphi })}& ={\lambda }_i(\mathit{\theta };\mathit{\varphi }),\hfill \\ \hfill {⟨N_i^2⟩}_{(\mathit{d}|\mathit{\theta },\mathit{\varphi })}& ={\lambda }_i(\mathit{\theta };\mathit{\varphi })+{\lambda }_i^2(\mathit{\theta };\mathit{\varphi })\mathrm{for}\forall i\in [1,N].\hfill \end{array}$$(C.6)

The first equality follows from the conditional independence of the Poisson realizations within each grid cell, which always holds once the underlying density field λ(θ, ϕ) is given; the second results from the expression of variance of the Poisson distribution. Therefore, the diagonal terms are:

$$\begin{array}{c}\hfill \mathcal{I}(\mathit{\theta }|\mathit{\varphi })={\displaystyle \sum _i}{\left(\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2\frac{1}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}.\end{array}$$(C.7)

From the above equation, we see that the contribution from a particular volume element to the Fisher information is given by:

$$\begin{array}{c}\hfill h_i^{\prime }(\mathit{\theta }|\mathit{\varphi })={\left(\frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2\frac{1}{{\lambda }_i(\mathit{\theta };\mathit{\varphi })}\equiv {\lambda }_i(\mathit{\theta },\mathit{\varphi }){\left(\frac{\partial ln{\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2.\end{array}$$(C.8)

This allows us to obtain a 3D Fisher map. Now, we also have to specify the Poisson intensity in terms of the output of our physics simulator. This relation can be written as:

$$\begin{array}{c}\hfill {\lambda }_i(\mathit{\theta };\mathit{\varphi })=B(\delta \equiv G(\mathit{\theta };\mathit{\varphi })),\end{array}$$(C.9)

where B(x) is an arbitrary non-negative bias function and G(θ, ϕ) is a physics simulator of the cosmic large-scale structures with output δ corresponding to the 3D matter density contrast amplitudes. For illustrative purposes, we assume B(x) = 1 + x. We note that any bias model monotonic in the density will change the quantitative results – but not the qualitative results, as can also be seen from Eq. (C.8). We then obtain:

$$\begin{array}{c}\hfill {\lambda }_i(\mathit{\theta };\mathit{\varphi })=1+G_i(\mathit{\theta };\mathit{\varphi }),\end{array}$$(C.10)

which leads to the following expression for the derivative with respect to cosmological parameters θ:

$$\begin{array}{c}\hfill \frac{\partial {\lambda }_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}=\frac{\partial G_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}.\end{array}$$(C.11)

We can therefore approximate the Fisher information elements through finite differencing as:

$$\begin{array}{cc}\hfill h_i^{\prime }(\mathit{\theta };\mathit{\varphi })& ={\left(\frac{\partial G_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2\frac{1}{1+G_i(\mathit{\theta };\mathit{\varphi })}\hfill \\ \hfill & \approx {\left(\frac{G_i({\mathit{\theta }}_0+\mathrm{\Delta }\mathit{\theta };\mathit{\varphi })-G_i({\mathit{\theta }}_0;\mathit{\varphi })}{\mathrm{\Delta }\mathit{\theta }}\right)}^2\frac{1}{1+G_i({\mathit{\theta }}_0;\mathit{\varphi })}.\hfill \end{array}$$(C.12)

We note, however, that the above calculation is done for a fixed realization of the phases ϕ. Since we do not know what particular realization is compatible with our Universe, we need to marginalize over them:

$$\begin{array}{cc}\hfill h_i(\mathit{\theta }|d)& ={\displaystyle \int }\mathcal{D}\mathit{\varphi }p(\mathit{\varphi }|d)h_i^{\prime }(\mathit{\theta }|\mathit{\varphi })\hfill \\ \hfill & \equiv {\displaystyle \int }\mathcal{D}\mathit{\varphi }p(\mathit{\varphi }|d){\left(\frac{\partial G_i(\mathit{\theta };\mathit{\varphi })}{\partial \mathit{\theta }}\right)}^2\frac{1}{1+G_i(\mathit{\theta };\mathit{\varphi })}.\hfill \end{array}$$(C.13)

This task is made possible by the BORG algorithm, as described in Appendix B. In this way, we are able to use the information content of the previously obtained constraints and update the Fisher information. Projecting this 3D Fisher map onto a HEALPIX grid and visualizing a particular spherical slice through it, we obtain an all-sky map which we refer to in our study as the "Fisher information map." A specific example is provided in Fig. 3.

Now, in order to calculate the G_i(θ₀ + Δθ; ϕ) term above, we must first perform a forward model evaluation by setting the cosmological parameters to their corresponding fiducial values: θ₀ ≡ {Ω_m = 0.3111, Ω_b = 0.049, Ω_k = 7 × 10⁻⁴, h = 0.6766, σ₈ = 0.8102, n_s = 0.9665}, according to the latest Planck best-fit ΛCDM cosmology (Planck Collaboration VI 2020), for all the MCMC realizations of the initial conditions from the BORG 2M++ analysis. We then compute the mean galaxy field marginalized over the forward model output realizations. We subsequently repeat this procedure for a new set of perturbed cosmological parameters, θ′≡θ₀ + Δθ, ensuring Δθ does not exceed values outside 1-sigma width of a Gaussian centered at the fiducial cosmology θ₀ as characterized by Planck best-fit values. Using the respective means of the fiducial and perturbed galaxy fields, we compute the gradient of the forward model output with respect to the cosmological parameters, Δ_G/Δ_θ, using a finite differencing scheme, with:

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_G& ={\displaystyle \int }\mathcal{D}\mathit{\varphi }p(\mathit{\varphi }|d)[G({\mathit{\theta }}_0+\mathrm{\Delta }\mathit{\theta };\mathit{\varphi })-G({\mathit{\theta }}_0;\mathit{\varphi })]\hfill \\ \hfill & \approx \frac{1}{N_{\mathit{\varphi }}}{\displaystyle \sum _i}[G({\mathit{\theta }}_0+\mathrm{\Delta }\mathit{\theta };{\mathit{\varphi }}_i)-G({\mathit{\theta }}_0;{\mathit{\varphi }}_i)]\text{with}{\mathit{\varphi }}_i\hookleftarrow p(\mathit{\varphi }|\mathit{d})\hfill \end{array}$$

and Δθ = θ′−θ₀. The distinct gradient components correspond to the cosmological sensitivity maps displayed in Fig. 2. The gradient is then squared and divided by the mean fiducial galaxy field and marginalized over phase realizations, as required by Eq. (C.12), to finally yield the desired Fisher information elements. The resulting Fisher information map for a particular spherical slice through the 3D field is illustrated in Fig. 3.

In this approach, we do not explicitly enforce that ∑_iΩ_i = 1 for θ′, since we are interested in infinitesimal and independent variations of the cosmological parameters as required by the Fisher information. Therefore, although this is marginally inconsistent from the cosmological perspective, it is perfectly consistent from the information theory perspective. Furthermore, enforcing the ∑_iΩ_i = 1 in our forward model introduces dependencies of the variations required for the evaluation of the Fisher information, which renders the interpretation more convoluted. Nonetheless, we also tested the impact of enforcing ∑_iΩ_i = 1 in our forward model and found that this has a negligible effect on the derived Fisher information map. In addition, we verified the robustness of our results, as expected, to the choice of fiducial cosmology.

Appendix D: Marginalizing phase realizations and targeted searches

The derivation presented in the previous section took care of explicitly keeping, where needed, the conditional dependence on the specific realization of the phases of the initial density field, also referred to as the initial white noise field. The reason being that in practice, the exact phase information of the initial density field is not available:

$$\begin{array}{c}\hfill \mathcal{I}(\mathit{\theta })={\displaystyle \int }\mathrm{d}\mathit{\varphi }\mathcal{P}(\mathit{\varphi })\mathcal{I}(\mathit{\theta }|\mathit{\varphi }),\end{array}$$(D.1)

where 𝒫(ϕ) is a Gaussian prior with zero mean and unit variance for the initial white noise field. We effectively use a Markov approximation of the integral by drawing random realizations from the white noise prior and evaluating the corresponding averaged Fisher information as:

$$\begin{array}{c}\hfill \mathcal{I}(\mathit{\theta })=\frac{1}{N_{\mathit{\varphi }}}{\displaystyle \sum _i}\mathcal{I}(\mathit{\theta }|{\mathit{\varphi }}_i).\end{array}$$(D.2)

where the ϕ_i are independent white noise realizations. Similarly, we can just decompose the averaged Fisher information into the individual components per volume element:

$$\begin{array}{c}\hfill h_j(\mathit{\theta })=\frac{1}{N_{\mathit{\varphi }}}{\displaystyle \sum _i}{h}_j(\mathit{\theta }|{\mathit{\varphi }}_i).\end{array}$$(D.3)

It is clear that since we marginalize over all possible white noise realizations, there cannot be any variation in the h_i(θ) since every configuration is equally likely. Thus, without any further prior information on the phases, every point in the Universe, on average, is expected to contribute exactly the same amount of information.

This is one of the logical reasons behind performing homogeneous cosmological surveys. However, we are no longer in a regime of cosmology where complete ignorance about the Universe prevails. We have access to data, which can inform us about the specific realization of our Universe. If we account for this fact, we can condition the average Fisher information on the existing data:

$$\begin{array}{c}\hfill \mathcal{I}(\mathit{\theta }|\mathit{d})={\displaystyle \int }\mathrm{d}\mathit{\varphi }\mathcal{P}(\mathit{\varphi }|\mathit{d})\mathcal{I}(\mathit{\theta }|\mathit{\varphi }).\end{array}$$(D.4)

Now, 𝒫(ϕ|d) is the data constrained posterior of the initial phase distribution. Obtaining this posterior distribution is a non-trivial task, but one that is solved by our BORG algorithm. In particular, BORG provides a Markov approximation to the high-dimensional posterior distribution 𝒫(ϕ|d), such that we can approximate the integral as:

$$\begin{array}{c}\hfill \mathcal{I}(\mathit{\theta }|\mathit{d})=\frac{1}{N_{\mathit{\varphi }}^{\mathit{borg}}}{\displaystyle \sum _i}\mathcal{I}(\mathit{\theta }|{\mathit{\varphi }}_i^{\mathit{borg}}),\end{array}$$(D.5)

with this inferred posterior of initial phase distributions being fairly robust to the details of the physical model adopted (Nguyen et al. 2021). Similarly, we can use the results of BORG to estimate the Fisher information per volume element as:

$$\begin{array}{c}\hfill h_j(\mathit{\theta }|\mathit{d})=\frac{1}{N_{\mathit{\varphi }}^{\mathit{borg}}}{\displaystyle \sum _i}{h}_j(\mathit{\theta }|{\mathit{\varphi }}_i^{\mathit{borg}}).\end{array}$$(D.6)

The key point to realize here is that the result does not depend on the white noise phases ϕ as they have been marginalized out. More precisely, this has been achieved by evaluating Eq. (D.6) with $N_{\mathit{\varphi }}^{\mathrm{B}\mathrm{O}\mathrm{R}\mathrm{G}}=2422$ samples characterizing the posterior 𝒫(ϕ|d).

All Figures

Fig. 1. Schematic of the Bayesian physical forward-modeling framework of BORG, which solves a large-scale Bayesian inverse problem by fitting a dynamical structure formation model to galaxy observations and subsequently inferring the primordial initial conditions (ICs) that lead to the formation of the presently observed cosmic structures via gravitational evolution. The BORG forward-modeling approach naturally marginalizes over unknown galaxy bias and accounts for all relevant physical effects, such as redshift space distortions (RSDs) resulting from the peculiar velocities of galaxies, as well as instrumental selection effects.

In the text

Fig. 2. Components of the gradient of the matter density field with respect to the cosmological parameters for a spherical slice of thickness ∼2.65 h−1 Mpc at a comoving distance of 100 h−1 Mpc from the observer, with the corresponding galaxies from the 2M++ catalog denoted via black dots. Visually, we find that the baryon density Ωb has the largest influence on the spatial distribution of the cosmic structures, with the regions surrounding the filamentary galaxy distribution being particularly sensitive to changes in the baryon density. Conversely, the least significant response emanates from the cosmic curvature Ωk, with only the vicinity of the dense galaxy clusters reacting to changes in the geometry of the Universe.

In the text

Fig. 3. Fisher information map for the same spherical slice as in Fig. 2. The observed galaxy distribution from the 2M++ catalog, lying in the corresponding spherical shell centered around the observer, is represented by the red dots. We find that the regions in the vicinity of the massive cosmic structures, as traced by the galaxy distribution, are the most informative according to the Fisher information map. These regions correspond to the regime of gravitational infall of the galaxy clusters.

In the text

Fig. 4. Fisher information map of the Coma cluster, with its corresponding mass density overlaid as a contour, for the central slice of thickness ∼5 h−1 Mpc through a 3D region centered on the cluster that extends over 40 h−1 Mpc. According to the Fisher information map, the central region containing the core of the cluster encodes fairly limited information gain, whilst the peripheral regions of the filaments and cluster core, where mass accretion occurs via gravitational infall, constitute the greatest proportion of information gain.

In the text

Fig. 5. Flowchart of the targeted search approach. The core idea is to use the existing knowledge as quantified by the posterior p(θ, ϕ|d) in order to calculate the Fisher map marginalized over this posterior, ⟨ℐ(θ|d0)⟩(ϕ|θ, d0), while keeping the observables of interest. This allows this Fisher map to provide us with regions of the sky with highest information gain potential for acquiring new data that are optimal for testing our model predictions, ⟨ℐ(θ|d0)⟩(ϕ|θ, d0): d0 → d. This data will in turn be used to update our knowledge about the model through updating the posterior and the procedure would repeat until the information content is fully depleted.

In the text