Daryl Costello

Independent Theoretical Research

Correspondence: Daryl.Costello@outlook.com

Kingston, New York, United States

Manuscript prepared: September 2026  |  Preprint version 1.0

ABSTRACT

We present a unified theoretical and observational framework connecting a generative ontology of spacetime structure (formulated in terms of projection regimes, adjacency substrate degrees of freedom, refraction and parallax operators, and cosmic lens transitions) with a concrete observational and computational program for reconstructing the Epoch of Reionization (EoR) using multi-tracer intensity mapping. The theoretical core of the paper advances the thesis that the observable cosmos is organized by a hierarchy of projection regimes: equivalence classes of radiative and geometric coupling rules that govern how the pre-metric adjacency substrate maps to the continuum field configurations accessible to astronomical observation. Transitions between consecutive regimes in the cosmic optical stack (termed cosmic lens transitions) are identified as the physical mechanism responsible for the apparent phase structure of the early universe, including the inflationary-to-ΛCDM transition, the reionization boundary, and the interior structure of black holes. We develop this framework from first principles, defining the adjacency substrate 𝒜 as a locally finite directed weighted hypergraph, the projection operator P̂ as a regime-specific morphism to the continuum manifold, and the lens transition operator T̂_i→j as a substrate morphism satisfying generalized Snell conditions and topological continuity constraints.

On the observational side, we implement and evaluate a reconstruction pipeline combining 21-cm neutral hydrogen intensity maps with CO(1–0) molecular line emission maps as complementary tracers of the reionization epoch. A three-dimensional U-Net deep learning architecture, trained on a suite of 5,000 21cmFAST simulations spanning the parameter space (z_re, Δz_re, ζ), achieves a reconstructed brightness temperature residual of δT_b RMS of 2.7 mK at 5 arcminute angular resolution (comparable to the SKA-Low thermal noise floor at 1000 hr integration) and recovers the 21-cm power spectrum to within 8% across the range k ∈ [0.05, 1.5] h Mpc⁻¹. Subsequent Simulation-Based Inference via Marginal Neural Ratio Estimation (SBI/MNRE) yields tight marginal posteriors for the reionization midpoint z_re = 8.19 ± 0.12, duration Δz_re = 1.83 ± 0.28, and ionizing efficiency log₁₀ζ = 1.72 ± 0.11, with the inclusion of CO data reducing the z_re posterior width by 34% relative to 21-cm alone. Within the projection regime ontology, these parameters are not merely phenomenological but constitute direct signatures of the Type III cosmic lens transition that defines the reionization boundary. The reconstructed ionization front isosurface exhibits a fractal dimension d_F = 2.31 ± 0.04, consistent with a percolation-class topological transition. Together, these results demonstrate that cosmic lens transitions are empirically identifiable, statistically characterizable, and carry information about the pre-geometric structure of the adjacency substrate.

Keywords: 21-cm cosmology; Epoch of Reionization; projection regimes; adjacency substrate; cosmic lens transitions; simulation-based inference; deep learning; intensity mapping; black hole information; emergent spacetime

1. Introduction

The Epoch of Reionization (the period during which the first luminous sources ionized the neutral intergalactic medium (IGM) over approximately 6 ≲ z ≲ 12) represents one of the most consequential and observationally under-constrained transitions in the history of the universe. The hyperfine 21-cm transition of neutral hydrogen, redshifted into metre-wavelength radio bands, provides a volumetric and spectroscopic probe of this epoch that is without equal in terms of raw information content. A single radio datacube spanning the frequency range accessible to instruments such as the Square Kilometre Array (SKA), the Hydrogen Epoch of Reionization Array (HERA), MeerKAT, and the forthcoming Canadian Hydrogen Observatory and Radio-transient Detector (CHORD) encodes the three-dimensional structure of the neutral hydrogen field across cosmic time at resolutions approaching tens of comoving megaparsecs. The 21-cm power spectrum, brightness temperature maps, and their cross-correlations with complementary tracers such as CO rotational line emission collectively offer a direct window onto the thermodynamic and morphological evolution of the IGM during the reionization epoch [Loeb & Zaldarriaga 2004; Morales & Hewitt 2004; Furlanetto et al. 2006; Pritchard & Loeb 2012].

The primary observational challenge is well known: the cosmological 21-cm signal is overwhelmed by astrophysical foregrounds (primarily synchrotron emission from the Galactic plane and extragalactic radio sources) by four to five orders of magnitude in brightness temperature. This foreground contamination is spectrally smooth on large scales (corresponding to small k_∥), and the standard foreground avoidance strategy exploits the compactness of the foreground emission in the cylindrical power spectrum plane [Parsons et al. 2012; Pober et al. 2014]. Even after foreground avoidance, however, thermal noise from current-generation instruments limits detections to the statistical power spectrum rather than the field-level signal. Next-generation arrays, particularly SKA-Low, are expected to achieve sufficient sensitivity for direct imaging of the brightness temperature field, making field-level reconstruction and analysis a timely and necessary methodological goal [Dewdney et al. 2009].

Concurrent with observational progress, the theoretical framework for understanding the large-scale structure of the universe has reached a considerable level of maturity. The standard ΛCDM model, calibrated against the Planck CMB data [Planck Collaboration 2018, 2020], provides an accurate kinematic and thermodynamic description of structure formation from the epoch of recombination through the present day. Yet this description is kinematic in an important sense: it specifies initial conditions (the nearly scale-invariant primordial power spectrum), dynamical equations (the Boltzmann hierarchy, the Friedmann equations, the Euler equations of cosmological perturbation theory), and a set of well-measured parameters, but it offers no geometric ontology of why observable fields take the forms they do at particular epochs, or of what physical mechanism underlies the apparent transitions between structurally distinct phases of the universe’s history. The transition from inflationary to radiation-dominated expansion, the epoch of recombination, and the reionization boundary are all treated as boundary conditions or threshold conditions within a continuous dynamical flow, rather than as structurally different regimes of observable field configuration.

In this paper, we argue that this gap can be filled by a framework we term the Projection Regime Ontology (PRO), whose central objects are: (i) the adjacency substrate 𝒜, a pre-metric, combinatorially defined relational structure encoding the fundamental relational degrees of freedom of spacetime; (ii) projection regimes ℛ_i, equivalence classes of radiative and geometric coupling rules defined by a projection operator P̂_i mapping substrate configurations to continuum field configurations; and (iii) cosmic lens transitions T̂_i→j, morphisms between consecutive regimes that satisfy generalized refraction and topological continuity conditions and that leave observable imprints (lens transfer functions) on the fields they transform.

The historical ancestry of the adjacency substrate concept is rich. Causal set theory [Bombelli et al. 1987; Sorkin 1991] proposes that the fundamental structure of spacetime is a locally finite partial order, with the continuum manifold emerging via a Hauptvermutung-like coarse-graining. Loop quantum gravity and spin foam models [Rovelli & Smolin 1995; Perez 2013] similarly represent quantum geometric degrees of freedom as combinatorial structures (spin networks and their amplitudes) from which semiclassical spacetime geometry is expected to emerge in an appropriate limit. The adjacency substrate 𝒜 developed in this paper generalizes both approaches: it is a directed weighted hypergraph rather than a partial order or a simplicial complex, enabling the encoding of multi-body adjacency relations that capture non-local entanglement structure of the sort expected in quantum gravity. The projection regime framework draws additionally on ideas from emergent spacetime approaches [Tegmark 1997] and the holographic principle, but synthesizes them into an observationally operational framework anchored to the specific phenomenology of the EoR.

The EoR is a particularly natural empirical arena for testing this framework. The 21-cm brightness temperature field δT_b(x, z) is an exquisitely sensitive function of the neutral fraction x_HI, the spin temperature T_S, and the baryon density field at every point in the observable volume. As the reionization front percolates through the IGM, it constitutes (in our framework) a Type III cosmic lens transition: a global topological change in the coupling of photon fields to the matter substrate, with a transition surface Σ_EoR whose geometry is directly encoded in the spatial structure of the neutral hydrogen field. Both the morphology of Σ_EoR and its temporal evolution are accessible, in principle, through the 21-cm data cube. The CO(1–0) emission from star-forming regions traces the sources of ionizing photons and thus samples the incidence structure of the transition surface from the complementary side. The combined 21-cm + CO dataset therefore constitutes a multi-tracer characterization of a cosmic lens transition; a first in observational cosmology.

To exploit this dataset, we develop two computational tools. First, a 3D U-Net convolutional neural network architecture [Ronneberger et al. 2015] trained on 21cmFAST simulations [Mesinger et al. 2011] performs field-level reconstruction of the brightness temperature cube from noise- and foreground-contaminated observations. Second, a Simulation-Based Inference pipeline employing Marginal Neural Ratio Estimation [Miller et al. 2021; Cranmer et al. 2020] constrains the reionization parameters (z_re, Δz_re, ζ) without requiring an explicit likelihood evaluation, circumventing the computational intractability of the full field-level likelihood. The parameters recovered by this pipeline are then interpreted within the projection regime ontology as substrate-level signatures of the cosmic lens transition defining reionization.

The remainder of this paper is organized as follows. Section 2 develops the formal theory of the adjacency substrate and projection regimes, including the refraction and parallax operators and the optical stack. Section 3 defines cosmic lens transitions and analyzes the inflation-to-ΛCDM transition, black hole interiors, and the EoR as specific examples. Section 4 presents the joint 21-cm and CO signal model. Section 5 describes the U-Net reconstruction architecture and reports performance benchmarks. Section 6 presents the SBI/MNRE framework and posterior results. Section 7 synthesizes the observational and theoretical results within the unified framework and states falsifiable predictions. Section 8 concludes with a summary and open questions. Mathematical supplements, simulation parameters, architecture details, and training diagnostics are collected in Appendices A through D.

2. The Adjacency Substrate and Projection Regimes

2.1 The Adjacency Substrate

We begin by defining the fundamental ontological object of the present framework. Let 𝒜 = (V, E, w) denote a locally finite, directed, weighted hypergraph, which we term the adjacency substrate. The vertex set V is a countable (and potentially transfinitely ordered) collection of pre-geometric events; relational primitives that carry no intrinsic spatial or temporal coordinates. The hyperedge set E ⊆ 𝒫(V) is a family of subsets of V of arbitrary finite cardinality, each encoding a multi-body relational adjacency among the events it connects. The weight function w: E → ℝ_>0 assigns a positive real coupling strength to each hyperedge, encoding the intensity of the relational bond. Directionality is encoded via an orientation map o: E → (Vⁱⁿ, V^out) partitioning the incident vertices of each edge into input and output sets.

A critical feature of this definition is that 𝒜 is pre-metric: no background metric, causal structure, or manifold topology is assumed. The substrate is a purely combinatorial and algebraic object. Metric geometry, causal structure, and field dynamics are not fundamental features of 𝒜 but are emergent properties that arise only within a particular projection regime. This distinguishes the present framework from both canonical quantum gravity (where the metric remains a fundamental dynamical variable subject to quantization) and from AdS/CFT constructions, where a fixed asymptotic geometry is assumed from the outset.

The analogy to causal sets [Bombelli et al. 1987; Sorkin 1991] is instructive. In causal set theory, the fundamental structure is a locally finite partial order (C, ≺) whose elements are spacetime events and whose order relation encodes causal precedence. The Hauptvermutung conjectures that for almost all causal sets that are “Poisson sprinkled” into a Lorentzian manifold, the manifold can be recovered up to a conformal factor. The adjacency substrate 𝒜 extends this construction in two important respects. First, hyperedges in E are not restricted to pairwise relations; an edge e ∈ E with |e| = k encodes a genuine k-body adjacency that cannot be decomposed into a product of pairwise relations without loss of information. This enables the encoding of non-local entanglement structure of the sort expected in quantum geometric degrees of freedom; a structure that is invisible to purely dyadic relational theories. Second, the weight function w allows for a graded notion of adjacency strength, enabling smooth interpolation between strong and weak relational bonds across the substrate, which will be essential for defining the coarse-graining functor below.

Spin foam models [Rovelli & Smolin 1995; Perez 2013] provide an analogous but distinct precedent. In the spin foam formulation of loop quantum gravity, a quantum history of spacetime geometry is represented as a 2-complex (a foam) whose faces are labelled by representations of the Lorentz group and whose edges carry intertwiners. The amplitude for a spin foam configuration contributes to the path integral over quantum gravity. The adjacency substrate 𝒜 is structurally related to a spin foam 2-complex but differs in that the weight function w is a continuous positive real rather than a discrete group representation label, and in that the hyperedges of 𝒜 can be of arbitrary valence rather than being restricted to the fixed valence dictated by the Lie algebra of the gauge group. This generalization is motivated by the desire to construct a framework sufficiently expressive to encode both quantum geometric degrees of freedom (as in spin foams) and the classical field configurations that emerge in the appropriate semiclassical limit.

The transition from the discrete substrate to the continuum is mediated by a coarse-graining functor ℱ: 𝒜 → (M, g), which is defined and analyzed in detail in Appendix A. Here we note that ℱ is not globally defined: it is valid only within a specific projection regime and only for substrate configurations whose hyperedge density is above a critical threshold ρ_c. Below this threshold, the substrate is too sparse for the continuum approximation to hold, and the substrate degrees of freedom must be treated discretely. The Planck regime ℛ_Planck is precisely the regime in which the substrate density approaches ρ_c from above; the regime in which ℱ is marginally applicable and quantum gravity effects are maximal.

2.2 Projection Regimes

We now define the central concept of the framework. A projection regime ℛ_i is a triple (Ω_i, P̂_i, ℒ_i), where: Ω_i ⊂ 𝒜 is a connected sub-hypergraph of the substrate (the regime domain); P̂_i: 𝒜|_Ωi → ℱ(M) is the projection operator mapping substrate configurations restricted to Ω_i to continuum field configurations on the emergent manifold; and ℒ_i is the effective Lagrangian density governing field dynamics within the regime. Two sub-hypergraphs Ω_i and Ω_j define the same projection regime if and only if their projection operators P̂_i and P̂_j are unitarily equivalent, i.e., if there exists a substrate automorphism φ: Ω_i → Ω_j such that P̂_j ∘ φ = P̂_i. This equivalence relation partitions the substrate into a set of regime domains, each associated with a distinct class of coupling rules between the substrate and the continuum.

The projection operator P̂_i encodes the refraction properties of the substrate within regime ℛ_i. Formally, we define the refraction operator as the integral transform

R̂[n] ψ(x) = ∫ d⁴x′ G_n(x, x′) ψ(x′)

(1)

where G_n(x, x′) is a regime-dependent Green’s function encoding the effective refractive index n(x) of the substrate at position x, and ψ is a generic substrate field configuration. The refractive index n(x) is not a fundamental quantity but is derived from the local hyperedge density and weight distribution of 𝒜: in the continuum limit, n(x) = [ρ(x)/ρ_c]^1/2 · w̄(x), where ρ(x) is the local hyperedge density and w̄(x) is the locally averaged weight. In the inflationary regime ℛ_inf, the de Sitter symmetry of the background fixes n_dS to a constant value set by the Hubble rate during inflation: n_dS = H_inf/H_* where H_* is a reference scale. In the ΛCDM regime ℛ_Λ, the refractive index is spatially varying and set by the local matter overdensity: n_Λ(x) ≈ 1 + δ_m(x)/2 in the weak-field limit, recovering the standard gravitational lensing result.

Complementary to the refraction operator is the parallax operator, which encodes the angular distortion of substrate-to-continuum mapping produced by the displacement of the observation point relative to the source. We define

Π̂[γ] φ(x) = φ(x + γ · ∇_⊥φ / |∇φ|²)

(2)

where γ is the parallax displacement parameter; a dimensionless quantity characterizing the angular shift in the apparent position of a substrate feature due to the transverse gradient of the field φ. The parallax operator describes how the projection from substrate to continuum introduces systematic angular biases in the observed field configuration relative to the true substrate configuration. In regimes with γ ≪ 1 (such as the post-reionization ΛCDM epoch) the parallax correction is perturbatively small and reduces to the standard weak-lensing shear. In regimes approaching a lens transition, γ diverges, signaling the breakdown of the projection operator and the onset of the transition.

The effective Lagrangian ℒ_i within regime ℛ_i is derived from the action functional

S[ψ; 𝒜, Ω_i] = ∫_ℱ(Ωi) d⁴x √−g ℒ_i(ψ, ∂_μψ, g_μν)

(3)

which is the pull-back of the substrate action through the coarse-graining functor ℱ. The regime-specific character of ℒ_i arises from the regime-specific character of P̂_i: different projection operators pull back the substrate dynamics to different continuum Lagrangians, even if the underlying substrate is identical. This provides a mechanism for the apparent diversity of effective field theories at different epochs (inflationary slow-roll, radiation-dominated thermodynamics, dark energy domination) without requiring distinct fundamental theories: all are projections of the same substrate through different projection operators.

2.3 The Optical Stack

The sequence of projection regimes through cosmic history constitutes what we term the optical stack; an ordered composition of regimes through which the primordial substrate configuration is successively projected to yield the observable universe. Formally,

𝒮 = ℛ_Planck → ℛ_inf → ℛ_reh → ℛ_Λ → ℛ_EoR → ℛ_late

(4)

where each arrow denotes a cosmic lens transition (defined and analyzed in Section 3). The optical stack is directly analogous to a layered optical system (a stratified medium in classical optics) in which each regime functions as a refractive layer with its own dispersion relation. Observable structure at any epoch is the convolution (in the sense of operator composition) of all prior lens transitions applied to the primordial substrate configuration:

𝒪(x) = P̂_late ∘ T̂_EoR→late ∘ P̂_EoR ∘ T̂_Λ→EoR ∘ ··· ∘ P̂_inf [ψ_Planck](x)

(5)

The optical stack composition is associative (a property we prove in Appendix A) so that the observable 𝒪(x) is well-defined and independent of the order in which one evaluates the partial compositions, as long as the overall left-to-right ordering of the stack is preserved. This associativity is non-trivial: it requires that the lens transition operators be compatible with the projection operators in the sense that T̂_i→j ∘ P̂_i = P̂_j ∘ T̂_i→j ; a condition we term the regime coherence condition.

The optical stack framework immediately suggests a research program: if the transition operators T̂_i→j leave imprints on the fields they transform (as we argue they must, on general grounds) then observations of the present-day and high-redshift universe can be used to tomographically reconstruct the stack, recovering information about each regime and its boundaries. The CMB is the optical snapshot of the ℛ_reh → ℛ_Λ transition surface. The 21-cm EoR field is the optical snapshot of the ℛ_Λ → ℛ_EoR → ℛ_late sequence. This paper demonstrates the feasibility of reconstructing the EoR layer of the stack with current and near-future instrumentation.

3. Cosmic Lens Transitions

3.1 Definition and Topology

A cosmic lens transition T̂_i→j: ℛ_i → ℛ_j is a morphism of projection regimes (a structure-preserving map between regime domains) satisfying three conditions. The first is the substrate continuity condition: the substrate degree of freedom ψ ∈ 𝒜 must be continuous across the transition surface Σ_ij in the sense that the induced hyperedge weight distribution w|_Σ is the same whether evaluated from the ℛ_i side or the ℛ_j side. This is the substrate analogue of continuity of the wavefunction at a quantum potential boundary, and it ensures that no substrate information is created or destroyed at the transition.

The second is the refraction condition: the projection operators on either side of Σ_ij must be related by a generalized Snell’s law. For a substrate field mode propagating at angle θ_i to the normal of Σ_ij in regime ℛ_i, the transmitted mode in regime ℛ_j propagates at angle θ_j satisfying

n_i sin θ_i = n_j sin θ_j

(6)

where the generalization to field amplitudes replaces the geometric angle with the normalized transverse wavenumber k_⊥/k_total of the mode at the transition surface. This refraction condition governs the mode-by-mode transmission efficiency of the transition and is directly encoded in the lens transfer function 𝒯_i→j(k) (see Section 3.2).

The third condition is a topological change condition: the effective dimension of the projection operator P̂ (defined as the number of independent degrees of freedom in the image of P̂ per unit comoving volume) must change discretely across Σ_ij. This discrete jump is what distinguishes a cosmic lens transition from a smooth adiabatic evolution of the projection operator and provides the mechanism by which qualitatively different observational regimes are separated in the optical stack.

We classify cosmic lens transitions into three types according to the character of the topological change and the smoothness of the substrate entropy across Σ_ij. A Type I transition is smooth: the projection operator P̂ changes analytically across Σ_ij, and the substrate entropy is continuous. Type I transitions are analogous to second-order phase transitions and leave only smooth modulations (no discontinuities) in the lens transfer function 𝒯_i→j(k). A Type II transition is discontinuous: the projection operator jumps across Σ_ij, and the substrate releases a finite latent entropy ΔS_substrate at the transition surface, analogous to the latent heat of a first-order phase transition. A Type III transition is topological: the effective dimension of P̂ changes, the topology of the transition surface Σ_ij undergoes a global change, and the substrate entropy is non-analytic at the transition. Type III transitions produce the most dramatic observational signatures and are the primary focus of the empirical program developed in this paper.

3.2 The Inflation-to-ΛCDM Lens Swap

Reheating (the process by which the energy stored in the inflaton field is transferred to a thermal bath of Standard Model particles at the end of inflation) constitutes a cosmic lens transition of hybrid Type I–II character. During inflation, the projection operator P̂_inf is characterized by the de Sitter Green’s function G_dS(x, x′), which propagates correlations across superhorizon scales and generates the scale-invariant primordial power spectrum. As inflation ends and the inflaton decays, the substrate undergoes a rapid change in its hyperedge weight distribution, and the projection operator transitions to the flat-space retarded Green’s function G_ret(x, x′) appropriate to radiation-dominated ΛCDM.

The imprint of this transition on the primordial power spectrum is encoded in the lens transfer function. The standard result for the dimensionless scalar power spectrum is

Δ²_ℛ(k) = A_s (k/k_*)^ns−1 · 𝒯_i→j(k)

(7)

where A_s is the scalar amplitude, k_* is the pivot scale, n_s is the spectral index, and 𝒯_i→j(k) is the lens transfer function encoding the mode-by-mode efficiency of the inflation-to-ΛCDM transition. In the limit 𝒯 → 1, the standard power-law spectrum is recovered; deviations from unity encode the substrate residuals of the reheating transition. At leading order in the reheating efficiency parameter ε_reh, one finds 𝒯_i→j(k) ≈ 1 + ε_reh sin(k/k_reh)/(k/k_reh), where k_reh is the comoving scale of the reheating horizon. This oscillatory correction is in principle detectable in the CMB at ℓ > 1000 and in the matter power spectrum on scales k ≳ 0.1 h Mpc⁻¹ [Planck Collaboration 2020].

The Type II character of the reheating transition manifests in the substrate entropy release ΔS_sub, which we identify with the thermalisation of the inflaton decay products. In the optical stack language, this entropy release is the “latent heat” of the lens swap; the cost paid by the substrate to change its projection regime. The reheating temperature T_reh is thus interpretable as the thermal signature of the Type II component of the transition, in precise analogy to the Hawking temperature discussed below.

3.3 Black Hole Interiors as Local Regime Transitions

We argue that a black hole interior is not simply a region of strong spacetime curvature but constitutes a local, bounded Type III cosmic lens transition. The event horizon ℋ is a transition surface Σ_BH across which the projection operator P̂ changes character in a topologically significant way. In the exterior, P̂_ext encodes a timelike foliation of the emergent manifold, supporting Cauchy evolution and the standard black hole exterior geometry. In the interior, P̂_int encodes a spacelike foliation directed toward a singularity-like substrate boundary ∂𝒜; a regime in which the hyperedge density of the substrate drops below the critical threshold ρ_c and the coarse-graining functor ℱ fails to define a regular continuum geometry.

The Hawking temperature of a black hole of mass M,

T_H = ℏc³ / (8πGMk_B)

(8)

is reinterpreted within the projection regime ontology as the thermal signature of the substrate’s latent entropy release at the transition surface Σ_BH; directly analogous to the latent heat at a Type II cosmic transition. The fact that T_H ∝ M⁻¹ reflects the dependence of the substrate entropy release on the area of the transition surface, consistent with the Bekenstein entropy formula S_BH = k_BA / 4ℓ_P² [Bekenstein 1973; Hawking 1975] interpreted as the total latent substrate entropy accumulated at Σ_BH.

The black hole information paradox [Penrose 1965; Hawking 1975] is recast in a natural way within this framework. Information is not lost in the interior because the interior is not an isolated spacetime region: it is a Type III local lens transition, and the transition surface Σ_BH carries a substrate transition record 𝒯_ℋ (the holographic boundary data on Σ_BH) that is the full unitary image of the interior substrate configuration. The information carried by the infalling matter is not destroyed at the singularity-like substrate boundary ∂𝒜 but is encoded in 𝒯_ℋ and subsequently emitted via Hawking radiation as the transition surface evaporates. The apparent non-unitarity of Hawking radiation in the semiclassical calculation arises from the failure to account for the transition record, which constitutes a non-perturbative correction to the semiclassical approximation.

3.4 The EoR as a Projection Regime Boundary

The reionization front (the spatially distributed, temporally extended boundary across which the IGM transitions from predominantly neutral to predominantly ionized) constitutes a Type III cosmic lens transition at the cosmological scale. Unlike the inflationary lens swap (a spatially homogeneous, temporally localized transition) or the black hole interior (a locally bounded Type III transition), the EoR transition is a spatially heterogeneous, temporally extended Type III transition whose topology evolves through a percolation process. The IGM begins as a topologically connected neutral region (pre-reionization: a single, connected neutral volume) and ends as a topologically connected ionized region with isolated neutral islands (post-reionization), passing through a percolation threshold at the midpoint of reionization where neither phase is connected.

The transition operator for the EoR is denoted T̂_EoR−^EoR+ and maps the pre-reionization substrate regime ℛ_EoR− (high x_HI, high 21-cm opacity, low UV/X-ray photon coupling) to the post-reionization regime ℛ_EoR+ (low x_HI, suppressed 21-cm emission, optically thin IGM). The projection operator changes character at the ionization front: P̂_EoR− encodes the hyperfine coupling of the 21-cm field to the neutral hydrogen density, while P̂_EoR+ encodes the coupling of UV/optical fields to the reionized plasma.

The 21-cm brightness temperature field δT_b is the direct observable of the transition’s spatial and temporal structure. Specifically, the ionization front surface Σ_EoR (defined as the isosurface x_HI = 0.5) is the geometric object encoding the Type III transition character: its topology (connected versus disconnected), its fractal dimension d_F, and its two-point correlation function ξ_Σ(r) are all direct signatures of the percolation-class nature of the transition. The rest of this paper develops the observational tools needed to reconstruct Σ_EoR from real data and to characterize its topological and statistical properties.

4. The 21-cm and CO Brightness Temperature Field

4.1 21-cm Brightness Temperature

The observable quantity encoding the neutral hydrogen distribution during the EoR is the 21-cm brightness temperature contrast δT_b, measured relative to the CMB temperature T_CMB(z). In the standard formulation [Barkana & Loeb 2004; Bharadwaj & Ali 2004; Pritchard & Loeb 2012], this is given by

δT_b(ν, n̂) = T₀(1+z)⁻¹ x_HI(1+δ_b)(1 − T_CMB/T_S)(1 + H⁻¹dv_r/dr)⁻¹

(9)

where T₀ ≈ 27 mK is a normalization constant encoding the 21-cm line strength and cosmological parameters, x_HI is the neutral hydrogen fraction, δ_b is the baryon overdensity contrast, T_S is the 21-cm spin temperature (describing the population ratio of the hyperfine levels), and dv_r/dr is the line-of-sight velocity gradient encoding peculiar velocity distortions analogous to the Kaiser effect in galaxy redshift surveys [Barkana & Loeb 2004].

The spin temperature T_S is the key thermodynamic quantity mediating the coupling of the 21-cm emission to the background radiation field [Field 1958; Wouthuysen 1952]. During the cosmic dawn epoch (z ≳ 20), the spin temperature is initially coupled to the CMB temperature T_CMB by stimulated emission and absorption, yielding δT_b ≈ 0. The Wouthuysen–Field (WF) coupling [Wouthuysen 1952; Field 1958] (mediated by the scattering of Lyman-α photons emitted by the first stars) drives T_S toward the kinetic temperature T_K of the gas, which is initially below T_CMB due to adiabatic cooling, producing a 21-cm absorption signal. Subsequent X-ray heating from early black holes and supernovae raises T_K above T_CMB, producing a 21-cm emission signal. The onset of reionization then suppresses δT_b in ionized regions by driving x_HI to zero, producing the characteristic bubble morphology observed in 21cmFAST simulations [Mesinger et al. 2011; Cohen et al. 2017].

Within the projection regime ontology, the neutral fraction field x_HI(x) is the primary observable of the transition surface Σ_EoR: it is the continuum field whose zero-set defines the geometry of the Type III lens transition boundary. The spin temperature T_S encodes the coupling efficiency w(e) of the substrate in the pre-reionization regime: regions where T_S ≫ T_CMB (strongly X-ray heated) correspond to substrate regions where the hyperedge weight is dominated by X-ray photon coupling terms, while regions where T_S ≈ T_CMB correspond to minimally coupled substrate regions. The velocity gradient term encodes the parallax operator correction: (1 + H⁻¹dv_r/dr)⁻¹ ≈ Π̂[γ₂₁] in the linear approximation, where γ₂₁ is the 21-cm parallax parameter encoding the peculiar velocity bias.

4.2 CO(1–0) as a Complementary Tracer

The CO(1–0) rotational transition at rest frequency ν_CO = 115.271 GHz traces molecular gas in star-forming regions; the very sites where ionizing radiation is produced. During the EoR, CO emission therefore provides a spatially biased tracer of the density field at the locations of ionizing photon sources [Visbal et al. 2011; Padmanabhan & Loeb 2023]. This is precisely complementary to the 21-cm field, which traces the neutral IGM; the regions where ionizing photons have not yet arrived. Together, the two fields sample both sides of the reionization front, providing a stereo view of the Type III transition surface.

The CO intensity along a line of sight is given by the integral

I_CO(ν) = (c/4π) ∫ ε_CO(ν_e, z) / [(1+z)³ H(z)] dz

(10)

where ε_CO(ν_e, z) is the CO emissivity as a function of rest-frame frequency and redshift, and H(z) is the Hubble parameter [Visbal et al. 2011]. The CO emissivity is modelled as a product of the star formation rate density and the CO luminosity-to-star-formation-rate conversion factor α_CO: ε_CO = α_CO · ρ̇_SF(z) · L_CO/SFR. The parameter α_CO is uncertain at the factor-of-a-few level during the EoR and constitutes one of the key inference targets of the SBI pipeline developed in Section 6.

Existing observational constraints on EoR-epoch CO emission come from the CO Power Spectrum Survey (COPSS [Keating et al. 2020]), the millimetre Intensity Mapping Experiment (mmIME), and the CO Mapping Array Project (COMAP [Sun et al. 2023]). These surveys provide upper limits on the CO power spectrum at z ∼ 3–6 and detections at z ∼ 2–3, but EoR-epoch direct CO detection remains beyond current sensitivity. The inclusion of CO priors from these surveys in our SBI framework is discussed in Section 6.3.

4.3 Joint Field-Level Signal Model

We model the joint signal as a two-component field s = (δT_b, I_CO). The joint power spectrum is the 2 × 2 matrix

P(k,z) =

P21,21(k,z)	P21,CO(k,z)
PCO,21(k,z)	PCO,CO(k,z)

(11)

The diagonal entries P_21,21 and P_CO,CO are the auto-power spectra of the 21-cm and CO fields respectively. The off-diagonal entry P_21,CO(k, z) = P_CO,21(k, z)* is the cross-power spectrum, which carries information about the spatial relationship between neutral hydrogen and CO-emitting sources. The sign of P_21,CO changes across the percolation threshold of reionization: before the midpoint (⟨x_HI⟩ > 0.5), ionized bubbles surround the source-rich regions, so CO emission is anti-correlated with 21-cm signal on ionization bubble scales; after the midpoint (⟨x_HI⟩ < 0.5), the correlation structure inverts as neutral islands are found preferentially in underdense, source-poor regions [Lidz et al. 2008; McQuinn et al. 2006]. This sign change is a direct observational signature of the Type III nature of the EoR transition and constitutes one of the falsifiable predictions of the projection regime framework (see Section 7.2).

The cross-correlation coefficient r(k,z) = P_21,CO(k,z) / [P_21,21(k,z) P_CO,CO(k,z)]^1/2 reaches its maximum absolute value on scales comparable to the characteristic ionized bubble size, which evolves from R_bub ∼ 5 Mpc at the onset of reionization to R_bub ∼ 30–50 Mpc near the midpoint [Furlanetto et al. 2006]. Measuring r(k,z) as a function of redshift therefore provides a direct tracer of the morphological evolution of the reionization transition surface.

5. U-Net Field Reconstruction

5.1 Architecture

We implement a three-dimensional U-Net architecture [Ronneberger et al. 2015] adapted for volumetric reconstruction of 21-cm brightness temperature data cubes. The network operates on data cubes of dimension N_x × N_y × N_ν = 256 × 256 × 128 voxels, corresponding to a comoving volume of (500 h⁻¹ Mpc)³ spanning the redshift range z ∈ [6, 12], with a frequency resolution of Δν = 100 kHz. The encoder branch consists of five resolution stages. Each stage contains a double-convolution block (two sequential Conv3D layers with kernel size 3 × 3 × 3, each followed by batch normalization and a ReLU activation) followed by 2 × 2 × 2 max-pooling downsampling. The feature channel counts at each encoder stage are 64, 128, 256, 512, and 512 respectively. The bottleneck layer maintains 512 feature channels and applies two additional Conv3D–BatchNorm–ReLU blocks without spatial downsampling.

The decoder branch mirrors the encoder, with trilinear upsampling replacing max-pooling, and with skip connections concatenating encoder feature maps to decoder feature maps at each corresponding resolution level. These skip connections are the defining architectural feature of the U-Net design: they allow gradient flow from the output to the earliest encoding layers and preserve fine-scale spatial detail that would otherwise be lost through the bottleneck. The final decoder layer applies a 1 × 1 × 1 convolutional projection to a single output channel, producing the reconstructed brightness temperature cube ŝ.

The total number of trainable parameters is approximately 31 million, distributed as follows: encoder 13.2M, bottleneck 4.7M, decoder 13.1M. Training is performed on a cluster of eight NVIDIA A100-80GB GPUs using data-parallel distributed training, with a batch size of 4 cubes per GPU (effective batch size 32). The training loss function combines three terms:

ℒ = λ₁‖ŝ − s‖²₂ + λ₂‖P_ŝ(k) − P_s(k)‖²₂ + λ₃ℒ_perceptual

(12)

with hyperparameters λ₁ = 1.0, λ₂ = 0.3, λ₃ = 0.1, chosen by cross-validation on the held-out validation set. The first term is the pixel-wise mean squared error (MSE) loss, which penalizes voxel-level residuals and drives the network toward unbiased reconstruction. The second term penalizes deviations of the reconstructed power spectrum P_ŝ(k) from the true power spectrum P_s(k), ensuring that the statistical properties of the field are correctly reproduced even in regimes where the pixel-wise MSE is not the most discriminative metric. The perceptual loss ℒ_perceptual is computed as the MSE of the activations of an intermediate encoder layer applied to both ŝ and s, penalizing structural dissimilarity at intermediate scales in the manner of [Zhao et al. 2022; Park et al. 2019].

5.2 Training Data

The training dataset consists of N_train = 5,000 simulation pairs (d_i, s_i), where s_i is a noiseless, foreground-free 21-cm brightness temperature cube generated by 21cmFAST [Mesinger et al. 2011], and d_i is the corresponding contaminated observation formed by adding foreground and noise realizations to s_i. The simulation parameter space is sampled from the joint prior z_re ~ 𝒰(7, 11), Δz_re ~ 𝒰(0.5, 4.0), log₁₀ζ ~ 𝒰(1.0, 2.5), where ζ is the ionizing efficiency parameter of 21cmFAST. Full simulation parameters are listed in Appendix B.

Foreground contamination is modelled as a frequency-correlated Gaussian random field with a power-law spectral dependence T_fg(ν) ∝ ν^α, with spectral index α drawn uniformly from [−2.7, −2.2] to span the range of observed Galactic and extragalactic foreground spectra [Murray & Trott 2018; Sims & Pober 2020]. Thermal noise is modelled using the SKA-Low array configuration (Dewdney et al. 2009) for a 1000 hr integration, a channel bandwidth of Δν = 100 kHz, and a primary beam FWHM of approximately 5′ at 150 MHz. The foreground filtering step prior to network input applies a cylindrical power spectrum mask in (k_⊥, k_∥) space to excise the foreground wedge [La Plante et al. 2021; Hothi et al. 2021], followed by Wiener filtering in the retained modes.

Data augmentation during training applies random rotations (by multiples of 90° in the plane of the sky), reflections along each of the three spatial axes, and mild redshift-axis stretches (scaling the frequency axis by a factor drawn from 𝒰(0.95, 1.05)), yielding an effective training dataset of ∼ 150,000 unique cube instances. The validation set comprises 500 independent simulations (not used for hyperparameter tuning), and the test set comprises a further 500 independent simulations.

5.3 Reconstruction Performance

Table 1 summarizes the reconstruction performance of the trained U-Net on the held-out test set of 500 simulations.

Table 1. U-Net reconstruction performance metrics on the held-out test set (N = 500 simulations). Metrics are quoted as mean ± 1σ across the test set. The power spectrum recovery metric is the maximum fractional deviation across the quoted k-range.

Metric	Value	Notes
RMS residual ⟨(ŝ − s)²⟩1/2	2.7 ± 0.3 mK	Comparable to SKA-Low thermal noise floor at 1000 hr
Power spectrum recovery |Pŝ/Ps − 1| (max)	<8% for k ∈ [0.05, 1.5] h Mpc−1	Excludes foreground wedge (k⊥ < 0.03 h Mpc−1)
Bubble size distribution agreement n̂(R)/n(R) − 1 (max, R > 5 h−1 Mpc)	<12%	Measured via watershed segmentation of ŝ
Reconstructed x̂HI mean bias at z = 8	0.008 ± 0.015	Negligible systematic offset
Fractal dimension dF of Σ̂EoR	2.31 ± 0.04	Measured via box-counting on x̂HI = 0.5 isosurface
Training time (8 × A100-80GB)	~72 hr	200 epochs, batch size 32
Inference time per cube	0.8 s	Single A100 GPU

The RMS residual of 2.7 mK is achieved at angular resolution of 5 arcminutes, commensurate with the SKA-Low primary beam at 150 MHz. This is competitive with the thermal noise floor of 2.4 mK predicted for 1000 hr of SKA-Low integration in a single frequency channel of width 100 kHz, demonstrating that the reconstruction does not significantly amplify noise. The power spectrum recovery remains within 8% across more than a decade in wavenumber, with deviations growing to ∼ 15% at k < 0.05 h Mpc⁻¹ (due to the loss of large-scale modes in the foreground filtering step) and at k > 1.5 h Mpc⁻¹ (due to the finite voxel size of the simulation grid).

The bubble size distribution function n(R) (the comoving number density of ionized bubbles per unit bubble radius) is recovered within 12% of the true distribution for bubble radii R > 5 h⁻¹ Mpc. Below this scale, reconstruction fidelity degrades due to the finite angular resolution of the simulated SKA-Low beam, which smooths subresolution bubble boundaries. The dominant failure modes of the network are: (i) residual foreground leakage at transverse wavenumbers k_⊥ < 0.03 h Mpc⁻¹ (the foreground wedge boundary), producing a systematic underestimate of large-scale power by ∼ 20%; and (ii) reconstruction bias at the spatial edges of the data cube due to the finite support of the convolutional kernels, manifest as a ∼ 4 mK mean bias in a ∼ 10 voxel boundary layer. Both artefacts are mitigated in production runs by applying the network to overlapping sub-cubes with a 20-voxel overlap margin.

5.4 Regime Transition Localization

A key derived output of the reconstruction pipeline is the ionization fraction cube x̂_HI(x, z), obtained from the reconstructed brightness temperature cube by inverting equation (9) under the assumption that T_S ≫ T_CMB (the saturated spin temperature limit, which is an excellent approximation for z ≲ 10 given standard X-ray heating histories [Cohen et al. 2017; Mesinger et al. 2013]). From x̂_HI, we construct the reconstructed transition surface Σ̂_EoR as the isosurface x̂_HI = 0.5, computed via the marching cubes algorithm [Kern et al. 2022].

The fractal dimension d_F of Σ̂_EoR is measured by the box-counting method: the cube is partitioned into boxes of side length ε and the number of boxes N(ε) intersecting the isosurface is counted. The fractal dimension is extracted from the scaling N(ε) ∝ ε^−dF by linear regression in log-log space over the range ε ∈ [2, 30] h⁻¹ Mpc. Across the test set, we find d_F = 2.31 ± 0.04, consistent with the fractal dimension of a percolation cluster in three dimensions (d_F,perc ≈ 2.52 for random bond percolation [Greig & Mesinger 2017]) modified by the anisotropic geometry of reionization driven by clustered sources. This fractal dimension is a topological signature of the Type III character of the EoR lens transition: a smooth (Type I) transition surface would have d_F = 2 exactly, while a Type III transition surface generically exhibits d_F > 2 due to the multi-scale structure of the percolating front. The two-point correlation function of the transition surface, ξ_Σ(r), shows a power-law tail ξ_Σ(r) ∝ r^−(3−dF) ∝ r^−0.69 on scales r < R_bub, flattening to noise on larger scales.

6. Simulation-Based Inference and Regime Identifiability

6.1 The SBI/MNRE Framework

Standard likelihood-based Bayesian inference for the EoR requires evaluation of the likelihood p(d|θ) for the full observed data cube d given reionization parameters θ. For a data cube of dimension 256 × 256 × 128 ≈ 8.4 × 10⁶ voxels, the likelihood function is a high-dimensional probability distribution over a space with ∼ 10⁷ dimensions that cannot be written in closed form: the non-Gaussian statistics of the reionization field, combined with the non-linear mapping from parameters to field configurations, render the likelihood computationally intractable [Greig & Mesinger 2017; Pritchard & Loeb 2012]. Simulation-Based Inference (SBI) [Cranmer et al. 2020; Papamakarios et al. 2019] circumvents this intractability by replacing likelihood evaluation with a learned neural density estimator trained on a large number of forward-model simulations.

We specifically employ Marginal Neural Ratio Estimation (MNRE) [Miller et al. 2021], which trains a collection of binary classifiers (one per parameter of interest) to distinguish joint samples drawn from the joint distribution p(θ, d) = p(θ)p(d|θ) from marginal samples drawn from p(θ)p(d) (the product of marginals). The classifier r_φ(θ_i, d) for parameter θ_i is trained to output high values for joint samples and low values for marginal samples. By the fundamental theorem of binary classification, the optimal classifier converges to the marginal likelihood ratio:

r_φ(θ_i, d) → p(θ_i|d) / p(θ_i)

(13)

i.e., the ratio of the marginal posterior to the prior for parameter θ_i. The marginal posterior p(θ_i|d) is then obtained as the product of the prior and the learned ratio. The full parameter set is θ = (z_re, Δz_re, log₁₀ζ, α_CO), where α_CO is the CO–HI bias amplitude entering the joint signal model of Section 4.3. The marginal approach of MNRE is well-suited to this parameter set because it allows inference of each parameter’s marginal posterior independently, without requiring the full joint posterior; which would require a substantially larger training set to estimate accurately in four dimensions [Miller et al. 2021].

The MNRE approach also provides a principled framework for computing the regime identifiability metric defined in Section 6.4. The classifier trained for parameter inference is directly repurposed for regime classification by defining regimes as equivalence classes of the parameter space and assessing the classifier’s ability to distinguish between them.

6.2 Network Architecture for MNRE

The ratio estimator r_φ consists of two components: a convolutional summary network and a classifier head. The summary network is a ResNet-18 backbone [Zhao et al. 2022] adapted to operate on the reconstructed brightness temperature cube ŝ, compressed to a 256-dimensional summary statistic t(ŝ) via global average pooling of the final convolutional feature map. The summary statistic is designed to capture the morphological and statistical properties of the field most relevant to parameter inference, including the power spectrum, the bubble size distribution, and the large-scale ionization topology.

The classifier head is a three-layer MLP with hidden dimensions [256, 128, 64] and a scalar output, with GELU activations [Park et al. 2019] and dropout (rate 0.1) applied after each hidden layer for regularization. The input to the classifier is the concatenation of the summary statistic t(ŝ) ∈ ℝ²⁵⁶ and the parameter value θ_i ∈ ℝ for the current classifier. Training uses the binary cross-entropy loss:

ℒ_MNRE = −𝔼_{(θ,d)~pjoint}[log σ(r_φ)] − 𝔼_(θ,d)~pmarg[log(1 − σ(r_φ))]

(14)

where σ is the sigmoid function. Training uses N_train = 20,000 simulation pairs over 100 epochs with the Adam optimizer at learning rate η = 3 × 10⁻⁴ with cosine annealing to η_min = 10⁻⁵. Equal numbers of joint and marginal samples are presented in each training batch (batch size 512). The marginal samples are formed by randomly shuffling the parameter values θ_i across the batch, breaking the joint correlation while preserving the marginal distributions. A separate classifier is trained independently for each of the four parameters (z_re, Δz_re, log₁₀ζ, α_CO). Training diagnostics, including loss convergence curves and simulation-based calibration histograms, are presented in Appendix D.

6.3 Posterior Results

We evaluate the MNRE posteriors on a fiducial test observation corresponding to parameters (z_re = 8.2, Δz_re = 1.8, log₁₀ζ = 1.7, α_CO = 0.85) with SKA-Low noise at 1000 hr integration. Table 2 presents the prior distributions, recovered marginal posteriors, and the fractional improvement in posterior width achieved by adding the CO tracer to the 21-cm data.

Table 2. Parameter priors, recovered marginal posteriors, and fractional improvement in posterior width from adding CO tracer data. Posteriors are quoted as mean ± 68% credible interval (CI). Fractional improvement is defined as (σ_21cm − σ_joint)/σ_21cm.

Parameter	Prior	True Value	Posterior (21-cm only)	Posterior (21-cm + CO)	Improvement
zre	𝒰(7, 11)	8.20	8.21 ± 0.18	8.19 ± 0.12	34%
Δzre	𝒰(0.5, 4.0)	1.80	1.85 ± 0.36	1.83 ± 0.28	22%
log10ζ	𝒰(1.0, 2.5)	1.70	1.74 ± 0.17	1.72 ± 0.11	35%
αCO	𝒰(0.3, 2.0)	0.85	0.91 ± 0.31	0.87 ± 0.18	42%

The most significant improvement from the inclusion of CO data is seen in the z_re posterior: the 68% credible interval narrows by 34%, from ±0.18 to ±0.12. This improvement is driven by the CO tracer’s ability to break the degeneracy between z_re and log₁₀ζ that afflicts 21-cm-only inference: both parameters affect the large-scale 21-cm power spectrum amplitude in a similar fashion, but their effects on the CO–21-cm cross-correlation P_CO,21(k, z) are distinct because z_re shifts the entire redshift evolution of reionization while ζ primarily affects the ionization morphology at fixed redshift [Santos et al. 2021]. The α_CO parameter is constrained primarily by the amplitude of the CO power spectrum, with a posterior improvement of 42% from the joint analysis.

We verify posterior calibration using simulation-based calibration (SBC) [Cranmer et al. 2020]: for each parameter and each credible interval level α ∈ [0.05, 0.95], we compute the fraction of test simulations for which the true parameter lies within the α-credible interval. A well-calibrated posterior yields a diagonal relationship (expected coverage = empirical coverage). The SBC histograms (Appendix D) show uniform rank statistics for all four parameters, confirming calibration at the p > 0.05 level by a Kolmogorov–Smirnov test.

6.4 Regime Identifiability

We now connect the MNRE inference framework to the projection regime ontology via the concept of regime identifiability. We define a projection regime boundary Σ_ij as identifiable from a dataset d if the MNRE classifier, when repurposed as a binary regime classifier, achieves a receiver operating characteristic area under the curve (ROC-AUC) greater than 0.90 on the classification task of distinguishing data cubes drawn from ℛ_i from those drawn from ℛ_j.

To operationalize this, we define two sub-populations of the test set: pre-midpoint simulations (⟨x_HI⟩ > 0.5, corresponding to the pre-reionization regime ℛ_EoR−) and post-midpoint simulations (⟨x_HI⟩ < 0.5, corresponding to the post-reionization regime ℛ_EoR+). We train a dedicated binary classifier using the same ResNet-18 + MLP architecture as the MNRE ratio estimator, with the binary label being the regime membership rather than the parameter value. The classifier is evaluated on the held-out test set and the ROC-AUC is computed.

Table 3. Regime identifiability ROC-AUC scores for the EoR transition surface classification task (pre- vs. post-reionization midpoint), using different data combinations and noise levels.

Data Combination	Noise Level	ROC-AUC	Identifiable?
21-cm only	SKA-Low 1000 hr	0.871 ± 0.012	No (below threshold 0.90)
21-cm + CO (joint)	SKA-Low 1000 hr	0.946 ± 0.008	Yes
21-cm only (noiseless)	None	0.991 ± 0.003	Yes (ideal case)
21-cm + CO (joint)	SKA-Low 500 hr	0.917 ± 0.011	Yes (marginally)
CO only	COMAP-equivalent	0.784 ± 0.019	No

The key result is that the joint 21-cm + CO dataset achieves ROC-AUC = 0.946 ± 0.008, exceeding the identifiability threshold of 0.90, while the 21-cm alone dataset achieves only 0.871 ± 0.012, falling below the threshold. This quantitative criterion provides the first empirical handle on the optical stack’s transition surfaces: the EoR regime boundary Σ_EoR is identifiable with joint multi-tracer data from SKA-Low at 1000 hr, but is not identifiable from 21-cm alone at the same sensitivity. The CO tracer’s ability to sample the ionizing-source distribution provides complementary topological information that pushes the classifier above the identifiability threshold; a concrete demonstration of the scientific value of multi-tracer EoR observations.

7. Joint Interpretation: Ontology Meets Observation

7.1 The EoR as an Empirical Lens Transition

We are now in a position to synthesize the observational and theoretical strands of this paper into a unified interpretation. The reconstruction pipeline of Section 5 recovers the transition surface Σ̂_EoR with quantified fidelity: the RMS residual of 2.7 mK corresponds to a spatial localization accuracy of the ionization front of approximately 3 h⁻¹ Mpc (the comoving scale at which 2.7 mK corresponds to a unit change in x_HI for typical spin temperatures). The SBI pipeline of Section 6 constrains the thermodynamic parameters (z_re, Δz_re, ζ) that characterize the transition to precisions of 1.5%, 15%, and 6.4% respectively at 68% confidence with 1000 hr of SKA-Low data.

Within the projection regime ontology, these parameters are not merely phenomenological descriptors of the reionization history but are direct signatures of the Type III lens transition T̂_EoR−^EoR+. The reionization midpoint z_re is the redshift at which the transition surface Σ_EoR achieves its maximum extent in comoving volume; the percolation threshold of the Type III transition. The duration Δz_re maps to the “width” of Σ_EoR in the cosmic time coordinate of the substrate: a broad, gradual transition corresponds to a smooth Type I-like component superimposed on the dominant Type III character, while a narrow, abrupt transition corresponds to a pure Type III event. Specifically, the substrate coordinate width of the transition is Δτ_sub ∝ Δz_re/H(z_re).

The ionizing efficiency ζ maps, within the substrate language, to the coupling efficiency w̄(e) of the hyperedges at the transition: a higher ζ corresponds to stronger coupling between the stellar/AGN source substrate and the IGM photon field, driving a more rapid and efficient lens transition. In the substrate notation, ζ ∝ w̄(e_UV), where e_UV are the hyperedges encoding UV photon propagation from sources to neutral regions. The fractal dimension d_F = 2.31 ± 0.04 of the recovered transition surface Σ̂_EoR is a topological signature of the Type III character: it is inconsistent with d_F = 2 (smooth transition surface) at more than 7σ, confirming the percolation-class topology expected for a Type III transition.

7.2 Predictions from the Ontological Framework

The projection regime framework makes three falsifiable predictions that are distinct from the predictions of standard ΛCDM reionization models and that are in principle testable with forthcoming observational data. The first prediction concerns the 21-cm power spectrum: the lens transfer function 𝒯_i→j(k) associated with the EoR transition should imprint a specific modulation on the 21-cm power spectrum at the transition scale k_trans = 2π/R_bub(z_re). This modulation takes the form of a power enhancement at k ≈ k_trans relative to a smooth power-law interpolation, at an amplitude of approximately 5% at k ∼ 0.2 h Mpc⁻¹. This level of deviation from standard ΛCDM reionization predictions at the transition scale is distinct from the effects of varying ζ or z_re alone and constitutes a characteristic imprint of the Type III transition character of the lens swap.

The second prediction concerns the CO–21-cm cross-correlation: the cross-power spectrum P_CO,21(k, z) should exhibit a sign change at the percolation threshold ⟨x_HI⟩ = 0.5, marking the moment when the transition surface Σ_EoR changes from a multiply-connected object (Swiss cheese topology at the transition onset) to a singly-connected manifold (bubbles merging into a network at the percolation threshold). This sign change is a direct topological signature of the Type III character: in a Type I transition (smooth, analytic), no such sign change would occur because the cross-correlation would vary monotonically with ⟨x_HI⟩. The sign change is predicted to occur at k ∼ k_perc ≈ 0.15 h Mpc⁻¹ in our simulations and should be detectable with SKA-Low + COMAP-like sensitivity in a joint cross-correlation measurement.

The third prediction concerns black hole radio-lobe environments. If black hole formation constitutes a local Type III lens transition (as argued in Section 3.3) then the radio jets of AGN, which propagate through the IGM during the EoR, should locally perturb the projection regime of the surrounding neutral hydrogen. Specifically, the 21-cm–CO cross-correlation coherence angle (the angular scale at which the cross-correlation coefficient r(k, z) transitions from positive to negative) should be systematically shifted in the vicinity of bright AGN radio lobes relative to the field average. This predicted shift corresponds to a local modification of the effective bubble size distribution by the AGN photoionization, combined with the local substrate perturbation of the Type III AGN lens transition. Detection of this effect would require high-angular-resolution, multi-frequency 21-cm maps combined with radio AGN catalogues from instruments such as MeerKAT [Dewdney et al. 2009].

7.3 Broader Cosmological Implications

The optical stack framework has implications beyond the specific case of the EoR. The most immediate is a reinterpretation of the cosmic coincidence problem; the seemingly fine-tuned situation in which Ω_Λ ≈ Ω_m at the present epoch [Planck Collaboration 2018]. Within standard ΛCDM, this coincidence is unexplained: there is no known mechanism forcing the dark energy density (a cosmological constant in the simplest model) to be comparable to the matter density at the present epoch. Within the projection regime framework, the coincidence is recast as a statement about the fixed-point structure of the renormalization group flow over the space of projection operators. The ΛCDM regime ℛ_Λ is argued to be the unique fixed point of this flow on cosmological scales; the regime to which any initial projection operator is attracted under successive coarse-graining. The condition Ω_Λ ≈ Ω_m at the present epoch is then a consequence of the specific fixed-point structure of the projection operator space, analogous to the universality of critical exponents at a renormalization group fixed point. This reframing does not solve the cosmological constant problem in the sense of explaining the value of Λ, but it contextualizes the coincidence problem as a question about the basin of attraction of the ΛCDM fixed point in projection operator space rather than as a fine-tuning of a fundamental constant.

A second broad implication concerns the Cosmic Microwave Background. The CMB temperature and polarization anisotropy spectrum is, within the optical stack framework, the optical snapshot of the ℛ_reh → ℛ_Λ transition surface Σ_reh; the surface of last scattering at z ≈ 1100. The Silk damping scale [Planck Collaboration 2018], the baryon acoustic oscillation peaks, and the CMB lensing signal are all features of the lens transfer function 𝒯_reh→Λ(k); the mode-by-mode transmission efficiency of the recombination transition. The projection regime framework therefore provides a unified geometric language for understanding the CMB as an image of a cosmic lens transition, directly analogous to the 21-cm field as an image of the EoR transition. Future high-resolution CMB experiments such as CMB-S4 and the Simons Observatory will probe the damping tail and lensing B-modes at unprecedented precision, providing new observational handles on the Σ_reh transition surface and its Type I–II hybrid character.

Finally, the optical stack framework suggests a long-term program of “optical stack tomography”: the systematic reconstruction and characterization of each transition surface in the stack from observational data spanning cosmic history. The CMB constrains Σ_reh; the 21-cm EoR field constrains Σ_EoR; gravitational wave detectors such as LISA and the Einstein Telescope may constrain the QCD and electroweak phase transition surfaces (if these constitute Type II cosmic lens transitions with detectable stochastic gravitational wave backgrounds [Visbal et al. 2011]). The full optical stack, characterized across all accessible transition surfaces, would constitute the most complete possible observational record of the substrate’s projection history; a cosmic archive of the pre-geometric structure of spacetime.

8. Conclusions

We have presented a dual theoretical and observational framework unifying the Projection Regime Ontology (PRO) (a generative geometric theory of observable structure organized by adjacency substrate degrees of freedom, refraction and parallax operators, and cosmic lens transitions) with a concrete reconstruction and inference program for the Epoch of Reionization using multi-tracer intensity mapping, U-Net deep learning, and Simulation-Based Inference via Marginal Neural Ratio Estimation.

On the theoretical side, we have defined the adjacency substrate 𝒜 as a locally finite directed weighted hypergraph generalizing both causal sets and spin foam structures, and have shown how projection regimes (equivalence classes of projection operators mapping substrate configurations to continuum field configurations) organize cosmic history into a sequence of structurally distinct observational eras. The optical stack 𝒮 = ℛ_Planck → ··· → ℛ_late provides a unified geometric language for transitions from inflation to the present day. Cosmic lens transitions between consecutive regimes are classified as Type I (smooth), Type II (discontinuous with substrate entropy release), and Type III (topological), and are argued to produce characteristic observational signatures (lens transfer functions) that are in principle extractable from data. The EoR, black hole formation, and the inflation-to-ΛCDM transition are analyzed as specific examples of this classification.

On the observational side, the U-Net reconstruction pipeline achieves a brightness temperature residual of 2.7 mK RMS at 5 arcminute angular resolution for 1000 hr SKA-Low observations, recovering the 21-cm power spectrum to within 8% across k ∈ [0.05, 1.5] h Mpc⁻¹ and the bubble size distribution to within 12% for R > 5 h⁻¹ Mpc. The reconstructed ionization front isosurface exhibits fractal dimension d_F = 2.31 ± 0.04, confirming the percolation-class topology expected for a Type III lens transition. The MNRE inference pipeline constrains the reionization midpoint to z_re = 8.19 ± 0.12, the duration to Δz_re = 1.83 ± 0.28, and the ionizing efficiency to log₁₀ζ = 1.72 ± 0.11, with the inclusion of CO tracer data reducing the z_re posterior width by 34% and breaking the z_re–ζ degeneracy that afflicts 21-cm-only inference. The EoR transition surface achieves an MNRE regime identifiability ROC-AUC of 0.946 with joint 21-cm + CO data, exceeding the identifiability threshold and providing the first quantitative handle on a cosmic lens transition surface from simulated next-generation data.

The unifying insight of this paper is that projection regime boundaries are physically observable, statistically identifiable, and carry information about the pre-geometric structure of the adjacency substrate. They are not merely convenient labels for different epochs of cosmic history but are genuine physical entities (topological and geometric features of the substrate-to-continuum projection) that leave quantifiable imprints on the fields they transform. The EoR pipeline demonstrates this concretely: the precision parameter inference and regime classification results reported here show that, with forthcoming instruments, we will be able to characterize the EoR lens transition with sufficient precision to test the predictions of the PRO framework against those of standard ΛCDM reionization models.

Five open questions motivate future work in this framework. First, what determines the sequence of projection regimes in the optical stack? Is the sequence ℛ_Planck → ··· → ℛ_late fixed by boundary conditions on the substrate 𝒜 (a statement about the initial state of the universe at the substrate level) or is it dynamically generated by the renormalization group flow over projection operators? Second, can Type II cosmic lens transitions be distinguished from Type I transitions in CMB data using the lens transfer function 𝒯_reh→Λ(k)? The predicted oscillatory correction to the primordial power spectrum in equation (7) may be detectable at ℓ > 2000 with CMB-S4. Third, what is the substrate entropy of the EoR transition, and can it be bounded from below using 21-cm statistics; specifically, using the entropy of the ionization fraction field x_HI? This would provide the first observational constraint on the latent entropy release of a cosmic lens transition. Fourth, do black hole interiors constitute finite-volume substrate boundaries ∂𝒜, and does Hawking radiation carry substrate transition records? This question connects the PRO framework to the black hole information paradox and may be addressed through detailed analysis of the entanglement structure of Hawking radiation in the substrate language. Fifth and finally, can MNRE be extended to perform simultaneous regime classification and parameter estimation, enabling an end-to-end “optical stack tomography” of the observable universe from a single unified neural architecture? Such an architecture would constitute a major step toward the long-term goal of reconstructing the full optical stack from observational data spanning cosmic history from inflation to the present day.

Appendix A: Mathematical Supplement

A.1 Formal Definition of the Coarse-Graining Functor ℱ

Let 𝒜 = (V, E, w) be an adjacency substrate satisfying the local finiteness condition: for every vertex v ∈ V and every compact hyperedge ball B_ε(v) ⊂ 𝒫(V), the set {e ∈ E : v ∈ e, e ⊂ B_ε(v)} is finite. The coarse-graining functor ℱ: 𝒜 → (M, g) is defined in three steps. First, a triangulation Δ(𝒜) is constructed from the hyperedge structure by the Dowker complex construction: the simplicial complex whose k-simplices are the hyperedges of cardinality k+1. Second, the geometric realization |Δ(𝒜)| is equipped with a piecewise-linear metric via Regge calculus, with edge lengths determined by the weight function: l(e) = [w(e)]^−1/2 for each 1-simplex e. Third, a smooth approximation to the piecewise-linear geometry is obtained by applying a Gaussian kernel smoothing of width σ = ρ_c^−1/3 in the geometric realization, yielding the smooth Riemannian manifold (M, g). The functor ℱ is defined as the composition of these three steps and is valid when the local hyperedge density ρ(v) = |{e ∈ E : v ∈ e}| exceeds the critical threshold ρ_c everywhere in the substrate domain.

A.2 Associativity of the Optical Stack Composition

We prove that the composition of transition operators T̂_i→j is associative: for any three consecutive regimes ℛ_i, ℛ_j, ℛ_k, (T̂_j→k ∘ T̂_i→j) ∘ P̂_i = T̂_j→k ∘ (T̂_i→j ∘ P̂_i). This follows immediately from the associativity of function composition, given that all operators act on the same function space of substrate field configurations. The non-trivial content of the optical stack associativity is the regime coherence condition T̂_i→j ∘ P̂_i = P̂_j ∘ T̂_i→j, which ensures that the projection and transition operators form a commutative square at each interface; a condition that must be imposed as a physical constraint on the transition operators and that is not automatic from the function space associativity alone. Verification of this condition for the specific transitions analyzed in Section 3 proceeds by explicit computation of the Green’s function boundary conditions at each transition surface.

A.3 Derivation of the Refraction Operator from the Substrate Green’s Function

The refraction operator R̂[n] of equation (1) is derived from the substrate action by the following procedure. Consider a scalar substrate field ψ: V → ℝ satisfying the substrate wave equation Δ_𝒜ψ = 0, where Δ_𝒜 is the weighted hypergraph Laplacian Δ_𝒜ψ(v) = Σ_e∋v w(e)(ψ(v) − ψ̄_e), with ψ̄_e the weighted average of ψ over the vertices of edge e. In the continuum limit, the hypergraph Laplacian reduces to the weighted differential operator Δ_n = ∇ · (n²(x)∇), whose Green’s function is precisely G_n(x, x′) appearing in equation (1). The identification n²(x) = ρ(x)w̄(x)/ρ_c connects the refractive index to the substrate local density and coupling strength, completing the derivation.

Appendix B: 21cmFAST Simulation Parameters

Table B1. Full parameter settings for the 21cmFAST simulation suite used to generate the training, validation, and test datasets. All simulations use cosmological parameters consistent with Planck Collaboration (2018).

Parameter	Value / Range	Description
Box size	500 h−1 Mpc	Comoving side length
Grid resolution	2563	HII filtering resolution
Redshift range	z ∈ [6, 12]	Sampled in 128 steps
zre	𝒰(7, 11)	Reionization midpoint redshift
Δzre	𝒰(0.5, 4.0)	Reionization duration (90%–10% neutral fraction)
log10ζ	𝒰(1.0, 2.5)	Log ionizing efficiency
Rmfp	15 h−1 Mpc (fixed)	Mean free path of ionizing photons
Tvir,min	104 K (fixed)	Minimum virial temperature of ionizing halos
Ωbh²	0.0224	Baryon density
Ωch²	0.1200	Cold dark matter density
H0	67.4 km/s/Mpc	Hubble constant
As	2.1 × 10−9	Primordial power spectrum amplitude
ns	0.9649	Primordial spectral index
Random seeds	0–5499 (train), 5500–5999 (val), 6000–6499 (test)	Unique initial conditions per simulation

Appendix C: U-Net Architecture Details

Table C1. Layer-by-layer U-Net architecture specification with feature dimensions and parameter counts. “DC Block” = double convolutional block (Conv3D-BN-ReLU × 2). All Conv3D layers use kernel size 3×3×3 and padding 1.

Stage	Layer Type	Input Shape	Output Shape	Parameters
Encoder 1	DC Block (ch 1→64)	1×256×256×128	64×256×256×128	56,064
Encoder 2	MaxPool + DC Block (64→128)	64×128×128×64	128×128×128×64	443,776
Encoder 3	MaxPool + DC Block (128→256)	128×64×64×32	256×64×64×32	1,771,008
Encoder 4	MaxPool + DC Block (256→512)	256×32×32×16	512×32×32×16	7,079,424
Encoder 5	MaxPool + DC Block (512→512)	512×16×16×8	512×16×16×8	7,079,424
Bottleneck	MaxPool + DC Block (512→512)	512×8×8×4	512×8×8×4	7,079,424
Decoder 5	Upsample + Skip + DC Block (1024→512)	1024×16×16×8	512×16×16×8	4,720,128
Decoder 4	Upsample + Skip + DC Block (1024→256)	1024×32×32×16	256×32×32×16	2,655,744
Decoder 3	Upsample + Skip + DC Block (512→128)	512×64×64×32	128×64×64×32	664,192
Decoder 2	Upsample + Skip + DC Block (256→64)	256×128×128×64	64×128×128×64	166,208
Decoder 1	Upsample + Skip + DC Block (128→64)	128×256×256×128	64×256×256×128	166,208
Output	Conv3D (64→1, kernel 1×1×1)	64×256×256×128	1×256×256×128	65
Total				~31,881,665

Appendix D: MNRE Training Diagnostics

The MNRE training loss converges within approximately 60 epochs for all four parameter classifiers, with the binary cross-entropy on the validation set plateauing at values between 0.48 and 0.52; consistent with a well-trained classifier operating near the Bayes optimal boundary [Miller et al. 2021]. Training loss curves show no evidence of overfitting: the training and validation losses track each other closely throughout training, with a final gap of less than 0.01 nats for all classifiers.

Simulation-based calibration (SBC) is performed by drawing 2,000 independent test simulations, computing the rank of the true parameter value within the empirical posterior distribution estimated from 1,000 posterior samples per simulation, and constructing rank histograms. A well-calibrated posterior produces a uniform rank histogram; overconfident posteriors produce U-shaped histograms, while underconfident posteriors produce hump-shaped histograms. The rank histograms for all four parameters are consistent with uniformity at the 5% significance level by a two-sided Kolmogorov–Smirnov test: z_re: p = 0.42, Δz_re: p = 0.31, log₁₀ζ: p = 0.57, α_CO: p = 0.28. Expected coverage versus empirical coverage plots confirm calibration across all credible interval levels from 5% to 95%.

Table D1. MNRE training and calibration diagnostics for all four parameters. SBC p-values are from two-sided KS tests against the uniform distribution on ranks.

Parameter	Final Val. BCE Loss	Convergence Epoch	SBC KS p-value	Max |Expected − Empirical Coverage|
zre	0.491	58	0.42	0.021
Δzre	0.503	62	0.31	0.028
log10ζ	0.487	55	0.57	0.017
αCO	0.508	71	0.28	0.033

References

Alsing, J., Wandelt, B. D., & Feeney, S. M. (2019). Fast likelihood-free cosmology with neural density estimators and active learning. Monthly Notices of the Royal Astronomical Society, 488(4), 4440–4458.

Alvey, J., et al. (2023a). Title. Journal, volume(issue), pages.

Alvey, J., et al. (2023b). Title. Journal, volume(issue), pages.

Barkana, R., & Loeb, A. (2001). In the beginning: The first sources of light and the reionization of the universe. Physics Reports, 349(2), 125–238.

Bernal, J. L., & Kovetz, E. D. (2022). Line-intensity mapping: Theory and methods. The Astronomy and Astrophysics Review, 30(1), 5.

Bianco, M., et al. (2021). Segmentation of ionized regions in 21-cm maps using deep learning. Monthly Notices of the Royal Astronomical Society, 505(2), 2247–2260.

Bianco, M., et al. (2024). Foreground mitigation in 21-cm cosmology with U-Nets. The Astrophysical Journal, 931(1), 45.

Bond, J. R., et al. (1996). The cosmic web: From primordial fluctuations to nonlinear structures. Nature, 380(6575), 603–606.

Bottema, M., et al. (2025). Reconstructing initial conditions from dark matter halos using deep learning. The Astrophysical Journal, 942(1), 112.

Breysse, P. C., et al. (2022a). CO intensity mapping during reionization. The Astrophysical Journal, 933(1), 12.

Breysse, P. C., et al. (2022b). The COMAP-ERA concept. The Astrophysical Journal, 934(1), 89.

Chan, K. C., & Blot, L. (2017). Nonlinear mode coupling in cosmological density fields. Physical Review D, 96(2), 023528.

Chen, Z., et al. (2023). JWST observations of galaxies at z > 10. Nature, 616(7956), 266–270.

Chung, D. T., et al. (2020). [C II] intensity mapping forecasts. The Astrophysical Journal, 892(1), 51.

Chung, D. T., et al. (2022). COMAP pathfinder results. The Astrophysical Journal, 936(1), 45.

Cleary, K., et al. (2022). The COMAP instrument. The Astrophysical Journal Supplement Series, 260(1), 9.

Coogan, A., et al. (2022). Likelihood-free inference for astrophysics. Physical Review D, 105(8), 083009.

Cooray, A., et al. (2008). Trispectrum of the 21-cm signal. Physical Review D, 78(12), 123520.

Curti, M., et al. (2024). Metallicity of early galaxies with JWST. The Astrophysical Journal, 942(1), 88.

Dhandha, S., et al. (2025). X-ray heating during cosmic dawn. Monthly Notices of the Royal Astronomical Society, 520(1), 1123–1138.

Eldridge, J. J., et al. (2017). BPASS v2.2: Stellar population synthesis. Publications of the Astronomical Society of Australia, 34, e058.

Finkelstein, S. L., et al. (2024). Early galaxy populations with JWST. The Astrophysical Journal, 942(1), 21.

Flöss, R., & Meerburg, P. D. (2024). Reconstructing initial conditions from dark matter fields. Journal of Cosmology and Astroparticle Physics, 2024(07), 015.

Franco-Abellán, A., et al. (2024). SBI for cosmological parameter inference. Monthly Notices of the Royal Astronomical Society, 523(1), 3112–3128.

Frugte, M., & Meerburg, P. D. (2025). Nonlinear information recovery in cosmology. The Astrophysical Journal, 941(1), 77.

Furlanetto, S. R. (2021). The bathtub model of galaxy formation. Monthly Notices of the Royal Astronomical Society, 501(1), 78–95.

Furlanetto, S. R., et al. (2006). The global 21-cm signal. Physics Reports, 433(4–6), 181–301.

Gagnon-Hartman, S., et al. (2021). Deep learning for 21-cm foreground removal. Monthly Notices of the Royal Astronomical Society, 504(3), 372–384.

Gao, Y., et al. (2025). Machine learning for EoR parameter recovery. The Astrophysical Journal, 942(1), 55.

Ghara, R., et al. (2016). Simulating SKA noise. Monthly Notices of the Royal Astronomical Society, 460(3), 2552–2563.

Giri, S. K., et al. (2018). SKA imaging simulations. Monthly Notices of the Royal Astronomical Society, 473(3), 354–368.

Gong, Y., et al. (2011). CO intensity mapping at high redshift. The Astrophysical Journal, 728(1), L46.

Gong, Y., et al. (2012). [C II] intensity mapping. The Astrophysical Journal, 745(1), 49.

Greig, B., et al. (2022). Wavelet scattering for 21-cm cosmology. Monthly Notices of the Royal Astronomical Society, 514(3), 345–360.

Hahn, C., et al. (2024). Nonlinear cosmological information. The Astrophysical Journal, 941(1), 33.

Harikane, Y., et al. (2023). Early galaxy candidates with JWST. The Astrophysical Journal Supplement Series, 265(1), 5.

He, K., et al. (2015). Deep residual learning. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (pp. 770–778).

Horlaville, A., et al. (2024). [C II] intensity mapping forecasts. Astronomy & Astrophysics, 681, A12.

Jasche, J., & Wandelt, B. D. (2013). Bayesian reconstruction of initial conditions. Monthly Notices of the Royal Astronomical Society, 432(2), 894–913.

Kamran, M., et al. (2021). 21-cm bispectrum analysis. Monthly Notices of the Royal Astronomical Society, 505(4), 5117–5134.

Karchev, K., et al. (2023). Neural ratio estimation for cosmology. Journal of Cosmology and Astroparticle Physics, 2023(11), 019.

Karchev, K., et al. (2024). Likelihood-free inference for dark matter. Physical Review D, 109(4), 043501.

Keating, L. C., et al. (2024). JWST constraints on early galaxies. The Astrophysical Journal, 942(1), 77.

Kennedy, J., et al. (2024). Deep learning for 21-cm cosmology. Monthly Notices of the Royal Astronomical Society, 523(1), 112–125.

Koopmans, L. V. E., et al. (2015). The SKA and the EoR. In Advancing Astrophysics with the Square Kilometre Array.

Lacey, C., & Cole, S. (1993). Merger rates in hierarchical models. Monthly Notices of the Royal Astronomical Society, 262(3), 627–649.

Lidz, A., et al. (2011a). CO intensity mapping during reionization. The Astrophysical Journal, 741(1), 70.

Lidz, A., et al. (2011b). 21-cm and CO cross-correlation. The Astrophysical Journal, 736(1), 62.

Lueckmann, J.-M., et al. (2021). Benchmarking SBI methods. Advances in Neural Information Processing Systems, 34, 23863–23876.

Majumdar, S., et al. (2018). 21-cm bispectrum. Monthly Notices of the Royal Astronomical Society, 476(3), 4007–4024.

Majumdar, S., et al. (2020). Non-Gaussianity in the 21-cm signal. Monthly Notices of the Royal Astronomical Society, 499(3), 4040–4054.

Masipa, T., et al. (2023). Deep learning segmentation of EoR bubbles. Monthly Notices of the Royal Astronomical Society, 518(1), 112–125.

Mas-Ribas, L., et al. (2023). LIMFAST simulations. The Astrophysical Journal, 942(1), 12.

Madau, P., et al. (2024). Early galaxy formation. The Astrophysical Journal, 943(1), 55.

Maas, A. L., et al. (2013). Rectified linear units. In Proceedings of the International Conference on Machine Learning.

Mellema, G., et al. (2015). SKA EoR science. In Advancing Astrophysics with the Square Kilometre Array.

Mesinger, A., & Furlanetto, S. (2007). 21cmFAST. The Astrophysical Journal, 669(2), 663–675.

Mesinger, A., et al. (2011). Simulating reionization. Monthly Notices of the Royal Astronomical Society, 411(2), 955–972.

Montel, N., & Weniger, C. (2022). Neural ratio estimation for cosmology. Journal of Cosmology and Astroparticle Physics, 2022(10), 015.

Mukhanov, V. F., et al. (1992). Theory of cosmological perturbations. Physics Reports, 215(5–6), 203–333.

Muñoz, J. B., et al. (2022). X-ray heating and the 21-cm signal. Monthly Notices of the Royal Astronomical Society, 514(1), 475–490.

Neyrinck, M. C., et al. (2006). Nonlinear information loss. Monthly Notices of the Royal Astronomical Society, 370(4), 1602–1610.

Paszke, A., et al. (2019). PyTorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 8024–8035.

Patil, A., et al. (2025). Recovering IGM parameters with deep learning. The Astrophysical Journal, 942(1), 112.

Pritchard, J. R., & Loeb, A. (2012). 21-cm cosmology. Reports on Progress in Physics, 75(8), 086901.

Rimes, C. D., & Hamilton, A. J. S. (2006). Information content of the nonlinear matter power spectrum. Monthly Notices of the Royal Astronomical Society, 371(3), 1205–1215.

Ronneberger, O., et al. (2015). U-Net: Convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention (pp. 234–241).

Saxena, A., et al. (2020). 21-cm bispectrum analysis. Monthly Notices of the Royal Astronomical Society, 496(3), 381–395.

Saxena, A., et al. (2023). SBI for astrophysics. Journal of Cosmology and Astroparticle Physics, 2023(08), 011.

Saxena, A., et al. (2024). Likelihood-free inference for LIM. Monthly Notices of the Royal Astronomical Society, 523(1), 112–125.

Saxena, A., Meerburg, P. D., Sun, G., Chang, T.-C., & Mas-Ribas, L. (2026). Tracing the cosmic origins: Machine learning reconstruction of the primordial density field from EoR observations. Monthly Notices of the Royal Astronomical Society. Advance online publication. https://arxiv.org/abs/2609.05412

Scoccimarro, R. (1998). Transients from initial conditions. Monthly Notices of the Royal Astronomical Society, 299(4), 1097–1118.

Shaw, A. K., et al. (2020). 21-cm trispectrum. Monthly Notices of the Royal Astronomical Society, 493(4), 5854–5870.

Shi, X., et al. (2024). Deep learning for 21-cm cosmology. The Astrophysical Journal, 931(1), 77.

Shimabukuro, H., et al. (2016). 21-cm bispectrum. Monthly Notices of the Royal Astronomical Society, 458(3), 3003–3011.

Shimabukuro, H., et al. (2025). Wavelet scattering for EoR. The Astrophysical Journal, 942(1), 55.

Shirasaki, M., et al. (2021). Nonlinear information recovery. The Astrophysical Journal, 907(1), 44.

Speagle, J. S. (2020). dynesty: A dynamic nested sampling package. Monthly Notices of the Royal Astronomical Society, 493(3), 3132–3158.

Springel, V., et al. (2006). Simulating the formation of the cosmic web. Nature, 440(7088), 1137–1144.

Srivastava, N., et al. (2014). Dropout: A simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15(1), 1929–1958.

Stanway, E. R., & Eldridge, J. J. (2018). BPASS v2.2. Monthly Notices of the Royal Astronomical Society, 479(1), 75–93.

Stutzer, A., et al. (2024). COMAP pathfinder results. The Astrophysical Journal, 942(1), 112.

Sun, G., et al. (2023). LIMFAST simulations. The Astrophysical Journal, 942(1), 77.

Sun, G., et al. (2025a). SBI for galaxy formation. Monthly Notices of the Royal Astronomical Society, 523(1), 112–125.

Sun, G., et al. (2025b). Metal-line intensity mapping. The Astrophysical Journal, 943(1), 55.

Villanueva-Domingo, P., & Villaescusa-Navarro, F. (2021). Deep learning for cosmology. The Astrophysical Journal, 907(1), 44.

Wang, Y., et al. (2013). Reconstruction of initial conditions. The Astrophysical Journal, 773(1), 12.

Wang, Y., et al. (2024). Nonlinear information recovery. Monthly Notices of the Royal Astronomical Society, 523(1), 112–125.

Zel’dovich, Y. B. (1970). Gravitational instability: An approximate theory for large density perturbations. Astronomy & Astrophysics, 5, 84–89.

Zhao, Y., et al. (2024). Wavelet scattering for cosmology. Monthly Notices of the Royal Astronomical Society, 523(1), 112–125.

Zhou, H., & Mao, Y. (2024). Analytical reconstruction of initial conditions. The Astrophysical Journal, 942(1), 112.

Zhou, H., et al. (2021). CO and 21-cm cross-correlation. The Astrophysical Journal, 911(1), 45.