Lyα at Cosmic Dawn with a Simulated Roman Grism Deep Field

The slitless grism on the Nancy Grace Roman Space Telescope will enable deep near-infrared spectroscopy over a wide field of view. We demonstrate Roman's capability to detect Lyα galaxies at z > 7 using a multiple position angle (PA) observational strategy. We simulate Roman grism data using a realistic foreground scene from the COSMOS field. We also input fake Lyα galaxies spanning redshift z = 7.5–10.5 and a line-flux range of interest. We show how a novel data-cube search technique—CUBGRISM—originally developed for the Galaxy Evolution Explorer can be applied to Roman grism data to produce a Lyα flux-limited sample without the need for continuum detections. We investigate the impact of altering the number of independent PAs and exposure time. A deep Roman grism survey with 25 PAs and a total exposure time of 70 hr can achieve Lyα line depths comparable to the deepest z = 7 narrowband surveys (L Lyα ≳ 1043 erg s−1). Assuming a null result, where the opacity of the intergalactic medium (IGM) remains unchanged from z ∼ 7, this level of sensitivity will detect ∼400 deg−2 Lyα emitters from z = 7.25 to 8.75. A decline from this expected number density is the signature of an increasing neutral hydrogen fraction and the onset of reionization. Our simulations indicate that a deep Roman grism survey has the ability to measure the timing and magnitude of this decline, allowing us to infer the ionization state of the IGM and helping us to distinguish between models of reionization.


INTRODUCTION
Lyα emission is a bright spectral feature in starforming galaxies that many observational surveys have targeted to efficiently obtain large samples of Lyα emitting galaxies over a wide redshift range z = 2 − 7 (e.g., Cowie & Hu 1998;Rhoads et al. 2000;Malhotra & Rhoads 2002;Ouchi et al. 2003;Gawiser et al. 2007;Gronwall et al. 2007;Blanc et al. 2011;Finkelstein et al. 2013;Konno et al. 2014;Matthee et al. 2015;Tilvi et al. 2020).With the arrival of JWST, Lyα continues to play an essential role in the study of high-redshift galaxies and the ionization state of the intergalactic medium (IGM) (Endsley et al. 2022;Tang et al. 2023;Bunker et al. 2023;Saxena et al. 2023).The most exciting JWST Lyα-observation to date is the discovery of an emitter at z=10.6 (Bunker et al. 2023).This object shows that Lyα can find a way to escape only 430 Myr after the Big Bang, where even extremely early reionization models (e.g., Finkelstein et al. 2019) predict a neutral hydrogen fraction of ∼ 80%.However, the current Lyα JWST observations are limited in volume and are pre-selected based on broadband criteria.We need a deep, blind search for Lyα over a large volume to realize the full potential of Lyα-based reionization constraints.
We know that reionization should fall within the 6 < z < 9 redshift range from the saturation of Lyα absorbers in z ∼ 6 quasar spectra (Fan et al. 2006) and by polarization measurements of the cosmic microwave background (Planck Collaboration et al. 2018).To go beyond this rough picture, many groups have used the redshift evolution of the number density of Lyα emitters (LAEs) to constrain the timing of reionization.Some groups have found a large drop in their density at z ∼ 7 and have attributed this decline to the increasing opacity of the IGM and the onset of the epoch of reionization (e.g., Stark et al. 2010;Konno et al. 2014;Hoag et al. 2019), while other groups have found a more modest evolution over this same redshift range consistent with a largely ionized universe (e.g., Tilvi et al. 2010;Itoh et al. 2018;Zheng et al. 2017;Hu et al. 2019;Wold et al. 2022).Some of the variation in the results may be due to spatially inhomogeneous reionization since cosmological reionization is thought to have proceeded through the growth and eventual overlap of ionized bubbles, driven by ultraviolet photon production in the earliest galax-ies.This produces a non-uniform reionization history that necessitates observational studies with large survey volumes to mitigate cosmic variance.
The current variation and uncertainty in the observational constraints allow for a diverse set of reionization histories.These reionization scenarios range from a rapid neutral-to-ionized transition dominated by relatively massive galaxies (Naidu et al. 2020) to a gradual transition dominated by more numerous low-mass galaxies (Finkelstein et al. 2019).The redshift where these reionization models differ the most is z ∼ 8.For example, at this redshift Naidu et al. predict an almost fully neutral universe, while the Finkelstein et al. model predicts a much lower neutral fraction of 35%.Moving to lower (higher) redshifts the models converge to a fully ionized (neutral) IGM, making z ∼ 8 an ideal redshift to distinguish between models and tell us about the sources of ionization radiation -whether large galaxies or small.
However, it is difficult to obtain z ∼ 8 Lyα-based observational constraints from ground facilities largely due to the increasing night sky background and the need for expensive spectroscopic followup of Lyα candidates to eliminate contamination.To effectively study reionization at z ∼ 8, we need to acquire greater statistical leverage by obtaining spectroscopically complete samples of LAEs over larger volumes without the interference of the Earth's atmosphere.
In this paper, we investigate Roman's potential to obtain a LAE sample at z 7 with a multi-positionangle grism deep field.Roman's Wide Field Instrument (WFI) has visible to near-infrared imaging capability and a slitless grism capability spanning a wavelength range of 1.00 − 1.93 µm with an ∼ 11 Å per pixel dispersion.The WFI has a field of view of 0.281 square degrees or over 100 times larger than similar instruments on HST and JWST.This large field-of-view (FOV) is tiled by 18 detectors in 3 rows and 6 columns with a 0.11 pixel scale.It is the combination of Roman's wide field of view and grism wavelength range that will clear the way for z > 7 LAE searches with unprecedented volumes -a volume of ∼ 3 × 10 6 Mpc 3 from z = 7.25 to 8.75 in a single pointing.
Slitless grism capabilities are well suited for spacebased telescopes where the sky background is low, and where slitless modes are often the simplest way to add spectroscopy to a mission.Furthermore, slitless observations are extremely efficient at spectroscopic surveys because all objects within the field-of-view have their light dispersed across the detector.This is particularly a desirable feature for finding strong emission line sources that are easily missed in slit-mask observations which often select their targets based on continuum brightness (Malhotra et al. 2005;Tilvi et al. 2016;Larson et al. 2018).
These strengths also present challenges.In particular, neighboring sources aligned in the dispersion direction can have overlapping spectra.This means that a highredshift slitless grism survey will need to contend with contamination caused by the foreground scene of stars and galaxies, and a simple exposure time calculation will fail to capture the true survey performance given this source confusion.For these reasons, we simulate a realistic extra-galactic foreground scene that captures the key characteristics of a grism deep field and then assess our ability to conduct a blind search for Lyα galaxies at cosmic dawn.We show how a novel data cube search technique -CUBGRISM -originally developed for GALEX grism data can be applied to Roman data allowing us to isolate z > 7 Lyα galaxies.

GRISM SIMULATIONS
We use HST imaging and best-fit SEDs in the COS-MOS field to define our simulated foreground stars and galaxies.We also construct fake high-redshift (7.5 < z < 10.5) Lyα emitters with a variety of intrinsic characteristics and then measure our ability to detect these emitters with different observational strategies.In the following sections, we describe the construction of our grism simulations in detail.

Morphology
We define the morphology of simulated objects with image postage stamps.We use the 3D-HST F160W mosaic to make these images for all of our simulated foreground sources.For each object, we extract a subimage from the science F160W image centered on the object's coordinates, where all pixels not assigned to the object via the segmentation map are zeroed out.To prepare these sub-images for input, we resample the F160W mosaic's drizzled pixel scale of 0.06 to match Figure 1.Roman WFI grism response functions for the in-focus spectral orders.For display purposes, we increased the 2-2 and 0-0 grism response functions by a factor of ten.All three in-focus spectral orders are included in our grism simulations.
Roman's pixel scale of 0.11 .Note, the native pixel scale of WFC3-IR is about 0.13 /pixel, so we're producing a final model image with a pixel scale very near the actual input data.The F160W mosaic has a median 5σ depth of 25.8 mag and an average FWHM of 0.19 (Skelton et al. 2014).Our foreground simulations include all F160W objects with magnitudes brighter than 25.5 mag.This corresponds to the expected 1σ per pixel sensitivity at 1.5µm for our deepest simulated Roman grism survey with a total t exp = 69.4hrs 1 .

Spectra
We use the best-fit EAZY SEDs from Skelton et al. (2014) to define the spectra of all the simulated foreground objects.These spectral templates were fit to 44 bandpasses with 13 spanning our simulated 1 to 2 µm wavelength range.Skelton et al. find excellent agreement with their photometric redshifts and the available spectroscopic redshifts, with a scatter of σ NMAD = 0.007 and an outlier percentage of 1%.This level of agreement suggests that their templates SEDs, which we use directly, are able to accurately represent the spectral data.
For our Roman grism simulations, we are most concerned with reproducing a realistic foreground scene for the 1 to 2 µm wavelength range, and we normalize the EAZY SEDs to the observed F160W photometry to help ensure that our simulated grism image has a realistic distribution of foreground flux.
For the Lyα spectral shape, the continuum has a spectral slope of −2 such that f λ ∝ λ −2 and is attenuated by the intervening IGM.The Lyα line is a Gaussian with rest-frame FWHM= 1 Å or 250 km s −1 which is consistent with observational constraints (Hu et al. 2010).We insert a total of 288 Lyα emitters into our foreground grism scene.To ensure that all position angles contain our simulated emission lines, we randomly place our Lyα emitters within a radius of 1 arcmin of the field center.

Simulation Construction
We use the aXeSIM (Kümmel et al. 2009) software package to simulate 25 independent position angles over a 2048×2048 pixel or 14.1 arcmin 2 field of view.aXeSIM was originally designed to simulate grism and direct images for the Hubble Space Telescope and a number of modifications were needed to simulate Roman data.We modified aXe's instrument parameter file to match Roman's trace and wavelength dispersion.We simulated a wavelength range of 9500 − 19000 Å and assumed a constant wavelength dispersion of 265.7, 10.8, and 5.6 Å pix −1 for orders 0-0, 1-1, and 2-2, respectively.We set the zero deviation wavelength of the 1-1 order to be 0.95µm, which differs from Roman's undeviated wavelength of ∼ 1.55µm.Given the field of view of the input scene relative to our smaller simulated region, we do not expect this difference to impact our measurement of recovery fractions for simulated Lyα emitters.
We used the morphological and spectral inputs described in the previous sections to run three separate aXe simulations to form grism images for Roman's 0-0, 1-1, and 2-2 spectral orders.To further facilitate parallelization and to allow for alterations to our Lyα emitter sample, we simulated the grism foreground scene separately from the simulated LAE grism images.: 0-0, 1-1, and 2-2 are simulated independently and then co-added to form our combined grism image.Finally, we add in the Poisson noise appropriate to our desired exposure time (see Figure 3).For the four panels above, we show a 3.75 × 3.75 field of view.
In Figure 1, we show the Roman grism response functions that we used in our aXe simulations.Compared to the science 1-1 order, the 0-0 and 2-2 orders have less-sensitive response functions.However, these orders can significantly complicate the foreground scene.The 0-0 order is compact with the full 1µm grism range only occupying tens of pixels.A compact 0-0 order from a bright foreground object may be mistaken as an emission line, and we wish to verify that our LAE search does not select these potential contaminants.The 2-2 order is spread out over ∼ 2000 pixels, about twice the length of the science order.This means that bright objects well outside of the simulated FOV can have their 2-2 light dispersed into the grism image and contribute to the foreground flux that we need to contend with when searching for high-redshift LAEs.
In Figure 2, we show an example of our simulated spectral orders and their summation for a position angle (PA) aligned with the x-axis or PA= 0 by aXe convention.Non-zero position angles are not naturally supported by the aXe simulation software.Thus, we rotated all objects' morphologies and positions about the field center to simulate position angles not equal to zero.
Once all three spectral order simulations were completed for both the foreground and LAE scene, we combined them and applied a Poisson noise level of 0.8 e − pix −1 s −1 which is the expected near-infrared background for the High Latitude Survey (see Figure 3).For our deepest simulated Roman grism survey, we explored an observational strategy consisting of 25 exposureseach with a unique PA -observed for 10 ks giving a total exposure time of 69.4 hours.The simulated PAs were selected to uniformly sample 360 degrees, while avoiding 180 degree reflections.This was done because PA 0 and PA 0 + 180 will have largely the same spectra that are aligned and overlapping, complicating our effort to deblend sources.
We also investigated two lesser surveys: one with 25 PAs and a total exposure time of 41.7 hours (25 × 6 ks) and another with a subset of 15 PAs and the same total exposure time of 41.7 hours (15 × 10 ks).For brevity, we refer to these three scenarios as 25 × 10 ks, 25 × 6 ks, 15 × 10 ks.

Simulation Caveats
Roman's wide field slitless spectroscopy with a wavelength dispersion of ∼ 11 Å per pixel over a 1 µm range is achieved with a compound grism with two diffractive surfaces.This grism assembly mitigates wavelengthdependent aberrations at the expense of producing additional unwanted spectral orders.Our simulations are designed to capture the grism characteristics most relevant for a deep LAE survey.In particular, we aim to assess Roman's ability to conduct a deep LAE survey given a realistic foreground scene.Accordingly, we simulate the spectral off-orders that have the highest surface brightness (0-0 and 2-2) and ignore out-of-focus spectral orders (1-0, 1-2, 0-1, 2-1).Additionally, we do not model field-dependent distortions of the spectral trace or dispersion which are expected to have a minor impact on the performance of a deep LAE survey.
Our simulations are designed to predict the needed ontarget time to perform a deep LAE survey, and we do not investigate different dithering/mosaicing strategies.For simplicity, all 25 simulated PAs are assigned the same field center maximizing the area that is covered by all PAs.

CUBGRISM
We use a novel data cube search technique -CUB-GRISM -originally developed for GALEX grism data to test our ability to recover high-redshift Roman LAEs.This technique does not require a broad-band detection and was constructed to produce a line-flux limited sample of emission-line objects.These features are welltailored to our goal of characterizing our ability to recover Lyα emitters given an assumed Roman grism deep field observing strategy.We have highlighted a foreground emission line galaxy (blue circle) and an artifact caused by the 0-0 order of a bright star (underlined in blue).See Section 3.1 for additional discussion of our continuum subtraction method.We demonstrate that our emission line search technique, CUB-GRISM, can effectively remove contamination from 0-0 spectra and recover a clean sample of emission line sources.

Continuum Subtraction
In preparation for the data cube construction and our emission line search, we subtracted continua from the grism images by applying a 1 × 21 pixel median filter with the 21 pixel dimension aligned with the spectral trace.The size of the median filter (21 pixels or 231 Å) is designed to not bias our search for LAEs which will have FHWMs a factor of 20 smaller.In Figure 4, we demonstrate our continuum subtraction method showing that the remaining sources are primarily emission lines (both simulated LAEs and foreground sources) and 0-0 spectra.

Cube Construction
In Barger et al. (2012); Wold et al. (2014Wold et al. ( , 2017)), we describe CUBGRISM (previously referred to as the data cube search) in detail and demonstrate our method for converting multiple slitless spectroscopic images into a three-dimensional data cube with two spatial axes and one wavelength axis.Here, we provide a brief overview of this process with emphasis on alterations needed for Roman grism images.
For each simulated grism exposure, we know the spectral trace and wavelength dispersion, and this allows us to extract a spectrum for each spatial pixel position, thus forming an initial data cube.A data cube constructed from a single slitless spectroscopic image may suffer from overlapping spectra caused by neighboring objects that are oriented in line with the dispersion direction.For our simulated Roman deep field, we have multiple position angles -one per exposure -and objects that overlap in one position angle are unlikely to overlap in another position angle.This allows us to disentangle overlaps by constructing data cubes for each exposure and then combining these initial data cubes, applying a 5σ cut to remove contamination from overlapping sources.In more detail, we take our library of data cubes (N = 25 for our deep survey) and for each spaxel location, we measure the variation in flux values, quantified by the normalized median absolute deviation.Any spaxel that deviates from the median by more than 5σ is eliminated and the remaining spaxels are used to compute the median flux value which is then used to populate the final data cube.
This procedure results in a final data cube that has a wavelength step of 10 Å and wavelength coverage of approximately 1.0-1.5 µm.We designed the wavelength slices to have a wavelength extent that matches the spec-Noiseless Direct Image 1.155 !" 1.277 !" 1.398 !" Continuum subtracted wavelength slices tral resolution of the Roman Space Telescope.We are interested in our ability to recover high-redshift Lyα emitters (z = 7.5-10.5),and accordingly, we artificially limit the wavelength range of the output data.

Line Search
The final Roman data cubes cover a wavelength range of 1.00 to 1.49 µm or a Lyα redshift range of z = 7.24−11.25.For each wavelength slice, we used SExtractor (SE; Bertin & Arnouts 1996) to identify > 7.5σ detections using SNR-optimised aperture photometry.We computed an optimal aperture size of D = 0.4 by determining the aperture that maximizes the median Lyα data-cube flux relative to the measured aperture noise (for details on SNR-optimal apertures, see Gawiser et al. 2006).As discussed in the next section, this level of significance (> 7.5σ) is needed to minimize false detections within the cube.

False detections
CUBGRISM enables a blind search for Lyα emitters -with no continuum detection required -at the expense of requiring a PA-rich observational strategy to accurately isolate emitter locations.We favor this continuum-independent reduction strategy because our science goals require us to compare Roman Lyα populations at cosmic dawn to blind Lyα narrow-band searches just outside the reionization epoch at z 7 (Ouchi et al. 2008(Ouchi et al. , 2010;;Shibuya et al. 2012;Konno et al. 2014Konno et al. , 2018;;Matthee et al. 2015;Santos et al. 2016;Ota et al. 2017;Itoh et al. 2018;Hu et al. 2019Hu et al. , 2021;;Tilvi et al. 2020;Wold et al. 2022).
However, for reduced PA observational strategies, a potential drawback of our technique is the inability to isolate source locations resulting in spurious detections.
CUBGRISM has been used successfully on GALEX grism deep fields observed with hundreds of PAs to identify populations of low-redshift Lyα emitters.These far-ultraviolet and near-ultraviolet surveys found a high optical-spectroscopic confirmation rate of 85% (Wold et al. 2014(Wold et al. , 2017)).In addition to a 4σ detection within data cube, these GALEX studies also visually inspected profile-weighted one-dimensional spectra to help eliminate spurious detections.For our Roman simulations, we did not implement a hard-to-quantify visual inspection step.Instead, we verify that our adopted CUB-GRISM extraction parameters result in a manageable contamination rate ( 15%) which would likely be improved with a visual inspection or machine learning step.
For a deep Roman grism survey, we aim to detect the Lyα number density decline caused by the increasing opacity of the IGM and the onset of reionization.As a baseline for determining the rate of false detections, we define the false detection fraction as F/(F + E) where F is the number of recovered spurious objects within the z Lyα = 7.25-8.75redshift range and E is the number of expected LAEs over the same redshift range assuming no number density evolution from the z ∼ 7 LAGER survey (e.g., Wold et al. 2022;Hu et al. 2019).Given our survey flux limits, we expect a 69.4 hr survey to detect N ∼ 400 LAEs while our 41.7 hr surveys should detect N ∼ 150 LAEs over a nominal deep grism survey area of 1 deg 2 (see Section 4 for details).
It is not possible to simulate a larger ∼ 1 deg 2 FOV because of the limited HST near-infrared imaging data a Survey flux limits are computed from the 50% completeness threshold for all simulated LAEs z = 7.5-10.5.b The false detection fraction is defined as F/(F + E) where F is the number of recovered spurious objects within the zLyα = 7.25-8.75redshift range and E is the number of expected LAEs over the same redshift range assuming no number density evolution from the z ∼ 7 LAGER survey.and our COSMOS observation-based simulations.Instead, we assume that our limited FOV is representative and perturb our simulations with different realizations of Poisson noise.For each realization, we find the number of false detections over a wavelength range of 1.00 to 1.19 µm, which corresponds to the wavelength range of observed z = 7.25-8.75Lyα emission.In total, we look for false detections over an effective survey area 350 arcmin 2 and consider all three survey scenarios: 25 × 10 ks, 25 × 6 ks, 15 × 10 ks.
For all three survey configurations, we find that a 7.5σ detection threshold can keep the false detection rate at 15%.In Table 1, we summarize these results with their 1σ Poisson errors.We caution that using N 15 PAs will result in reduced performance of our CUB-GRISM reduction strategy.Using the same procedure as employed with the 25 × 10 ks, 25 × 6 ks, and 15 × 10 ks surveys, we estimate that a 10 × 15 ks PA survey would result in a ∼ 70% false detection rate.Reduced performance with a PA-poor survey is expected because of CUBGRISM's method of disentangling overlapping spectra (described in Section 3.2) which requires many PAs to measure the uncontaminated flux at each x, y, λ position.
We note that our false detection results do not account for foreground emission-line contamination.We assume that a deep Roman grism survey will be located in a field that contains deep ancillary imaging data that are capable of efficiently identifying foreground objects.As performed in the GALEX grism and narrow-band Lyα surveys, the false detection rate can be further investigated by independently confirming a representative subset of identified Roman grism Lyα emitters.

RESULTS
Given our generated realistic foreground scene and our reduction strategy, we now determine the recovery fraction of simulated Lyα emitters as a function of line-flux, redshift, and number of position angles.In Figure 5, we demonstrate this Lyα search by highlighting a small 20 × 20 region of our deepest simulation that contains 5 LAEs with different line fluxes and redshifts.For each wavelength slice, we search for LAEs and record the Lyα recovery fraction as a function of line flux.
As discussed in Section 2.3, we focused our investigation on three different deep-field scenarios.
1. 25 PAs with a total exposure time of 69.4 hours (25 × 10 ks) 2. 25 PAs with a total exposure time of 41.7 hours (25 × 6 ks) 3. 15 PAs with a total exposure time of 41.7 hours (15 × 10 ks) For each scenario, we show the fraction of recovered LAEs as a function of line-flux in Figure 6.Following Negrello et al. (2013), we fit our recovery fractions with an error function where logf and γ are free parameters that control the location of the zero-to-one transition and the rate of this transition, respectively.We find that the γ parameter is very similar for all curves, and we fix its value to the 25 × 10 ks best-fit value of γ = 0.52.Based on these best-fit functions, we find that our three scenarios are 50% complete at a Lyα line fluxes of 8.2, 10.8, and 10.4 × 10 −18 erg s −1 cm −2 , respectively.This is consistent with the survey depth increasing proportionally to the square root of the total exposure time.In Table 1, we list our surveys' flux completeness with 1σ errors computed with a Monte Carlo (MC) procedure that perturbs our completeness data-points by Poisson random deviates.This MC procedure is used to create N = 10, 000 perturbed completeness curves and the measured variation in the survey completeness is used to estimate the 1σ error.We investigate Roman's Lyα survey depth as a function of redshift for our deepest 25 × 10 ks survey configuration.Binning by line flux and redshift reduce the number of simulated LAEs per bin by a factor of four when compared to binning by line flux alone.To reduce the Poisson error in our recovery fraction results, we have produced four Lyα samples while keeping our foreground scene the same.For each Lyα sample, we independently select Lyα parameters as described in Section 2.2.This includes randomly selecting their positions.For each run, the Lyα sample is added to the foreground scene and Poisson noise is added.We combine Lyα recovery results from all four runs to compute the final completeness curves.
In Figure 7, we show these redshift-dependent Lyα recovery fractions.From z = 7.5 to 10.5, we find that our survey is 50% complete at Lyα line fluxes of 13.6, 8.3, 7.5, and 6.9 × 10 −18 erg s −1 cm −2 , respectively.The decline in survey depth at z = 7.5 is due to the grism's response function (see Figure 1).From a wavelength of 1.16 to 1.03 µm or a redshift range of z Lyα = 8.5 to 7.5, the grism response declines by a factor of 1.6.
In the right panel of Figure 7, we show that the reduced cosmological dimming at z = 7.5 makes up for the decline in the grism response function allowing for the recovery of L Lyα > 10 43 erg s −1 LAEs across all redshift bins.In terms of Lyα luminosity, our simulated Roman survey is most sensitive at z = 8.5 with a 50% survey completeness at 7.4 × 10 42 erg s −1 .We summarize the 25 × 10 ks survey's redshift performance in Table 2.
Considering a nominal deep survey area of 1 deg 2 with a z = 7.25-8.75redshift range and assuming no number density evolution from the z ∼ 7 LAGER survey (Wold et al. 2022), we estimate a survey sample size of 395 LAEs.To compute this quantity, we consider three redshift slices: z = 7.5, 8.0, and 8.5 each with a z ± 0.25 extent.For each slice, we compute the Lyα number density by integrating the z ∼ 7 LAGER LF down to our redshift dependent 50% completeness limits, where the z = 8.0 completeness is interpolated from our z = 7.5 and 8.5 measurements.Given the slice volumes (where V Total = 1.1 × 10 7 Mpc 3 ), the total number of LAEs is found to be 395 within an area of 1 deg 2 and a z = 7.25-8.75redshift range.Any significant decline in this sample size is a signature of the increasing opacity of the IGM within the reionization epoch.To give an estimate of this decline, we note that from a redshift of z = 7 to 8 the number density of star-forming galaxies declines by a factor of ∼ 2 (as traced by the evolution of the UV LFs Bouwens et al. 2015), indicating that a 2 Lyα number density decline over the same redshift range may be attributed to increasing IGM opacity during the reionization epoch.
We note that there are examples of spectroscopically confirmed z 8 LAEs within our detection limits (e.g., Oesch et al. 2015;Zitrin et al. 2015;Larson et al. 2018;Tilvi et al. 2020;Jung et al. 2020), indicating that these sources exist and can be recovered with a large-volume survey.With JWST observations the number of spectroscopically confirmed z 8 sources are increasing rapidly with the highest redshift LAE currently at a redshift of z = 10.6 (Bunker et al. 2023).However, these existing z 8 LAE samples have a relatively small survey volume or are pre-selected based on photometric redshifts.
For example, Tilvi et al. (2020) found one L Lyα > 10 43 erg s −1 LAE at z = 7.7 over a ∼ 4.5 × 10 4 Mpc 3 narrowband survey.Naively scaling the implied z = 7.7 LAE number density by the volume of our nominal deep Roman grism survey, implies 1) a Roman sample size of 230 LAEs, 2) a high IGM Lyα transmission, and 3) a small neutral hydrogen fraction.Clearly, redshift evolution, cosmic variance, and Poisson fluctuations make this estimate highly uncertain.Considering just the 1σ Poisson uncertainty (Gehrels 1986), allows for Roman sample sizes ranging from 40 to 770 LAEs assuming a 1 deg 2 FOV and a z = 7.25-8.75redshift range.To overcome these uncertainties, we need a deep Roman survey to conduct a blind search for strong Lyα emit- ters over a large volume.This will allow us to make a direct comparison to the number density evolution prior (z ∼ 5.7 via ground-based NB Lyα surveys) and during the reionization epoch (z 7 via a deep Roman grism survey).
At a redshift of z = 8, models of reionization predict significantly different volume-averaged neutral hydrogen fractions, x HI .For example, at this redshift Naidu et al. (2020) predict an almost fully neutral universe with x HI = 87%, while the Finkelstein et al. (2019) model predicts a much lower neutral fraction of x HI = 35%.One of the main differences in these models is whether massive galaxies or more numerous low-mass galaxies account for the majority of ionizing flux.With a deep Roman grism survey, we will obtain observational constraints on the evolution of the z ∼ 8 Lyα LF and be able to distinguish between a rapid neutral-to-ionized transition dominated by relatively massive galaxies or a gradual transition dominated by more numerous lowmass galaxies.
A deep Roman survey will also allow us to study the shape of the Lyα luminosity function.Our completeness simulations indicate that a deep Roman grism survey can achieve Lyα line depths comparable to the deep-est z 7 NB surveys (e.g., Clément et al. 2012;Hibon et al. 2012;Shibuya et al. 2012;Konno et al. 2014;Itoh et al. 2018;Ota et al. 2017;Hu et al. 2019;Tilvi et al. 2020;Wold et al. 2022).Based on these z ∼ 7 Lyα surveys, this depth will allow us to measure the shape of the Lyα luminosity function at the reionization epoch.Previous studies have suggested the bright-end of the luminosity function may evolve less rapidly due to bright LAEs preferentially residing within large ionized bubbles (Matthee et al. 2015;Hu et al. 2021;Taylor et al. 2020Taylor et al. , 2021;;Ning et al. 2022) and a deep Roman grism survey will allow us to study the bright-end of the LF at z 8 for the first time.

SUMMARY
We investigated Roman's potential to conduct a blind search for Lyα galaxies at cosmic dawn using a multiposition-angle observational strategy.We produced a realistic Roman WFI grism foreground scene based on observational constraints from the COSMOS field.We also simulated Lyα galaxies spanning our redshift and line-flux range of interest.We showed how a novel data cube search technique -CUBGRISM -originally developed for GALEX grism data can be applied to Roman grism data to produce a Lyα flux-limited sample without the need for a continuum detection.Given our adopted reduction technique, we investigated the impact of altering the number of independent position angles and total exposure time.Our results indicate that a proposed deep Roman grism survey can achieve Lyα line depths comparable to the deepest z = 7 NB surveys, allowing us to study the evolution of Lyα populations and infer the ionization state of the intergalactic medium at cosmic dawn.

Figure 2 .
Figure 2. Construction of the noiseless simulated grism images.All the in-focus spectral orders: 0-0, 1-1, and 2-2are simulated independently and then co-added to form our combined grism image.Finally, we add in the Poisson noise appropriate to our desired exposure time (see Figure3).For the four panels above, we show a 3.75 × 3.75 field of view.

Figure 3 .
Figure 3. Left: Observed HST F160W image.Middle: Our simulated F160W image.Right: Our simulated Roman grism image with a position angle aligned with the x-axis.For our deep field simulations, we consider up to 25 independent position angles and up to 69.4 hours of total integration time.We use a noise level of 0.8 e − pix −1 s −1 which is the expected near-infrared background for the High Latitude Survey.

Figure 4 .
Figure4.We show a small 40 × 40 region of our simulated texp = 10ks grism image before and after continuum subtraction.We have highlighted a foreground emission line galaxy (blue circle) and an artifact caused by the 0-0 order of a bright star (underlined in blue).See Section 3.1 for additional discussion of our continuum subtraction method.We demonstrate that our emission line search technique, CUB-GRISM, can effectively remove contamination from 0-0 spectra and recover a clean sample of emission line sources.

Figure 5 .
Figure 5.A 20 × 20 region of our simulated Roman direct image and three selected wavelength slices from our data cube constructed from 25 Roman grism images.From left to right, we show the Roman direct image, a 1.155 µm data cube slice containing 1 simulated z = 8.5 LAE, a 1.277 µm data cube slice containing 2 simulated z = 9.5 LAEs, and a 1.398 µm data cube slice containing 2 simulated z = 10.5 LAEs.The simulated LAEs are circled in red.The number associated with each LAE indicates the line-flux of the source, where the Lyα flux is in units of 10 −17 erg s −1 cm −2 .By design continuum sources are not recovered in our data cube due to our background subtraction procedure (see Section 3 for details).

Figure 6 .
Figure 6.The fraction of recovered Lyα emitters as a function of line-flux.Error-bars show the 1σ Poisson error.The red, green, and blue curves show examples of how our survey completeness varies with different observational strategies assuming a 0.8 e − pix −1 s −1 background noise level typical of the High Latitude Survey (see Section 4 for details).

Table 1 .
Roman Grism Survey Performance

Table 2 .
Deep Roman Grism Survey Redshift Performance Same as the left figure except completeness is plotted versus Lyα luminosity.Across all redshift bins, we find that our survey can recover LLyα > 10 43 erg s −1 LAEs at the 50% completeness level (see Table2).