We report systematic issues commonly observed in the Hydrogen Epoch of Reionization Array observations that cause abrupt changes in power over time, resulting in temporal discontinuities. To identify these effects, we applied the Temporal Discontinuity Index and Spectral Discontinuity Index metrics, which complement existing diagnostic tools by detecting discontinuities and revealing their potential relationships to bandpass and power-related issues. Analysis of 30 data sets, each corresponding to a single observing night, shows that such systematics appear consistently in many antennas and across multiple nights. Based on these metrics, we identify three main types of discontinuity-related issues: periodic broad-band discontinuities linked to Inter-Integrated Circuit (I2C) polling, broad-band discontinuities associated with degraded post-amplifier modules (PAMs), and sudden power changes related to the digital backend of the antenna system. The second type, associated with degraded PAMs, occurs less frequently within a single night compared to the periodic polling effect, though both are broadband. The third type affects the largest number of antennas and is most prevalent in the H6C data. The width of the discontinuities related to the digital backend corresponds to the frequency range handled by a single correlator box, scaled by the number of boxes involved. Although 30 data sets were examined, this paper presents a focused analysis of four representative nights to highlight the main discontinuity issues. We also present correlations between flagging due to discontinuities and flagging caused by other failure modes to investigate how these issues are interrelated.

A crucial physical quantity in determining the 21-cm signal during cosmic dawn is the inhomogeneous background of Ly alpha photons originating from the first galaxies. As these photons travel through the intergalactic medium (IGM), their scattering cross section is often approximated as a delta function at resonance due to computational cost. That is, photons with emitted wavelengths between Ly alpha and Ly beta are assumed to travel in straight lines until they redshift into the Ly alpha resonance. However, due to the damping wing in the Ly alpha cross section, this approximation fails as the frequency of the photon approaches the resonant frequency, resulting in multiple scatterings events that could be separated by non-negligible distances. These multiple scattering events effectively modify the intrinsic Ly alpha emissivity from galaxies. Some previous works studied this effect of Ly alpha multiple scattering by running computationally heavy radiative-transfer simulations. However, robustly interpreting the cosmic 21 cm signal requires exploring a large parameter space of astrophysical uncertainties, motivating more computationally efficient approaches. Here we incorporate Ly alpha multiple scatterings in the public, seminumerical simulation 21cmFAST. To do so, we employ Monte Carlo simulations to study the trajectories of Ly alpha photons on different scales. We find that the distance distributions of Ly alpha photons with respect to the absorption point can be modeled as analytical functions that are governed by a single parameter. Upon implementing the distance distributions in 21cmFAST, we find that the multiple scattering effect is important (about 50% difference in the 21-cm power spectrum) only at high redshifts before the spin temperature is fully coupled to the kinetic temperature. Furthermore, we find that Ly alpha multiple scattering does not enhance Ly alpha heating, and that the combined effect is negligible, especially under realistic x-ray heating scenarios.

Synergies with other instruments will be essential in making, verifying, and interpreting a detection of the cosmic 21-cm signal from the Epoch of Reionization (EoR) and Cosmic Dawn (CD) with the Square Kilometer Array (SKA) telescope. Such synergies can (i) provide prior information about galaxies and the intergalactic medium (IGM) during the EoR/CD; (ii) pave the road to a first 21cm detection by mitigating foregrounds and systematics through cross-correlations; and (iii) give complimentary physical insights into the galaxy – IGM connection. Here we review the current state of synergies and discuss what observations will best compliment SKA-low EoR/CD observations.

The redshifted 21-cm line is a promising probe of the Epoch of Reionization, during which the first generation of stars ionized the InterGalactic Medium (IGM). We aimed to constrain the IGM neutral Hydrogen fraction and spin temperature via redshifted 21-cm absorption against high-redshift radio sources. We analysed an 18-hour observation of the radio-loud quasar PSO J352.4034-15.3373 (z = 5.84) in the 203-222.5 MHz band. We obtained a continuum image with a 0.66 mJy/beam rms noise and a spectrum with 3.6 mJy/beam rms noise per 390 kHz-wide channel. We fit the quasar spectrum with a power-law continuum modified by an intervening 21-cm absorption, testing two scenarios: an island model, a residual cold neutral Hydrogen patch surviving at the end of reionization, and a global model, approximating the overall decline of the IGM HI fraction and spin temperature with redshift. Combining archival data spanning 150 MHz to 3.0 GHz with our measurements, we determined a quasar flux density of 86.8 ± 1.1 mJy at 200 MHz and a spectral index of -0.88 ± 0.02 across the full frequency range. We found no evidence for the 21-cm absorption from intervening neutral Hydrogen at 5.38<z<5.84. Assuming the island model, we set a 95

We report the first upper limits on the power spectrum of 21 cm fluctuations during the Epoch of Reionization and Cosmic Dawn from Phase II of the Hydrogen Epoch of Reionization Array (HERA) experiment. HERA Phase II constitutes several significant improvements in the signal chain compared to Phase I, most notably resulting in expanded frequency bandwidth, from 50–250 MHz. In these first upper limits, we investigate a small two-week subset of the available Phase II observations, with a focus on identifying new systematic characteristics of the instrument, and establishing an analysis pipeline to account for them. We report 2 σ upper limits in eight spectral bands, spanning 5.6 ≤ z ≤ 24.4 that are consistent with thermal noise at the 2 σ level for k ≳ 0.6–0.9 h Mpc ^−1 (band dependent). Our tightest limit during Cosmic Dawn ( z > 12) is 1.13 × 10 ^6 mK ^2 at ( k = 0.55 h Mpc ^−1 , z = 16.78), and during the EoR (5.5 < z < 12), it is 1.78 × 10 ^3 mK ^2 at ( k = 0.70 h Mpc ^−1 , z = 7.05). We find that mutual coupling has become our dominant systematic, leaking foreground power that strongly contaminates the low- k modes, resulting in the loss of modes from k = 0.35 to 0.55 compared to Phase I data.

Reconstructing the initial conditions (ICs) of the matter field for an observed volume gives a complete picture of that region's temporal evolution. IC reconstruction is common at low redshifts but largely unexplored during the epoch of reionization (EoR; z≳ 5), partly due to difficulties modeling inhomogeneous cosmic radiation fields. Yet the EoR spans half the observable Universe, offering unmatched potential for astrophysics and cosmology. Here we quantify how well upcoming galaxy and 21cm observations can constrain ICs during the EoR. We develop TILING (Tomographic Inference of Linear ICs via Network Grafting): a hybrid machine learning pipeline that first produces a point estimate of the ICs, which then improves a score-based diffusion network for generating posterior samples. We train TILING on mock galaxy maps at varying UV magnitude limits, plus corresponding 21cm maps at varying noise levels and foreground "wedge" contamination. For fiducial survey choices, we achieve accurate IC reconstruction (posterior mean cross-correlation coefficients >0.8 and power spectrum errors under a few percent) down to k≲ 0.2 cMpc^-1. While 21cm interferometry aids recovery of IC power spectra, most constraining power, especially for IC Fourier phases, comes from galaxy maps, underscoring the need for complementary observations when interpreting the 21cm signal. We show how TILING can recover 21cm power spectrum modes excised by foreground contamination. Our framework can also: (i) guide follow-up observations of sub-volumes of interest; (ii) reconstruct galaxy evolution and reionization morphology for specific volumes; and (iii) isolate the contribution to reionization from the vast majority of galaxies unobservable by optical/IR telescopes like JWST. Our code is publicly available.

The fraction of 'dark pixels' in the Ly alpha and other Lyman-series forests at z similar to 5-6 provides a powerful constraint on the end of the reionization process. Any spectral region showing transmission must be highly ionized, while dark regions could be ionized or neutral, thus the dark pixel fraction provides a (nearly) model independent upper limit to the volume-filling fraction of the neutral intergalactic medium, modulo choices in binning scale and dark pixel definition. Here, we provide updated measurements of the 3.3 comoving Mpc dark pixel fraction at z = 4.85-6.25 in the Ly alpha, Ly beta, and Ly gamma forests of 34 deep 5.8 less than or similar to z less than or similar to 6.6 quasar spectra from the (enlarged) XQR-30 sample. Using the negative pixel method to measure the dark pixel fraction, we derive fiducial 1 sigma upper limits on the volume-average neutral hydrogen fraction of < x(HI)><= {0.030 + 0.048, 0.095 + 0.037, 0.191 + 0.056, 0.199 + 0.087} at z = {5.481, 5.654, 5.831, 6.043} from the optimally sensitive combination of the Ly beta and Ly gamma forests. We further demonstrate an alternative method that treats the forest flux as a mixture of dark and transparent regions, where the latter are modelled using a physically motivated parametric form for the intrinsic opacity distribution. The resulting modeldependent upper limits on < x(HI)> are similar to those derived from our fiducial model-independent analysis. We confirm that the bulk of reionization must be finished at z > 6, while leaving room for an extended 'soft landing' to the reionization history down to z similar to 5.4 suggested by Ly alpha forest opacity fluctuations.

The Lyman α (Lyα) line from high-redshift galaxies is a powerful probe of the Epoch of Reionization (EoR). Neutral hydrogen in the intergalactic medium (IGM) can significantly attenuate the emergent Lyα line, even in the damping wing of the cross-section. However, interpreting this damping wing imprint relies on our prior knowledge of the spectrum that escapes from the galaxy and its environs into the IGM. This emergent spectrum is highly sensitive to the composition and geometry of the interstellar and circumgalactic media, and so exhibits a large galaxy to galaxy scatter. Characterizing this scatter is further complicated by non-trivial selection effects introduced by observational surveys. Here we build a flexible, empirical model for the emergent Lyα spectra. Our model characterizes the emergent Lyα luminosity, the velocity offset of the Lyα line with respect to the systemic redshift, and the Hα luminosity, with multivariate probability distributions conditioned on the UV magnitude. We constrain these distributions using z∼5-6 galaxy observations with VLT MUSE and JWST NIRCam, forward-modeling observational selection functions together with galaxy parameters. Our model results in Lyα equivalent width distributions that are a better match to (independent) Subaru observations than previous empirical models. The extended distributions of Lyα equivalent widths and velocity offsets we obtain could facilitate Lyα transmission during the early stages of the EoR. We also illustrate how our model can be used to identify GN-z11-like outliers, potentially originating from merging systems. We publish fitting functions and make our model publicly available.

We are witnessing a surge in observations of the cosmic dawn (CD) and epoch of reionisation (EoR), driving an increasing demand for fast and robust theoretical interpretation frameworks. In response, machine learning (ML), and emulation in particular, has emerged as a powerful approach to accelerate and enhance inference pipelines. In this work, we present 21cmEMUv3, an emulator trained on 21cmFASTv3 simulations that model both atomically and molecularly cooling galaxies. 21cmEMUv3 is conditioned on σ_8 and ten astrophysical parameters to produce seven summary observables: (i) the cylindrical 21cm power spectrum (PS), emulated for the first time at such high resolution and accuracy across a wide redshift range of z ∼ 6–30; (ii) the spherically-averaged 21cm PS; (iii) the mean neutral fraction of the intergalactic medium (IGM); (iv) the mean 21cm spin temperature; (v) the global 21cm signal; (vi) the ultraviolet (UV) luminosity functions (LFs); and (vii) the Thomson scattering optical depth. Notably, the cylindrical 21cm PS is emulated via score-based diffusion, while the remaining six summaries are emulated via long-short term memory (LSTM) networks, all achieving sub-percent median accuracy. We use the emulator to reinterpret current 21cm PS upper limits from HERA, for the first time using state-of-the-art hydrodynamical simulations to inform priors on star formation inside molecularly cooling galaxies. We find that our inferred soft-band X-ray luminosity per unit star formation rate is consistent with extrapolations of high-mass X-ray binaries to the low-metallicity regimes expected in the first galaxies, excluding values below 10^39.2 erg s^-1M^-1_⊙yr at 95% confidence. Finally, we produce forecasts for the detection of the cosmic 21cm PS with the Square Kilometre Array for different array configurations. The 21cmEMU package is publicly available.

State-of-the-art simulations of re-ionisation-era 21 cm signal have limited volumes, generally orders of magnitude smaller than observations. Consequently, the Fourier modes in common between simulation and observation have limited overlap, especially in cylindrical (2D) k-space that is natural for 21 cm interferometry. This makes sample variance (i.e. the deviation of the simulated sample from the population mean due to finite box size) a potential issue when interpreting upcoming 21 cm observations. Here, we introduce 21cmPSDenoiser, a score-based diffusion model that can be applied to a single, forward-modelled realisation of the 21 cm 2D power spectrum (PS), predicting the corresponding population mean on-the-fly during Bayesian inference. Individual samples of 2D Fourier amplitudes of wave modes relevant to current 21 cm observations can deviate from the mean by over 50% for 300 cmpc simulations, even when only considering stochasticity due to the sampling of Gaussian initial conditions. 21cmPSDenoiser reduces this deviation by an order of magnitude, outperforming current state-of-the-art sample variance mitigation techniques such as fixing and pairing by a factor of a few at almost no additional computational cost (similar to 2 s per PS). Unlike emulators, 21cmPSDenoiser is not tied to a particular model or simulator since its input is a (model-agnostic) realisation of the 2D 21 cm PS. Indeed, we confirm that 21cmPSDenoiser generalises to PSs produced with a different 21 cm simulator than those on which it was trained. To quantify the improvement in parameter recovery, we simulated a 21 cm PS detection by the Hydrogen Epoch of Reionization Arrays (HERA) and ran different inference pipelines corresponding to commonly used approximations. We find that using 21cmPSDenoiser in the inference pipeline outperforms other approaches, yielding an unbiased posterior that is 50% narrower in most inferred parameters.

The precise characterization and mitigation of systematic effects is one of the biggest roadblocks impeding the detection of the fluctuations of cosmological 21 cm signals. Missing data in radio cosmological experiments, often due to radio frequency interference (RFI), pose a particular challenge to power spectrum analysis as this could lead to the ringing of bright foreground modes in the Fourier space, heavily contaminating the cosmological signals. Here we show that the problem of missing data becomes even more arduous in the presence of systematic effects. Using a realistic numerical simulation, we demonstrate that partially flagged data combined with systematic effects can introduce significant foreground ringing. We show that such an effect can be mitigated through inpainting the missing data. We present a rigorous statistical framework that incorporates the process of inpainting missing data into a quadratic estimator of the 21 cm power spectrum. Under this framework, the uncertainties associated with our inpainting method and its impact on power spectrum statistics can be understood. These results are applied to the latest Phase II observations taken by the Hydrogen Epoch of Reionization Array, forming a crucial component in power spectrum analyses as we move toward detecting 21 cm signals in the ever more noisy RFI environment.

The cosmic 21-cm signal promises to revolutionize studies of the Epoch of Reionization (EoR). Radio interferometers are aiming for a preliminary, low signal-to-noise (S/N) detection of the 21-cm power spectrum. Cross-correlating 21-cm with galaxies will be especially useful in these efforts, providing both a sanity check for initial 21-cm detection claims and potentially increasing the S/N due to uncorrelated residual systematics. Here we self-consistently simulate large-scale (1 Gpc(3)) galaxy and 21-cm fields, computing their cross-power spectra for various choices of instruments and survey properties. We use 1080 h observations with SKA-low AA* and HERA-350 as our benchmark 21-cm observations. We create mock Lyman-alpha narrow-band, slitless and slit spectroscopic surveys, using benchmarks from instruments such as Subaru HyperSupremeCam, Roman grism, VLT MOONS, ELT MOSAIC, and JWST NIRCam. We forecast the resulting S/N of the galaxy-21-cm cross-power spectrum, varying for each pair of instruments the galaxy survey area, depth, and the 21-cm foreground contaminated region of Fourier space. We find that the highest S/N is achievable through slitless, wide-area spectroscopic surveys, with the proposed Roman HLS survey resulting in a similar to 55 sigma (similar to 13 sigma) detection of the cross-power with 21-cm as observed with SKA-low AA* (HERA-350), for our fiducial model and assuming similar to 500 sq. deg. of overlap. Narrow-band dropout surveys are unlikely to result in a detectable cross-power, due to their poor redshift localization. Slit spectroscopy can provide a high S/N detection of the cross-power for SKA-low AA* observations. Specifically, the planned MOONRISE survey with MOONS on the VLT can result in a similar to 3 sigma detection, while a survey of comparable observing time using MOSAIC on the ELT can result in a similar to 4 sigma detection. Our results can be used to guide survey strategies, facilitating the detection of the galaxy-21-cm cross-power spectrum.

A precise measurement of photometric redshifts (photo-z) is crucial for the success of modern photometric galaxy surveys. Machine learning (ML) methods show great promise in this context, but suffer from covariate shift in training sets due to selection bias where interesting sources, e.g., high redshift objects, are underrepresented, and the corresponding ML models exhibit poor generalisation properties. We present an application of the StratLearn method to the estimation of photo-z (StratLearn-z), validating against simulations where we enforce the presence of covariate shift to different degrees. StratLearn is a statistically principled approach which relies on splitting the combined source and target datasets into strata, based on estimated propensity scores. The latter is the probability for an object in the dataset to be in the source set, given its observed covariates. After stratification, two conditional density estimators are fit separately within each stratum, and then combined via a weighted average. We benchmark our results against the GPz algorithm, quantifying the performance of the two algorithms with a set of metrics. Our results show that the StratLearn-z metrics are only marginally affected by the presence of covariate shift, while GPz shows a significant degradation of performance, specifically concerning the photo-z prediction for fainter objects for which there is little training data. In particular, for the strongest covariate shift scenario considered, StratLearn-z yields a reduced fraction of catastrophic errors, a factor of 2 improvement for the RMSE as well as one order of magnitude improvement on the bias. We also assess the quality of the predicted conditional redshift estimates using the probability integral transform (PIT) and the continuous rank probability score (CRPS). The PIT for StratLearn-z indicates that predictions are well-centered around the true redshift value, if conservative in their variance; the CRPS shows marked improvement at high redshifts when compared with GPz. Our julia implementation of the method, StratLearn-z, is publicly available at .

Inflationary models that involve bursts of particle production generate bump-like features in the primordial power spectrum of density perturbations. These features influence the evolution of density fluctuations, leaving their unique signatures in cosmological observations. A detailed investigation of such signatures would help constrain physical processes during inflation. With this motivation, the goal of this paper is two-fold. First, we conduct a detailed analysis of the effects of bump-like primordial features on the sky-averaged 21 cm signal. Using semi-numerical simulations, we demonstrate that the primordial features can significantly alter the ionization history and the global 21 cm profile, making them a promising probe of inflationary models. We found a special scale (namely, the turnover wavenumber, k^ turn) at which the effect of primordial bump-like features on the global 21 cm profile vanishes. Also, we found that the behaviour of the primordial features on the global profile and ionization history are quite opposite for k > k^ turn and k < k^ turn. We trace the root cause of these behaviours to the effects of primordial features on the halo mass function at high redshifts. Furthermore, we discuss the degeneracy between the astrophysical parameters and the primordial features in detail. Secondly, for a fixed set of astrophysical parameters, we derive upper limits on the amplitude of bump-like features in the range 10^-1 < k [ Mpc^-1] < 10^2 using current limits on optical depth to reionization from CMB data by Planck.

Understanding the epochs of cosmic dawn and reionisation requires us to leverage multi-wavelength and multi-tracer observations, with each dataset providing a complementary piece of the puzzle. To interpret these data, we updated the public simulation code, 21cmFASTv4, to include a discrete source model based on stochastic sampling of conditional mass functions and semi-empirical galaxy relations. We demonstrate that our new galaxy model, which parametrises the means and scatters of well-established scaling relations, is flexible enough to characterise a range of predictions from different hydrodynamic cosmological simulations of high-redshift galaxies. Combining a discrete galaxy population with approximate, efficient radiative transfer allows us to self-consistently forward-model galaxy surveys, line intensity maps (LIMs), and observations of the intergalactic medium (IGM). Not only does each observable probe different scales and physical processes, but their cross-correlation will maximise the information gained from each measurement by probing the galaxy-IGM connection at high redshift. In this work, we found that a stochastic source field produces significant shot-noise in 21cm and LIM power spectra. Scatter in galaxy properties can be constrained using ultraviolet (UV) luminosity functions and/or 21cm power spectra, especially if astrophysical scatter is higher than expected (as might be needed to explain recent JWST observations). Our modelling pipeline is both flexible and computationally efficient, thereby facilitating high-dimensional, multi-tracer, field-level Bayesian inference of cosmology and astrophysics over the first billion years.

Mapping out the first billion years using the 21-cm line with the Square Kilometer Array (SKA) will revolutionize our understanding of the cosmic dawn, reionization and the galaxies that drove these milestones. However, synergies with other telescopes in the form of cross correlations will be fundamental in making and confirming initial, low signal-to-noise claims of a detection. Participants of the 2023 European Southern Observatory (ESO) - SKA Observatory workshop discussed such synergies for Epoch of Reionization (EoR) and Cosmic Dawn (CD) science. Here we highlight some of the most promising candidates for cross-correlating SKA EoR/CD observations with ESO instruments such as the Multi-Object Optical and Near-infrared Spectrograph (MOONS), the MOSAIC multi-object spectrograph, and the ArmazoNes high Dispersion Echelle Spectrograph (ANDES).

The Lyman alpha (Ly alpha) forest in the spectra of z>5 quasars provides a powerful probe of the late stages of the epoch of reionisation (EoR). With the recent advent of exquisite datasets such as XQR-30, many models have struggled to reproduce the observed large-scale fluctuations in the Ly alpha opacity. Here we introduce a Bayesian analysis framework that forward-models large-scale lightcones of intergalactic medium (IGM) properties and accounts for unresolved sub-structure in the Ly alpha opacity by calibrating to higher-resolution hydrodynamic simulations. Our models directly connect physically intuitive galaxy properties with the corresponding IGM evolution, without having to tune 'effective' parameters or calibrate out the mean transmission. The forest data, in combination with UV luminosity functions and the CMB optical depth, are able to constrain global IGM properties at percent level precision in our fiducial model. Unlike many other works, we recover the forest observations without invoking a rapid drop in the ionising emissivity from z similar to 7 to 5.5, which we attribute to our sub-grid model for recombinations. In this fiducial model, reionisation ends at z=5.44 +/- 0.02 and the EoR mid-point is at z=7.7 +/- 0.1. The ionising escape fraction increases towards faint galaxies, showing a mild redshift evolution at fixed UV magnitude, M-UV. Half of the ionising photons are provided by galaxies fainter than M-UV similar to-12, well below direct detection limits of optical/NIR instruments including JWST. We also show results from an alternative galaxy model that does not allow for a redshift evolution in the ionising escape fraction. Despite being decisively disfavoured by the Bayesian evidence, the posterior of this model is in qualitative agreement with that from our fiducial model. We caution, however, that our conclusions regarding the early stages of the EoR and which sources reionised the Universe are more model-dependent.

In 2018 the EDGES experiment claimed the first detection of the global cosmic 21 cm signal, which featured an absorption trough centered around z similar to 17 with a depth of approximately -500 mK. This amplitude is deeper than the standard prediction (in which the radio background is determined by the cosmic microwave background) by a factor of two and potentially hints at the existence of a radio background excess. While this result was obtained by fitting the data with a phenomenological flattened-Gaussian shape for the cosmological signal, here we develop a physical model for the inhomogeneous radio background sourced by the first galaxies hosting population III stars. Star formation in these galaxies is quenched at lower redshifts due to various feedback mechanisms, so they serve as a natural candidate for the excess radio background indicated by EDGES without violating present-day measurements by ARCADE2. We forward-model the EDGES sky temperature data, jointly sampling our physical model for the cosmic signal, a foreground model, and residual calibration errors. We compared the Bayesian evidence obtained by varying the complexity and prior ranges for the systematics. We find that the data are best explained by a model with seven log-polynomial foreground terms and a component accounting for calibration residuals. Interestingly, the presence of a cosmic 21 cm signal with a non-standard depth is decisively disfavored. This result is contrary to previous EDGES analyses in the context of extra radio background models, thus serving as a caution against using a "pseudo-likelihood" built on a model (flattened Gaussian) that is different from the one being used for inference. We make our simulation code and associated emulator publicly available.

Ionized bubble sizes during reionization trace physical properties of the first galaxies. JWST's ability to spectroscopically confirm and measure Lyman-alpha (Lyα) emission in sub-L* galaxies opens the door to mapping ionized bubbles in 3D. However, existing Lya-based bubble measurement strategies rely on constraints from single galaxies, which are limited by the large variability in intrinsic Lyα emission. As a first step, we present two bubble size estimation methods using Lya spectroscopy of ensembles of galaxies, enabling us to map ionized structures and marginalize over Lyα emission variability. We test our methods using Gpc-scale reionization simulations of the intergalactic medium (IGM). To map bubbles in the plane of the sky, we develop an edge detection method based on the asymmetry of Lyα transmission as a function of spatial position. To map bubbles along the line-of-sight, we develop an algorithm using the tight relation between Lyα transmission and the line-of-sight distance from galaxies to the nearest neutral IGM patch. Both methods can robustly recover bubbles with radius ≳10 comoving Mpc, sufficient for mapping bubbles even in the early phases of reionization, when the IGM is ∼70-90% neutral. These methods require ≳0.002-0.004 galaxies/cMpc^3, a 5σ Lyα equivalent width upper limit of ≲30Å for the faintest targets, and redshift precision Δ z ≲ 0.015, feasible with JWST spectroscopy. Shallower observations will provide robust lower limits on bubble sizes. Additional constraints on IGM transmission from Lyα escape fractions and line profiles will further refine these methods, paving the way to our first direct understanding of ionized bubble growth.