On the possibility of Baryon Acoustic Oscillation measurements at redshift $z>7.6$ with the Roman Space Telescope

BlueTides is a hydrodynamic simulation which was designed to cover the high redshift Universe (Feng et al., 2015, 2016). The model uses the Smoothed Particle Hydrodynamics code MP-Gadget to evolve a cubical volume with $2\times 7040^{3}$ particles. These particles are used to model dark matter, gas, stars and supermassive black holes. This simulation evolved a $(400/h\approx 577)^{3}$ cMpc³ cube to redshift $z=7.3$. The resolution, volume and number of galaxy halos in BlueTides are well matched to depth and area of current observational surveys targeting the high-redshift Universe. Additional details of the BlueTides simulation can be found in (Feng et al., 2015, 2016; Wilkins et al., 2016a; Wilkins et al., 2016b; Wilkins et al., 2017).

A Friends-of-Friends (FoF) algorithm was used to identify galaxy halos. At the high redshifts relevant here these halos are sufficiently accurate representations of the galaxies that would be found in the simulations if observational determinations were used. This was seen directly in Feng et al. (2016) by comparing FoF halos with galaxy halos extracted from BlueTides mock images and identified with the Source Extractor algorithm (Bertin & Arnouts, 1996).

Feng et al. (2016) used results from the BlueTides simulation to find good correspondence between the predicted galaxy luminosity function and observations in the observable range $-18\leq M_{\rm UV}\leq-22.5$ with some dust extinction required to match the abundance of brighter objects. Additionally, Feng et al. (2016) also find good agreement between the star formation rate density in BlueTides and current observations in the redshift range $8\leq z\leq 10$.

Wilkins et al. (2017) use galaxy halos from BlueTides simulation to show that intrinsic and observed luminosity functions at $z\in\{8,9,10\}$ exhibit the fast expected build up of the galaxy population at high-redshifts. The work presented in Wilkins et al. (2017) reveals that the intrinsic luminosity function is broadly similar to the observed luminosity function at faint luminosities ($M>-20$). At brighter luminosities, the luminosity function is shown to have a stronger steepening which indicates the increasing strength of dust attenuation. In their work, Wilkins et al. (2017) also demonstrate that a Schechter function provides a good overall fit to the shape of a the dust attenuated far-UV luminosity function for $z\in\{8,9,10\}$.

In Waters et al. (2016), it was shown that the number density of galaxies that will be observable by the RST High Latitude Survey should be high enough that BAO measurement at redshift $z\sim 8$ could be possible. The BlueTides volume being too small to make direct mock HLS observations and measurements of the BAO, it was necessary to follow up with different techniques, and this is the role of the current paper.

The BlueTides simulation is used in this paper in two roles: First, to make some example RST grism fields to check the degree of spectral overlap between high redshift galaxies and others at similar redshifts. Second, we use the BlueTides galaxies to populate the dark matter halos in a much larger dark matter only simulation. For the dark matter only simultion, we use the Outer Rim catalogue. Details of the Outer Rim simulation are described in  2.2.

The Outer Rim simulation (Heitmann et al., 2019; Li et al., 2019) is a large-scale dark matter only simulation which we use to measure large-scale clustering and BAO. It was run using the Hardware/Hybrid Accelerated Cosmology Code (HACC; Habib et al., 2013; Habib et al., 2016), on Mira, a BG/Q system at the Argonne Leadership Computing Facility. The cosmology used is a $\Lambda$CDM model close to the best-fit model from WMAP-7 (Komatsu et al., 2011), with parameters given in Table 1. The comoving box size of the simulation is $3000\ h^{-1}$Mpc, yielding a volume 400 times that of BlueTides. 1.07 trillion particles were evolved, leading to a particle mass of $m_{p}=1.85\times 10^{9}\ h^{-1}M_{\sun}$. The model was run to redshift $z=0$, but in this paper we use only outputs at $z=7$ and above.

Halos are identified with a Friends-of-Friends (Davis et al., 1985) halo finder with a linking length of $b=0.168$.

At the redshifts $z\sim 8$ which concern us here, we will simulate the Lyman-$\alpha$ line luminosity of galaxies, which will lie within the RST grism spectral range (see section 3.1 below). An example of Lyman-$\alpha$ line detection from a galaxy at these redshifts is given in the work of Lehnert et al. (2010). The redshift of the galaxy UDFy-381355395 is $z=8.6$, and Lehnert et al. (2010) measured a Lyman-$\alpha$ line flux of $(6.1\pm 1.0)\times 10^{-18}$ erg/s/cm², corresponding to a luminosity of $5.5\pm 1.0\pm 1.8\times 10^{42}$ erg/s.

In BlueTides, the star formation rate (SFR) for each galaxy is available to us, and we translate this quantity directly into Lyman-$\alpha$ luminosity using the following scaling equation (Cassata, P. et al., 2011).

In equation 1, $L_{{\rm Ly}\alpha}$ denotes the Lyman-$\alpha$ luminosity. This conversion factor is based on a stellar population with a Salpeter initial mass function (IMF) and is accurate to within a factor of a few for a range of population ages, high mass cutoffs of stars, and metallicities (Leitherer et al., 1999). Most relevant in our case is that there is no correction for effects like dust, escape fraction, and intergalactic or circumgalactic absorption. This means that using the raw star formation rate measurements from the simulation represents a best case for the detectability of galaxies and this should be borne in mind when interpreting the results. We will return to this in Section 5. Due to the paucity of observational data at these redshifts, the magnitude of these obscuring effects is not well known, but they are likely to be significant. In the example of the galaxy mentioned in Lehnert et al. (2010), the continuum SFR was estimated to be $2-4\ M_{\sun}\ {\rm yr}^{-1}$. At the same time, Lehnert et al. (2010) evaluated the SFR from the Lyman-$\alpha$ line to be $0.3-2.1\ M_{\sun}\ {\rm yr}^{-1}$, which leads to Lyman-$\alpha$ overall escape fraction from $8\%-100\%$.

The BlueTides simulation volume is not sufficiently large to make mock RST survey catalogues that can be used to measure large-scale structure. We use the halo catalogues from the Outer Rim simulation to do this, painting SFRs from BlueTides galaxies onto correspondingly massive Outer Rim halos. Fig. 1 shows the relationship between mass and SFR for BlueTides at redshift $z=8$, We can see that there is an approximately linear relationship between $\log(M)$ and $\log({\rm SFR})$. Our SFR assignment scheme is as follows: for each Outer Rim halo, we randomly pick a BlueTides halo with total mass within $5\times 10^{10}$ $h^{-1}M_{\sun}$, and assign the SFR for that halo to the Outer Rim object. In this way, we account for some aspects of the scatter between SFR and halo mass, although not any spatially correlated structure. For large mass galaxies we increase the range to make sure there are enough sample halos in BlueTides.

In one of the RST design documents, the $7\sigma$ line flux sensitivity for the High Latitude Survey is $10^{-16}\ {\rm erg/s/cm}^{2}$. ²²2 https://wfirst.gsfc.nasa.gov/science/sdt$\_$public/WFIRST-AFTA$\_$SDT$\_$Report$\_$150310$\_$Final.pdf We discuss RST and the flux limits on our mock surveys further in Section 3.1.1.

Fig. 2 shows the SFR distribution of BlueTides galaxies at $z=8$ with mass from $2\times 10^{11}\ h^{-1}M_{\sun}$ to $5\times 10^{11}\ h^{-1}M_{\sun}$. We can see that when the line sensitivity is $7\times 10^{-17}\ {\rm erg/s/cm}^{2}$, then the minimum detectable luminosity is $5\times 10^{43}\ {\rm erg/s}$ at $z=8$, which corresponds to $40\ M_{\sun}\ {\rm yr}^{-1}$. Only the most highly star-forming galaxies are therefore likely to be included in any BAO survey.

Working with grism spectroscopy is complex (McCarthy et al., 1999; Fossati et al., 2017; Colbert et al., 2018), with relatively low resolution possible, mixing of spatial and spectral structure in background and stray light, and overlap of sources. We leave full simulation of mock RST grism surveys (including foregrounds and interlopers) to other work (Gehrels et al., 2015; Colbert et al., 2018; Malhotra, 2019). We will however make some extremely simple grism images using BlueTides simulation data, in order to assess the level of spectral overlap from galaxies at similar redshifts.

The data we use from the Outer Rim simulation covers snapshots from $\text{z}=7.4$ to $\text{z}=10$. Our mock RST high-latitude survey (HLS) data, will cover four $3000\times 3000\times 700\ h^{-3}$Mpc³ cuboids drawn from the $3000^{3}\ h^{-3}$Mpc³ simulation volume. This corresponds to a sky area of $4\times 570=2280$ sq. deg, approximately equal to the total sky coverage of the RST HLS. The redshift coverage in the line of sight dimension is $\text{z}=7.4$ to $\text{z}=10$, with the survey volume covering $4\times 6.3=25.2$ comoving (Gpc/$h)^{3}$. We use the flat sky approximation for simplicity, and the conversion from distance to angle above refers to the edge of the volume, at $z=10$. Because the whole of the Outer Rim simulation is necessary to produce a single mock RST HLS, we will analyse the four separate cuboids drawn from Outer Rim separately, as quarter-mocks. The scatter between the results from the different quarter-mocks will give us a crude estimate of cosmic variance.

In order to include redshift evolution in the mock surveys, we choose a cuboid with fixed position in different snapshots, and calculate the comoving radial distance of each redshift. We then extract the corresponding slice in the cuboid from different snapshots according to the relevant comoving radial distance. Stacking all the slices in order of redshift yields the mock survey. The redshift evolution is discrete rather than a true lightcone - we use redshifts at z=7.4, 7.7, 7.9, 8.4, 8.6, 8.8, 8.0, 9.2, 9.4, 9.6, 9.8 and 10.0. We include the effect of redshift distortions in the mock surveys using the peculiar velocities of the simulated galaxies.

As the Outer Rim simulation covers a $3000\times 3000\times 3000\ h^{-3}$Mpc³ cube, we are able to create four mock cuboids from the Outer Rim data, using different parts of the volume. Each of these covers $\sim 1/4$ of the RST HLS area, so that there are therefore four approximately independent quarter-mock surveys to use for baryon acoustic peak fitting.

As stated above, we assign BlueTides galaxy luminosities to Outer Rim halos according to their mass. We make several sets of mock surveys corresponding to different galaxy detection thresholds (different line sensitivities):

(a) L7 quarter-mocks. The 2015 WFIRST technical report (Casertano et al., 2015) estimates that with current design parameters, the RST HLS grism survey will have a $7\sigma$ H$\alpha$ line detection sensitivity of $6-12\times 10^{-17}\ {\rm erg/s/cm}^{2}$ for sources with effective radius $<0.3$ arcsecs. This is for H$\alpha$ spectral line ($\lambda_{\rm rest}=6563\ \mbox{\AA}$) in the redshift range $z=1.05$ to $z=1.88$. In Spergel et al. (2015), a value of $1\times 10^{-17}\ {\rm erg/s/cm}^{2}$ at $7\sigma$ is quoted, which is in the same range. In our analysis, we use this value, but assume a 5 $\sigma$ detection limit, which results in a lower limit to the detectable line flux of $7\times 10^{-17}\ {\rm erg/s/cm}^{2}$. As the RST survey strategies are as yet unfinalized, we assume that there could be some flexibility, both in the selection of the fiducial sensitivity, and also possibilities of deeper integrations over part of the survey area.

We assume that this sensitivity can be applied to the Lyman-$\alpha$ line over the same observed range of wavelengths. The line sensitivity for resolved sources is approximately a factor of two worse, but we expect the galaxies to be unresolved. For example, Feng et al. (2015) find half light radii of $0.5\sim 1$ kpc for the most massive galaxies in BlueTides at redshift $z=8$, which corresponds to $0.1-0.2$ arcsec. For our first mocks, those closest to the RST design sensitivity, the mimimum Lyman-$\alpha$ line flux is $7\times 10^{-17}\ {\rm erg/s/cm}^{2}$. We name these the L7 quarter-mocks in what follows.
(b) L3.5 quarter-mocks. Unfortunately, even with the above parameter choices, which are in themselves relatively optimistic regarding the capabilities of RST, we shall see below that the possibility of detecting BAO is marginal. We therefore make another set of quarter-mocks with a 50% less restrictive line detection sensitivity limit, $3.5\times 10^{-17}\ {\rm erg/s/cm}^{2}$. This sensitivity limit could in principle be reached by increasing the integration time and lowering the sky area covered by the HLS. We name these the L3.5 quarter-mocks.
(c) L1 quarter-mocks. The final set of mocks we use for analysis have a line detection sensitivity of $1\times 10^{-17}\ {\rm erg/s/cm}^{2}$ 7 times better than the planned RST design. We create these (L1) quarter-mocks in order to explore BAO measurement when there is extremely good sampling.

In Fig. 2 we show galaxies with SFR more than $44.20\ {\rm yr}^{-1}M_{\sun}$ account for 15.9% of the total galaxy counts. We assume that if a galaxy mass falls between $2\times 10^{11}\ h^{-1}M_{\sun}$ and $5\times 10^{11}\ h^{-1}M_{\sun}$, the probability of it being able to be detected is 15.9%. We calculate this probability for different mass ranges by using data from BlueTides, and match galaxies from simulated cuboid to decide which galaxies will be detected in the mock survey. Fig. 6 shows the difference after we apply this selection. Only 2% galaxies is kept in the right figure compared with left one. And, Fig. 7 shows the relationship between probability and mass. The maximum mass in redshift 9 and 10 is much smaller than other redshift, but considering the maximum mass of the Outer Rim is $1.30\times 10^{12}\ h^{-1}M_{\sun}$ in redshift 9.2 and $6.7\times 10^{11}\ h^{-1}M_{\sun}$ in redshift 10. Fig. 7 is enough to cover the mass range in the Outer Rim. And Fig. 8 shows the distribution of galaxies along the z-axis under different line sensitivity. There are several sudden changes in the distribution, the reason is that we use snapshots from discrete redshift to construct mock data so there is discontinuity between two redshift slices.

We use the Landy Szalay estimator (Landy & Szalay, 1993) to compute anisotropic two point correlation functions $\hat{\xi}(s)$ of galaxies in our datasets. We prefer the Landy Szalay estimator, since its performance has been proved to be better than other comparable two point correlation functions (Davis & Peebles, 1983; Hamilton, 1993) at large scales (Pons-Bordería et al., 1999; Kerscher et al., 2000). Equation 2 gives the formula for the calculation of the Landy Szalay estimator.

In equation 2, $\mu=cos\ \theta$ (where $\theta$ is the angle between the line of sight and the pair separation vector for any given pair), $DD(s,\mu)$ depicts the pair count of galaxies with pair separation distance $s$ and orientation $\mu$. The symbol $DR(s,\mu)$ represents the cross-pair counts between galaxies and randoms which have pair separation $s$ and angle $\theta$, and $RR(s,\mu)$ is the number of pairs for a random distribution. We have computed the correlation functions as a function of $s$ and $\mu$ as an additional check on our calculations, to see that the clustering is behaving as expected. We do not present or use the two dimensional results further because the use of the full two-dimensional two point correlation functions in analyses is computationally expensive. One can use isotropised correlation functions ($\tilde{\xi}_{f}(s)$) to circumvent this. To obtain isotropised correlation functions, one uses specific kernels to marginalize out the $\mu$ coordinate in the two point correlation function. In Hamilton (1993) it was shown that one can use Legendre polynomials as a kernel ($f(\mu)=P_{\ell}(\mu)$ to find the isotropic correlation functions $\tilde{\xi}_{\ell}(s)$. Different values of $\ell$ will lead to isotropic correlation functions of different orders. In our research, we use isotropic correlation functions with Legendre polynomials and $\ell=0$ to obtain isotropized monopoles: ($\tilde{\xi}_{0}(s)$). For completeness, Equation 3 shows the process by which anisotropic correlationfunction ${\xi}_{\ell}(s)$ is obtained from the marginalization of anisotropic two correlation function ($\hat{\xi}(s,\mu)$).

In equation 3.3, we show how we can compute $\tilde{\xi}_{0}(s)$ from equation 3.

In the work presented in this paper, we have used 100 bins in $\mu$ in the computation of all two dimensional two-point correlation functions. Also, all the multipoles ($\xi_{0,2}(s)$) that we compute for BlueTides and Outer Rim galaxies have evenly spaced bins of width $3\ h^{-1}$Mpc in $s$. Apart from their use as an internal check here, we reserve the presentation of results for multipoles higher than the monopole for future work.

For fitting of the correlation functions, we use techniques which are outlined in Anderson et al. (2012); Xu et al. (2012, 2013); Vargas-Magaña et al. (2014); Ansarinejad & Shanks (2018); Sarpa et al. (2019). We briefly describe the methods that we have used in this section.

We use the following fiducial form to fit space correlation monopoles $\xi_{0}(s)$ over the range $33\ h^{-1}{\rm Mpc}<s<138\ h^{-1}{\rm Mpc}$ with bin sizes of $s=3\ h^{-1}$Mpc.

Here, the symbol ‘$s$’ represents the pair separation between galaxies in redshift space. The quantity $B$ is a multiplicative constant bias which adjusts the amplitude of the correlation functions ($\xi$), and lets the model account for any unknown large-scale bias. The term

is included so that the fitting model can marginalize over broad-band effects due to scale dependent bias, redshift space distortions and errors during our assumption of the fiducial cosmology. The parameters $a_{1},\ a_{2},\ a_{3}$ in equation 6 are nuisance parameters. The choice of number of parameters in $A(s)$ in our work is inspired by Anderson et al. (2012); Xu et al. (2012, 2013); Vargas-Magaña et al. (2014); Ansarinejad & Shanks (2018); Sarpa et al. (2019). The inclusion of $A(s)$ can help to mitigate the effects of assumption of a wrong cosmological model. At the same time, $A(s)$ can be also be thought of as an effectual description of mode coupling not affecting the BAO scale but biasing its measurement if not considered properly (Crocce & Scoccimarro, 2008).

The scale dilation parameter $\alpha$ measures the position of the BAO peak of the data relative to the model. One can also think of the parameter $\alpha$ as the extent to which the acoustic peak in the data is shifted relative to the model. This isotropic shift in the positions of the BAO peak in the data and the model occurs due to non-linear structure growth. A numerical definition of $\alpha$ can be obtained in equation 7.

In equation 7 $r_{d}$ represents the sound horizon at the drag epoch. The superscript ‘fid’ denotes the fiducial cosmology (given in Table 1.) A value of $\alpha<1$ would represent a shift towards larger scales, whereas a value of $\alpha>1$ denotes a shift towards smaller distance scales. The volume averaged distance $D_{\rm V}$ is related to the Hubble constant $H$, the angular diameter distance $D_{\rm A}$, the speed of light $c$ and the redshift $z$ by the following equation:

The template power spectrum which we use is defined in the below mentioned relation.

In equation 9, $P_{\rm lin}(k)$ represents the actual linear power spectrum at $z=0$. The term $P_{\rm smooth}(k)$ denotes the de-wiggled limit where the BAO feature is removed. More details of $P_{\rm smooth}(k)$ can be found in Eisenstein & Hu (1998). The term $\Sigma_{{\rm NL}}$ helps to damp the BAO features in the linear power spectrum $P_{\rm lin}(k)$. This parameter also helps to account for the broadening and shift of the BAO feature due to non linear structure evolution. We use $\Sigma_{{\rm NL}}=0\ h^{-1}$Mpc in our work, i.e. $P_{{\rm m}}(k)=P_{\rm lin}(k)$. Choice of this value of $\Sigma_{{\rm NL}}$ is justified given the high redshifts that we consider ($z_{\rm OuterRim}=9.4$) and the correspondingly low value of the growth factor, $g=0.123$. More details of this are given in Section 3.5.

We use CAMB (Lewis & Bridle, 2002) to compute the power spectrum $P_{{\rm m}}(k)$ using the same cosmological model that is adopted in the Outer Rim simulation.

The model correlation function $\xi_{\rm m}$ is the Fourier transform of the the template power spectrum given in equation 9.

Xu et al. (2013) showed that variations in the bias like term $B$ do not affect the position of the BAO peak in monopoles. At the same time, the parameter $\alpha$, which gives us a quantification of the acoustic scale, affects the position of the BAO peak. As such, $\alpha$ is the parameter that we are most interested in obtaining from our fits. The models that we use are non-linear in $\alpha$, so we can nest a linear least-squares fitter inside a non-linear fitting routine. The best fit value of $\alpha$ is obtained by minimizing the $\chi^{2}$ goodness-of-fit indicator (non-linear fitter). To obtain an optimum value of $\alpha$, we compute $\chi^{2}$ and shift models in the range $0.8<\alpha<1.2$ with intervals of $\alpha=0.01$. The choice of range within which values of $\alpha$ are examined, viz. $\left[0.8,1.2\right]$, is inspired by Vargas-Magaña et al. (2014) and Ansarinejad & Shanks (2018).

The best fit value of $\alpha$ is obtained by minimizing the $\chi^{2}$ goodness-of-fit indicator given in equation 10.

In equation 10, $\xi^{\rm fit}(\alpha)$ represents a model vector at a given value of $\alpha$, $\xi^{\rm mock}$ denotes the correlation function measured from mock cuboids corresponding to a given line detection sensitivity, and $\mathbf{\Sigma}$ represents a diagonal covariance matrix which is obtained from all mock cuboids in the considered line detection sensitivity. In an ideal case, one would have prefered to use a variance-covariance matrix from a large number (preferably $\sim$100 or higher) of mock catalogues. But, in our work, we are dealing with high redshifts ($z>7.6$) and we only have four mock cuboids for each of the three line detection sensitivities. We use these cuboids to construct separate covariance matrices for each line detection sensitivity. Hereafter, we shall use the symbol $\alpha^{\ast}$ to denote the best fit values of the scale dilation parameters $\alpha$. To obtain robust results for best fit values of scale dilation parameters and best fit model vectors ($\xi^{\rm fit}(\alpha^{\ast})$), we construct diagonal covariance matrices from mocks in two steps. In the first step, we compute a diagonal covariance matrix (say $\mathbf{\Sigma}_{0}$) from cuboids which correspond to the same line detection sensitivity. (Separate $\mathbf{\Sigma}_{0}$ are constructed for catalogues which correspond to without patchy reionization and with patchy reionization.) In the second step, we construct a new diagonal covariance matrix (say $\mathbf{\Sigma}$) by using a moving average (with a window size of three) over the diagonal elements of $\mathbf{\Sigma}_{0}$. That is to say,

To compute $\mathbf{\Sigma}^{\{0,0\}}$, we use equation 12

Similarly, to find the last diagonal element of $\mathbf{\Sigma}$, i.e. $\mathbf{\Sigma}^{\{N,N\}}$ we use equation 13. Here $N$ denotes the index of the largest distance scale ($138\ h^{-1}{\rm Mpc}$) in our fitting range.

The approach highlighted in equations 11 , 12 and 13 help in getting better agreement in positions of BAO peaks in mock and best fit model correlation functions.

Detailed description of the BAO fitting methodologies can also be found in Anderson et al. (2012); Xu et al. (2012, 2013); Vargas-Magaña et al. (2014).

We follow the semi-analytic model outlined in Battaglia et al. (2013) to construct reionization fields from our density fields. The density fluctuation field is supposed to be at the midpoint of reionization, which is $z=9.5$ for BlueTides. So, we use data from Outer Rim simulation which is closest to $z=9.5$ (i.e. $z_{\rm OuterRim}=9.4$). The fluctuation field for halos from the Outer Rim catalogue at $z_{\rm OuterRim}=9.4$ is defined as:

Here, $\bar{\delta}_{\rm halo}$ is the mean density of halos. The required linear bias, $b$ is obtained by using equation 15.

In equation 15, $g$ is the growth factor between $z=9.4$ and $z=0.0$, i.e. $g=0.123$. The symbol $\xi_{0}^{z=9.4\ ;\ {\rm OuterRim}}$ denotes the correlation function monopole computed from Outer Rim halos at $z=9.4$. The symbol $\xi_{0}^{z=0.0\ ;\ {\rm CDM}}$ represents the monopole from cold dark matter (CDM) at $z=0.0$. We use a value of $\sigma_{8}=0.8$. After the aforementioned steps, one can obtain the matter overdensity ($\delta_{\rm m}$) from equation 16

Before we outline the steps for the computation of the reionization redshift field, $z_{\rm RE}(\mathbf{x})$, it is important to define the fluctuation field $\delta_{z}(\mathbf{x})$.

In equation 17, $\bar{z}$ is the mean value of the $z_{\rm RE}$ field. At the same time, $\bar{z}$ sets the midpoint of reionization which is $z=9.4$ for Outer Rim data.

From here, we can get the reionization redshift field, $z_{\rm RE}(\mathbf{x})$, by following the below mentioned steps.

1. 1.

   We Fourier transform the density fluctuation field $\delta_{\rm m}$ from Outer Rim simulation at the midpoint of reionization (i.e. $z=\bar{z}=9.4$) and obtain $\tilde{\delta}_{\rm m}(k)$.

2. 2.

   We multiply the field $\tilde{\delta}_{\rm m}(k)$ with the bias factor given in equation 18 to obtain $\tilde{\delta}_{\rm z}(k)$.

   $\displaystyle b_{\rm mz}(k)=\frac{b_{0}}{\left(1+k/k_{0}\right)^{\alpha}}$

   (18)

   The above mentioned bias factor uses three parameters, viz. $b_{0}$, which is the bias amplitude on the largest scales, $k_{0}$, which is the scale threshold and $\alpha$, which is the asymptotic exponent. For the three parameters $b_{0}$, $k_{0}$, $\alpha$, we use the best-fit values which were used in Battaglia et al. (2013). More specifically, we use values of $b_{0}=0.593$ , $k_{0}=0.185\ {\rm Mpc}^{-1}h$ and $\alpha=0.564$ in our work. Since we don’t focus much on small distance scales in our work, we ignore the top hat filters that are used in Battaglia et al. (2013).

3. 3.

   We inverse Fourier transform $\tilde{\delta}_{\rm z}(k)$ to real space to get the $\delta_{\rm z}$ field.

4. 4.

   Now, we can convert the $\delta_{\rm z}$ field to $z_{\rm RE}$ using equation 17.

Once we have the $z_{\rm RE}(\mathbf{x})$ field, we can use that to modulate the galaxy distribution. Basically, we can use the same grid for all snapshots, and if a galaxy is in a cell that has reionized by that snapshot redshift, we treat it as being visible. If a galaxy is in a cell which hasn’t reionized (these will be the lower density cells, based on how the patchy reionization model works), then we say that it won’t be visible, and we remove it from the sample. As with the previous mock catalogues, we include the effect of redshift-space distortions in the galaxy positions: we use real space position to locate the galaxy in the reionization field, and then position the visible galaxies in redshift space.

We follow the following prescription to distinguish between lower and higher density cells. We use a Gaussian filter with $\sigma=10\ h^{-1}$Mpc to smooth the $z_{\rm RE}(\mathbf{x})$ field. In the smoothed $z_{\rm RE}(\mathbf{x})$ field, any cell whose value is less than the mean value of the smoothed $z_{\rm RE}(\mathbf{x})$ field is considered as a low density cell. Similarly, any cell in the the smoothed $z_{\rm RE}(\mathbf{x})$ field whose value is bigger than the mean is considered as a higher density cell. Fig. 9 shows a slice (${\rm r}_{\rm z}=1500h^{-1}$ Mpc) of the three dimensional grid of cells with ionized and non-ionized regions which is obtained from the aforesaid method. We compare all snapshots (redshifts) with this grid of high and low density cells and construct a new catalogue of galaxies which correspond to scenarios after patchy reionization. In Fig. 10, we see slices (at ${\rm r}_{\rm z}=1500h^{-1}$ Mpc) from Outer Rim data at $z=7.4$ with and without patchy reionization.

We measure correlation function monopoles in these galaxy catalogues which correspond the cases with and without patchy reionization. In the non-patchy case, we assume that the intervening intergalactic medium does not affect the visibility of the Lyman-$\alpha$ line from galaxies (i.e. the universe is optically thin at the time of Lyman-$\alpha$ emission).Results of monopoles and the ensuing BAO fitting are presented in Section 4.

The top four panels in Fig. 11 show plots of variation of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L7 quarter-mock catalogue which correspond to distribution of halos without patchy reionization. The plots in the bottom four panels of Fig 11 show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for the aforementioned mock cuboids. For the cuboids which correspond to distribution of halos without patchy reionization, the minimum values of $\chi^{2}$ are 27.41, 31.32, 28.17, 45.86. For these cuboids, the best fit vaues of $\alpha$ are 1.08, 1.16, 0.86, 1.08 respectively. Consequently, for the cuboids from L7 quarter-mocks the mean value of minimum $\chi^{2}$ is 33.19, and the mean value of $\alpha^{\ast}$ is 1.04. We estimate the fractional error on the distance to $z=7.8$ (the mean galaxy redshift) by computing the sample standard deviation of the $\alpha^{\ast}$ results for the four quarter-mocks, and then dividing by $\sqrt{4}$ to give the error on the mean. The fractional error computed this way is 0.064, i.e., a distance error accurate to $6.4\%$. We can see from the bottom panel of Fig. 11 however that one of the BAO peaks is not visibly detected.

One thing that is worth noting as this point is the fact that for a given mock catalogue (L7 or L3.5 or L1), plots at corresponding positions in the figures which illustrate results with and without patchy reionization are results gotten from the same cuboid. That is, if the leftmost and topmost subplot of Fig. 11 corresponds to Cuboid 1, then leftmost and topmost subplot of Fig. 12 also corresponds to Cuboid 1. The same holds of all mock catalogues (L7 or L3.5 or L1).

The left four panels in Fig. 12 illustrate the variations of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L7 quarter-mock catalogue which correspond to distribution of halos after applying patchy reionization. In the right four panels in Fig. 12, we show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for these four mock cuboids. The cuboids from L7 quarter-mock catalogue which correspond to distribution of halos with patchy reionization have the following minimum $\chi^{2}$ values: 24.41, 26.68, 34.96, 26.59. These cuboids have best fit vaues of $\alpha$ as 1.04, 0.81, 1.05, 1.19 respectively. As such, the mean value of minimum $\chi^{2}$ for these cuboids is 28.16, and the mean value of $\alpha^{\ast}$ is 1.02. Again, as for the previous figure, one of the BAO peaks (in the top right panel) is not detected (although there is a chi-squared minimum and an estimate of $\alpha^{\ast}$ ). From the scatter of $\alpha^{\ast}$ values between quarter-mocks, we find the fractional error on the mean distance to be 0.079.

The left four panels of Fig. 13 show plots of variation of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L3.5 quarter-mock catalogue which correspond to distribution of halos without patchy reionization. In the right four panels of Fig 13, we show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for the aforementioned mock cuboids. For the cuboids from L3.5 quarter-mock which correspond to the distribution of halos without patchy reionization, the values of minimum $\chi^{2}$ are 21.59, 20.03, 16.20, 26.42, and the values of $\alpha^{\ast}$ are 0.99, 1.08, 1.00, 1.12 respectively. Hence, the mean values of minimum $\chi^{2}$ and $\alpha^{\ast}$ are 21.06 and 1.05 respectively. The error on the BAO distance measurement is 0.032.

The plots in the left four panels of Fig. 14 illustrate the variations of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L3.5 quarter-mock catalogue which correspond to distribution of halos after applying patchy reionization. In the bottom two rows of Fig. 14, we show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for these four mock cuboids. These cuboids have minimum values of $\chi^{2}$ of 38.80, 22.52, 43.14, 32.26 and best fit values of $\alpha$ as 0.95, 0.99, 1.02, 0.98. As such, the mean values of minimum $\chi^{2}$ and $\alpha^{\ast}$ are 34.18 and 0.98 respectively. The error on the BAO distance measurement is 0.014.

The left four panels of Fig. 15 show plots of variations of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L1 quarter-mock catalogue which correspond to distribution of galaxies visible without patchy reionization. In the right four panels of Fig. 15, we show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for the aforementioned mock cuboids. These cuboids have the following minimum values of $\chi^{2}$: 19.49, 18.58, 11.79, 13.34. They have the following values of $\alpha^{\ast}$: 1.01, 1.03, 0.99, 0.95. Consequently the mean values of minimum $\chi^{2}$ and $\alpha^{\ast}$ for these mock cuboids are 15.80 and 1.00 respectively. The fractional distance error is 0.015.

The left four panels of Fig. 16 illustrate the variations of $\chi^{2}$ with $\alpha$ for the four mock cuboids in L1 quarter-mock catalogue which correspond to distribution of galaxies visible with patchy reionization. In the right four panels of Fig 16, we show monopoles computed from Landy Szalay estimators and best fit monopoles ($\xi(\alpha^{\ast})$) for these four mock cuboids. The cuboids L1 quarter-mock catalogue which correspond to distribution of galaxies after patchy reionization have the following values of minimum $\chi^{2}$: 12.33, 14.26, 15.67, 14.26. The values of $\alpha^{\ast}$ for these cuboids are as follows: 0.96, 1.06, 1.01, 0.96. Hence, the mean values of minimum $\chi^{2}$ and $\alpha^{\ast}$ for these mock cuboids are 14.13 and 1.00. The fractional distance error is 0.024.

Table 3, gives summaries of our best fit results for all mocks in all three quarter-mock catalogues.

We now carry out an analysis of the BAO distance errors and their scatter, in order to check that they are statistically consistent. For the L7 quarter mock catalogues, we found the following values of best fit $\alpha$ ($\alpha^{\ast}$) from the four cuboids that represent the distribution of halos without patchy reionization: 1.08, 1.16, 0.86, 1.08. These values of $\alpha^{\ast}$ lead to an overall fractional error of 0.064. At the same time, we observed the following values of $\alpha^{\ast}$ in the four cuboids in the L7 quarter mock catalogues with patchy reionization : 1.04, 0.81, 1.05, 1.19. These values of $\alpha^{\ast}$ yield an error of 0.079 on BAO distance measurement. To reconcile the scatter in the values of $\alpha^{\ast}$ that we see for cuboids with and without patchy reionization, we carry out Fisher’s exact test (Fisher, 1922; Agresti, 1992). For each category of cuboids (with reionization or without), we are working with only four cuboids. This is a classic case of data with small sample sizes. Fisher’s exact test enables us to conduct accurate statistical significance tests on small samples. In this test, our null hypothesis is that there is no significant deviation in the distributions of the values of $\alpha^{\ast}$ for cases with and without patchy reionization. A two-tailed Fisher exact hypothesis test for comparison of the distributions of the aforesaid values, gives us a p-value which is approximately 0.99. This p value is very high and it strongly suggests that we cannot ignore the null hypothesis.

Similarly, for L3.5 mocks, we have the following values of $\alpha^{\ast}$ for halos without patchy reionization: 0.99, 1.08, 1.00, 1.12. For the four cuboids with patchy reionization, we have: 0.95, 0.99, 1.02, 0.98. We again conduct Fisher’s exact test. Our null hypothesis in this test is that there is no significant deviation in distributions of values of $\alpha^{\ast}$. From a two-tailed Fisher’s exact test, this time we get a p-value of approximately 1. This strongly suggests that there is significant evidence in support of our null hypothesis.

Finally, for L1 mocks, we get the following values of $\alpha^{\ast}$ for halos without patchy reionization: 1.01, 1.03, 0.99, 0.95. When cuboids have patchy reionization, we get the following values of $\alpha^{\ast}$: 0.96, 1.06, 1.01, 0.96. A two-tailed Fisher’s exact test gives us a p-value of approximately 1, again strong evidence of the fact that there is no significant deviation in distributions of values of $\alpha^{\ast}$ for cases with and without patchy reionization.

The results that we get from Fisher’s exact tests carried out on the mocks also suggest that we cannot make significant inferences from comparisons of scatter in values of $\alpha^{\ast}$ in cuboids with patchy reionization and without patchy reionization.

The results that we see in Figures 12, 11, 14, 13, 16 and 15 suggest that the clustering amplitude changes somewhat significantly with and without reionization and also depends on the line strength. Equation 19 below defines the linear bias due to reionization compared to no reionization for a given cuboid. It is defined using a ratio of two galaxy correlation functions and is of course different from the usual linear bias of galaxies with respect to mass.

We tabulate the values of $b_{\rm reion}$ for cuboids in the L7, L3.5 and L1 mocks in Table 4. For the L7 mock catalogue, we obtain a mean $b_{\rm reion}$ of 1.32. We estimate the fractional errors on $b_{\rm reion}$ by computing the sample standard deviation of the $b_{\rm reion}$ results for the four quarter-mocks, and then dividing by $\sqrt{4}$ to give the error on the mean. The fractional error computed this way is 0.136. In the case of the L3.5 mock catalogue, we get a mean $b_{\rm reion}$ of 1.24 and a fractional error of 0.089. Finally, for the four quarter-mocks in L1 mock catalogue, we obtain a mean $b_{\rm reion}$ of 2.41 and fractional error of 0.514. It is as we might expect: in all cases the effects of patchy reionization lead to higher amplitude clustering for the observable galaxies. The fainter galaxies (the L1 mocks) are also relatively much more strongly affected than the brighter galaxies.

Using mock catalogues, we have investigated the prospect of utilizing the BAO feature at redshifts $z>7.6$ as a standard ruler. We have utilized four mock catalogues, each representing approximately one quarter of the volume that will be surveyed by RST with grism spectroscopy.

We see evidence for BAO peaks in three of the four mocks when we use the planned RST HLS Lyman-$\alpha$ line detection sensitivity of $L=7\times 10^{-17}\ {\rm erg/s/cm}^{2}$. Hence, the BAO feature is likely detectable for this line detection sensitivity (measured in three out of four survey quadrants analysed independently). When we make mocks with better line detection sensitivities, we see evidence of BAO features in all mocks. For the standard line detection sensitivity, $L=7\times 10^{-17}\ {\rm erg/s/cm}^{2}$, we find a root mean square fractional distance error of 0.068 after comparison between mocks and fitted correlation function monopoles. For a factor of two fainter line detection sensitivity, $L=3.5\times 10^{-17}\ {\rm erg/s/cm}^{2}$, we find a distance error of 0.013. We find a distance error of 0.021 for mocks with line detection sensitivity $L=1\times 10^{-17}\ {\rm erg/s/cm}^{2}$.

Our analysis of the use of BAO at high redshifts has included some of the relevant features and difficulties, such as non-linear clustering, redshift distortions, intrinsic galaxy line luminosities, and the effect of reionization bubbles. The results are promising, as we have summarised in Section 5.1 above. We have however left much for future work, including the visibility of the Lyman-$\alpha$ line at these redshifts (including the effects of dust and the Lyman-$\alpha$ damping wing), and the effects of interlopers and low redshift foregrounds. These factors could easily derail the possibility of measuring large-scale galaxy clustering at redshift $z>7.6$, and should be investigated in detail. On the other hand, there are examples of galaxies being detected in Lyman-$\alpha$ emission at these redshifts, for example the $z=7.51$ object by Tilvi et al. (2016), and $z=7.452$ from Larson et al. (2018). Even if the Lyman-$\alpha$ emission line itself cannot be seen, the Lyman-$\alpha$ break can also be used to spectroscopically identify galaxies, as was done using the HST grism at $z=11.1$ by Oesch et al. (2016). This could also be investigated to see whether it could enable clustering measurements at the BAO scale.

We have seen that the detectability of the BAO peak is just at the level of being possible with the fiducial flux limit of the RST HLS. To do this, we have already reduced the nominal 7$\sigma$ detection limit to 5$\sigma$ ($7\times 10^{-17}\ {\rm erg/s/cm}^{2}$). The exact survey design choices made when the RST mission is finalized will therefore be extremely relevant to whether the BAO measurement will be feasible. We have found that with a sensitivity a factor of two better the BAO can be seen clearly in four survey quadrants. This means that it could be advantageous to make use of smaller and deeper areas.

The analysis presented in this work was based on only 4 quarter-mocks for each of the tested line detection sensitivities. In the future, large number of mock catalogues should be generated to allow for the build up of a full covariance matrix for the correlation function. Because of the small number of mocks here, we were only able to use the diagonal elements of the covariance matrix in model fitting. As a result, the distance error estimates, which were derived from the scatter between mocks, were necessarily conservative. In a future work with large numbers of mocks, the extra information available from the off diagonal covariance matrix elements will mean that the fits to the BAO peaks will be improved and the distance errors will be somewhat reduced. The uncertainty in the distance error will also be better estimated. In the present case, the error on the error on the mean from only four quarter mocks will be large, and this is why the distance scale error on the L1 mocks increases to 2.4% from 1.4% for the L3.5 mocks, in spite of a factor of 3 increase in sensitivity.

In our treatment of reionization, we have assumed that galaxies which lie within ionized bubbles will be visible in Lyman-$\alpha$ (consistent with the interpretation of Bagley et al. (2017)). We have compared the effects of including this patchy reionization, or not including it (all galaxies visible), and have found that the effect on BAO peak detection is quite small. For example, for the fiducial L7 sensitivity case, the effect of patchy reionization only increases the distance error from 6.4% to 7.9 %. This is likely to be because in the patchy model for reionization employed, overdense regions reionize first, and this is where most of the galaxy clustering signal lies. The reionization bubbles are also on smaller scales, than the BAO signal, and have a range of sizes, and therefore there is no masking of the BAO peak by bubbles. We note however that we have carried out no Lyman-$\alpha$ radiative transfer on the scales of individual galaxies, and that the visibility or not of the line in different environments is likely to be much more complex than our simple bubble model.

One significant difference between the BAO measurements at redshifts $z>7$ and those for $z<1$ is the lack of non-linearities in gravitational clustering. For example we can see that in the right panels of Fig. 15 there is no obvious smoothing of the BAO peak compared to the linear CDM fit. This is an advantage of working at higher redshifts, as no reconstruction techniques (Sherwin & White, 2019) are needed to recover the maximum benefit. The high redshift galaxies do have a high bias (Bhowmick et al., 2019) and so their clustering is detectable where clustering of the density field would not be.

Measuring BAO in a new regime would at the very least be useful as a consistency check on measurements that have been made at low redshifts and in the CMB. It is important to fill in the large gap between the Lyman-$\alpha$ forest measurements at $z\sim 2.2$ (Aubourg et al., 2015) and the features seen at $z\sim 1100$ by Planck and WMAP (Planck Collaboration et al., 2018b; Bennett et al., 2013). Such measurements might shed some light on the tension between BAO and other Hubble parameter measurements (Cuceu et al., 2019), and whether early dark energy may be acting before the recombination epoch (Verde et al., 2019; Poulin et al., 2019).

One of the targets of cosmological analysis with the RST satellite is large scale structure measurements from galaxies at redshifts $z\sim 1-3$. This will extend the successful redshift space distortion, BAO, and other techniques to significantly earlier epochs than current large scale structure measurements from ground based surveys. Although much attention has rightly been focussed on this redshift range, we have seen in this paper that it may be possible to make an even greater leap, and use the large scale structure traced out by the very first massive galaxies as a cosmological probe. Although many potential pitfalls remain to be investigated, it is worth bearing this in mind when considering the RST survey strategies and parameters.