EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 709 (May 2026) A&A, 709 (2026) A110 Full HTML
Open Access
Issue
A&A
Volume 709, May 2026
Article Number A110
Number of page(s) 16
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202555817
Published online 12 May 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1. Introduction

The primordial density fluctuations provide the initial conditions for structure formation and encode statistical signatures of the early Universe. These early inhomogeneities, seeded during the initial moments after the Big Bang, evolved under gravity to form the vast cosmic-web structure we observe today. Reconstructing the initial conditions from present-day observations of luminous tracers – such as galaxies, quasars, Lyman-α forests, intensity maps, or from the 21 cm line – offers a powerful probe of fundamental physics, including inflation, dark matter, dark energy, and gravity itself. The scientific community is making significant efforts to map the matter distribution in the Universe, as exemplified by current large-scale spectroscopic surveys such as the Dark Energy Spectroscopic Instrument (DESI) (Levi et al. 2013; DESI Collaboration 2025), Euclid (Amendola et al. 2016; Euclid Collaboration: Aussel et al. 2025); and future ones such as the Prime Focus Spectrograph (PFS) (Takada et al. 2014), MUltiplexed Survey Telescope (MUST) (Zhao et al. 2024), and the Nancy Grace Roman Space Telescope (Roman) (Wang et al. 2022).

The first attempts to recover the primordial density fluctuations from the galaxy distribution were based on inverse reconstruction methods (e.g. Bertschinger & Dekel 1989; Nusser & Dekel 1992; Monaco & Efstathiou 1999). These approaches proved to be highly effective for sharpening features such as baryon acoustic oscillations (BAOs; Eisenstein et al. 2007), and were even proposed as tools to constrain primordial non-Gaussianities (Shirasaki et al. 2021). However, the accuracy of these inverse techniques is fundamentally limited by the non-invertibility of the gravitational evolution in the non-linear regime in the absence of full phase-space information, particularly due to shell-crossing. Once multiple streams of matter overlap, information about the initial conditions is irreversibly lost, as the full phase-space dynamics – crucial for a unique reconstruction – are no longer accessible through the evolved tracer distribution alone.

To overcome the limitations of inverse methods, forward modelling was introduced as a principled alternative. Instead of inverting the non-linear structure formation process, forward approaches start from initial conditions and evolve them through a physical model of gravitational dynamics to generate predictions for present-day observables (Kitaura 2013; Jasche & Wandelt 2013; Wang et al. 2013). The first application of a forward-modelling reconstruction (including peculiar motions) to observational data was presented by Kitaura et al. (2012b). All these works demonstrated the viability of sampling high-dimensional posterior distributions consistent with both data and gravitational physics, laying the foundation for later advances in bias modelling and inference methods. Subsequent developments considered improved mask handling (Jasche et al. 2015; Lavaux & Jasche 2016; Ata et al. 2021), improved gravity solvers using particle-mesh codes (Wang et al. 2014; Jasche & Lavaux 2019; Horowitz et al. 2019b), and emulators (Doeser et al. 2024). Further extensions incorporated light-cone evolution (Kitaura et al. 2021; Lavaux et al. 2019), primordial non-Gaussianities (Andrews et al. 2023), effective field theory (EFT) approaches (Stadler et al. 2023), and hydrodynamics (Horowitz & Lukic 2025). These techniques have also been applied to alternative tracers, such as the Lyman-α forest (e.g. Kitaura et al. 2012c; Horowitz et al. 2019a; Porqueres et al. 2019; Horowitz et al. 2021, 2022). A benchmark of computational efficiency in this context is presented in Bayer et al. (2023).

Bayesian statistics offers the ideal framework for field-level inference, enabling clear specification of prior knowledge and observational uncertainties (Kitaura & Enßlin 2008). In this setting, reconstruction becomes a matter of sampling from the posterior distribution of initial conditions, informed by both the underlying physics and the observed distribution of tracers.

This approach, however, faces a significant computational barrier. High-fidelity dark-matter simulations capable of resolving the halos that host galaxies require substantial computational resources – often hundreds of thousands to millions of CPU hours – rendering brute-force posterior sampling impractical (see state-of-the-art N-body simulations, e.g. Garrison et al. 2018; Chuang et al. 2019; Ishiyama et al. 2021).

To mitigate this challenge, Lagrangian perturbation theory (LPT) offers a computationally efficient alternative to full gravity solvers (Bernardeau et al. 2002), such as N-body simulations (Angulo & Hahn 2022). While being approximate, LPT accurately captures the gravitational evolution of matter on large, quasi-linear scales and serves as an effective forward model. Additional advancements have further extended and regularised these approximations, enabling reliable predictions down to megaparsec – or even sub-megaparsec – scales (Kitaura & Heß 2013; Kitaura et al. 2024).

This strategy is particularly effective when the forward model is constrained to its domain of validity and coupled to the observed tracers through an effective field-level bias prescription. The core challenge then lies in constructing a robust and physically motivated bias model. A critical first step is to ensure that large-scale bias is accurately captured, which necessitates the inclusion of non-linear and non-local contributions up to the third order (McDonald & Roy 2009).

A longstanding limitation of perturbation-theory-based bias models is that truncation at a fixed order often introduces unphysical behaviour (McDonald & Roy 2009; Schmittfull et al. 2019; Werner & Porciani 2020). Specifically, these models can yield non-positive densities and exhibit oscillatory or noisy behaviour when extended into the highly non-linear regime. However, recognising that non-local-bias terms at different orders correspond to distinct morphological features of the cosmic web opens the door to a more stable and flexible framework (Kitaura et al. 2022). For each morphological feature, we can then assume a local non-linear-bias model tailored to the specific properties of that structure. In this way, the classification of the cosmic web into distinct patterns effectively acts as a diagonalising operation, isolating the dominant local contributions to the tracer-matter relationship within each structural environment. This enables the systematic inclusion of higher order terms while preserving physical consistency and numerical stability down to small scales, making the framework well-suited for precision field-level inference. This ansatz provides a general and flexible framework applicable to a wide range of tracer populations and has demonstrated unprecedented accuracy in field-level bias modelling of point-like tracers such as dark-matter halos (Balaguera-Antolínez et al. 2023; Coloma-Nadal et al. 2024) and the baryonic components resolved in cosmological hydrodynamical simulations – including ionised gas, neutral hydrogen, and the Lyman-α forest (Sinigaglia et al. 2021, 2022, 2024b,a).

In this paper, we introduce BRIDGE1, a differentiable, GPU-accelerated framework for field-level inference written in JAX (Bradbury et al. 2018). In Sect. 2, we detail the technical implementation of the structure formation and bias models. Section 3 presents a series of computational tests demonstrating that the code successfully recovers both the primordial density field and the bias parameters of a long-range non-local-bias model at resolutions of 5 and 10 h−1 Mpc. Finally, we describe how we applied BRIDGE using a combined long- and short-range bias model, showing that the framework can robustly recover the primordial density field even in the presence of highly complex tracer bias.

2. Methods

In this section we present the technical details of the BRIDGE code. A flowchart showing the architecture of BRIDGE is presented in Fig. 1.

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Schematic overview of the BRIDGE pipeline. The framework combines a scalable, differentiable structure-formation model with a flexible field-level bias prescription, all implemented in a GPU-accelerated JAX environment. This enables efficient Bayesian inference of the primordial density field from tracer observations, with support for complex, non-local-bias models and multi-resolution analysis. The modular design allows for the seamless integration of physical models while maintaining end-to-end differentiability and high computational performance.

2.1. Field-level inference

We first consider a galaxy distribution, mapped onto a comoving volume by converting redshifts to distances, as the observed tracer field. With n ∗ ∈ N N 3 , Mathematical equation: $ \boldsymbol{n}^*\in \mathbb{N} ^{N^{3}}, $ we denote the observed tracer intensity on a cubic, comoving grid consisting of N3 cells. A differentiable forward model, ℳ, maps a white-noise realisation, ν, which defines the initial density field through a linear power spectrum, together with tracer bias parameters b and cosmological parameters Ω to the expected tracer field n ¯ Mathematical equation: $ \bar{\boldsymbol{n}} $:

M : ( ν , b , Ω ) ↦ n ¯ . Mathematical equation: $$ \begin{aligned} \mathcal{M} \!:\!(\boldsymbol{\nu },\boldsymbol{b}, \boldsymbol{\Omega })\mapsto {\bar{\boldsymbol{n}}} \ . \end{aligned} $$(1)

The forward model, ℳ, is constructed as a composition of three sequential maps, ℳ = ℳb ° ℳΨ ° ℳδ, where ℳδ is a linear transformation that maps a white-noise realisation, ν, to the initial overdensity field, δ0, at high redshift using the linear matter power spectrum (see Appendix A for details on the input field parametrisation); ℳΨ is the gravity solver that evolves the initial field, δ0, into the non-linear dark-matter over-density field, δ, at a specified redshift using an efficient gravity solver; and ℳb maps the evolved matter field, δ, together with the bias parameters, b, to the mean tracer field, n ¯ Mathematical equation: $ \bar{\boldsymbol{n}} $, employing a cosmic-web-dependent, non-local-bias model.

For each cell, we assume the data arise from an independent discrete distribution with mean n ¯ i Mathematical equation: $ \bar{n}_i $, which may include additional parameters, p, modelling survey completeness or super-Poissonian dispersion. The likelihood, which incorporates the forward model, reads

L ( ν , b , p , Ω ) ≡ P ( n ∗ | ν , b , p , Ω ) , Mathematical equation: $$ \begin{aligned} \mathcal{L} (\boldsymbol{\nu }, \boldsymbol{b}, \boldsymbol{p},\boldsymbol{\Omega }) \equiv \mathcal{P} (\boldsymbol{n}^* | \boldsymbol{\nu }, \boldsymbol{b}, \boldsymbol{p},\boldsymbol{\Omega }) \, , \end{aligned} $$(2)

where 𝒫 denotes the probability of observing n* given the model parameters.

Performing Bayesian inference at the field level amounts to sampling from the posterior distribution of the latent initial conditions and model parameters, given the observed tracer field n*. This posterior is given by

P ( ν , b , p , Ω ∣ n ∗ ) ∝ P ( ν ) P ( b , p , Ω ) L ( ν , b , p , Ω ) , Mathematical equation: $$ \begin{aligned} \mathcal{P} (\boldsymbol{\nu }, \boldsymbol{b}, \boldsymbol{p}, \boldsymbol{\Omega } \mid \boldsymbol{n}^*) \propto \mathcal{P} (\boldsymbol{\nu })\, \mathcal{P} (\boldsymbol{b}, \boldsymbol{p}, \boldsymbol{\Omega })\, \mathcal{L} (\boldsymbol{\nu }, \boldsymbol{b}, \boldsymbol{p}, \boldsymbol{\Omega }) \,, \end{aligned} $$(3)

where 𝒫(ν) encodes prior information on the initial conditions, and 𝒫(b, p, Ω) represents prior knowledge on the bias and cosmological parameters. The forward model was constructed such that the latent field, ν, is whitened, and therefore 𝒫 (ν) is an N3-dimensional multi-variate Gaussian with zero mean and identity covariance. The bias and additional likelihood model parameters, b and p, which in our case were strictly positive, were assigned independent log-normal priors. The probability density function of the log-normal distribution for the bias parameters implemented in the BRIDGE code is given by

f ( b ; μ , σ ) = 1 b σ 2 π exp [ − ( ln b − μ ) 2 2 σ 2 ] , b > 0 , Mathematical equation: $$ \begin{aligned} { f(b;\, \mu , \sigma ) = \frac{1}{b \sigma \sqrt{2\pi }} \exp \left[ -\frac{(\ln b - \mu )^2}{2\sigma ^2} \right] \, , \quad b > 0 }\, , \end{aligned} $$(4)

where μ and σ are fixed hyperparameters that define the mean and standard deviation of ln b. Throughout the remainder of this work, we fixed the cosmological parameters, Ω, to their fiducial values and omitted their explicit dependence for notational simplicity.

2.2. Posterior sampling

We exploited the differentiability of the forward model to efficiently sample the high-dimensional posterior using gradient-based Hamiltonian Monte Carlo (HMC) sampling (Duane et al. 1987). HMC was first applied to Bayesian large-scale-structure inference by Jasche & Kitaura (2010) and subsequently to observational data by Jasche et al. (2010).

We implemented within the BRIDGE code various posterior sampling methods, such as second- and fourth-order HMC sampling and the Blanes et al. (2014) HMC variant. The numerical tests shown in this work were finally performed with the publicly available No-U-Turn Sampler (NUTS) in NumPyro (Hoffman & Gelman 2011; Phan et al. 2019). During the adaptation (burn-in) phase, the algorithm tunes the step size, trajectory length, and the mass matrix, which is kept diagonal due to memory constraints. During sampling, we fixed the adapted mass matrix and step size and integrated the trajectories with a NUTS tree depth equal to 10.

To accelerate convergence, each chain was initialised from a deterministic state obtained by Wiener-filtering the observed tracer density contrast (Zaroubi et al. 1995), δ tr ∗ ≡ n ∗ N 3 / N tr − 1 Mathematical equation: $ \boldsymbol{\delta }^*_\text{tr} \equiv \boldsymbol{n}^* N^3 / N_\text{tr} - 1 $, back to an approximate initial density field, δ 0 ( 0 ) Mathematical equation: $ \boldsymbol{\delta }^{(0)}_{0} $, which was transformed into white noise, ν(0); Ntr is the total number of tracers. In this first step, we assumed the linear power spectrum to be known and ignored peculiar motions.

2.3. Gravity solver

We modelled large-scale-structure (LSS) formation with Augmented Lagrangian perturbation theory (ALPT; Kitaura & Heß 2013), which merges second-order LPT (2LPT; Zel’dovich 1970; Buchert 1994; Bouchet et al. 1995; Catelan 1995) for long-range tidal displacements and a spherical-collapse (SC) solution (Bernardeau 1994; Neyrinck 2016) for short-range dynamics. This was achieved by decomposing the particle displacement field as implemented in the BRIDGE code:

Ψ = Ψ L + Ψ S , Mathematical equation: $$ \begin{aligned} {\boldsymbol{\Psi } = {\boldsymbol{\Psi }}_{\rm L}+{\boldsymbol{\Psi }}_{\rm S}}\,, \end{aligned} $$(5)

with ΨL = 𝒦(q, rs)°Ψ2LPT,ΨS = (1 − 𝒦)°ΨSC, and 𝒦 a smoothing kernel of fixed radius, rs = 4 h−1 Mpc, following Kitaura & Heß (2013). In future work, we plan to sample this parameter within the Bayesian framework and consider nLPT for the long-range interaction.

This hybrid formulation retains analytic dependence on cosmological parameters while suppressing shell crossing. Coupled with empirical bias schemes such as the Bias Assignment Method (BAM). (Balaguera-Antolínez et al. 2018, 2020), ALPT reproduces N-body halo statistics down to k ≲ 0.4 h Mpc−1, the scale at which most cosmological information resides (Hahn & Villaescusa-Navarro 2021). An even higher accuracy was achieved with the HICOBIAN bias model considered in this work (see Sect. 2.7).

We note that ALPT is substantially more efficient in terms of both computational speed and memory usage compared to fast particle mesh (PM) methods (e.g. Tassev et al. 2013; Feng et al. 2016; Klypin & Prada 2018), as shown in the comparison work by Blot et al. (2019), making it especially advantageous for scalable Bayesian inference frameworks. Moreover, ALPT can be extended via iteration (eALPT) to reach sub-megaparsec accuracy (Kitaura et al. 2024) by indirectly emulating the effect of a viscous stress tensor in the equations of motion, thereby capturing aspects of non-linear dynamics such as vorticity. Recent GPU-based PM methods have shown noteworthy advancements (e.g. Modi et al. 2021; Li et al. 2024). Alternative promising approaches to approximate structure formation are offered by emulators trained on large ensembles of simulations (Kodi Ramanah et al. 2020; Conceição et al. 2024). However, for the time being, we opt to rely on analytic solutions with explicit cosmological dependence.

2.4. Bias model framework

When employing an approximate gravity solver, the evolved density field should not be directly trusted at the particle level. Instead, it is only reliable on a coarse mesh with a resolution of a few megaparsecs, where the approximation remains accurate. As a result, the bias relation between matter and tracers is no longer deterministic – as in approaches that apply a halo finder directly to particle distributions – it must be modelled as an effective field-level bias (Schmittfull et al. 2019). In this framework, the luminous tracers of LSS exhibit a complex, non-linear, non-local, and stochastic relationship with the underlying matter field (McDonald & Roy 2009; Desjacques et al. 2018). Accurately capturing this relationship is essential for unbiased inference and requires flexible and physically motivated modelling approaches that operate at the field level.

We took advantage of a positive, non-local-bias modelling framework that employs local-bias expansions within distinct morphological regions of the cosmic web (Kitaura et al. 2022; Coloma-Nadal et al. 2024; Sinigaglia et al. 2024b).

2.5. Deterministic local-bias model

The bias relation between tracers and the underlying matter field has long been known to be non-linear. One of the earliest approaches involved a local perturbative expansion of the tracer overdensity, introduced by Fry & Gaztanaga (1993):

δ tr = ∑ m = 0 ∞ b m m ! δ m , Mathematical equation: $$ \begin{aligned} \delta _{\rm tr} = \sum _{m = 0}^\infty \frac{b_m}{m!} \delta ^m\,, \end{aligned} $$(6)

where δtr denotes the tracer overdensity, δ the matter overdensity, and bm the bias coefficients. While conceptually straightforward, this expansion suffers from the drawback that it may yield unphysical negative tracer densities when truncated at any order.

To address this issue, an alternative logarithmic formulation was proposed around the same time by Cen & Ostriker (1993):

log ( 1 + δ tr ) = ∑ m = 0 ∞ c m m ! [ log ( 1 + δ ) ] m , Mathematical equation: $$ \begin{aligned} \log (1 + \delta _{\rm tr}) = \sum _{m = 0}^\infty \frac{c_m}{m!} \left[\log (1 + \delta )\right]^m\,, \end{aligned} $$(7)

which ensures positivity of the tracer field by construction. To a linear order, this expansion corresponds to a power-law bias model, which has been particularly effective in modelling the Lyman-α forest. In this context, it leads to the fluctuating Gunn–Peterson approximation, where the optical depth, τ, is related to the matter density field by τ ∝ (1 + δ)α (Bi & Davidsen 1997). An analogous relationship was found for the functional dependency between the molecular gas temperature and dark-matter density (Hui & Gnedin 1997).

Interestingly, a power-law bias model in Eulerian space can be related to a linear bias in Lagrangian space – i.e. in the coordinate system of the initial conditions – according to the continuity equation, which leads to a log-normal evolution of the density field (Coles & Jones 1991). In this picture, a linear Lagrangian bias naturally evolves into a power-law bias at later times.

In the context of galaxy bias, the polynomial perturbative expansion has historically been preferred. A foundational physical interpretation of bias was introduced even earlier by Kaiser (1984) through the excursion set model, where galaxies preferentially form in high-density peaks of the matter distribution. Building on this idea, Kitaura et al. (2014) proposed a hybrid model that combines the power-law bias of Cen & Ostriker (1993) with the threshold bias of Kaiser (1984) to model galaxy number counts as a function of the underlying dark-matter field. This framework was later refined by Neyrinck et al. (2014), the authors of which replaced the hard threshold with an exponential cutoff. The resulting model captures the average bias behaviour of halos with remarkable accuracy, particularly in low-density environments.

Although the power-law and threshold bias models are degenerate with respect to two-point statistics – both amplifying power on large scales – this degeneracy can be broken using three-point statistics, which provide a means to distinguish between different biasing mechanisms (see also Kitaura et al. 2015 for the definitions of three-point statistics calculations used in this work). The hybrid model has proven especially useful in this regard and was employed to generate thousands of mock catalogues of luminous red galaxies for the final BOSS data release (see the PATCHY mocks in Kitaura et al. 2016b).

To date, the suppression of tracer densities has primarily focused on low-density regions, following the picture introduced by Kaiser. However, another important effect in the context of halo and galaxy formation is halo exclusion (Baldauf et al. 2013), which accounts for the saturation of tracer objects in highly overdense environments. This mechanism introduces a natural suppression of tracer abundance in high-density regions, complementing the low-density suppression from the original threshold bias picture (Coloma-Nadal et al. 2024). A similar concept has also been applied in modelling the Lyman-α forest, where saturation effects in flux absorption become significant in dense regions (Sinigaglia et al. 2024b).

Motivated by these considerations, we modelled and implemented the expected tracer number density in the BRIDGE code using a power-law component modulated by both a low-pass and a high-pass filter:

n ¯ i = C ( 1 + δ i ) α exp [ − ( 1 + δ i ρ ) ϵ ] exp [ − ( 1 + δ i ρ ′ ) − ϵ ′ ] , Mathematical equation: $$ \begin{aligned} {\bar{n}_i = C \left(1 + \delta _i\right)^{\alpha } \exp \left[-\left(\frac{1 + \delta _i}{\rho } \right)^{ \epsilon } \right] \exp \left[- \left( \frac{1 + \delta _i}{\rho ^{\prime }} \right)^{ -\epsilon ^{\prime }} \right]} \ , \end{aligned} $$(8)

where C is a normalisation constant defined such that ∑ i n ¯ i = ∑ i n i ∗ ≡ N tr Mathematical equation: $ \sum _i \bar{n}_i = \sum _i n^*_i \equiv N_\text{tr} $, where Ntr is the total number of tracers in the data. All parameters are positive. We assume that the observed tracer count in each cell is a stochastic realisation drawn from a discrete probability distribution, P i tr Mathematical equation: $ \mathcal{P} ^{\text{ tr}}_i $, characterised by a mean expected value, n ¯ i Mathematical equation: $ \bar{n}_i $, provided by the forward model.

2.6. Stochastic bias model

The coarse resolution inherent to the field-level framework necessitates the introduction of a stochastic component to the bias model (Dekel & Lahav 1999; Sheth & Lemson 1999), which is not explicitly present in deterministic high-resolution N-body simulations. Nonetheless, the lack of full phase-space information after shell-crossing also introduces an uncertainty that requires a probabilistic treatment in the likelihood.

The most basic discretisation assumes a Poisson likelihood, where the variance equals the mean tracer count. The complete integration of this likelihood beyond the second moment within the Bayesian inference framework has been studied for a long time (Kitaura & Enßlin 2008; Kitaura et al. 2010). However, this assumption holds only in very limited regimes – typically for low tracer densities and in relatively homogeneous environments. As originally predicted by Peebles (1980), (anti-)correlations between tracers below the mesh resolution introduce deviations from pure Poisson statistics. These unresolved sub-grid correlations give rise to either sub- or super-Poissonian noise characteristics. To account for this, various probability distribution functions (PDFs) have been proposed in the literature (Saslaw & Hamilton 1984; Sheth 1995), most of which introduce additional degrees of freedom to model the over-dispersion commonly observed in galaxy and halo distributions relative to the Poisson expectation (Somerville et al. 2001; Casas-Miranda et al. 2002). Pellejero-Ibañez et al. (2020) demonstrated that the negative binomial distribution can be parametrised to reproduce a wide range of functional dependencies in the over-dispersed deviations from Poisson statistics.

We implemented the negative binomial distribution to model P i tr Mathematical equation: $ \mathcal{P} ^{\text{ tr}}_i $, thereby accounting for potential over-dispersion in the tracer counts. This distribution was first considered in the context of field-level bias in Kitaura et al. (2014) and later implemented in the Bayesian framework in Ata et al. (2015). Our implementation of the negative binomial distribution in this work relies on the gamma–Poisson mixture modelling, in which the Poisson rate, λi, is treated as a random variable drawn from a gamma distribution. The corresponding hierarchical formulation implemented in the BRIDGE code is given by

λ i ∼ Gamma ( β , n ¯ i / β ) Mathematical equation: $$ \begin{aligned} {\lambda _i \sim \text{ Gamma}(\beta , \bar{n}_i / \beta )}&\end{aligned} $$(9)

n i ∗ ∼ Poisson ( λ i ) , Mathematical equation: $$ \begin{aligned} {n^*_i \sim \text{ Poisson} (\lambda _i)}&\,, \end{aligned} $$(10)

where n ¯ i Mathematical equation: $ \bar{n}_i $ denotes the expected tracer count in cell i, and β is the over-dispersion parameter. In the limit β → ∞, the gamma distribution becomes sharply peaked, recovering the standard Poisson likelihood. In the gamma distribution, the first argument is the shape parameter and the second argument is the scale parameter. With this model, the (negative log) likelihood reads

− log L ( ν , b , β ) ≡ ∑ i = 1 N 3 − log P i tr ( n i ∗ | M ( ν , b ) , β ) , Mathematical equation: $$ \begin{aligned} - \log \mathcal{L} (\boldsymbol{\nu }, \boldsymbol{b}, \beta ) \equiv \sum _{i = 1}^{N^3} - \log \mathcal{P} ^{\text{ tr}}_i(n^*_i | \mathcal{M} (\boldsymbol{\nu }, \boldsymbol{b}), \beta ) \ , \end{aligned} $$(11)

with b = (α, ρ, ϵ, ρ′,ϵ′) from Eq. (8). While sub-Poissonian behaviour is indirectly captured by the high-density damping exponential factor, an explicit treatment can be implemented using the Conway–Maxwell–Poisson distribution (Daly & Gaunt 2016) or extensions of the binomial distribution (Vos-Ginés et al. 2024).

2.7. Non-local bias and cosmic-web dependence

In the previous sections, we introduced a local non-linear-bias model to describe the average bias relation and subsequently incorporated scatter to account for stochasticity. However, this stochastic component should be understood as a feature of the coarse-grained nature of the effective bias approach, without superseding deterministic contributions from non-local-bias effects.

Galaxy bias is shaped not only by the local matter density, but also by the surrounding large-scale tidal environment and the small-scale curvature of the density field. These non-local dependencies stem from the intrinsic nature of gravitational collapse and galaxy formation. A general framework for incorporating such effects into bias expansions was introduced by McDonald & Roy (2009). In another particularly relevant study, Chan et al. (2012) connected long-range tidal tensor invariants to higher order bias terms. Kitaura et al. (2022) later related the invariants of the Hessian of a field to the morphological classifications of the cosmic web.

These developments motivated the extension of bias models to explicitly include cosmic-web morphology (see Appendix B). Galaxies inhabiting different large-scale structures, including voids, sheets, filaments, and knots, exhibit systematically different biasing behaviours, as shown in simulations and observations (e.g. Nuza et al. 2014; Filho et al. 2015).

Several approaches exist for incorporating cosmic-web morphology. One method involves binning tidal field tensor invariants, which can result in a high-dimensional classification if multiple invariants and their dependence on local density are included (Kitaura et al. 2022). Alternatively, one can compress this information into four characteristic environments, sacrificing some modelling granularity for interpretability – this leads to the so-called Φ-web classification, based on the large-scale tidal field (Hahn et al. 2007).

To recover finer biasing distinctions, we introduced a hierarchical structure by embedding a second classification based on the Hessian of the density field (the δ-web) within each Φ-web region. This short-range non-local bias captures small-scale tidal effects rooted in tidal torque theory (Heavens & Peacock 1988). The result is the HICOBIAN model, which partitions the space into 16 regions (Coloma-Nadal et al. 2024). For each region, it applies a different local-bias relation, which is determined by the specific bias parameters of model Eq. (8). This approach gains physical insight with computational efficiency and outperforms fine binning of the tidal field in terms of both precision and interpretability, as both short- and long-range non-local biases are incorporated.

Given that gravitational evolution is inherently non-local, the shortcomings of ALPT relative to full gravity solvers can be effectively absorbed into non-local-bias terms. This explains the excellent performance of more recent studies employing ALPT (e.g. Balaguera-Antolínez et al. 2023; Coloma-Nadal et al. 2024). This insight provides a foundation for extending the bias framework to capture deviations arising from alternative gravity models (García-Farieta et al. 2024), where characteristic bias parameters and higher order statistics including Legendre multipole expansions would serve as indicators of such modified dynamics.

In the context of this work, it is important to note that most cosmic-web classifiers rely on hard, non-differentiable thresholds (e.g. eigenvalue sign changes), which are incompatible with gradient-based inference techniques such as HMC sampling. To address this, we adopted a fuzzy, differentiable classification scheme using sigmoid-based transitions. This enables smooth, physically motivated boundaries between web types and integrates seamlessly with our differentiable, field-level inference framework.

2.8. Differentiable fuzzy cosmic-web classification

In the top-down structure-formation scenario proposed by Zel’dovich (1970), cosmic structures form via anisotropic collapse governed by the tidal field tensor. The process unfolds along the tensor’s eigenvectors, with collapse first occurring along the direction associated with the largest eigenvalue – a phenomenon known as pancake formation. The term cosmic web was coined by Bond et al. (1996) to describe the large-scale structure of the Universe, characterised by the emergence of distinct morphological components such as filaments, walls, and nodes arising from anisotropic gravitational collapse. This interpretation was supported by shape parameters such as ellipticity and prolateness, derived from the eigenvalues of the deformation tensor and linked to the Zel’dovich approximation. A more systematic and operational classification of the cosmic web was later proposed by Hahn et al. (2007), the authors of which used the eigenvalues of the tidal (equal to the velocity shear within the Zel’dovich approximation) tensor to categorise regions of space into knots, filaments, sheets, and voids, depending on the number of collapsing directions. This method provided a quantitative tool to identify web environments in cosmological simulations.

While widely used, cosmic-web classification schemes remained largely phenomenological, lacking a firm quantitative foundation. A more rigorous interpretation began to emerge with the connection to the framework of non-local bias, which provided a physically motivated context for understanding how large-scale structure influences halo and galaxy formation – an approach we adopted in this work (see Appendix B).

We now consider a cosmic-web classification scheme based on the tidal field tensor, ∂j ∂k Φ, where Φ = ∇−2δ is the gravitational potential. Let λ i ( 1 ) < λ i ( 2 ) < λ i ( 3 ) Mathematical equation: $ \lambda _i^{(1)} < \lambda _i^{(2)} < \lambda _i^{(3)} $ be the eigenvalues of the tidal field tensor at the ith cell. The cosmic-web classification at the ith cell can be compactly expressed as

W i = ∑ n = 1 3 h ( λ i ( n ) − λ th ) , Mathematical equation: $$ \begin{aligned} W_i = \sum _{n = 1}^{3} h (\lambda _i^{(n)} - \lambda _\text{th}) \ , \end{aligned} $$(12)

where h is the Heaviside step function and λth is a free parameter. In this work, we set λth = 0.05 (see Forero-Romero et al. 2009 for the introduction of a threshold in the cosmic-web classification). For recent studies regarding the choice of λth, see Coloma-Nadal et al. (2024) and Olex et al. (2025). Wi counts the number of eigenvalues above the threshold; the four resulting values of 0, 1, 2, and 3 correspond to voids, sheets, filaments, and knots, respectively. This cosmic-web definition is problematic if introduced in a forward model that aims to be differentiable with respect to ν. This is because the cosmic-web classification mesh, W, depends on ν through the potential of δ = ℳΨ ° ℳδ (ν), but it does so in a manifestly non-differentiable way due to the presence of the Heaviside step function. As a consequence, the gradients ∂Wi/∂qj are ill-defined, which precludes the use of gradient-based sampling methods. To address this, we relied within the BRIDGE code on a fuzzy cosmic-web classification defined through the sigmoid weights, with

w i ( n ) = σ [ s ( λ i ( n ) − λ th ) ] , Mathematical equation: $$ \begin{aligned} {w^{(n)}_i = \sigma [s (\lambda ^{(n)}_i - \lambda _\text{th})]}\,, \end{aligned} $$(13)

where σ ( x ) = ( 1 + e − x ) − 1 Mathematical equation: $ \sigma (x) = (1+\mathrm{e} ^{-x})^{-1} $, and s is the steepness parameter controlling how sharply w i ( n ) Mathematical equation: $ w^{(n)}_i $ transitions around λth. Because σ is strictly increasing and the eigenvalues are sorted, the weights inherit the order w(1) > w(2) > w(3). We converted these cumulative exceedance probabilities into mutually exclusive ‘fuzzy’ memberships for knots (K), filaments (F), sheets (S), and voids (V):

p i ( K ) = w i ( 3 ) , Mathematical equation: $$ \begin{aligned} p^{(\text{ K})}_i&= w^{(3)}_i \ , \end{aligned} $$(14)

p i ( F ) = w i ( 2 ) − w i ( 3 ) , Mathematical equation: $$ \begin{aligned} p^{(\text{ F})}_i&= w^{(2)}_i - w^{(3)}_i\ , \end{aligned} $$(15)

p i ( S ) = w i ( 1 ) − w i ( 2 ) , Mathematical equation: $$ \begin{aligned} p^{(\text{ S})}_i&= w^{(1)}_i - w^{(2)}_i\ , \end{aligned} $$(16)

p i ( V ) = 1 − w i ( 1 ) , Mathematical equation: $$ \begin{aligned} p^{(\text{ V})}_i&= 1- w^{(1)}_i \ , \end{aligned} $$(17)

so that they represent the probabilities of having exactly three, two, one, or zero eigenvalues above the threshold, and, by construction, p i ( V ) + p i ( S ) + p i ( F ) + p i ( K ) = 1 Mathematical equation: $ p^{(\mathrm{V} )}_i+p^{(\mathrm{S} )}_i+p^{(\mathrm{F} )}_i+p^{(\mathrm{K} )}_i = 1 $. After introducing this smooth classification scheme, the forward model remained differentiable. Regarding the steepness parameter s, a larger value makes the fuzzy classification closer to the original classification, but it also risks causing numerical instabilities in the gradients. We therefore chose the largest s that preserves stable gradient computations. In Fig. 2, we show a comparison of the ‘sharp’ and smooth cosmic-web classification.

Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Example of fuzzy cosmic-web classification with N = 256 and ΔL = 1.7 h−1 Mpc. Top left: Evolved dark-matter density contrast field. Bottom left: ‘Hard’ Φ-web classification. Right panels: Fuzzy membership weights, pi(V), pi(S), pi(F), and pi(K), for voids, sheets, filaments, and knots, respectively, as defined in Eq. (17).

We now consider a cosmic-web-dependent bias model where each of the bias parameters from the local-bias model depend on the cosmic-web classification of the voxel. That is, given a bias parameter of b ∈ {α, ρ, ϵ, ρ⌂,ϵ⌂}, it now takes different values in each of the cosmic-web regions: b(K), b(F), b(S), b(V). To ensure differentiability, we used the fuzzy cosmic-web classification, and the bias parameter at voxel i takes the value

b i = b ( K ) p i ( K ) + b ( F ) p i ( F ) + b ( S ) p i ( S ) + b ( V ) p i ( V ) . Mathematical equation: $$ \begin{aligned} {b_i = b^{(\text{ K})} p^{(\text{ K})}_i + b^{(\text{ F})} p^{(\text{ F})}_i + b^{(\text{ S})} p^{(\text{ S})}_i + b^{(\text{ V})} p^{(\text{ V})}_i} \ . \end{aligned} $$(18)

2.9. Redshift-space-distortion modelling

Redshift-space distortions (RSDs) arise when the observed position of a tracer deviates from its Hubble-flow distance due to peculiar velocities, making them sensitive probes of the growth rate of structure formation (see Kaiser 1987; Hamilton 1998). Redshift-space distortions are governed by two distinct types of peculiar motions: coherent large-scale flows, which arise from gravitational infall towards overdense regions, and small-scale random motions within quasi-virialised structures. For a variety of ways of modelling coherent flows, we refer the reader to Kitaura et al. (2012a) and references therein (see also Wang et al. 2012). There are several strategies for incorporating RSDs within a Bayesian framework (we chose the third), which we list below.

  1. Iterative sampling of the redshift-to-real-space mapping. This approach iteratively reconstructs the real-space distribution from redshift-space observations using Gibbs sampling (Kitaura & Enßlin 2008; Kitaura et al. 2012c, 2016a). It is inspired by earlier iterative RSD correction methods (see Yahil et al. 1991; Monaco & Efstathiou 1999).

  2. Forward modelling of RSDs at the tracer level. Originally introduced by Kitaura et al. (2012b), this method evolves tracers from sampled initial conditions to redshift space, where their final positions are compared with observed galaxy distributions to constrain the primordial density field. Small-scale virial motions can be addressed in two ways: either by collapsing elliptical groups of galaxies – as proposed by Tegmark et al. (2004) – and modelling only coherent flows, as in the original study, or by retaining fingers-of-God features and modelling them with a stochastic velocity dispersion component, as done by Heß et al. (2013).

  3. Forward modelling of RSDs at the dark-matter-field level. In this work, we followed the approach of mapping the matter field directly to redshift space in a fully differentiable Bayesian inference framework. We set the observer at the center of the cubical volume and computed the RSD effects accordingly along lines of sight. See Bos et al. (2019) for a first implementation of such a method and calculation details.

We note that the modelling of RSDs for biased tracers involving non-linear transformations – such as cosmic voids and Lyman-α forests – requires the inclusion of velocity bias. For the Lyman-α forest, RSD affects the optical depth, τ, but the observable is the transmitted flux F = exp(−τ), introducing a non-trivial mapping. Accurate modelling in these cases must account for velocity-bias effects (McDonald et al. 2000; Seljak 2012; Sinigaglia et al. 2022, 2024b). Similarly, cosmic voids are identified from the galaxy distribution in redshift space, which effectively applies a non-linear and non-local transformation (Chuang et al. 2017). For the time being, we restrict the implementation tests within the BRIDGE code to unbiased peculiar motions (bv = 1); i.e. the dark-matter-density contrast in redshift space is given by

δ ( s ) = δ ( r + b v v · r a H r 2 r ) , Mathematical equation: $$ \begin{aligned} {\delta (\boldsymbol{s}) = \delta \left(\boldsymbol{r}+b_v\frac{\boldsymbol{v}\cdot \boldsymbol{r}}{aHr^2}\boldsymbol{r}\right)} \,, \end{aligned} $$(19)

with s and r being the redshift- and real-space coordinates, respectively; v the peculiar velocity; a the scale factor; and H the Hubble constant for a given Hubble distance.

3. Validation of field-level inference with fuzzy cosmic-web bias

To validate the BRIDGE code, we considered three numerical tests with different resolutions and bias models: TEST1, TEST2, and TEST3 (see Table 1). Results of TEST2 and TEST3 are presented in Appendix C. The numerical tests consist of the following steps.

  1. Generate a ground-truth number-count catalogue of biased tracers from an initial Gaussian density field with the same forward model employed during inference. All runs are performed to redshift z = 0 with a particular random seed and specific bias parameters as listed in Tables 2 and 3.

  2. Perform field-level inference of the posterior describing the white-noise and bias parameters for TEST1 and TEST2, applying the BRIDGE code.

  3. Perform quality assessment of the reconstructions by evaluating the recovery of both the initial and final conditions using two-point and three-point statistical analyses.

Regarding the hyperparameters of the bias parameter prior of Eq. (4), the power-law (α) and negative-binomial (β) hyperparameters are μα = 0.27, σα = 0.30, μβ = 2.12, and σβ = 0.20.

Table 1.

Setting the different numerical cases.

Table 2.

Values of bias parameters by cosmic-web region of the reference for the TEST1 and TEST2 cases (8 bias parameters).

Table 3.

Values of bias parameters by combined cosmic-web regions for the TEST3 case (64 bias parameters).

3.1. Motivation of the numerical tests

The numerical tests were chosen to demonstrate a number of scientific relevant cases. We list these below.

  • TEST1 and TEST2: The non-local-bias description used here is based on the Φ-web classification. Balaguera-Antolínez et al. (2018) demonstrated that this model can accurately reproduce key statistical properties of the halo distribution, including one-point PDFs, two-point power spectra at percent-level accuracy, and bispectra with reasonable agreement, generally within 15% (see also Balaguera-Antolínez et al. 2020, Pellejero-Ibañez et al. 2020, and Balaguera-Antolínez et al. 2023, who confirm these results for different halo resolutions). These studies confirm the Φ web as a physically meaningful and effective framework for modelling non-local bias within a reasonable degree of accuracy through non-parametric approaches. However, for field-level inference applications – such as the one developed here – a parametric bias model is required to enable efficient sampling and marginalisation within a Bayesian framework. The same non-local-bias treatment incorporating an equivalent non-linear-bias prescription to the one considered here (see Eq. 8) was employed by Sinigaglia et al. (2024b,a) to model the Lyman-α forest. Their work demonstrated excellent agreement with high-resolution hydrodynamical simulations, achieving power-spectrum accuracy within 5% up to k ∼ 1 h Mpc−1, along with consistent bispectrum predictions. In both cases, i.e. haloes and Lyman-α forests, the underlying structure-formation model was ALPT. For both the TEST1 and TEST2 cases, we jointly sampled the bias parameters starting with an initial guess randomly sampled from the prior distribution (see Eq. 4).

  • TEST1 focuses on a resolution of 10 h−1 Mpc, corresponding to a Nyquist frequency of k ≃ 0.3 h Mpc−1. This resolution surpasses the typical scale used in most cosmological analyses, which are commonly limited to k < 0.2 h Mpc−1 (see e.g. Ivanov et al. 2025). Also, the ideal resolution at which BAO reconstruction techniques are applied is slightly lower, i.e.15 h−1 Mpc (see e.g. Paillas et al. 2025).

  • TEST2 focuses on a resolution of 5 h−1 Mpc with a Nyquist frequency of k ≃ 0.6 h Mpc−1. This resolution has been reported to be high enough to produce accurate galaxy catalogues combined with a sub-grid model based on ALPT (see Forero Sánchez et al. 2024).

  • TEST3: This case study incorporates the full HICOBIAN non-local bias model, which has been shown to significantly enhance accuracy in both power spectra and bispectra, achieving full compatibility with N-body-based catalogues across a wide range of scales. While Coloma-Nadal et al. (2024) demonstrated that this model remains valid at mesh resolutions below 4 h−1 Mpc, in this study we restricted our analysis to a coarser resolution of 8 h−1 Mpc. This choice is motivated by the fact that the corresponding Nyquist frequency, k ≃ 0.4 h Mpc−1, marks the scale beyond which cosmological information is effectively saturated (Hahn & Villaescusa-Navarro 2021). For computational reasons, we kept the bias parameters fixed in this case, assuming they are known. This assumption is reasonable if the parameters can be extracted from high-fidelity reference catalogues, as was done for the PATCHY BOSS mocks, where the MULTIDARK reference catalogue was designed to reproduce the observed galaxy distribution with halo abundance matching (Rodríguez-Torres et al. 2016). Even more sophisticated reference catalogues are now being developed to reproduce the clustering statistics of the observed DESI catalogues (DESI Collaboration 2025), from which the parameters of the HICOBIAN model can be accurately extracted (Favole et al. 2025).

We considered ρ′ = 0 and ϵ′ = 1 for simplicity, interpreting this by continuity for δ > −1, so the high-density damping factor, exp [ − ( 1 + δ i ρ ′ ) − ϵ ′ ] Mathematical equation: $ \exp \left[- \left( \frac{1 + \delta_i}{\rho\prime} \right)^{ -\epsilon\prime} \right] $, was set to one in all three test cases.

3.2. Numerical results of TEST1

The first numerical test was performed on a grid with 1283 voxels within a comoving volume of 1280 h−1 Mpc per side, corresponding to a spatial resolution of 10 h−1 Mpc. We used a Φ-web-dependent bias modelled with a power-law and negative-binomial likelihood for each region, for a total of eight bias parameters. Using this model and a fixed random seed, we constructed a mock catalogue of object number counts, n*, in redshift space, which we define as the ground truth. The priors of the power-law (α) and negative-binomial (β) bias parameters were set to be log-normally distributed (see Eq 4) with fixed hyperparameters: μα = 0.27, σα = 0.30, μβ = 2.12, σβ = 0.20.

The HMC chain was initialised with an state (ν(0), b(0)), with ν(0) being the white-noise representation of the density contrast coming from Wiener-filtering, n*, and b(0) being bias parameters randomly sampled from their priors. During burn-in, chain convergence was achieved after approximately 250 samples. The adapted step size was ∼ 10−3, and the NUTS trajectory length adaptation saturated at 1024 leapfrog steps (corresponding to a maximum tree depth of 10).

After convergence, we assessed the correlation length in the white-noise samples, based on their autocorrelation at lag m:

ξ i ( m ) = 1 σ i 2 ( M − m ) ∑ j = 0 M − m ( ν i ( j ) − ν ¯ i ) ( ν i ( j + m ) − ν ¯ i ) , Mathematical equation: $$ \begin{aligned} \xi _i(m) = \frac{1}{\sigma ^2_i {\left(M-m\right)}} \sum _{j = 0}^{M-m} \left(\nu _i^{(j)} - \bar{\nu }_i\right) \left(\nu _i^{(j+m)} - \bar{\nu }_i\right) \ , \end{aligned} $$(20)

where νi(j) is the ith voxel of the jth sample, M the total number of samples of a given chain run, ν ¯ i Mathematical equation: $ \bar{\nu }_i $ the within-chain mean, and σ i 2 Mathematical equation: $ \sigma _i^2 $ the within-chain variance. For a chain of M = 1000 samples, we computed ξi(m) for lags m = 0 to 100 across 104 randomly selected voxels. Figure 3 shows the mean auto-correlation over this subset.

Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

TEST1: Mean auto-correlation ξ(n) as a function of lag, n, for a chain of length of M ≈ 1000, evaluated over 104 randomly selected voxels. The red curve shows the ensemble average, ⟨ξ(n)⟩,  across parameters, with the shaded band indicating the ±1σ dispersion. The horizontal dashed line marks the threshold, ξ = 0.1, used to assess mixing efficiency.

We defined the correlation length as the smallest lag, n, such that ⟨ ξ (n) ⟩ ≤ 0.1, yielding an estimate of ∼40 samples. Given this auto-correlation, we drew 15 000 post-burn-in states, providing about 500 effectively independent posterior samples for the subsequent statistical analysis.

We assessed the computing times of single-gradient calculations across different mesh sizes and bias configurations finding that they are below one second (see Fig. 4).

Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Computing times for single-gradient evaluations of the forward model employing ALPT evolution and the HICOBIAN bias model for different mesh sizes, run on a single NVIDIA A100-SXM4 GPU equipped with 40 GB of on-board HBM2 memory. For reference, the full burn-in and sampling procedure for the 1283 case in this study required approximately 100 GPU hours.

In addition, we performed the Gelman–Rubin convergence diagnostic, which confirms that the drawn samples have converged to the same target distribution (see Fig. 5).

Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Gelman–Rubin convergence diagnostic computed for three independently initialised chains of 250 samples after convergence. This statistic is calculated for all cells in Fourier space, and we show their distribution with scale, k. Values close to one indicate that the drawn samples have converged to the same distribution.

As a first step, we conducted a visual assessment of the true and reconstructed maps of the initial and final density fields (see Fig. 6). The reconstructed samples exhibit a high degree of visual similarity to the ground truth for both the initial and final conditions. On a more quantitative level, we find that the standard deviations of the reconstructed fields display structured patterns at amplitudes approximately an order of magnitude lower than the true density field at the initial conditions. In contrast, the final conditions exhibit higher variance, reflecting the increased non-linearity and complexity of the evolved structures.

Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

TEST1: Spatial summary of 500 independent reconstructions (N = 128). Slices of ΔL = 10 h−1 Mpc (one voxel width). Top row: Reconstructed initial density contrast. Bottom row: Reconstructed tracer field. Columns from left to right: (1) Pixel-wise standard deviation across the 500 independent reconstructions; (2) mean reconstructed field; (3) one representative reconstruction sample; and (4) true reference field.

To quantitatively assess the quality of the reconstruction, we computed both two- and three-point statistics in Fourier space for the initial (Fig. 7) and final conditions (Fig. 8). Specifically, we evaluated the monopole and quadrupole moments of the power spectrum, along with the reduced bispectrum for a configuration particularly sensitive to non-linear effects (k1 = 0.1 and k2 = 0.2 h Mpc−1). In addition, we computed the propagator, defined as the cross-correlation between the mock tracer density contrast and the initial density reference. The motivation for investigating the quadrupole and bispectrum of a the reconstructed Gaussian field was to assess how well the samples recover deviations introduced by cosmic variance in these statistics.

Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

TEST1: Comparison of summary statistics of the initial density contrast field for the mock reference (dashed red line) and 500 independent reconstructions (solid blue line). Columns from left to right: (1) Monopole power spectrum; (2) quadrupole; (3) reduced bispectrum as a function of triangle opening angle, θ; and (4) cross-correlation coefficient C(k) between the reference and reconstruction (solid blue line) and the propagator (dashed red line), defined as the cross-correlation between the mock tracer density contrast and the initial density reference. In columns (1)–(3), the lower sub-panels display the ratios with respect to the reference. Shaded bands denote the 1σ intervals from the posterior samples. Vertical dashed lines indicate the isotropic Nyquist frequency, k nyq / 3 Mathematical equation: $ k_{\mathrm{nyq}}/\sqrt{3} $.

In Figs. 7 and 8, it is remarkable to observe that the monopoles, quadrupoles, and bispectra are accurately recovered not only in the final conditions, but also in the initial conditions, even reproducing the specific features inherent to the particular cosmic variance realisation. Although local-bias models fail to reproduce three-point statistics without short- and long-range non-local terms (Coloma-Nadal et al. 2024), we show that only neglecting long-range non-local bias also leads to significant differences in the bispectrum (see the third panel in Fig. 8).

Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

TEST1: Comparison of summary statistics of the mock reference and 500 independent reconstructions evolved through the forward model. Panel definitions are identical to those in Fig. 7, except that in the rightmost panel the propagator is replaced by the maximal tracer field cross-correlation and that we show the bispectrum of the dark-matter field and a tracer field using a local bias (without cosmic web) with bias parameters set to those of knots.

The propagator ensures a significant gain in information, demonstrating the effectiveness of the reconstruction. To evaluate how close the result is to an ideal scenario, we computed the optimal cross at final conditions by generating a synthetic-tracer number-count field from the evolved ground-truth dark-matter field using the exact same bias parameters and stochastic prescription as in the model, but a different seed. The cross-correlation between this ideal tracer field and the ground truth serves as a benchmark, represented by the dashed-dotted red line in the figure. This allowed us to assess whether the reconstruction approaches the theoretical limit given the assumed bias model. As we see in Fig. 8, the red and blue lines of the right panel are on top of each other, indicating that the optimal cross-correlation limited by shot noise has been reached.

In Fig. 9, we show a corner plot of the posterior distributions of the eight parameters of the bias model compared to their values in the reference catalogue. The power-law parameters, α, were recovered with sub-percent accuracy, while the negative binomial β-parameter was recovered with precision of a few percent. The absence of significant correlations between α bias parameters in distinct cosmic-web regions supports the validity and stability of our inference framework.

Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

TEST1: Posterior distributions for the 500 independent samples of the α and β bias parameters in each cosmic-web region of the Φ web. Red dots and lines show the values of the parameters for the ground-truth catalogue. The posteriors are smoothed with a 2D Gaussian filter with standard deviation equal to 0.75 bin widths. TEST2 yields equivalent results.

Overall, the consistency of all diagnostics – maps, summary statistics, and bias-parameter posteriors – demonstrates that the complete field-level inference pipeline reliably recovers both the underlying density field and the cosmic-web-dependent bias parameters at the targeted resolution. The numerical studies TEST2 and TEST3 further confirm the validation of the BRIDGE code (details can be found in Appendix C). However, it is remarkable that we did not recover the original bias parameters when we allowed them to vary freely in case TEST3.

4. Conclusions

In this work we developed and validated a novel Bayesian field-level inference framework that incorporates a physically motivated, differentiable, and non-local bias model – HICOBIAN – based on the hierarchical cosmic-web. This model allows for smooth transitions between different cosmic environments via a fuzzy classification approach, making it well suited for machine-learning integration and GPU acceleration.

We demonstrated that, within a self-specified setting where inference is performed on data generated by the exact same forward model, our approach accurately reconstructs the initial density field, achieving high fidelity even in the presence of redshift-space distortions and complex, non-linear, non-local, and stochastic bias. The reconstruction quality was evaluated using two- and three-point statistics, including the power-spectrum monopole and quadrupole, the reduced bispectrum, and the propagator, showing excellent agreement with the ground truth and reaching the optimal limit set by shot noise.

We further tested the robustness of the model under various conditions, including varying spatial resolutions and complex bias scenarios involving four to 16 cosmic-web regions, each characterised by distinct bias parameters. These included power-law bias, exponential threshold bias, and deviations from Poisson statistics, totalling eight to 64 bias parameters. Across all settings, the method maintained its high accuracy, confirming the generality and adaptability of our approach. However, for the complex bias models, the original bias parameters were not recovered when left to vary freely. This hints at the complexity of bias modelling and the potential limitations of field-level inference approaches that ignore intra-halo dynamics. As demonstrated in recent work (e.g. Favole et al. 2025), a realistic halo-to-galaxy connection based on high-resolution N-body simulations requires a large number of parameters to capture processes such as halo exclusion, satellite dynamics, quenching, and intra-halo flows.

In future work, we plan to investigate the performance of the BRIDGE code with high-fidelity synthetic reference catalogues and directly apply it to redshift survey data across various tracers, including DESI-tracers, such as bright galaxies, luminous red galaxies, emission-line galaxies, quasars, and Lyman-α forests. This will require the incorporation of light-cone evolution and selection effects and a more detailed analysis of redshift-space distortions in the highly non-linear regime. Moreover, our results suggest that achieving full consistency may require a joint Bayesian inference framework that incorporates intra-halo dynamics informed by detailed N-body simulations to break the degeneracy in the bias parameters.

By implementing this framework in the GPU-accelerated BRIDGE code using JAX, we enable efficient and scalable inference on large datasets. This work thus marks a significant step towards a unified, interpretable, and fully differentiable framework for cosmological data analysis, with direct applicability to tracers of the large-scale structure.

Acknowledgments

The authors thank the Spanish Ministry of Science and Innovation for funding the Big Data of the Cosmic-Web project (PID2020-120612GB-I00/AEI/10.13039/501100011033, PI:FSK), under which this work was conceived and carried out. We are also grateful to the Instituto de Astrofísica de Canarias (IAC) for its continued support through the Cosmology with LSS Probes project (PI:FSK). The authors thank Jesús Sanz Serna and José María Coloma-Nadal for discussions. PR acknowledges support from the IAC Research Summer Grant 2024 for the project “Bayesian Inference of the LSS with ALPT and HICOBIAN bias using JAX”. This grant enabled his contribution to this work prior to the start of his master’s thesis on the same topic, as reflected in this manuscript. He also thanks Gustavo Yepes, the Universidad Autónoma de Madrid (UAM), and the Department of Theoretical Physics for their hospitality and for providing access to computing facilities during his visit – resources that were instrumental to the development of this work.

References

  1. Amendola, L., Appleby, S., Avgoustidis, A., et al. 2016, arXiv e-prints [arXiv:1606.00180] [Google Scholar]
  2. Andrews, A., Jasche, J., Lavaux, G., & Schmidt, F. 2023, MNRAS, 520, 5746 [Google Scholar]
  3. Angulo, R. E., & Hahn, O. 2022, LRCA, 8, 1 [Google Scholar]
  4. Ata, M., Kitaura, F.-S., & Müller, V. 2015, MNRAS, 446, 4250 [Google Scholar]
  5. Ata, M., Kitaura, F.-S., Lee, K.-G., et al. 2021, MNRAS, 500, 3194 [Google Scholar]
  6. Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., Zhao, C., & Abel, T. 2018, MNRAS, 483, L58 [Google Scholar]
  7. Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibáñez, M., et al. 2020, MNRAS, 491, 2565 [Google Scholar]
  8. Balaguera-Antolínez, A., Kitaura, F.-S., Alam, S., et al. 2023, A&A, 673, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Baldauf, T., Seljak, U., Smith, R. E., Hamaus, N., & Desjacques, V. 2013, Phys. Rev. D, 88, 083507 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bayer, A. E., Seljak, U., & Modi, C. 2023, arXiv e-prints [arXiv:2307.09504] [Google Scholar]
  11. Bernardeau, F. 1994, ApJ, 427, 51 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1 [Google Scholar]
  13. Bertschinger, E., & Dekel, A. 1989, ApJ, 336, L5 [NASA ADS] [CrossRef] [Google Scholar]
  14. Bi, H., & Davidsen, A. F. 1997, ApJ, 479, 523 [Google Scholar]
  15. Blanes, S., Casas, F., & Sanz-Serna, J. M. 2014, SIAM J. Sci. Comput., 36, A1556 [Google Scholar]
  16. Blot, L., Crocce, M., Sefusatti, E., et al. 2019, MNRAS, 485, 2806 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603 [NASA ADS] [CrossRef] [Google Scholar]
  18. Bos, E. G. P., Kitaura, F.-S., & van de Weygaert, R. 2019, MNRAS, 488, 2573 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575 [NASA ADS] [Google Scholar]
  20. Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, {JAX}: composable transformations of {P}ython+{N}um{P}y programs, http://github.com/jax-ml/jax [Google Scholar]
  21. Buchert, T. 1994, MNRAS, 267, 811 [NASA ADS] [Google Scholar]
  22. Casas-Miranda, R., Mo, H. J., Sheth, R. K., & Boerner, G. 2002, MNRAS, 333, 730 [NASA ADS] [CrossRef] [Google Scholar]
  23. Catelan, P. 1995, MNRAS, 276, 115 [NASA ADS] [Google Scholar]
  24. Cen, R., & Ostriker, J. P. 1993, ApJ, 417, 415 [NASA ADS] [CrossRef] [Google Scholar]
  25. Chan, K. C., Scoccimarro, R., & Sheth, R. K. 2012, Phys. Rev. D, 85, 083509 [CrossRef] [Google Scholar]
  26. Chuang, C.-H., Kitaura, F.-S., Liang, Y., et al. 2017, Phys. Rev. D, 95, 063528 [NASA ADS] [CrossRef] [Google Scholar]
  27. Chuang, C.-H., Yepes, G., Kitaura, F.-S., et al. 2019, MNRAS, 487, 48 [NASA ADS] [CrossRef] [Google Scholar]
  28. Coles, P., & Jones, B. 1991, MNRAS, 248, 1 [NASA ADS] [CrossRef] [Google Scholar]
  29. Coloma-Nadal, J. M., Kitaura, F.-S., García-Farieta, J. E., et al. 2024, JCAP, 2024, 083 [CrossRef] [Google Scholar]
  30. Conceição, M., Krone-Martins, A., da Silva, A., & Moliné, Á. 2024, A&A, 681, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Daly, F., & Gaunt, R. E. 2016, ALEA: Latin Am. J. Probability Math. Stat., 13, 635 [Google Scholar]
  32. Dekel, A., & Lahav, O. 1999, ApJ, 520, 24 [Google Scholar]
  33. DESI Collaboration (Abdul-Karim, M., et al.) 2025, arXiv e-prints [arXiv:2503.14745] [Google Scholar]
  34. Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1 [Google Scholar]
  35. Doeser, L., Jamieson, D., Stopyra, S., et al. 2024, MNRAS, 535, 1258 [NASA ADS] [CrossRef] [Google Scholar]
  36. Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D. 1987, Phys. Lett. B, 195, 216 [CrossRef] [Google Scholar]
  37. Eisenstein, D. J., Seo, H.-J., Sirko, E., & Spergel, D. N. 2007, ApJ, 664, 675 [NASA ADS] [CrossRef] [Google Scholar]
  38. Euclid Collaboration (Aussel, H., et al.) 2026, A&A, in press, https://doi.org/10.1051/0004-6361/202554610 [Google Scholar]
  39. Favole, G., Kitaura, F. S., Hadzhiyska, B., et al. 2025, ApJ, submitted [arXiv:2512.04362] [Google Scholar]
  40. Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273 [NASA ADS] [CrossRef] [Google Scholar]
  41. Filho, M. E., Sánchez Almeida, J., Muñoz-Tuñón, C., et al. 2015, ApJ, 802, 82 [NASA ADS] [CrossRef] [Google Scholar]
  42. Forero Sánchez, D., Kitaura, F. S., Sinigaglia, F., Coloma-Nadal, J. M., & Kneib, J. P. 2024, JCAP, 2024, 001 [Google Scholar]
  43. Forero-Romero, J. E., Hoffman, Y., Gottlöber, S., Klypin, A., & Yepes, G. 2009, MNRAS, 396, 1815 [Google Scholar]
  44. Fry, J. N., & Gaztanaga, E. 1993, ApJ, 413, 447 [NASA ADS] [CrossRef] [Google Scholar]
  45. García-Farieta, J. E., Balaguera-Antolínez, A., & Kitaura, F.-S. 2024, A&A, 690, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Garrison, L. H., Eisenstein, D. J., Ferrer, D., et al. 2018, ApJS, 236, 43 [NASA ADS] [CrossRef] [Google Scholar]
  47. Hahn, C., & Villaescusa-Navarro, F. 2021, JCAP, 2021, 029 [Google Scholar]
  48. Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hamilton, A. J. S. 1998, Astrophys. Space Sci. Libr., 231, 185 [NASA ADS] [CrossRef] [Google Scholar]
  50. Heavens, A., & Peacock, J. 1988, MNRAS, 232, 339 [NASA ADS] [CrossRef] [Google Scholar]
  51. Heß, S., Kitaura, F.-S., & Gottlöber, S. 2013, MNRAS, 435, 2065 [CrossRef] [Google Scholar]
  52. Hoffman, M. D., & Gelman, A. 2011, arXiv e-prints [arXiv:1111.4246] [Google Scholar]
  53. Horowitz, B., & Lukic, Z. 2025, arXiv e-prints [arXiv:2502.02294] [Google Scholar]
  54. Horowitz, B., Lee, K.-G., White, M., Krolewski, A., & Ata, M. 2019a, ApJ, 887, 61 [NASA ADS] [CrossRef] [Google Scholar]
  55. Horowitz, B., Seljak, U., & Aslanyan, G. 2019b, JCAP, 2019, 035 [CrossRef] [Google Scholar]
  56. Horowitz, B., Zhang, B., Lee, K.-G., & Kooistra, R. 2021, ApJ, 906, 110 [Google Scholar]
  57. Horowitz, B., Dornfest, M., Lukić, Z., & Harrington, P. 2022, ApJ, 941, 42 [NASA ADS] [CrossRef] [Google Scholar]
  58. Hui, L., & Gnedin, N. Y. 1997, MNRAS, 292, 27 [Google Scholar]
  59. Ishiyama, T., Prada, F., Klypin, A. A., et al. 2021, MNRAS, 506, 4210 [NASA ADS] [CrossRef] [Google Scholar]
  60. Ivanov, M. M., Obuljen, A., Cuesta-Lazaro, C., & Toomey, M. W. 2025, Phys. Rev. D, 111, 063548 [Google Scholar]
  61. Jasche, J., & Kitaura, F.-S. 2010, MNRAS, 407, 29 [NASA ADS] [CrossRef] [Google Scholar]
  62. Jasche, J., & Lavaux, G. 2019, A&A, 625, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Jasche, J., & Wandelt, B. D. 2013, MNRAS, 432, 894 [Google Scholar]
  64. Jasche, J., Kitaura, F. S., Li, C., & Enßlin, T. A. 2010, MNRAS, 409, 355 [Google Scholar]
  65. Jasche, J., Leclercq, F., & Wandelt, B. 2015, JCAP, 2015, 036 [Google Scholar]
  66. Kaiser, N. 1984, ApJ, 284, L9 [NASA ADS] [CrossRef] [Google Scholar]
  67. Kaiser, N. 1987, MNRAS, 227, 1 [Google Scholar]
  68. Kitaura, F.-S. 2013, MNRAS, 429, L84 [NASA ADS] [CrossRef] [Google Scholar]
  69. Kitaura, F.-S., & Enßlin, T. A. 2008, MNRAS, 389, 497 [NASA ADS] [CrossRef] [Google Scholar]
  70. Kitaura, F.-S., & Heß, S. 2013, MNRAS, 435, L78 [Google Scholar]
  71. Kitaura, F.-S., Jasche, J., & Metcalf, R. B. 2010, MNRAS, 403, 589 [Google Scholar]
  72. Kitaura, F.-S., Angulo, R. E., Hoffman, Y., & Gottlöber, S. 2012a, MNRAS, 425, 2422 [Google Scholar]
  73. Kitaura, F.-S., Erdoǧdu, P., Nuza, S. E., et al. 2012b, MNRAS, 427, L35 [NASA ADS] [Google Scholar]
  74. Kitaura, F.-S., Gallerani, S., & Ferrara, A. 2012c, MNRAS, 420, 61 [Google Scholar]
  75. Kitaura, F.-S., Yepes, G., & Prada, F. 2014, MNRAS, 439, L21 [NASA ADS] [CrossRef] [Google Scholar]
  76. Kitaura, F.-S., Gil-Marín, H., Scóccola, C. G., et al. 2015, MNRAS, 450, 1836 [NASA ADS] [CrossRef] [Google Scholar]
  77. Kitaura, F.-S., Ata, M., Angulo, R. E., et al. 2016a, MNRAS, 457, L113 [Google Scholar]
  78. Kitaura, F.-S., Rodríguez-Torres, S., Chuang, C.-H., et al. 2016b, MNRAS, 456, 4156 [NASA ADS] [CrossRef] [Google Scholar]
  79. Kitaura, F.-S., Ata, M., Rodríguez-Torres, S. A., et al. 2021, MNRAS, 502, 3456 [NASA ADS] [CrossRef] [Google Scholar]
  80. Kitaura, F.-S., Balaguera-Antolínez, A., Sinigaglia, F., & Pellejero-Ibáñez, M. 2022, MNRAS, 512, 2245 [NASA ADS] [CrossRef] [Google Scholar]
  81. Kitaura, F.-S., Sinigaglia, F., Balaguera-Antolínez, A., & Favole, G. 2024, A&A, 683, A215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Klypin, A., & Prada, F. 2018, MNRAS, 478, 4602 [NASA ADS] [CrossRef] [Google Scholar]
  83. Kodi Ramanah, D., Charnock, T., Villaescusa-Navarro, F., & Wandelt, B. D. 2020, MNRAS, 495, 4227 [NASA ADS] [CrossRef] [Google Scholar]
  84. Lavaux, G., & Jasche, J. 2016, MNRAS, 455, 3169 [Google Scholar]
  85. Lavaux, G., Jasche, J., & Leclercq, F. 2019, Systematic-Free Inference of the Cosmic Matter Density Field from SDSS3-BOSS Data [Google Scholar]
  86. Levi, M., Bebek, C., Beers, T., et al. 2013, arXiv e-prints [arXiv:1308.0847] [Google Scholar]
  87. Li, Y., Modi, C., Jamieson, D., et al. 2024, ApJS, 270, 36 [Google Scholar]
  88. McDonald, P., & Roy, A. 2009, JCAP, 2009, 020 [CrossRef] [Google Scholar]
  89. McDonald, P., Miralda-Escudé, J., Rauch, M., et al. 2000, ApJ, 543, 1 [NASA ADS] [CrossRef] [Google Scholar]
  90. Modi, C., Lanusse, F., & Seljak, U. 2021, Astron. Comput., 37, 100505 [NASA ADS] [CrossRef] [Google Scholar]
  91. Monaco, P., & Efstathiou, G. 1999, MNRAS, 308, 763 [NASA ADS] [CrossRef] [Google Scholar]
  92. Neyrinck, M. C. 2016, MNRAS, 455, L11 [Google Scholar]
  93. Neyrinck, M. C., Aragón-Calvo, M. A., Jeong, D., & Wang, X. 2014, MNRAS, 441, 646 [Google Scholar]
  94. Nusser, A., & Dekel, A. 1992, ApJ, 391, 443 [Google Scholar]
  95. Nuza, S. E., Kitaura, F.-S., Heß, S., Libeskind, N. I., & Müller, V. 2014, MNRAS, 445, 988 [NASA ADS] [CrossRef] [Google Scholar]
  96. Olex, E., Hellwing, W. A., & Knebe, A. 2025, A&A, 696, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Paillas, E., Ding, Z., Chen, X., et al. 2025, JCAP, 2025, 142 [Google Scholar]
  98. Peebles, P. J. E. 1980, The large-scale structure of the universe (Princeton University Press) [Google Scholar]
  99. Pellejero-Ibañez, M., Balaguera-Antolínez, A., Kitaura, F.-S., et al. 2020, MNRAS, 493, 586 [CrossRef] [Google Scholar]
  100. Phan, D., Pradhan, N., & Jankowiak, M. 2019, arXiv e-prints [arXiv:1912.11554] [Google Scholar]
  101. Porqueres, N., Jasche, J., Lavaux, G., & Enßlin, T. 2019, A&A, 630, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Rodríguez-Torres, S. A., Chuang, C.-H., Prada, F., et al. 2016, MNRAS, 460, 1173 [CrossRef] [Google Scholar]
  103. Saslaw, W. C., & Hamilton, A. J. S. 1984, ApJ, 276, 13 [Google Scholar]
  104. Schmittfull, M., Simonović, M., Assassi, V., & Zaldarriaga, M. 2019, Phys. Rev. D, 100, 043514 [NASA ADS] [CrossRef] [Google Scholar]
  105. Seljak, U. 2012, JCAP, 2012, 004 [Google Scholar]
  106. Sheth, R. K. 1995, MNRAS, 274, 213 [Google Scholar]
  107. Sheth, R. K., & Lemson, G. 1999, MNRAS, 304, 767 [CrossRef] [Google Scholar]
  108. Shirasaki, M., Sugiyama, N. S., Takahashi, R., & Kitaura, F.-S. 2021, Phys. Rev. D, 103, 023506 [Google Scholar]
  109. Sinigaglia, F., Kitaura, F.-S., Balaguera-Antolínez, A., et al. 2021, ApJ, 921, 66 [NASA ADS] [CrossRef] [Google Scholar]
  110. Sinigaglia, F., Kitaura, F.-S., Balaguera-Antolínez, A., et al. 2022, ApJ, 927, 230 [NASA ADS] [CrossRef] [Google Scholar]
  111. Sinigaglia, F., Kitaura, F.-S., Nagamine, K., & Oku, Y. 2024a, ApJ, 971, L22 [Google Scholar]
  112. Sinigaglia, F., Kitaura, F.-S., Nagamine, K., Oku, Y., & Balaguera-Antolínez, A. 2024b, A&A, 682, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Somerville, R. S., Lemson, G., Sigad, Y., et al. 2001, MNRAS, 320, 289 [Google Scholar]
  114. Stadler, J., Schmidt, F., & Reinecke, M. 2023, JCAP, 2023, 069 [Google Scholar]
  115. Takada, M., Ellis, R. S., Chiba, M., et al. 2014, PASJ, 66, R1 [Google Scholar]
  116. Tassev, S., Zaldarriaga, M., & Eisenstein, D. J. 2013, JCAP, 6, 036 [CrossRef] [Google Scholar]
  117. Tegmark, M., Blanton, M. R., Strauss, M. A., et al. 2004, ApJ, 606, 702 [Google Scholar]
  118. Vos-Ginés, B., Avila, S., Gonzalez-Perez, V., & Yepes, G. 2024, MNRAS, 530, 3458 [Google Scholar]
  119. Wang, H., Mo, H. J., Yang, X., & van den Bosch, F. C. 2012, MNRAS, 420, 1809 [NASA ADS] [CrossRef] [Google Scholar]
  120. Wang, H., Mo, H. J., Yang, X., & Van den Bosch, F. C. 2013, ApJ, 772, 63 [NASA ADS] [CrossRef] [Google Scholar]
  121. Wang, H., Mo, H. J., Yang, X., Jing, Y. P., & Lin, W. P. 2014, ApJ, 794, 94 [NASA ADS] [CrossRef] [Google Scholar]
  122. Wang, Y., Zhai, Z., Alavi, A., et al. 2022, ApJ, 928, 1 [Google Scholar]
  123. Werner, K. F., & Porciani, C. 2020, MNRAS, 492, 1614 [NASA ADS] [CrossRef] [Google Scholar]
  124. Yahil, A., Strauss, M. A., Davis, M., & Huchra, J. P. 1991, ApJ, 372, 380 [NASA ADS] [CrossRef] [Google Scholar]
  125. Zaroubi, S., Hoffman, Y., Fisher, K. B., & Lahav, O. 1995, ApJ, 449, 446 [Google Scholar]
  126. Zel’dovich, Y. B. 1970, A&A, 5, 84 [NASA ADS] [Google Scholar]
  127. Zhao, C., Huang, S., He, M., et al. 2024, arXiv e-prints [arXiv:2411.07970] [Google Scholar]
1

Bayesian Reconstruction and Inference of Data-driven Generative Environments: https://github.com/pererossello/bridge-repo.

Appendix A: White-noise parametrization

In the forward model we have adopted a white-noise parametrization, ν, of i.i.d. standard normal variates, which we transform into the linear density field through the linear operator

M δ ( ν ) ≡ δ 0 = ( Δ x ) − 3 F − 1 P F ν , Mathematical equation: $$ \begin{aligned} \mathcal{M} _\delta (\boldsymbol{\nu }) \equiv {\boldsymbol{\delta }}_0 = (\Delta x)^{-3}\, \mathsf F ^{-1}\sqrt{\mathsf{P }}\;\mathsf F \,\boldsymbol{\nu }, \end{aligned} $$(A.1)

where F is the discrete Fourier-transform (DFT) linear operator on an N3 mesh of side length L, given by

F nm = ( Δ x ) 3 exp ( − i k n · x m ) , F nm − 1 = L − 3 exp ( i k m · x n ) . Mathematical equation: $$ \begin{aligned} \mathsf F _{nm}&= (\Delta x)^3 \exp \left(-i \boldsymbol{k}_n \cdot \boldsymbol{x}_m\right), \\ \mathsf F ^{-1}_{nm}&= L^{-3} \exp \left(i \boldsymbol{k}_m \cdot \boldsymbol{x}_n\right) \, . \end{aligned} $$(A.2)

Here xn represents the Cartesian direction vector in commoving space at a voxel indexed by n. The same holds for the wave-vector kn in the Fourier-transformed space. We have defined the diagonal linear matter power spectrum matrix P nm ≡ δ nm K P lin ( | k n | ) Mathematical equation: $ \mathsf P _{nm}\equiv \delta ^{\mathrm{K} }_{nm}\,P_{\mathrm{lin} }(|{\boldsymbol{k}}_{n}|) $. The density contrast, δ0, is then Gaussian-distributed with zero mean and covariance:

C = ( Δ x ) − 3 F − 1 P F . Mathematical equation: $$ \begin{aligned} \mathsf C = (\Delta x)^{-3} \mathsf F ^{-1} \mathsf P \mathsf F \ . \end{aligned} $$(A.3)

In principle, one can choose to parametrize the field either in terms of δ0 directly or through the underlying white noise ν. Sampling δ0 directly removes the linear transformation from the forward model but the fast Fourier transforms must still be performed at every Markov chain Monte Carlo step to evaluate the negative log prior and to draw new field realisations. We have tested both parametrizations and did not observe significant differences in computational cost or sampling behaviour. Nonetheless, we favour the white-noise parametrization ν for its simpler implementation.

Another option is to parametrize the field directly in Fourier space. Defining

u ̂ = 2 N 3 L 6 F ν , Mathematical equation: $$ \begin{aligned} \hat{\boldsymbol{u}} = \sqrt{\frac{2 N^3}{L^6}} \, \mathsf F \boldsymbol{\nu }, \end{aligned} $$(A.4)

the components of u ̂ Mathematical equation: $ \hat{\boldsymbol{u}} $ are i.i.d standard normal variates for both real and imaginary parts. The density contrast is then given by

δ 0 = L 3 2 F − 1 P u ̂ . Mathematical equation: $$ \begin{aligned} {\boldsymbol{\delta }}_0 = \sqrt{\frac{L^3}{2}} \, \mathsf F ^{-1} \sqrt{\mathsf{P }} \, \hat{\boldsymbol{u}} \ . \end{aligned} $$(A.5)

With this parametrization every parameter corresponds to a definite scale/wavenumber. We expected a different chain sampling behaviour, especially when using diagonal mass matrix adaptation during burn-in, but we saw no significant differences with respect to the white noise parametrization. We note that this parametrization requires the careful construction of a Hermitian-packed Fourier (N,N,N/2+1)-shaped array produced by a 3D real-to-complex fast Fourier transform out of N3 random real numbers, with explicit handling of the array structure and Nyquist planes.

A third option we considered is to separate each Fourier mode into amplitude and phase, δ ̂ 0 = | δ 0 | exp ( i φ ) . Mathematical equation: $ \hat{\boldsymbol{\delta }}_{0} = |{\boldsymbol{\delta }}_{0}|\exp (i\boldsymbol{\varphi }) \ . $ The parameters to be sampled are then split into:

  • N3/2−4 phases uniformly distributed in [0, 2π)

  • N3/2−4 amplitudes following a Rayleigh distribution.

  • 8 Nyquist modes following a Gaussian distribution.

This parametrization may be advantageous if we want to fix the amplitudes to the theoretical linear matter power spectrum, effectively reducing the parameter space in half. But again, for our work, we did not come across any obvious advantages of using this approach over the more straightforward white-noise parametrization.

Appendix B: Nonlocal bias and cosmic-web relation

Let us consider a scalar field η, which could be, for example, the matter overdensity field η = δ or the gravitational potential η = Φ. We begin by constructing the Hessian matrix of this field, defined as

H jk = ∂ 2 η ∂ x j ∂ x k = η , j k , Mathematical equation: $$ \begin{aligned} \mathcal{H} _{jk} = \frac{\partial ^2 \eta }{\partial x_j \partial x_k} = \eta _{,jk}\,, \end{aligned} $$(B.1)

which encodes the second derivatives of η with respect to spatial coordinates.

We then compute the eigenvalues of the Hessian, denoted by λ(1) ≥ λ(2) ≥ λ(3). These eigenvalues characterize the local curvature of the field in each principal direction. From them, we define the three rotationally invariant scalar quantities, or invariants, of the Hessian:

  • I1 = λ(1) + λ(2) + λ(3), the trace of the Hessian,

  • I2 = λ(1)λ(2) + λ(1)λ(3) + λ(2)λ(3), the sum of principal minors,

  • I3 = λ(1)λ(2)λ(3), the determinant of the Hessian.

According to Kitaura et al. (2022), the invariants of the Hessian matrix can be directly linked to different regions of the cosmic-web. For simplicity, and without loss of generality, this classification can be expressed for a threshold λth = 0 as follows:

  • Knots: I3 > 0, I2 > 0, I1 > λ(1)

  • Filaments: I3 < 0, I2 < 0, or I3 < 0, I2 > 0, λ ( 3 ) < I 1 < λ ( 3 ) − λ ( 2 ) λ ( 3 ) λ ( 1 ) Mathematical equation: $ \lambda ^{(3)} < I_1 < \lambda ^{(3)} - \frac{\lambda ^{(2)} \lambda ^{(3)}}{\lambda ^{(1)}} $

  • Sheets: I3 > 0, I2 < 0, or I3 < 0, I2 > 0, λ ( 1 ) − λ ( 2 ) λ ( 3 ) λ ( 1 ) < I 1 < λ ( 1 ) Mathematical equation: $ \lambda ^{(1)} - \frac{\lambda ^{(2)} \lambda ^{(3)}}{\lambda ^{(1)}} < I_1 < \lambda ^{(1)} $

  • Voids: I3 < 0, I2 > 0, I1 < λ(1).

Moreover, these invariants can be directly related to the nonlocal bias operators commonly used in perturbation theory. Specifically for η = Φ, i.e., long-range nonlocal bias:

  • δ= I1, the local density contrast,

  • s 2 = 2 3 I 1 2 − 2 I 2 Mathematical equation: $ s^2 = \frac{2}{3}I_1^2 - 2I_2 $, the tidal shear squared,

  • s 3 = − I 1 I 2 + 3 I 3 + 2 9 I 1 3 Mathematical equation: $ s^3 = -I_1 I_2 + 3 I_3 + \frac{2}{9}I_1^3 $, the cubic tidal bias term.

This formulation thus bridges cosmic-web morphology with perturbative bias theory, showing that both short-range (η = δ) and long-range (η = Φ) nonlocal bias terms can be systematically described in terms of the same invariant-based framework. We refer to the cosmic-web classification based on η = δ as the δ-web, and when based on η = Φ, as the Φ-web.

Appendix C: Numerical results of TEST2 and TEST3

Here we present the summary statistics resulting from TEST2 and TEST3. The case of TEST2 – when the resolution is doubled to 5 h−1 Mpc with respect to TEST1 – yields similarly accurate results across all summary statistics (see Figs. C.1, C.2, and C.3).

Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

TEST2: Analogous to Fig. 6.

Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

TEST2: Initial conditions. Analogous to Fig. 7.

Thumbnail: Fig. C.3. Refer to the following caption and surrounding text. Fig. C.3.

TEST2: Final conditions. Analogous to Fig. 8.

The additional numerical case TEST3 in which the bias parameters are assumed to be known, allows us to implement the most complex bias model as in Coloma-Nadal et al. (2024). Specifically, we define distinct bias prescriptions across the hierarchical cosmic-web, using a classification based on both the tidal field tensor and the Hessian of the density field. This results in 16 regions, each characterised by its own power-law, exponential threshold, and stochastic bias parameters. As shown in Figs. C.4C.5, and C.6, the reconstruction maintains the same level of precision under this more sophisticated setting, as in TEST1 and TEST2, assuming, however, the bias parameters to be known.

Thumbnail: Fig. C.4. Refer to the following caption and surrounding text. Fig. C.4.

TEST3: Analogous to Fig. 6.

Thumbnail: Fig. C.5. Refer to the following caption and surrounding text. Fig. C.5.

TEST3: Initial conditions. Analogous to Fig. 7.

Thumbnail: Fig. C.6. Refer to the following caption and surrounding text. Fig. C.6.

TEST3: Final conditions. Analogous to Fig. 8.

All Tables

Table 1.

Setting the different numerical cases.

In the text
Table 2.

Values of bias parameters by cosmic-web region of the reference for the TEST1 and TEST2 cases (8 bias parameters).

In the text
Table 3.

Values of bias parameters by combined cosmic-web regions for the TEST3 case (64 bias parameters).

In the text

All Figures

Thumbnail: Fig. 1. Refer to the following caption and surrounding text. Fig. 1.

Schematic overview of the BRIDGE pipeline. The framework combines a scalable, differentiable structure-formation model with a flexible field-level bias prescription, all implemented in a GPU-accelerated JAX environment. This enables efficient Bayesian inference of the primordial density field from tracer observations, with support for complex, non-local-bias models and multi-resolution analysis. The modular design allows for the seamless integration of physical models while maintaining end-to-end differentiability and high computational performance.
In the text
Thumbnail: Fig. 2. Refer to the following caption and surrounding text. Fig. 2.

Example of fuzzy cosmic-web classification with N = 256 and ΔL = 1.7 h−1 Mpc. Top left: Evolved dark-matter density contrast field. Bottom left: ‘Hard’ Φ-web classification. Right panels: Fuzzy membership weights, pi(V), pi(S), pi(F), and pi(K), for voids, sheets, filaments, and knots, respectively, as defined in Eq. (17).
In the text
Thumbnail: Fig. 3. Refer to the following caption and surrounding text. Fig. 3.

TEST1: Mean auto-correlation ξ(n) as a function of lag, n, for a chain of length of M ≈ 1000, evaluated over 104 randomly selected voxels. The red curve shows the ensemble average, ⟨ξ(n)⟩,  across parameters, with the shaded band indicating the ±1σ dispersion. The horizontal dashed line marks the threshold, ξ = 0.1, used to assess mixing efficiency.
In the text
Thumbnail: Fig. 4. Refer to the following caption and surrounding text. Fig. 4.

Computing times for single-gradient evaluations of the forward model employing ALPT evolution and the HICOBIAN bias model for different mesh sizes, run on a single NVIDIA A100-SXM4 GPU equipped with 40 GB of on-board HBM2 memory. For reference, the full burn-in and sampling procedure for the 1283 case in this study required approximately 100 GPU hours.
In the text
Thumbnail: Fig. 5. Refer to the following caption and surrounding text. Fig. 5.

Gelman–Rubin convergence diagnostic computed for three independently initialised chains of 250 samples after convergence. This statistic is calculated for all cells in Fourier space, and we show their distribution with scale, k. Values close to one indicate that the drawn samples have converged to the same distribution.
In the text
Thumbnail: Fig. 6. Refer to the following caption and surrounding text. Fig. 6.

TEST1: Spatial summary of 500 independent reconstructions (N = 128). Slices of ΔL = 10 h−1 Mpc (one voxel width). Top row: Reconstructed initial density contrast. Bottom row: Reconstructed tracer field. Columns from left to right: (1) Pixel-wise standard deviation across the 500 independent reconstructions; (2) mean reconstructed field; (3) one representative reconstruction sample; and (4) true reference field.
In the text
Thumbnail: Fig. 7. Refer to the following caption and surrounding text. Fig. 7.

TEST1: Comparison of summary statistics of the initial density contrast field for the mock reference (dashed red line) and 500 independent reconstructions (solid blue line). Columns from left to right: (1) Monopole power spectrum; (2) quadrupole; (3) reduced bispectrum as a function of triangle opening angle, θ; and (4) cross-correlation coefficient C(k) between the reference and reconstruction (solid blue line) and the propagator (dashed red line), defined as the cross-correlation between the mock tracer density contrast and the initial density reference. In columns (1)–(3), the lower sub-panels display the ratios with respect to the reference. Shaded bands denote the 1σ intervals from the posterior samples. Vertical dashed lines indicate the isotropic Nyquist frequency, k nyq / 3 Mathematical equation: $ k_{\mathrm{nyq}}/\sqrt{3} $.
In the text
Thumbnail: Fig. 8. Refer to the following caption and surrounding text. Fig. 8.

TEST1: Comparison of summary statistics of the mock reference and 500 independent reconstructions evolved through the forward model. Panel definitions are identical to those in Fig. 7, except that in the rightmost panel the propagator is replaced by the maximal tracer field cross-correlation and that we show the bispectrum of the dark-matter field and a tracer field using a local bias (without cosmic web) with bias parameters set to those of knots.
In the text
Thumbnail: Fig. 9. Refer to the following caption and surrounding text. Fig. 9.

TEST1: Posterior distributions for the 500 independent samples of the α and β bias parameters in each cosmic-web region of the Φ web. Red dots and lines show the values of the parameters for the ground-truth catalogue. The posteriors are smoothed with a 2D Gaussian filter with standard deviation equal to 0.75 bin widths. TEST2 yields equivalent results.
In the text
Thumbnail: Fig. C.1. Refer to the following caption and surrounding text. Fig. C.1.

TEST2: Analogous to Fig. 6.
In the text
Thumbnail: Fig. C.2. Refer to the following caption and surrounding text. Fig. C.2.

TEST2: Initial conditions. Analogous to Fig. 7.
In the text
Thumbnail: Fig. C.3. Refer to the following caption and surrounding text. Fig. C.3.

TEST2: Final conditions. Analogous to Fig. 8.
In the text
Thumbnail: Fig. C.4. Refer to the following caption and surrounding text. Fig. C.4.

TEST3: Analogous to Fig. 6.
In the text
Thumbnail: Fig. C.5. Refer to the following caption and surrounding text. Fig. C.5.

TEST3: Initial conditions. Analogous to Fig. 7.
In the text
Thumbnail: Fig. C.6. Refer to the following caption and surrounding text. Fig. C.6.

TEST3: Final conditions. Analogous to Fig. 8.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.