The tumultuous effects of ultraviolet photons that source cosmic reionization the subsequent compression and shockheating of lowdensity regions and the modulation of baryons in shallow potential wells induced by the passage of ionization fronts collectively introduce perturbations to the evolution of the intergalactic medium in the postreionization era. These enduring fluctuations persist deep into the postreionization era casting a challenge upon precision cosmology endeavors targeting tracers in this cosmic era. Simultaneously these relics from reionization also present a unique opportunity to glean insights into the astrophysics that govern the epoch of reionization. In this work we propose a first study of the crosscorrelation of lya forest and 21 cm intensity mapping accounting for the repercussions of inhomogeneous reionization in the postreionization era. We investigate the ability of SKA times DESIlike SKA times MUSTlike and PUMA times MUSTlike instrumental setups to achieve a high signaltonoise ratio SNR in the redshift range 3.5 leq z leq 4. Moreover we assess how alterations in integration time survey area and reionization scenarios impact the SNR. Furthermore we forecast the crosscorrelations potential to constrain cosmological parameters under varying assumptions considering or disregarding reionization relics marginalizing over reionization astrophysics and assuming perfect knowledge of reionization. Notably our findings underscore the remarkable capability of a futuristic PUMA times MUSTlike setup with a modest 100hour integration time over a 100 sq. deg. survey to constrain the ionization efficiency error to sigmazeta  3.42 .
['Astrophysics - Cosmology and Nongalactic Astrophysics']
Reionization relics in the crosscorrelation between the Lyalpha forest and 21 cm intensity mapping in the postreionization era
{'Reionization relics in the cross-correlation between the Ly 𝛼 forest and 21 cm intensity mapping in the post-reionization era': 'Paulo Montero-Camacho, 1 , 2 ★ Catalina Morales-Gutiérrez, 3 , 4 Yao Zhang, 2 Heyang Long, 5 , 6 and Yi Mao 2  \n- 1 Department of Mathematics and Theory, Peng Cheng Laboratory, Shenzhen, Guangdong 518066, China\n- 2 Department of Astronomy, Tsinghua University, Beĳing 100084, China\n- 3 Department of Physics, University of Costa Rica, 11501 San José, Costa Rica.\n- 4 Space Research Center (CINESPA), University of Costa Rica, 11501 San José, Costa Rica.\n- 5 Department of Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA\n- 6 Center for Cosmology and AstroParticle Physics (CCAPP), The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': "The tumultuous effects of ultraviolet photons that source cosmic reionization, the subsequent compression and shock-heating of low-density regions, and the modulation of baryons in shallow potential wells induced by the passage of ionization fronts, collectively introduce perturbations to the evolution of the intergalactic medium in the post-reionization era. These enduring fluctuations persist deep into the post-reionization era, casting a challenge upon precision cosmology endeavors targeting tracers in this cosmic era. Simultaneously, these relics from reionization also present a unique opportunity to glean insights into the astrophysics that govern the epoch of reionization. In this work, we propose a first study of the cross-correlation of Ly 𝛼 forest and 21 cm intensity mapping, accounting for the repercussions of inhomogeneous reionization in the post-reionization era. We investigate the ability of SKA × DESI-like, SKA × MUST-like, and PUMA × MUST-like instrumental setups to achieve a high signal-to-noise ratio (SNR) in the redshift range 3 . 5 ≤ 𝑧 ≤ 4. Moreover, we assess how alterations in integration time, survey area, and reionization scenarios impact the SNR. Furthermore, we forecast the cross-correlation's potential to constrain cosmological parameters under varying assumptions: considering or disregarding reionization relics, marginalizing over reionization astrophysics, and assuming perfect knowledge of reionization. Notably, our findings underscore the remarkable capability of a futuristic PUMA × MUST-like setup, with a modest 100-hour integration time over a 100 sq. deg. survey, to constrain the ionization efficiency error to 𝜎 𝜁 = 3 . 42.  \nKey words: intergalactic medium - dark ages, reionization, first stars", '1 INTRODUCTION': "Cosmic reionization, the transformative phase during which our Universe shifted from predominantly neutral to highly ionized, is estimated to have occurred approximately around redshift 𝑧 ∼ 8 (e.g. Planck Collaboration et al. 2020). While the overarching mechanisms driving cosmic reionization are broadly understood (Mesinger 2016, 2019), the lack of direct observations introduces significant uncertainties, particularly regarding the timeline of reionization (Keating et al. 2020; Gnedin 2022; Jin et al. 2023; Roth et al. 2023). With the advent of the James Webb Space Telescope (JWST; Gardner et al. 2006) our knowledge of the sources of ultraviolet photons responsible for reionization is likely to increase dramatically. Moreover, reionization unfolds as a markedly inhomogeneous inside-out process (Lee et al. 2008; Choudhury et al. 2009) wherein denser regions undergo ionization first since sources of ultraviolet photons are expected to predominantly be situated within these dense environments.  \nThe post-reionization era is a treasure trove of cosmological information. The relatively high redshifts characteristic of this pe-  \nriod not only facilitate the exploration of dark matter candidates (Palanque-Delabrouille et al. 2020; Puchwein et al. 2022) but also offer a unique opportunity to study cosmic evolution and structure formation (Wyithe & Loeb 2009; Visbal et al. 2009) before the onset of nonlinearities imposes significant constraints on the observational landscape.  \nThe Lyman𝛼 (Ly 𝛼 ) forest, the absorption features observed in the spectra of background quasars, stands as a pivotal tool for exploring the Universe in the post-reionization era. Its applications extend to the study of H /i.pc and He /i.pc/i.pc reionization (Cen et al. 2009; Upton Sanderbeck & Bird 2020; Montero-Camacho & Mao 2020), the investigation of the matter power spectrum on scales beyond the reach of galaxy surveys (Weinberg et al. 2003), constrain the evolution of the Universe (du Mas des Bourboux et al. 2020), and the inference of cosmological parameters (Chabanier et al. 2019). Moreover, the forest also offers a unique window into the impact of neutrino masses (Yèche et al. 2017). However, obtaining reliable cosmological measurements of the Ly 𝛼 forest at high redshifts proves challenging due to the sparse sampling of quasars (Yèche et al. 2020; Montero-Camacho & Mao  \n2021; Chaussidon et al. 2023), particularly regarding potential measurements of the 3D flux power spectrum 1 .  \nIntensity mapping (IM) involves a trade-off of sacrificing angular resolution to concentrate on the integrated emission from unresolved sources. Analogous to the Ly 𝛼 forest, IM utilizing the 21 cm hyperfine transition of hydrogen presents a versatile probe capable of delving into various aspects of cosmology in the post-reionization era. It can probe H /i.pc reionization (Long et al. 2023), explore early universe cosmological parameters (e.g. primordial non-gaussianity), and constrain the evolution of the universe and structure formation (Castorina & White 2019; Bull et al. 2015; Wyithe & Loeb 2009; Visbal et al. 2009). Additionally, it can investigate the nature of dark matter (Carucci et al. 2015) and more (Square Kilometre Array Cosmology Science Working Group et al. 2020; Cosmic Visions 21 cm Collaboration et al. 2018).  \nNevertheless, the full potential of this observable is impeded by foregrounds, such as galactic synchrotron emission, which surpass the amplitude of the cosmological signal by several orders of magnitude (Wolz et al. 2015). Despite dedicated efforts to mitigate foreground effects (see e.g. Zuo et al. 2023; Diao et al. 2024), a conservative perspective suggests that confirming the cosmic origin of the signal might necessitate cross-correlation with another probe 2 . This cross-correlation typically considers high𝑧 galaxies (Furlanetto & Lidz 2007; La Plante et al. 2023; Hutter et al. 2023) as the additional tracer.  \nThe intense heating of the intergalactic medium (IGM) during cosmic reionization triggers a substantial increase in the IGM temperature, reaching a few times 10 4 K. As ionization fronts propagate, inducing a rise in temperature, low-density regions experience shocks that both heat and compress the gas. These shocks, which originate in denser regions due to the increase in Jeans mass with temperature, will push gas in minivoids to higher adiabats compared to gas in denser environments (Hirata 2018) - such as minihalos. This highentropy mean-density gas is then responsible for the long-lasting impact of reionization in the Ly 𝛼 forest (Montero-Camacho et al. 2019; Montero-Camacho & Mao 2020), essentially constituting the long-lasting memory of reionization at 𝑧 ∼ 2 in the forest. Simultaneously, reionization exerts influence over the baryon abundance within a given halo (Long et al. 2022). The thermal kick resulting from the passage of an ionization front expels some baryons from the halo, particularly for halos with shallow potential wells. This modulation is responsible for long-lasting reionization relics in 21 cm IM, persisting up to 𝑧 ∼ 3 (Long et al. 2023).  \nTherefore, both the Ly 𝛼 forest and H /i.pc 21 cm intensity mapping exhibit broad-spectrum contamination originating from the imprints of cosmic reionization. This contamination is especially pronounced at large scales and high redshifts in the post-reionization era, where the lingering effects of reionization are more pronounced due to less time to dissipate the additional injected energy (Montero-Camacho & Mao 2020). Managing this systematic is not merely a challenge essential for obtaining cosmological information free from bias in  \nthe post-reionization era, but it also represents a novel opportunity to gain insights into the intricate processes governing the reionization history of the Universe (Montero-Camacho & Mao 2021).  \nThe observational programs of these two probes are at different developmental stages. The Ly 𝛼 forest has achieved commendable SNR, e.g. eBOSS (Castorina & White 2019; du Mas des Bourboux et al. 2020) and DESI early (Ramírez-Pérez et al. 2023; Gordon et al. 2023; Ravoux et al. 2023) and year-1 results (DESI Collaboration et al. 2024). In stark contrast, the 2-point functions of the cosmological 21 cm signal at higher redshifts ( 𝑧 > 2) are currently constrained only by upper limits (see e.g., Munshi et al. 2023; Abdurashidova et al. 2022; Mertens et al. 2020). Despite this, because of the inherent challenges in measuring the 21 cm auto-power spectrum, there is considerable interest in the cross-correlation of these two distinct tracers as a promising avenue for cosmological studies.  \nThe Ly 𝛼 forest × H/i.pc 21 cm IM cross-correlation was initially proposed in Guha Sarkar et al. (2011) as a robust and independent probe of the post-reionization IGM. Carucci et al. (2017) underscored the reduced sensitivity of this cross-correlation to foreground contamination, emphasizing its potential to untangle degeneracies inherent in modeling parameters. Furthermore, the Ly 𝛼 × 21 cm crosscorrelation exhibits promising prospects for constraining dark energy (Dash & Guha Sarkar 2021) and 𝑓 ( 𝑟 ) gravity models, particularly when coupled with cross-correlations involving cosmic microwave background (CMB) lensing (Dash et al. 2023). Moreover, this observable has also emerged as a potential avenue for constraining the nature of dark matter (Sarkar et al. 2019).  \nRegarding the observability of this cosmological probe, the Owens Valley Widefield Array (OWFA) (Ali & Bharadwaj 2014; Bharadwaj et al. 2015) and a spectroscopic instrument such as BOSS/eBOSS were anticipated to achieve a robust 6 𝜎 detection at 𝑧 = 3 . 35 with 200 hours of integration (Sarkar et al. 2018). Meanwhile, the Square Kilometre Array Phaew 1 Mid-frequency (SKA1-Mid) in conjunction with an eBOSS-like survey could reach a peak SNR of 15 at 𝑧 = 2 . 5 (Guha Sarkar & Datta 2015).  \nHowever, the existing studies/forecasts have yet to incorporate the enduring impact of inhomogeneous reionization on the postreionization IGM. In this work, we delve into the repercussions of these reionization remnants in the Ly 𝛼 forest × H/i.pc 21 cm IM crosscorrelation within the redshift range 3 . 5 ≤ 𝑧 ≤ 4. We specifically target this range because of the heightened abundance of Ly 𝛼 spectra compared to higher redshifts and the pronounced strength of reionization relics in the forest during this period (Montero-Camacho & Mao2020). Nevertheless, it is crucial to note that the impact of reionization imprints will remain significant for this cross-correlation in the broader redshift range of 3 ≲ 𝑧 < 6.  \nThe rest of this work is organized as follows. In section 2, we outline our model for capturing the impact of reionization in the postreionization IGM. We present the simulations required to compute the reionization relics in section 3. The instrumental configurations considered in this study are introduced in section 4. Section 5 presents our findings regarding the effectiveness of the instrumental setups in detecting the Ly 𝛼 × 21 cm cross-correlation. Additionally, we consider deviations in survey design and strategy, providing insights into possible gains. In section 6, we demonstrate the importance of accounting for this novel effect through three distinct Fisher forecasts. Finally, we summarized our findings in 7.  \nThroughout this work, we use ℎ = 0 . 6774, Ω 𝑏 ℎ 2 = 0 . 0223, Ω 𝑐 ℎ 2 = 0 . 1188, 𝐴 𝑠 = 2 . 148 × 10 -9 , and 𝑛 𝑠 = 0 . 9667. In agreement with the ' TT + TE + EE + lowP + lensing + ext ' cosmology from Planck Collaboration et al. (2016). Furthermore, we use the values of 𝑏 H/i.pc and Ω H/i.pc reported in Crighton et al. (2015). In con-  \nSD parameters for the forest from Arinyo-i-Prats et al. (2015).", '2 MODELING THE CROSS-CORRELATION OF Ly 𝛼 AND 21 CM IM WITH REIONIZATION RELICS': "In this work, we consider the inclusion - at first order - of the memory of reionization, that is the long-lasting impact of inhomogeneous reionization, in both Ly 𝛼 flux fluctuations and H /i.pc 21 cm fluctuations. We examine the repercussions of this inclusion when embedded into the traditional cross-correlation of the Ly 𝛼 forest and 21 cm intensity mapping. Mathematically,  \n𝑃 21 ,𝐹 ( 𝒌 , 𝑧 ) = 𝑃 Fid . 21 ,𝐹 ( 𝒌 , 𝑧 ) + 𝑃 Mem . 𝐹 ( 𝒌 , 𝑧 ) + 𝑃 Mem . 21 ( 𝒌 , 𝑧 ) , (1)  \nwhere the first term on the right-hand side corresponds to the conventional cross-correlation of the Ly 𝛼 forest and 21cm IM without reionization relics. The second term accounts for the memory of reionization present in the Ly 𝛼 forest, originating in the underdense regions. The final term represents the memory of reionization in H /i.pc 21 cm IM, i.e. the memory of reionization sourced by overdense regions. The three terms are given by (Carucci et al. 2017; MonteroCamacho et al. 2019; Long et al. 2023)  \n𝑃 Fid . 21 ,𝐹 = 𝑏 𝐹 ( 𝑧 ) 𝑏 H/i.pc ( 𝑧 ) [ 1 + 𝛽 𝐹 ( 𝑧 ) 𝜇 2 ] [ 1 + 𝛽 H/i.pc ( 𝑧 ) 𝜇 2 ] 𝑃 𝑚 ( 𝑘, 𝑧 ) ,  \n𝑃 Mem . 21 = 𝑏 𝐹 ( 𝑧 ) [ 1 + 𝛽 𝐹 ( 𝑧 ) 𝜇 2 ] 𝑃 𝑚, Ξ ( 𝑘, 𝑧 ) ,  \n(2) 𝑃 Mem . 𝐹 = 𝑏 H/i.pc ( 𝑧 ) [ 1 + 𝛽 H/i.pc ( 𝑧 ) 𝜇 2 ] 𝑏 Γ ( 𝑧 ) 𝑃 𝑚,𝜓 ( 𝑘, 𝑧 ) , (3) (4)  \nwhere 𝑏 𝐹 and 𝛽 𝐹 are the usual flux bias and redshift space distortion parameter, respectively - and similarly for 21 cm quantities. Moreover, 𝜇 is the angle with respect to the line of sight, 𝑏 Γ is the radiation bias defined in Arinyo-i-Prats et al. (2015); Hirata (2018), and 𝑃 𝑚 is the linear matter power spectrum.  \nFurthermore, 𝑃 𝑚,𝜓 and 𝑃 𝑚, Ξ are the cross-power spectrum of matter and change of transparency of the IGM due to the impact of reionization, and the cross-power spectrum of matter and change of neutral hydrogen density induced by the passage of reionization fronts, respectively. These cross-power spectra are defined as follows (Montero-Camacho et al. 2019; Long et al. 2023)  \n𝑃 𝑚,𝜓 ( 𝑘, 𝑧 obs ) = -∫ d 𝑧 𝜕𝜓 ( 𝑧, 𝑧 obs ) 𝜕𝑧 𝑃 𝑚,𝑥 HI ( 𝑘, 𝑧 ) 𝐷 ( 𝑧 obs ) 𝐷 ( 𝑧 ) , (5) 𝑃 𝑚, Ξ ( 𝑘, 𝑧 obs ) = -∫ d 𝑧 𝜕 Ξ ( 𝑧, 𝑧 obs ) 𝜕𝑧 𝑃 𝑚,𝑥 HI ( 𝑘, 𝑧 ) 𝐷 ( 𝑧 obs ) 𝐷 ( 𝑧 ) , (6)  \nwhere the integration covers the epoch of reionization, 𝑃 𝑚,𝑥 HI is our proxy for the correlation of matter and ionized bubble spatial distribution, which accounts for the patchy nature of reionization. 𝐷 is the growth factor, while 𝜓 denotes the transparency of the IGM, and Ξ represents the modulation of the neutral hydrogen due to the passage of ionization fronts. They are given by  \n𝜓 ( 𝑧 re , 𝑧 obs | 𝑧 re ) = Δ ln 𝜏 1 = ln GLYPH<20> 𝜏 1 ( 𝑧 re , 𝑧 obs ) 𝜏 1 ( 𝑧 re , 𝑧 obs ) GLYPH<21> , (7)  \nΞ ( 𝑧 re , 𝑧 obs | 𝑧 re ) = Δ ln 𝜌 HI = ln GLYPH<20> 𝜌 HI ( 𝑧 re , 𝑧 obs ) 𝜌 HI ( 𝑧 re , 𝑧 obs ) GLYPH<21> , (8)  \nwhere 3 . 5 ≤ 𝑧 obs ≤ 4, 𝑧 re is the local redshift of reionization, and 𝑧 re serves as a reference redshift of reionization. The parameter 𝜏 1 is the optical depth that must be assigned in simulations to a patch of gas with mean density and temperature 𝑇 = 10 4 K in order for the mean transmitted flux to match observations, i.e. 𝜏 1 is a normalization factor that guarantees that an optical depth cube return sensible results by matching to the observations of Kim et al. (2007). The transparency  \nis defined as a relative measure, benchmarked against a fiducial scenario with 𝑧 re = 8. Conversely, Ξ characterizes the response of halos with shallow potential wells to the passage of ionization fronts, modeling the perturbations in neutral hydrogen density induced by the response relative to those that occur in a fiducial local reionization scenario. We plot some of the ingredients of the reionization relics in Figure B1 for reference.  \nSimilarly, the Ly 𝛼 flux and the 21 cm power spectra are respectively given by  \n𝑃 3D 𝐹 ( 𝒌 , 𝑧 ) = 𝑏 2 𝐹 ( 1 + 𝛽 𝐹 𝜇 2 ) 2 𝑃 𝑚 + 2 𝑏 𝐹 𝑏 Γ ( 1 + 𝛽 𝐹 𝜇 2 ) 𝑃 𝑚,𝜓 (9) 𝑃 21 ( 𝒌 , 𝑧 ) = 𝑏 2 H/i.pc ( 1 + 𝛽 H/i.pc 𝜇 2 ) 2 𝑃 𝑚 + 2 𝑏 H/i.pc ( 1 + 𝛽 H/i.pc 𝜇 2 ) 𝑃 𝑚, Ξ (10)  \nwhere we have neglected terms with higher order in Ξ and 𝜓 . Likewise, we disregard non-linear corrections, justified by our focus on large scales.  \nEven though we extract 𝑏 𝐹 and 𝛽 𝐹 from Arinyo-i-Prats et al. (2015), the maximum redshift in their tables is 3. Hence, for our 𝑧 > 3 calculations, we include a redshift evolution factor [( 1 + 𝑧 )/( 1 + 𝑧 𝑝 )] 3 . 55 (Palanque-Delabrouille et al. 2013b) with pivot 𝑧 𝑝 = 3 to account for the evolution of the flux bias and redshift space distortion parameter (see Eq.(10) and surrounding text of Montero-Camacho et al. 2024b).  \nWhile both Eq. (5, 6) emerge from the IGM's response to inhomogeneous reionization, they represent distinct response mechanisms. Eq. (7) covers the response to the local reionization process by underdense gas (Hirata 2018). Physically, it originates in minivoids, as underdense gas undergoes ultraviolet heating, shock heating, and compression during the reionization process. In contrast, Eq. (8) captures the response to local reionization by denser regions, predominantly arising from minihalos. Essentially, this is the modulating effect on the number of baryons allowed inside a given halo following the thermal kick delivered by an ionization front (Long et al. 2022). Note that both mechanisms differ in their origin on the small scales but are influenced equivalently by the patchy nature of reionization - parametrized here as 𝑃 𝑚,𝑥 HI and present in both Eq. (5, 6).  \nTo facilitate comparison between the different components of Eq. (1), we have chosen to utilize the spherically-averaged power spectrum 3 , i.e.  \n𝑃 Sph . 21 ,𝐹 ( 𝑘, 𝑧 ) = 1 2 𝜋 ∫ 2 𝜋 0 𝑑𝜙 ∫ 1 0 𝑑𝜇𝑃 21 ,𝐹 ( 𝒌 , 𝑧 ) , (11)  \nwhere the normalization factor already accounts for the extra factor of 2 due to the 𝜇 -symmetry.  \nWeplot Eq. (11) in Figure 1 at a few selected redshifts. The impact of reionization is more pronounced at higher redshifts, introducing a competition between the 21 cm effect, arising from the modulation of baryons due to reionization, and the Ly 𝛼 forest effect, originating in underdense regions. The former aims to diminish the signal's amplitude while the latter enhances it. Furthermore, the more substantial (in absolute value) impact of the Ly 𝛼 sourced term compared to the H /i.pc 21 cm term is anticipated due to its greater influence on voids, evident in the bimodal temperature-density relation (see Figure 3 in Montero-Camacho & Mao 2020 versus bayron modulation displayed in Figure 8 of Long et al. 2023). While our focus is on the 3 . 5 < 𝑧 < 4 . redshift range, we expect this qualitative trend to  \n<!-- image -->  \n<!-- image -->  \nFigure 1. The magnitude of the spherically-averaged cross-power spectrum for the Lyman-alpha forest × H/i.pc 21 cm IM cross-correlation is illustrated, encompassing its components, including the imprint of reionization in both probes. The blue dash-dotted line corresponds to the conventional signal without any relics from cosmic reionization. The purple dotted and green dashed curves represent the memory of reionization in the Lyman𝛼 forest and 21 cm power, respectively. We highlight that the memory of reionization in the 21 cm signal is negative, leading to a competition in which the Ly 𝛼 forest, which is positive, dominates. The influence of reionization becomes more pronounced at higher redshifts and larger scales.  \n<!-- image -->  \npersist at higher redshifts especially given that the impact of selfshielding in minihalos will also enhance the Ly 𝛼 power spectrum at large scales (Park et al. 2023). Regarding the wavenumber trend, the larger deviation at large scales compared to the conventional term is explained by the coupling to the reionization bubble scale in the cross-correlation of the matter and neutral hydrogen fraction field, 𝑃 𝑚,𝑥 H/i.pc .", '3 SIMULATIONS': 'To accurately model the reionization process, it is imperative to probe under the Jeans length prior to the passage of an ionization front ( ∼ 100 kpc). Failure to do so would result in the loss of the ability to track the response by the small-scale structure due to the significant wiping out of these structures caused by the increasing Jeans scale throughout reionization (Hirata 2018). Besides, reionization occurs in an inhomogeneous way. Thus, to secure enough statistical power, there is a need to simulate large comoving volumes (a few hundred Mpc) with substantial ionized bubbles each with a radius of a few Mpc. However, the dynamical range required for these considerations is too large to accomplish with reasonable computational resources using a single simulation.  \nWeadopt the methodology introduced by Montero-Camacho et al. (2019); Long et al. (2023) and employ a hybrid approach to compute Eq. (5,6). This approach relies on seminumerical simulations using 21/c.pc/m.pcFAST (Mesinger et al. 2011; Murray et al. 2020), which track the patchy nature of reionization within a 400 Mpc box. In addition, we leverage high-mass resolution small-box ( 𝐿 = 1152 ℎ -1 ckpc) simulations run using a modified version of G/a.pc/d.pc/g.pc/e.pc/t.pc-2 (Springel 2005; Hirata 2018) to monitor the response to ionization fronts in both dense (crucial for H /i.pc 21 cm) and underdense (Ly 𝛼 forest) regions. This hybrid simulation strategy allows us to capture the impact of reionization across a range of scales and redshifts with minimal trade-offs. Naturally, this hybrid methodology will not be sensitive  \nto higher-order correlations between the reionization field and the small-scale density field (Montero-Camacho et al. 2019).  \nThe small-box simulations include adiabatic expansion (including Hubble expansion), shock heating, Compton heating, and cooling for neutral gas (accounting for residual ionization). In addition, we incorporate the long-lasting impact of X-ray preheating in the post-reionization IGM (Montero-Camacho et al. 2024b). For ionized gas, we also include Compton cooling, He /i.pc/i.pc cooling, recombination cooling, photoionization heating, and free-free cooling. These simulations have been described - and tested - in detail in Hirata (2018); Long et al. (2022).  \nTo surpass the limitation imposed by the size of our 21/c.pc/m.pcFAST boxes, we utilize a simple liner bias to estimate the effect of reionization at 𝑘 ≲ 0 . 06 Mpc -1 as follows  \n𝑃 𝑚,𝑋 ( 𝑘, 𝑧 ) = 𝑃 𝑚,𝑋 ( 𝑘 cut ≈ 0 . 06 , 𝑧 ) 𝑃 𝑚 ( 𝑘 cut , 𝑧 ) 𝑃 𝑚 ( 𝑘, 𝑧 ) , (12)  \nwhere 𝑋 = { 𝜓, Ξ } . This approximation should be sufficient for smaller scales than the ionized bubble scale. Note that this restriction is due to the box side of our large-box simulations and it is a self-imposed restriction based on available computational resources.', '4 TELESCOPES': "To measure the cross-correlation between 21 cm and the Ly 𝛼 forest, we require a radio interferometer and a spectroscopic telescope with some non-negligible overlap in sky coverage and redshift. Unfortunately, our choice for estimating the SNR is not straightforward due to the timeline of the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2022) versus the timeline of the Square Kilometer array (SKA; Braun et al. 2019). While DESI may potentially extend its operational life as DESI-II, allowing for direct operational overlap with SKA, the exact timeline overlap remains uncertain.  \nConcurrently, upcoming spectroscopic instruments such as  \nz  \n=3.94  \nMegaMapper (Schlegel et al. 2019), the MUltiplexed Survey Telescope (MUST; Zhang et al. 2024), and MaunaKea Spectroscopic Explorer (MSE; Percival et al. 2019), categorized as Stage V spectroscopic instruments, are anticipated to be operational around the time when 21 cm radio interferometers come online. However, we underscore that the detailed designs of these instruments are currently in the developmental phase. Navigating the evolving landscape of these instrument timelines and designs will be vital to optimizing an observational program capable of measuring the cross-correlation between 21 cm and the Ly 𝛼 forest.  \nGiven these 'instrumental' constraints, we have opted to explore several distinct instrumentation scenarios.  \n- · SKA1-Low × DESI-like spectroscopic telescope (referred to as SKA in relevant figures).\n- · SKA1-Low × DESI-MUST-like hybrid (referred to as MUST).\n- · PUMA-like × DESI-MUST-like hybrid (referred to as PUMA).  \nAcross all scenarios, we maintain a redshift coverage of 3 . 5 ≤ 𝑧 ≤ 4 and consider two sky coverage overlaps - 100 and 1000 square degrees. Likewise, we vary the integration times, considering both 100 and 1000 hours, which will mainly affect the radio observations as described below.  \nOur base reference scenario involves the SKA1-Low × DESI-like pair, with 100 hours of integration time and a 100-square-degree overlap between the instruments. This reference scenario assumes a Planck-like reionization timeline, which may be underestimating the impact of reionization since Planck uses a hyperbolic tangent as a sigmoid to model reionization that proves a poor fit of astrophysical constraints on the timeline of reionization (see Figure 4 of MonteroCamacho et al. 2024a).", '4.1 SKA1-Low': "SKA1-Low, our baseline 21 cm instrument, operates as an aperture array radio interferometer. We consider only the dense core of the array, which is most sensitive to the 21 cm power spectrum. Do note that the outer stations in the array are crucial for calibration and foreground removal purposes. The SKA1-Low, upon completion, offers versatility in supporting various observing modes. Previous works have often focused on observing mode 1, characterized by a dense core comprising 224 stations of 40 meters in diameter, each containing 256 dipole antennas (Square Kilometre Array Cosmology Science Working Group et al. 2020; Zhao et al. 2022). These core stations are distributed over a diameter of approximately one kilometer.  \nNevertheless, as pointed out in SKA1-LOW Configuration - Constraints & Performance Analysis §5.4.3, power spectrum measurements are better served by 'substations' with 10 m of diameter. Here we chose to consider only the observing mode 4, as detailed in Table 2 of SKA1-LOW Configuration - Constraints & Performance Analysis. This mode, supported by the instrument's correlator, boasts six times more correlatable elements (substations) at the expense of a reduced bandwidth (8.4 MHz compared to 300 for mode 1).  \nTo quantify the unique baselines covered by the dense core, we utilize the number density of baselines 𝑛 𝑏 from Villaescusa-Navarro et al. (2015) and normalize it to recover the total amount of distinct baselines covered by the dense core, as outlined in §2 of Long et al. (2023).  \nThe system temperature - the sum of the sky temperature, the  \nground reflections, and the temperature of the receiver - is given by  \n𝑇 sys ( 𝜈 ) = 60 GLYPH<18> 300 MHz 𝜈 GLYPH<19> 2 . 55 × 1 . 1 + 40 [K] , (13)  \nwhile the effective collecting area per station 𝐴 𝑒 is given by (Square Kilometre Array Cosmology Science Working Group et al. 2020)  \n𝐴 𝑒 ( 𝜈 ) = 𝐴 𝑒, crit × ( GLYPH<0> 𝜈 crit 𝜈 GLYPH<1> 2 , 𝜈 > 𝜈 crit 1 , 𝜈 ≤ 𝜈 crit , (14)  \nwhere 𝜈 crit = 110 MHz and 𝐴 𝑒, crit is the collecting area in m 2 for the 256 dipole antennas of 3.2 m 2 each for observing mode 1. Thus, we approximate the effective area for observing mode 4 by decreasing Eq. (14) by a factor of 6.  \nMoreover, the field of view is given by  \nFoV ( 𝜆 ) = 𝜆 √ 0 . 7 𝐷 phys ! 2 , (15)  \nwhere 𝐷 phys is the diameter of the correlatable element and 𝜆 is the observing wavelength.", '4.2 PUMA': 'ThePacked Ultra-wideband Mapping Array design (PUMA; Cosmic Visions 21 cm Collaboration et al. 2018; Slosar et al. 2019) currently - consists of 32000 antennas distributed in a hexagonal-close packing in a compact circle of roughly 1.5 km diameter. Each antenna has a diameter of 6 meters. We take the baseline number density 𝑛 b directly from Appendix D of Cosmic Visions 21 cm Collaboration et al. (2018).  \nFor PUMA, the system temperature is given by  \n𝑇 sys ( 𝜈 ) = 25 GLYPH<18> 400 MHz 𝜈 GLYPH<19> 2 . 75 + 2 . 7 + 300 9 + 50 0 . 81 [K] . (16)  \nThe effective collecting area per antenna is simply the effective area of the antenna dish, i.e. 𝜋 ( √ 0 . 7 𝐷 phys / 2 ) 2 . Simultaneously, the field of view is given by Eq. (15) but with 𝐷 phys = 6 m. We set the bandwidth to 8.4 MHz (same as SKA1-Low observing mode 4).  \nWeconsider PUMA as an optimistic/futuristic 21 cm telescope for our purposes, essentially our best benchmark.', '4.3 DESI-like': 'We assume a DESI-like spectrograph as our baseline Ly 𝛼 forest instrument, which is likely a conservative estimate given the projected timeline of SKA. To gauge the performance of the final - 5 years - DESI data, we use the Quasar Luminosity Function (QLF) to quantify the expected number of Ly 𝛼 quasars observable with DESI (Palanque-Delabrouille et al. 2013a; Yèche et al. 2020) and the spectrograph performance of DESI - see Eq. (21). Detailed information about the DESI instrument is available in DESI Collaboration et al. (2022), the QSO and Ly 𝛼 QSO target selection can be found in Chaussidon et al. (2023), and the spectroscopic pipeline is described in Guy et al. (2023).', '4.4 DESI-MUST-like hybrid': "Werefine our initially conservative DESI-like scenario by enhancing its performance to align more closely with the standards expected for Stage V spectroscopic instruments. Pragmatically, this improvement involves adjusting the aliasing term in the covariance computation,  \neffectively augmenting the effective density of lines of sight beyond what is anticipated with our original DESI specs. In summary, we implement an optimistic factor of 3 reduction in both the second and third terms of Eq. (21). Given that MUST will likely be the World's first Stage V spectroscopic instrument 4 , we decided to refer to this scenario as MUST.", '5 RESULTS': 'In the redshift range of interest, the cross-correlation of Ly 𝛼 × 21 cm IMpresents a pragmatically more accessible measurement compared to the auto-correlation analyses of both the Ly 𝛼 forest or of the H /i.pc 21 cm field. This disparity arises from distinct challenges associated with each observable. For the Ly 𝛼 forest, the QLF peaks around 𝑧 ∼ 2, resulting in a substantial reduction in the number density of quasars at 𝑧 ∼ 4. Conversely, for 21 cm intensity mapping, the signal strength is significantly weaker than the foregrounds by several orders of magnitude. Notably, these challenges are uncorrelated, rendering the cross-correlation less susceptible to their impacts (Carucci et al. 2017). Since this is the first study of the effects of the relics from reionization in this observable, we will not consider foreground contamination nor several Ly 𝛼 forest systematics, like continuum fitting (Sun et al. 2023), spectra with broad absorption lines (Filbert et al. 2023), UV clustering (Long & Hirata 2023), spectra with damped Lyman𝛼 systems (Wang et al. 2022), and several other astrophysical and instrumental systematics. Note that this simpler approach is often used in other studies of this cross-correlation (e.g. Guha Sarkar et al. 2011; Sarkar et al. 2019; Carucci et al. 2017) at other redshifts.', '5.1 Impact on SNR': "Under these considerations, we can compute the errors on the observables following Carucci et al. (2017); Villaescusa-Navarro et al. (2015)  \n𝜎 2 [ 𝑃 21 ( 𝒌 , 𝑧 )] = GLYPH<16> 𝑃 T 21 ( 𝒌 , 𝑧 ) GLYPH<17> 2 , (17)  \n𝜎 2 [ 𝑃 𝐹 ( 𝒌 , 𝑧 )] = GLYPH<16> 𝑃 T 𝐹 ( 𝒌 , 𝑧 ) GLYPH<17> 2 and (18)  \n𝜎 2 [ 𝑃 21 ,𝐹 ( 𝒌 , 𝑧 )] = 1 2 GLYPH<16> 𝑃 2 21 ,𝐹 + 𝜎 [ 𝑃 21 ] 𝜎 [ 𝑃 𝐹 ] GLYPH<17> . (19)  \nFor 21 cm, 𝑃 Tot . has contributions from cosmic variance and thermal noise, i.e.  \n𝑃 T 21 ( 𝒌 , 𝑧 ) = 𝑃 H/i.pc + 𝑇 2 sys ( 𝑧 ) 𝜒 2 ( 𝑧 ) 𝜆 ( 𝑧 ) 1 + 𝑧 𝐻 ( 𝑧 ) GLYPH<18> 𝜆 2 ( 𝑧 ) 𝐴 𝑒 GLYPH<19> 2 GLYPH<18> 𝑆 area FoV ( 𝑧 ) GLYPH<19> × GLYPH<18> 1 𝑁 pol 𝑡 int 𝑛 b ( 𝑢 = 𝑘 ⊥ 𝜒 ( 𝑧 )/ 2 𝜋 ) GLYPH<19> , (20)  \nwhere 𝜒 is the comoving distance, 𝐻 is the Hubble parameter, and 𝜆 = 21cm /( 1 + 𝑧 ) . Moreover, 𝑇 sys stands for the system temperature of the radio interferometer, 𝑁 pol = 2 is the number of polarizations, 𝐴 𝑒 the effective area of the correlatable unit (dishes or stations), 𝑆 area is the survey area, FoV is the field of view, 𝑛 b quantifies the number density of baselines in the 𝑢𝑣 -plane as a function of wavenumber, and 𝑡 int corresponds to the integration time. The shot noise is subdominant for the redshift range of interest here (Castorina & Villaescusa-Navarro 2017), thus we do not include it in Eq. (20).  \nIn contrast, for the Ly 𝛼 forest we have (McDonald & Eisenstein 2007; Font-Ribera et al. 2014; Montero-Camacho & Mao 2021)  \n𝑃 T 𝐹 ( 𝒌 , 𝑧 ) = 𝑃 3D 𝐹 ( 𝒌 , 𝑧 ) + 𝑃 1D 𝐹 ( 𝑘 ∥ , 𝑧 ) 𝑃 2D 𝑤 ( 𝑧 ) + 𝑃 eff 𝑁 ( 𝑧 ) . (21)  \nThe first term corresponds to cosmic variance, i.e. it comes from ⟨ 𝛿 𝐹 𝛿 ∗ 𝐹 ⟩ , while the second term represents the aliasing noise due to the sparse sampling of quasars. The third term describes the effective noise due to the spectrograph performance. Note that at high redshifts ( 𝑧 ≳ 3 . 2), the aliasing term tends to dominate Eq. (21) since the limited number of quasars becomes more restrictive (see Figure 6 of Montero-Camacho & Mao 2021). Naturally, this trend reduces the expected signal-to-noise ratio (SNR) of Ly 𝛼 forest surveys attempting a measurement of the 3D flux power spectrum at high redshifts.  \nGiven Eq. (17, 18, or 19), the SNR for 𝑃 𝑖 can be written as  \nSNR 2 𝑖 = 𝑁 𝑘 𝑃 2 𝑖 𝜎 2 [ 𝑃 𝑖 ] = 𝑉 Survey 𝑘 3 𝜖𝑑𝜇 4 𝜋 2 𝑃 2 𝑖 𝜎 2 [ 𝑃 𝑖 ] , (22)  \nwhere 𝜖 = 𝑑𝑘 / 𝑘 and 𝑉 Survey is the - overlap - volume between the radio interferometer and the spectrograph. We consider three redshift bins centered 5 at 𝑧 c = [ 3 . 61 , 3 . 77 , 3 . 94 ] with widths defined by the bandwidth of the radio interferometer, i.e. Δ 𝜈 = 8 MHz. Consequently, our survey volume for a redshift bin centered at 𝑧 c is given by  \n𝑉 Survey ( 𝑧 c ) = 4 𝜋 3 𝑓 sky [ 𝜒 ( 𝑧 max ) 3 -𝜒 ( 𝑧 min ) 3 )] (23)  \nwith 𝜒 as the comoving distance. Note that 𝑓 sky , the sky coverage, depends on the area of overlap.  \nIn what follows, we focus on the spherical-averaged power spectrum, i.e. Eq. (11) and analogs, to compute the appropriate SNR for Ly 𝛼 × 21 cm, Ly 𝛼 forest, and 21 cm intensity mapping.  \nWe illustrate the observability of the cross-correlation compared to the auto-correlation of Ly 𝛼 and 21 cm in Figure 2 for our base scenario (SKA1-Low × DESI-like, 𝑡 int = 100 h, 𝐴 sky = 100 sq. deg. and Planck's reionization timeline). We highlight that all of the curves here include the long-lasting relics from cosmic reionization, i.e. Eqs. (5, 6). Unsurprisingly, the SNR for the Ly 𝛼 forest is less than unity because of the drop in available line-of-sight at these redshifts and because of the survey volume. In principle, it is necessary to observe more quasars to reduce the mean separation between forests to make the 3D flux power spectrum observation feasible. Although Figure 2 demonstrates that the 21 cm signal has a larger SNR, this is only true in the absence of foregrounds, which will severely bury the auto-power spectrum. In contrast, the cross-correlation is more robust against foreground contamination (Furlanetto & Lidz 2007; Sarkar et al. 2018; Zhou et al. 2021) and it is already competitive with the 21 cm auto-correlation with the baseline instrumental setup.  \nFigure 2 provides a compelling rationale for the significance of the Ly 𝛼 × 21 cm cross-correlation, particularly in the presence of foregrounds. Detection in the cross-correlation can guarantee the cosmological nature of the 21 cm signal. Besides, the SNR is similar to the auto-correlation at the ionized bubble scale (and for our first redshift bin). This is consistent with the findings of Carucci et al. (2017) at a lower redshift ( 𝑧 = 2 . 4). However, we underscore that in the absence of reionization relics, the cross-correlation will not have a similar SNR to that of the auto 21 cm SNR in the lowest redshift bin considered here. The boost provided by the memory of reionization is crucial to be competitive, yet it is not sufficient at higher redshift  \nFigure 3 showcases the dependence of the SNR for different integration times and telescope pairs. Unsurprisingly, there is a consistent enhancement in SNR across all redshifts and configurations. In particular, the SKA1-Low configurations benefit the most from this strategy, manifesting improvements of ≈ 2 . 6 in overall SNR (see Table 1). Intriguingly, the PUMA-like × DESI-MUST-like hybrid exhibits negligible improvement at large scales.  \n<!-- image -->  \n<!-- image -->  \nFigure 2. The SNR for the baseline configuration of the Ly 𝛼 forest (purple squares), 21 cm intensity mapping (green pentagons), and Ly 𝛼 × 21 cm (orange diamonds) is presented as a function of wavenumber in our baseline scenario (see text for details). The three panels correspond to different redshift bins. The cross-correlation is, naturally, the more pragmatic measurement since the 21 cm signal would be severely impacted by foregrounds. Notably, the SNR for the forest is somewhat impeded by the chosen survey volume, although the drop in the number density of quasars leads to difficulties at high redshifts regardless of sensible choices for survey volume. The gray-shaded region corresponds to SNR ≤ 0.  \n<!-- image -->  \nbins because of a significant rise in aliasing noise 𝑃 2D 𝑤 in Eq. (21) - at those redshifts.  \nFurthermore, akin to the findings in Carucci et al. (2017), the cross-correlation can surpass the SNR of the 21 cm auto-correlation at smaller scales. However, this trend is obscured at high redshifts due to the poor sampling of quasars.  \nHaving established the interest in the cross-correlation of the Ly 𝛼 forest and 21 cm IM at the redshift range of interest (3 . 5 ≤ 𝑧 ≤ 4), we now ponder the impact of the survey strategy.", '5.2 Importance of survey design': "Our main objective here is to underscore the importance of a nonnegligible overlap between spectroscopic instruments and radio interferometers. Given that SKA1-Low is under development in the southern hemisphere and DESI - also likely DESI II - is in the northern hemisphere, we investigate the potential advantages and insights that could be derived from this cross-correlation with diverse instrumental configurations.  \nIn order to reduce the noise in the radio interferometer, a widely employed strategy is to consider an increase in the integration time since the thermal noise is inversely proportional to 𝑡 int . Hence, we explore the impact of a tenfold increase in 𝑡 int across various instrumental configurations - specifically SKA1-Low × DESI-like, SKA1-Low × DESI-MUST-like hybrid, and PUMA-like × DESIMUST-like hybrid.  \nWe assume that the Ly 𝛼 forest observations will be completed regardless of the integration time; nevertheless, we reward the longer observation time with a factor of two reduction in the spectrograph effective noise 𝑃 eff 𝑁 . In other words, this implies that the extended observation time allows for more exposures. Note that we do not modify the aliasing error - second term in Eq. (21) that significantly dominates the error budget - even when considering additional observation time.  \nFrom Eq. (20), the thermal noise of the radio interferometer is inversely proportional to the total integration time. The reason Figure 3 shows a larger impact for the pairs using SKA1-Low is because the radio telescope dominates the error budget for those setups, i.e. it is larger than the cosmic variance contribution and 𝜎 [ 𝑃 H/i.pc ] > 𝜎 [ 𝑃 𝐹 ] . This is also why the change from DESI-like to DESI-MUST-hybrid does not result in a major improvement and largely conserves the increase in SNR that was present in the 100-hour baseline scenario. In contrast, for the PUMA × MUSTsetup, the Ly 𝛼 forest survey now dominates the error, gains are then possible at larger wavenumbers given the large drop at large scales (as can be seen in Figure 2). We also highlight that the number density of baselines also restricts the shape of the SNR at these scales, which is why the SKA1-Low × MUST can be unexpectedly competitive with PUMA × MUST's SNR at some wavenumbers for a 1000 hours of integration time.  \nIt is unclear whether major projects like SKA or a Stage V spectroscopic instrument could allocate 1000 hours/many exposures to observe the same overlapping area, given the promising science programs they aim to address. Nevertheless, it is evident that with such a time investment the community could expedite the progress in measurements of the cross-correlation by a decade, achieving observations competitive with what may be available with Stage II 21 cm radio telescopes - at least at large scales. We acknowledge, however, that this is a significant time commitment.  \nNaively, one might expect that expanding the overlapping area would perhaps be a more straightforward enterprise, especially for spectroscopic instruments located in the southern hemisphere like  \n<!-- image -->  \n<!-- image -->  \nFigure 3. SNR for the Ly 𝛼 forest × 21 cm intensity mapping accounting for variations in the integration time of the radio telescopes. We assume that the Ly 𝛼 forest survey is also completed during that integration time and reward the additional integration time (see text for details). Shown are the SNR for: PUMA × MUST 100 hours (blue circles) and 1000 hours (cyan squares); SKA1-Low × DESI 100 hours (green pentagons) and 1000 hours (orange octagons); SKA1-Low × MUST 100 hours (purple crosses) and 1000 hours (magenta diamonds). Increasing the total integration time allows for a better measurement in all configurations.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 4. Similar to Figure 3, but modifying the survey area. Notably, increasing the survey area does little to improve the SNR for the configurations with SKA1-Low since the error is dominated by the interferometric instrument and the thermal noise increases with survey area - see Eq. (20).  \n<!-- image -->  \nthe 4-metre Multi-Object Spectroscopic Telescope (4MOST; Richard et al. 2019). However, it is important to note that increasing the survey area produces a larger thermal noise - see Eq. (20). Conversely, the number of modes 𝑁 𝑘 is directly proportional to the survey area. Thus, as long as the thermal noise is manageable, an increase in the survey area will improve SNR due to the additional modes. For the PUMA setup, where Ly 𝛼 dominates the error budget at large scales, one could anticipate that the increase in modes may give rise to an enhancement in SNR up to a large 𝑘 value. In contrast, the impact of this choice is expected to generate limited gains for the  \nSKA and MUST configurations due to the dominance of the 21 cm error, particularly with a tenfold increase in thermal noise.  \nWe confirm these expectations in Figure 4 by augmenting the area of overlap by a factor of ten from the baseline scenario (from 100 to 1000 square degrees).  \nThe interpretation of this finding is convoluted. Initially, it may be more appealing to allocate 100 hours of a telescope's time to an observation spanning different pointings. Unfortunately, this approach does not lead to a significant enhancement in observations for the late 2020s configurations. In terms of the overall SNR, aggregated across  \nthe 𝑘 -bins, there is only a very modest increase of approximately 1 . 05 compared to the baseline scenario.  \nMoreover, if radio telescopes happen to be located in the global south and spectroscopic instruments keep the trend of being located in the northern hemisphere, then the overlap may be small to begin with and hence preference would be given to securing a small overlap in the survey footprints. To our knowledge the possible location of PUMA is not decided yet, consequently, our results advocate for the consideration of overlap with Stage V spectroscopic instruments. Furthermore, our results also showcase that a larger overlap can produce significant gains for a PUMA cross-correlation at large scales. For PUMA, the overall SNR augments by a factor of 1 . 7 when enlarging the survey area by a factor of 10. Besides, this enhancement is limited by the assumption of Stage V instrument performance. In principle, PUMA may instead share the skies with potential Stage VI spectroscopic instruments.  \nWhile we have only scratched the surface regarding the implications of survey strategy for the Ly 𝛼 × H/i.pc 21 cm cross-correlation at 3 . 5 ≤ 𝑧 ≤ 4, we deem further exploration of the instrumental setup beyond the scope of this work. This study focuses on demonstrating the impact of cosmic reionization in the cross-correlation during the post-reionization era. Therefore, it would be remiss of us not to address the dependence of the signal on the uncertain timeline of reionization (Montero-Camacho & Mao 2021; Jin et al. 2023).", '5.3 Dependence on reionization history': "The timeline of reionization remains uncertain, although abundant direct detection of the ionizing sources is becoming possible thanks to the James Webb Space Telescope (JWST) (Gardner et al. 2006). In fact, JWST has already started to revolutionize our understanding of galaxy formation at the epoch of reionization (EoR) redshifts (e.g. Bradley et al. 2022; Donnan et al. 2023; Atek et al. 2023; Adams et al. 2023). Nevertheless, even though perfect knowledge of the reionization timeline would help to model the EoR, it would still be non-trivial to translate these constraints into a single mapping of reionization astrophysics.  \nConsequently, it is intriguing to ponder the impact of different reionization scenarios in our SNR. We limit our study of reionization scenarios to three distinct reionization timelines. Our fiducial case aligns with Planck's timeline of reionization 𝑧 mid = 𝑧 ( 𝑥 H/i.pc = 0 . 5 ) ≈ 7 . 7, besides its duration, defined as Δ 𝑧 = 𝑧 ( 𝑥 H/i.pc = 0 . 10 ) -𝑧 ( 𝑥 H/i.pc = 0 . 90 ) is 3 . 83. Furthermore, we complement this model with scenarios representing later and earlier reionization scenarios, having midpoints at 𝑧 mid = 6 . 88 and 8 . 41, and durations of Δ 𝑧 = 4 . 15 and 3 . 59, respectively. Note that these additional models are roughly consistent with the 1 𝜎 error reported by Planck Collaboration et al. (2020) but we caution that the best fit from Planck is not able to reproduce astrophysical constraints in a satisfactory way - see Figure 4 of Montero-Camacho et al. (2024a) - due to the use of the hyperbolic tangent reionization model (Lewis 2008). A more adequate choice would be to consider the range of reionization profiles allowed in the Gompertz model of reionization. Its best fit to CMB data, in conjunction with astrophysical constraints, indicates a midpoint of reionization at ≈ 7 (Montero-Camacho et al. 2024a), more closely aligned with our late reionization scenario.  \nWe obtain our three reionization models by modifying the ionization efficiency 𝜁 in 21/c.pc/m.pcFAST. This parameter governs the overall timeline of reionization, although there are degeneracies (Pober et al. 2014; Park et al. 2019; Montero-Camacho & Mao 2021). Physically, it quantifies the ability of photons to escape their parent galaxies and reach the intergalactic medium (Dayal & Ferrara 2018).  \nTable 1. Summary of the total SNR at 𝑧 = 3 . 61 for the different scenarios considered in this work. The Fiducial scenario corresponds to 100 hours, 100 square degrees of survey area, and a Planck-like reionization scenario. We highlight that increasing the survey area by tenfold results in a marginal increase to the SNR for the baseline instrumental setup (see §5.2).  \nIn Figure 5, we illustrate the effect of including the memory of reionization in 𝑃 21 ,𝐹 compared to neglecting its existence. Naturally, the influence of the remnants of reionization is more pronounced in the highest redshift bin and for the delayed reionization scenario. As time progresses, gas dissipates the additional injected energy during reionization, resulting in a lesser effect in the Ly 𝛼 forest. Meanwhile, the modulation of the baryons in shallow potential wells becomes more subdued due to the growth of affected galaxies, as a result, the significance of the memory of reionization in H /i.pc 21 cm gradually diminishes as well.  \nAs expected, the primary gains in SNR occur at wavenumbers associated with the reionization-bubble scales, and they are minor at small scales. Note that the maximum importance overlaps with the best window in terms of SNR, as shown in Figure 2. Hence, the crosscorrelation of Ly 𝛼 and 21 cm intensity mapping holds significant promise for shedding light on the astrophysics of reionization and it could be used to indirectly constrain the timeline of reionization through its impact in the post-reionization intergalactic medium. Future work will investigate what could be gained from such efforts (see §6 for a preliminary forecast).  \nFurthermore, Figure 5 demonstrates that the impact of inhomogeneous reionization on the post-reionization IGM will increase the SNR for all timelines of reionization and redshifts considered here. In fact, the enhancement can reach 30% at its peak. However, we underscore that there is a fundamental competition in 𝑃 21 ,𝐹 between the response of underdense regions (forest) and that of denser regions (intensity mapping) to the reionization process. It is plausible that the overall effect would be a reduction of the signal-to-noise at redshifts close to the tail end of reionization. Fortunately, further scrutiny of this hypothesis will not require high-mass resolution simulations like the ones used throughout this work since the role of the small-scale structure should be subdued at these redshifts. Nevertheless, the feasibility of such a measurement is at best uncertain given the current capabilities of spectroscopic instruments and the sparse density of quasars at such high redshifts.  \nWe summarize our findings regarding observation strategy and dependence with reionization history in terms of signal-to-noise ratio in Table 1.  \nFigure 6. Forecast for the ability of SKA1-Low × DESI-like to constrain 𝜎 8 and 𝑛 𝑠 using 100 sq. degrees and 100 hours of integration time. Shown are the confidence ellipses - 1 𝜎 and 2 𝜎 - for different strategies where we consider the case of no reionization (I), i.e. the conventional signal; impact of reionization with known reionization (II), and impact of reionization with uncertain reionization (III). Notably this configuration, which uses only three redshift bins, is not very powerful compared to other cosmological programs expected to be operational by the end of the 2020s; however, it still showcases how reionization relics will jeopardize our ability to constrain cosmology from this observable.  \n<!-- image -->  \n<!-- image -->  \nFigure 5. Impact of reionization history in the signal-to-noise ratio of the cross-correlation assuming the SKA1-Low × DESI baseline configuration, i.e. 𝑡 int = 100 h and 𝑆 survey = 100 deg 2 . Shown are the SNR for a late reionization (purple squares), a Planck-like reionization (orange diamonds), and an early reionization scenario (green pentagons). The behavior at small 𝑘 is driven by Eq. (12).  \n<!-- image -->", '6 A SIMPLE FORECAST': "Having explored the potential gains in SNR due to different strategies and instrumental setups, we now turn our attention to quantifying the potential cosmological gains using the Fisher Matrix formalism (Heavens 2009).  \nSince no other work has looked at the cross-correlation of Ly 𝛼 × 21 cm at this redshift range, we will forecast based on three different scenarios: (i) Conventional signal, i.e. no memory of reionization, (ii) memory of reionization with perfect reionization knowledge, and (iii) marginalizing over the memory of reionization due to uncertainty regarding its timeline/modeling. Hence, we can illustrate the potential advantage gained by the reionization relics but simultaneously demonstrate the dangers of ignoring this effect. Furthermore, note that our strategies do not include any sort of priors from other datasets 6 , and hence future forecasts are likely to improve in these projections.  \nThe Fisher matrix is given by  \n𝐹 𝛼𝛽 = 𝑧 -bins ∑︁ 𝑖 𝑘 -bins ∑︁ 𝑗 𝜎 -2 𝑧 𝑖 𝑘 𝑗 [ 𝑃 21 ,𝐹 ] 𝜕𝑃 21 ,𝐹 𝜕𝜃 𝛼 ( 𝑧 𝑖 , 𝑘 𝑗 ) 𝜕𝑃 21 ,𝐹 𝜕𝜃 𝛽 ( 𝑧 𝑖 , 𝑘 𝑗 ) , (24)  \nwhere we have assumed uncorrelated errors between the different bins and that the posterior distribution can be reasonably welldescribed by a Gaussian. For simplicity, we consider only two cosmological parameters to forecast 𝜃 𝑖 = { 𝜎 8 , 𝑛 𝑠 } in (i) and (ii) . For (iii) , the conservative point of view on the progress of the field, we add the ionization efficiency 𝜁 to 𝜽 .  \nThe uncertainty in Eq. (24) is obtained from Eq. (19). Besides, We use the same three redshift bins centered at 3.61, 3.77, and 3.94 used throughout the rest of this paper. For the wavenumbers, we consider twenty logarithmic bins from 0.01 to 1 Mpc -1 .  \n<!-- image -->  \nn  \nTo compute the derivatives in Eq. (24), we run additional simulations that cover deviations of 3% in 𝜎 8 and 𝑛 𝑠 around their fiducial values of 0.8159 and 0.9667, respectively. Meanwhile, for 𝜁 we also consider 3% variation around 𝜁 = 30 but we check the convergence of our choice for the ionization efficiency in Appendix A. Furthermore, we illustrate some of the tendencies for different properties of this simulation suite in Appendix B.  \nWe plot the results of the SKA × DESI-like (PUMA × MUSTlike) forecast in Figure 6 (Figure 7) for the three different strategies  \nTable 2. Projected errors for the Fisher forecast of §6. Shown are the 1 𝜎 errors obtained for two instrumental configurations (SKA1-Low × DESI-like and PUMA × MUST-like) and for three different scenarios that correspond to no reionization (I), perfect knowledge of reionization/no marginalization over EoR astrophysics (II), and marginalization over ionization efficiency (III), respectively. Both instrumental setups have 100 hours of integration time and 100 sq. deg. of sky coverage.  \nPUMA × MUST  \n8  \nFigure 7. Similar to Figure 6, but for PUMA × MUST-like. In contrast to Figure 6, this setup places competitive constraints on the cosmological parameters. Therefore, our results emphasize the importance of considering the impact of reionization to prevent miscalculations of errors and the potential introduction of biases. Note that the combination with other datasets and the use of more than three redshift bins will significantly enhance the real constraining power expected from these instruments.  \n<!-- image -->  \nconsidered in this work. The confidence ellipses correspond to the 1 and 2 𝜎 contours for the 𝑛 𝑠 × 𝜎 8 plane. In addition, we tabulate the projected errors in Table 2.  \nBoth instrumental configurations exhibit a trend of greater ability to constrain the tilt of the primordial power spectrum when the impact of reionization is accounted for and known but when marginalization over the ionization efficiency is required the constraint power for the tilt degrades for the PUMA × MUST-like configuration. This trend of better 𝑛 𝑠 constraints was also present in the forecast for Ly 𝛼 forest done by Montero-Camacho & Mao (2021). As seen in their Figures 9 and 10, which include the memory of reionization in the 3D flux power spectrum from the Ly 𝛼 forest, the tilt plays a more important role than the amplitude. We attribute this to the increase (or decrease) of faint galaxies that would happen by modifying the value of 𝑛 𝑠 while a similar increase in 𝜎 8 will affect the environment more uniformly, hence influencing the ionizing sources to a lesser degree than the tilt. Note that reionization is driven by faint galaxies  \nin our 21/c.pc/m.pcFAST simulations. The larger significance of the impact of the tilt is also shown in Figure B1 where the impact of changing 𝑛 𝑠 and 𝜎 8 in the reionization history is illustrated.  \nIn the absence of reionization relics, the cross-correlation - without any external information from other cosmological probes - seems to constrain 𝜎 8 more than 𝑛 𝑠 . Interestingly, this is the opposite behavior of the results for H /i.pc 21 cm intensity mapping auto power spectrum albeit with 𝐴 𝑠 instead (see the 1 𝜎 error in Table 4 of Long et al. 2023). Disregarding the difference of 𝜎 8 with 𝐴 𝑠 , this trend could be due to the use of Planck priors in the auto-correlation, the additional parameters (like Ω 's that will likely be degenerate with 𝐴 𝑠 and 𝜎 8 ), or perhaps due to the anti-correlated nature of the crosscorrelation since the presence of extra structure would inevitably lead to more absorption in the Ly 𝛼 forest. This extra absorption or lack of could play a significant role in the range where the SNR is not dominated by the 21 cm error.  \nThe Ly 𝛼 × 21 cm cross-correlation will eventually become a competitive probe of the astrophysics that governs cosmic reionization. However, the constraining power for 100 hours of integration time is too weak to offer real insights into the reionization process. This research direction is likely to become promising with a more ambitious cross-correlation program, say 1000 hours in a SKA × MUST-like setup. Unexpectedly, even a conservative 100-hour cross-correlation survey will be highly competitive once PUMA (Cosmic Visions 21 cm Collaboration et al. 2018) starts taking data. For reference, 100 hours of integration time would greatly improve on the projected error for a full DESI (5 years) constraint using the 3D flux power spectrum ( 𝜎 𝜁 = 11 . 6, Montero-Camacho & Mao 2021) by roughly a factor of 3 . Its constraining power would be of similar strength to that obtained by demoting the optical depth to reionization to derivedparameter using symbolic regression and CMB data in conjunction with astrophysical data (Montero-Camacho et al. 2024a).", '6.1 Constraint on the timeline of reionization': "Here, we use the results of our PUMA × MUST-like Strategy III forecast - marginalize over reionization astrophysics - to constrain the timeline of reionization, assuming our Fiducial model correctly represents the Universe's neutral hydrogen fraction. In essence, we assume 𝜁 = 30 and we construct the range of allowed 𝑥 H/i.pc based on our forecasted error 𝜎 𝜁 .  \nFor clarity, we decided to only consider the PUMA configuration, which is significantly stronger than the SKA1-Low × DESI-like constraint (the error is smaller by almost a factor of 10). This stronger constraining power also justifies our choice of only considering 𝜁 due to the very small errors for 𝜎 8 and 𝑛 𝑠 expected from this instrumental configuration.  \nWe plot our findings in Figure 8 along with the current state of direct and indirect observations (i.e. optical depth constraint from the cosmic microwave background).  \nWe emphasize that the width of our inferred constraint is the true value added by our forecast. Our fiducial value, based on the default in 21/c.pc/m.pcFAST/v.pc3 and aligned with Planck's inferred reionization history(Planck Collaboration et al. 2020), provides a reasonable description of the Universe's reionization timeline. However, recent advancements in the more reliable quasar damping wing observations (see e.g. Ďurovčíková et al. 2024; Greig et al. 2024) suggest a slight tension with both our fiducial model and the Planck constraint. Notably, these developments also support the Gompertz reionization scenario (Montero-Camacho et al. 2024a), which successfully fits both CMB data and astrophysical constraints on reionization. Future  \nFigure 8. Inferred constraint on the reionization history by our PUMA × MUST setup using 100 hours of integration and assuming a 100 sq. deg. of overlap. Also shown are current observational constraints on the timeline of reionization including dark pixel fraction (McGreer et al. 2015; Jin et al. 2023), high-redshift galaxies through their clustering, luminosity evolution (LF) and equivalent width (EW) (Ouchi et al. 2010; Sobacchi & Mesinger 2015; Mesinger et al. 2015; Mason et al. 2018; Hoag et al. 2019; Mason et al. 2019; Morales et al. 2021), and high𝑧 quasars damping wings (DW) (Greig et al. 2022, 2024; Spina et al. 2024; Ďurovčíková et al. 2024). We also include the Planck indirect constraints on the timeline of reionization (Planck Collaboration et al. 2020), which relies on a hyperbolic tangent to parameterize the reionization process. In addition, we include the best fit using Gompertz reionization (Montero-Camacho et al. 2024a), which uses Planck data but does not rely on the hyperbolic tangent parametrization and provides a better fit to astrophysical observations. Note that we have anchored the constraint around the fiducial model, thus the key feature is the width of the constraint, rather than its exact location.  \n<!-- image -->  \nwork should consider a fiducial model that better aligns with the Gompertz timeline.", '7 SUMMARY': 'The cross-correlation between the Ly 𝛼 forest and 21 cm intensity mapping in the post-reionization era is a promising cosmological probe of the relatively high-redshift intergalactic medium. In particular, it can probe smaller scales before they become fully non-linear (compared to that of traditional galaxy surveys) and it is pragmatically an easier measurement than the auto-correlation of any of those fields (at 𝑧 ∼ 4) due to foregrounds or available quasars line of sight. However, just as the Ly 𝛼 forest (Montero-Camacho et al. 2019) and 21 cm intensity mapping (Long et al. 2023) are sensitive to relics from cosmic reionization, their cross-correlation will also be biased unless appropriate care is taken to handle this broadband effect.  \nRegardless of the impact of reionization, our results demonstrate the importance of overlap between radio interferometers and Stage V spectroscopic instruments. Since some of these telescopes are currently in the early planning/design stage, we underscore the significance of guaranteeing a small degree of overlap between the instruments. As shown in Figure 4, even 100 square degrees of overlap could lead to a detection and consequently, it would enhance observational programs aimed at the post-reionzation era.  \nFurthermore, we found that the gain in signal-to-noise is small for increased survey area in instrumental setups that use SKA1Low. This trend is caused by the increase in the SKA1-Low thermal noise. Nonetheless, the PUMA × MUST setup does exhibit significant gains, particularly at large scales with increased survey area,  \nsince the error is not dominated by thermal noise. In contrast, a longer integration time leads to considerably better measurements across the board, particularly for SKA configurations. Thus, we conclude that our baseline scenario, which we consider our cheap option with 100 hours of integration time and 100 square degrees of overlap, is a wellsuited observational setup although the measurements can likely be improved by further adjusting the survey strategy. For instance, a SKA1-Low × Stage V spectroscopic instrument is likely to perform comparatively to a Stage II radio interferometer at large scales with 1000 hours of integration time (see Figure 3). We emphasize that this result highlights the importance of progressing further than Stage V spectroscopic instruments in the 2040s.  \nThe inclusion of reionization relics in the cross-correlation increases the strength of the signal, particularly at high𝑧 and large scales. Interestingly, late reionization produces an enhancement of up to 30% over the predicted level without reionization imprints at 𝑧 ≈ 4. Separating this novel effect from the cosmological information - for instance using physics-inspired templates (Montero-Camacho et al. 2023) - will allow for unbiased inference of cosmological parameters and it would unseal an original methodology to investigate the astrophysics of reionization. Future work will focus on mitigation strategies for this observable.  \nWe have demonstrated the expected impact of reionization in the inference of cosmological parameters using a Fisher forecast. Cosmological parameters would be biased if one neglects the reionization relics and there would also be a significant underestimation of the error. However, if our knowledge of the reionization timeline improves significantly, the inclusion of the memory of reionization can result in a stronger constraint (as was the case for our PUMA × MUST forecast in Figure 7). Furthermore, the cross-correlation of the Ly 𝛼 forest and H /i.pc 21 cm IM will be a promising probe of reionization astrophysics in the next decades.  \nOur findings should also be interpreted as a cautionary tale - but simultaneously an exciting opportunity - to other post-reionization era high𝑧 ( 𝑧 ≳ 3) tracers, like Ly 𝛼 emission (e.g. Croft et al. 2016; Renard et al. 2021, 2024) and the CO rotational transitions (Breysse et al. 2022; Chung et al. 2022). Future work will assess the impact of reionization relics in other cosmological tracers of the post-reionization era.', 'ACKNOWLEDGEMENTS': "We are grateful to Chris Hirata for his helpful suggestions and comments. This work was supported by The Major Key Project of PCL. We acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University and PCL's Cloud Brain for providing computational and data storage resources that have contributed to the research results reported within this paper. This work made extensive use of the NASA Astrophysics DataSystem and the following open-source python libraries/packages: matplotlib (Hunter 2007), numpy (Harris et al. 2020), and scipy (Virtanen et al. 2020).", 'DATA AVAILABILITY': 'The data underlying this article will be shared on reasonable request to the corresponding authors.', 'APPENDIX A: DEPENDENCE OF FORECAST ON 𝛿𝜁': 'The ionization efficiency 𝜁 controls the timing of reionization in 21/c.pc/m.pcFAST by parametrizing the difficulty for an ionizing photon to escape into the IGM. Consequently, this parameter is not as well constrained or studied as cosmological parameters like 𝜎 8 or 𝑛 𝑠 . Here we confirm that the step size of this uncertain parameter does not have a significant impact on our results. We follow Facchinetti et al. (2023) strategy of checking the impact on the forecasted error instead of the performance of the numerical derivative.  \nFigure A1 highlights the small difference between the forecasted 𝜁 errors and justifies our choice of 3% variation in ionization efficiency for our Fisher forecast.', 'APPENDIX B: TRENDS FOR THE FISHER SIMULATIONS': 'Here we show some reionization-related properties of the simulations used in §6 to clarify their hierarchy.  \nFor instance, it is clear that 𝑛 𝑠 will be the most constrained parameter by the inclusion of reionization relics since it dominates most of the different metrics in Figure B1. In contrast, a weaker constraint in ionization efficiency is also expected, particularly since ionization efficiency does not affect the first term in Eq. (1).  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \n<!-- image -->  \nFigure A1. Convergence of the ionization efficiency error with step size. Shown are the forecasted errors for a variation of 3% (blue pentagons), 10% (purple diamonds), and 20% (green triangles) of the fiducial 𝜁 . Given the small difference between the forecasted errors, we opt for the 3% variation in ionization efficiency.  \n<!-- image -->  \nFigure B1. Properties of the simulations used in the Fisher forecast of section 6. (Top left) Reionization history. (Top right) Cross-correlation of matter and neutral hydrogen field as a function of redshift and evaluated at a 𝑘 representative of the ionized bubble scale. (Middle left) Relative transparency of the IGM 𝜓 as a function of local reionization redshift, Eq. (5). (Middle right) Relative neutral hydrogen density Ξ , Eq. (6) (Bottom left) The cross-power spectrum of matter and transparency of the model divided by that of the fiducial simulation, i.e. Eq. (7) divided by the reference scenario. (Bottom right) Ratio of the cross-power spectrum of matter and relative neutral hydrogen density, i.e. Eq. (8) normalized by the reference scenario.  \n<!-- image -->'}