Generating extremely large-volume reionisation simulations

The Trott et al. (2020) MWA upper-limits were determined from observations spanning 137–197 MHz in two slightly overlapping 30.72 MHz bands. Over these bands, the MWA frequency resolution is $\sim 80$ kHz, which sets the minimum spatial resolution for our large-volume simulations and the maximal line-of-sight length required (384 channels). At 150 MHz ($z\sim 8.5$), the MWA field of view is $\sim 25^{2}$ ${\rm deg.}^{2}$, corresponding to a comoving side-length of $\sim 4$ Gpc. To exceed these requirements, we generate a single simulation volume with a comoving transverse side-length of 7.5 Gpc over 6400 voxels and 600 voxels ($703.1$ Mpc) along the line-of-sight direction³³3Such a simulation required $\sim 6$ days, using 32 CPUs and $\sim 10$ GB of memory. The majority of this time is spent reading/writing to file, which could be considerably improved using solid state hard-drives.. This simulation corresponds to a cell resolution of $\sim 1.17$ cMpc. To our knowledge, this is the largest volume 21-cm simulation in the literature⁴⁴4As a verification of our method we also performed a simulation with a transverse side-length of 9.6 Gpc and 8092 voxels.. Using this simulation, we construct two 21-cm light-cones spanning the observing bands ($z=6.2-7.5$ and $z=7.5-9.4$) by stitching together 2D slices of the 21-cm signal from our simulation volume by linearly interpolating in cosmic time (e,g, Datta et al. 2012, 2014; La Plante et al. 2014; Ghara et al. 2015; Mondal et al. 2018).

For this simulation, we adopt $f_{\ast,10}=0.05$, $\alpha_{\ast}=0.5$, $f_{\rm esc,10}=0.08$, $\alpha_{\rm esc}=-0.5$, $t_{\ast}=0.5$ and $M_{\rm turn}=10^{8.7}$, consistent with the fiducial model adopted in Park et al. (2019), which was shown to match all recent observational constraints on the reionisation epoch. Additionally, in neglecting inhomogeneous recombinations we are required to set an effective mean free path for the ionising photons inside the ionised regions, $R_{\rm mfp}$. We adopt $R_{\rm mfp}=15$ comoving Mpc motivated by the fiducial model in Greig & Mesinger (2015).

In Figure 1 we present a 2D slice of the IGM ionisation fraction, $x_{\mathrm{H\,{\scriptscriptstyle I}}{}}$, at $z=8$ for which we obtain an IGM neutral fraction of $\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}\sim 0.55$. As a visual demonstration of the physical extent of these simulations, we overlay the respective fields of view for each of the MWA (red), HERA (blue), LOFAR (purple) and SKA (teal). Note that unlike the MWA, LOFAR and SKA which use phase tracking to follow a single patch of the sky, HERA is a drift-scan experiment (sky rotates over the fixed zenith pointing dishes), thus, here we show its instantaneous field of view. In Table 1 we provide the collecting area of a single station/element and the corresponding field of view for each 21-cm interferometer experiment.

Typically, simulations tailored for a 21-cm interferometer experiment will be designed to match the physical extent of the instruments primary beam. Although problematic for full tests of end-to-end pipelines, for most usage cases this assumption should be sufficient. However, by only simulating a volume equivalent to the primary beam, these simulations forgo modes on much longer scales which would otherwise be present in the observed 21-cm signal. In principle, the absence of these large-scale modes could result in an under-estimate of the true large-scale power.

We can test the impact of neglecting modes beyond the physical extent of the primary beam by using our 7.5 Gpc simulation as a reference simulation to represent the true Universe. We then compare the 3D spherically averaged 21-cm power spectrum (PS) obtained from our 7.5 Gpc simulation to the PS obtained from the smaller, independent realisations of the 21-cm signal tailored to match each interferometer experiment. Throughout, we define the 21-cm PS as $\Delta^{2}_{21}(k,z)\equiv k^{3}/(2\pi^{2}V)\,\delta\bar{T}^{2}_{\rm b}(z)\,\langle|\delta_{21}(\boldsymbol{k},z)|^{2}\rangle_{k}$ where $\delta_{21}(\boldsymbol{x},z)\equiv\delta T_{\rm b}(\boldsymbol{x},z)/\bar{\delta T_{\rm b}}(z)-1$. Further, when calculating the 21-cm PS we neglect all $k_{\perp}=0$ modes. That is, we remove all pure line-of-sight modes as these are not visible to an interferometer⁵⁵5Note, this is somewhat of a crude approximation as interferometers do in fact have a response to these $k_{\perp}=0$ modes, although it is suppressed due to the tapered response. While this does not impact the results of this work, realistically one should include the true instrument response rather than taking this approximation when comparing simulations to observations. (e.g. Datta et al., 2012).

In Figure 2 we compare the 21-cm PS obtained from our 7.5 Gpc simulation (black curve) to the mean and 68th percentile scatter obtained from 10 independent realisations of the 21-cm signal. In each panel we compare the 21-cm PS obtained from the primary field of view of the MWA (red), HERA (blue), LOFAR (purple) and SKA (teal). In all cases, we compute the 21-cm PS over a simulation volume with the same observing bandwidth of $\sim 30$MHz (transverse distance is equivalent to the field of view).

Immediately from Figure 2 we see no obvious deviations at large scales. The mean 21-cm PS determined over the 10 independent realisations for each survey footprint overlaps with the 21-cm PS from the 7.5 Gpc simulation. On the largest scales for each individual instrument, we observe some scatter in the 21-cm PS owing to cosmic variance, where the binned 21-cm power has few modes (i.e. Poisson sampling error). Thus, provided the simulation volume exceeds the field of view, this should not be a concern. We can then conclude that omitting large-scales modes beyond the primary field of view does not alter the recovered statistics. This is consistent with the hybrid, large-volume ($>1$ Gpc) reionisation simulations generated by Kim et al. (2016), whereby the amplitude of the large-scale 21-cm PS converges for increasing simulation volumes beyond 500 Mpc.

Note, this investigation slightly differs from previous works by Iliev et al. (2014), La Plante et al. (2014) and Kaur et al. (2020). Here, in each case the authors explore sub-volumes from their largest reionisation simulation, which in the case of Iliev et al. (2014) and Kaur et al. (2020) is used to determine a minimum volume for reionisation simulations to achieve convergence of the large-scale radiative effects. Since this is a sub-division of a larger volume, these smaller simulations contain modes beyond the scale of the sub-volume. In our work, we generate new independent realisations of our smaller volume simulations that do not include modes beyond the scale of the simulation.

The 3D spherically averaged 21-cm PS is calculated by averaging over all Fourier modes that fall within spherical shells. Naturally, the variance of this statistic is then given by the Poisson sampling error based on the total number of Fourier modes in each spherical shell,

Here, $\Delta^{2}_{21}(k,z)$ is the dimensional 21-cm PS measured within a volume, $V$, with spherical shells separated by $\Delta k$.

However, this assumes that the 21-cm PS fully encodes all the information about the 21-cm signal (i.e. that the 21-cm signal is Gaussian). Instead it is well known that the 21-cm signal is highly non-Gaussian (e.g. Furlanetto et al., 2004a, b; Morales & Hewitt, 2004; Cooray, 2005), in which case the higher-order moments of the 21-cm signal will be non-zero. As a result of this Mondal et al. (2016) have shown that the cosmic variance uncertainty on the 21-cm PS is instead,

where $T(k,k)$ is the Trispectrum component which arises from the four-point correlation function of the brightness temperature fluctuations.

To explore this non-Gaussianity in the 21-cm PS cosmic variance, we extract a large number of sub-volumes from our large 7.5 Gpc simulation. In particular, we extract sub-volumes equivalent in spatial extent to the field-of-view for each 21-cm interferometer experiment (as listed in Table 1), each with the same 30 MHz bandwidth along the line-of-sight. To create these sub-volumes, we first split the full $7.5\times 7.5$ Gpc² sky-area into our requisite smaller volumes without overlapping regions. Then, to increase our number statistics, we generate an inset sky-area offset from the simulation edge by half the size of the respective instrument field-of-view. We then additionally split this inset area into non-overlapping sub-volumes. In doing so, we obtain a total of 5, 41, 221 and 365 different sub-volumes from our single, large volume simulation for the MWA, HERA, LOFAR and the SKA, respectively. By considering sub-volumes equivalent to each instruments field-of-view, we consider volumes $\sim 2-10\times$ larger than those considered by Mondal et al. (2016) and Mondal et al. (2017). Further, by splitting up our large volume simulation we can obtain up to $\sim 7$ times as many realisations, increasing our statistical sampling of the numerical cosmic variance uncertainty, albeit at the expense of notably less physics rich simulations.

In Figure 3 we show the resultant $1\sigma$ uncertainty on the 21-cm PS obtained from our distribution of sub-volumes (solid curves) for each 21-cm interferometer experiment. In contrast, we also show the expected theoretical estimate of the uncertainty from Equation 10 if we assumed Gaussian statistics (dashed curves). The dotted curve corresponds to the non-Gaussian (Trispectrum) component of the cosmic-variance as highlighted in Equation 11.

On relatively large scales ($k<0.5$ Mpc^-1), we find that the true cosmic variance uncertainty is well approximated by a Gaussian component for the theoretically expected error. However, progressing to intermediate and smaller scales, the amplitude of the non-Gaussianity increases significantly, with the Trispectrum component dominating by $\geq 10\times$ that of the Gaussian term. This transition scale of $k>0.5$ Mpc^-1 for a non-Gaussian dominated cosmic variance uncertainty on the 21-cm PS is consistent with that of Mondal et al. (2016) and Mondal et al. (2017) for a similar stage of reionisation ($\bar{x}_{\mathrm{H\,{\scriptscriptstyle I}}{}}\sim 0.55$). On these scales, as the EoR progresses, the growth and percolation of the ionised regions drives the non-Gaussianity in the cosmic 21-cm signal. Importantly, we find the relative amplitude (the amplitude simply scales as $\propto V^{-1}$) and shape of the full cosmic variance uncertainty to be consistent across the differently sized sub-volumes. This indicates no additional sources of cosmic variance uncertainty due to the finite-size of the simulation volume considered.

Note however, that while the non-Gaussian term for the 21-cm PS cosmic variance begins to dominate on scales larger than $k>0.5$ Mpc^-1, this need not be too concerning for 21-cm experiments. Typically, on these scales, the 21-cm PS will be dominated by the instrument thermal noise (see figure 2 of Greig et al. 2020) making these two terms comparable. Nevertheless, care must be taken when estimating the 21-cm PS cosmic variance error.

The cosmic 21-cm signal observed by radio interferometer experiments is measured in visibility space ($uv$ coverage), which is the 2D Fourier transform of the brightness signal in the image (sky) plane. Measuring the power spectrum then requires gridding these visibilities. However, these visibilities are frequency dependent. Line-of-sight power in Fourier space (frequency dependent) can then leak into the transverse (frequency independent) Fourier modes. This results in the well known ‘wedge’ feature in cylindrical 2D Fourier space, where measured power within this regime is contaminated (Datta et al., 2010; Vedantham et al., 2012; Morales et al., 2012; Parsons et al., 2012; Trott et al., 2012; Thyagarajan et al., 2013; Liu et al., 2014a, b; Thyagarajan et al., 2015a, b; Pober et al., 2016; Murray & Trott, 2018). While it is plausible to mitigate or ‘clean’ these contaminated modes (see e.g. Chapman & Jelić 2019 for a comprehensive review), it is typically easier to simply avoid this contaminated region entirely, measuring the power spectrum in the relatively pristine region above the ‘wedge’.

The spatial location of this ‘wedge’ in 2D cylindrical space is determined by,

where $k_{\parallel}$ and $k_{\perp}$ are the line-of-sight and transverse Fourier modes, respectively. The constant $b$ corresponds to an additive buffer region above the ‘wedge’ of $\Delta k_{\parallel}=0.1\,h$ Mpc^-1 extending beyond the horizon limit while,

Here, $D_{\rm C}$ is the comoving distance, $H_{0}$ is the Hubble constant, $E(z)=\sqrt{\Omega_{\rm m}(1+z)^{3}+\Omega_{\Lambda}}$ and $\theta$ is the angular radius of the field of view which we conservatively take as $\theta=\pi/2$ (i.e. observing down to the horizon).

Previously, we numerically estimated the cosmic variance in the 21-cm PS measured over the full 30 MHz bandwidth simulation volume for each interferometer experiment. To highlight the impact of ‘wedge’ avoidance on estimating the 21-cm PS, in particular on the amplitude of the statistical uncertainty, we perform the same analysis as in Section 3.3 instead calculating the 21-cm PS using only Fourier modes above the ‘wedge’.

In Figure 4 we compare the 1$\sigma$ cosmic variance uncertainty after ‘wedge’ avoidance (dashed curves) relative to the 21-cm PS estimates from the full simulation volume (solid curves). As expected, in avoiding the ‘wedge’ region, we limit the total number of Fourier modes available in each spherical shell to measure the 21-cm PS. Thus, quite simply, our Poisson uncertainty increases owing to the smaller number of available modes driving the increasing uncertainty error. Further, as the number of available modes quickly diminishes as we move to larger spatial scales (smaller $k$) the cosmic variance uncertainty rapidly increases. For smaller scales, the cosmic variance uncertainty still becomes dominated by the Trispectrum component (as indicated by similar upticks in amplitude at $k>1$ Mpc^-1 as highlighted in Figure 3), though on slightly smaller scales than found previously. Nevertheless, the non-Gaussianity in the cosmic variance uncertainty is still important when performing ‘wedge’ avoidance.

Importantly, the differences between the two cases are amplified by our survey geometry. In general, our simulation volumes are considerably larger along the transverse direction relative to the line-of-sight. Thus, previously, when measuring the 21-cm PS over the full volume, our number statistics for the $k_{\perp}$ modes were considerably larger than for the $k_{\parallel}$ modes. Now, when measuring the 21-cm PS using foreground avoidance, most of these $k_{\perp}$ modes are excluded since they fall within the ‘wedge’.

In following this avoidance strategy, we also significantly limit the largest scale modes accessible to our 21-cm PS measurement. In all cases, the 21-cm PS is only measurable down to $k\sim 0.1$ Mpc^-1, with no measurable modes beyond this limit. This limit is simply set by the adopted size of our buffer ($b=0.1\,h$ Mpc^-1) from Equation 12 which is somewhat arbitrary. In principle, it can be reduced by considering either larger bandwidths or with an improved understanding of our systematics. Note however, that the smallest accessible scale for each instrument is not exactly $k\sim 0.1$ Mpc^-1, but rather slightly larger. As we additionally ignore all $k_{\perp}=0$ modes, from Equation 12 the minimum scale is simply $b$ plus $m$ times the smallest non-zero $k_{\perp}$ in our simulation volume. Thus, the minimum $k$ for the MWA is smaller than the SKA due to the smaller $k_{\perp}$ (larger transverse scale).

In Figure 5 we present the fractional increase in cosmic variance when removing the foreground contaminated wedge modes. We find that the relative increase is comparable across all scales irrespective of the total simulation volume, highlighting that the increase in the cosmic variance is primarily driven by the fraction of excised Fourier modes (i.e. increased Poisson noise in each bin due to having fewer modes free of foregrounds).

Previously, we explored the impact of foreground avoidance on the resultant cosmic variance uncertainty on a 21-cm PS measurement. Next, we consider whether the 21-cm PS measured following foreground avoidance can be biased relative to the 21-cm PS measured from a full simulation volume. This is particularly pertinent given that recent upper limits on the 21-cm PS measured by LOFAR (Mertens et al., 2020), MWA (Trott et al., 2020) and HERA (The HERA Collaboration et al., 2022b) have for the first time allowed Bayesian astrophysical parameter inference to be employed to understand the physics of reionisation (Ghara et al., 2020; Ghara et al., 2021; Greig et al., 2021a; Greig et al., 2021b; Mondal et al., 2020; The HERA Collaboration et al., 2022a). However, when these inference pipelines are applied to these upper-limits, the simulated 21-cm PS is generated from the full simulation volume, rather than the characteristics of the observation, such as foreground avoidance.

While likely irrelevant for the present upper-limits since these are not very aggressive and only disfavour relatively extreme models of reionisation. Nevertheless, as we approach a first detection, these inference pipelines will be vital to building our understanding of the galactic physics driving reionisation. Thus, it is important to explore the relative impact of any biases on the modelled 21-cm PS when not accounting for the characteristics of the observation.

Firstly, the observed 21-cm PS is asymmetric owing to redshift-space distortions from the peculiar motions of the gas along the line-of-sight. Thus unlike in the earlier sections, we must include the effects of redshift-space distortions on our estimate of the 21-cm brightness temperature signal (note, the inference pipelines above include these effects). These redshift-space distortions increase the relative power along the line-of-sight direction, which could amplify any biases when measuring the 21-cm PS following foreground avoidance. This is because modes are preferentially excised transverse to the line-of-sight, resulting in a higher amplitude signal after binning over all possible modes within spherical shells. Thus here, we compute the 21-cm brightness temperature signal for the full simulation following the multiplication of Equation 9 by the pre-factor $\left(1+\frac{{\rm d}v_{\rm r}}{H{\rm d}r}\right)^{-1}$ and computing the anisotropic 21-cm PS.

In Figure 6 we present the mean (solid curve) and the 68th and 95th percentiles (dark and light, respectively) for the ratio of the 21-cm PS with the ‘wedge’ removed (avoidance) relative to the 21-cm PS from the full simulation volume. Here, these are determined using the same sub-volumes as obtained previously for each individual experiment from the full 7.5 Gpc simulation. It is clearly evident that for all survey areas, the 21-cm PS measured following foreground avoidance returns a biased estimate of the true 21-cm PS obtained from the full simulation volume. On average, we find the relative amplitude of this bias to be relatively modest, at $\sim 10-20$ per cent over all measured $k$-scales. On the largest scales (small $k$), when foreground avoidance has limited sampling of the Fourier modes (see Figure 4) this can result in increases/decreases in the 21-cm PS amplitude of larger than $40$ per cent.

Although fairly modest in amplitude, it is important to note that these results are model dependent. That is, this bias corresponds specifically to our particular source parameterisation and frequency band ($z\sim 7.5-9.4$). Therefore, it is plausible that the amplitude of this bias could increase under different astrophysical models or stages of reionisation. Hence, this could have important consequences for parameter inference studies if the statistics of the simulated models are not consistent with the characteristics of the observation, potentially biasing the inferred astrophysical parameters. Although the model dependency of the work is important to consider here, in previous sections it is less relevant. When investigating the cosmic variance uncertainty the relative amplitude of the uncertainty is dependent on the number of modes within each spherical shell when estimating the 21-cm PS, not the 21-cm PS amplitude itself.

To more clearly illustrate the relative increase in amplitude in the 21-cm PS following foreground avoidance, in Figure 7 we show in the same panel the mean ratio of the 21-cm PS with/without foreground avoidance for all interferometer experiments. This clearly demonstrates that the relative increase in the 21-cm PS amplitude following foreground avoidance is consistent at $10-20$ per cent across all the different survey footprints.

This $10-20$ per cent bias is broadly consistent with a similar investigation performed by Jensen et al. (2016). Here, this bias was determined for two astrophysical models over the full reionisation history. For their fiducial model, which more closely matches our model, at a similar stage of reionisation, $\bar{x_{\mathrm{H\,{\scriptscriptstyle I}}{}}}\sim 0.5$, they find a similar bias of $10-20$ per cent, however, as a decrement rather than an amplification. Likely, these differences arise due to the specific modelling of the ionising sources, but more likely due to the simplifications present in our modelling (e.g. linear theory) relative to their $N$-body simulations. Importantly, though, they only consider a single realisation for this exploration. For sub-volumes equivalent in size to their single simulation (corresponding to the LOFAR field-of-view) in the lower left panel of Figure 6 we can accommodate individual models with a negative bias. Factoring in the modelling differences mentioned above, we are confident that these two works are broadly consistent.

Note, there also appears to be a slight gradient in the relative bias for the different experiments. The SKA, with the smallest survey footprint has the lowest amplitude whereas the MWA, with the largest footprint has the highest amplitude. However, this is simply due to the geometry of the survey volume. Foreground avoidance with the MWA removes on average more Fourier modes in each spherical shell below the ‘wedge’ than that of the SKA relative to the 21-cm PS computed from the full simulation volume. This is because the MWA has a considerably larger transverse size (higher sampling of transverse modes) which results in a higher fraction of modes cut at the same fixed $k$ than the SKA. Equally, the uptick at the largest scales occurs earliest for the SKA compared to the other, larger area surveys. This occurs for the same reason as discussed earlier, the removal of the $k_{\perp}=0$ modes coupled with the $k_{\parallel}$ needing to satisfy Equation 12.

Typically, when it is infeasible to generate a single, large volume simulation, the requisite size and volume can be achieved by stitching together many smaller simulations generated using periodic boxes (see e.g. Nasirudin et al., 2020). Although sufficient volumes are achieved in this manner, the statistics of the simulations will only be accurate on the scales modelled by the smaller simulation. That is, it is impossible to inject information on scales beyond the largest scale of the smallest simulation. Using our 7.5 Gpc simulation, we can explore the relative impact of tiling simulations on the large-scale information in the 21-cm PS.

The tiling of periodic simulations additionally leads to repeating structures on scales characteristic of the simulation volume. When considering discrete objects such as galaxies in a redshift survey, these repeating structures can be minimised through rotating the smaller simulations prior to tiling. Unfortunately, since the 21-cm signal arises from the neutral gas, it is a continuous rather than discrete quantity preventing the usage of rotations. Rotations would lead to edge effects and disjoint ionised regions on the boundaries of the smaller simulations leading to strange artefacts in the statistics. Thus, for 21-cm studies, we can only consider tiled simulations.

To investigate the impact on the large-scale power due to tiling, we consider three smaller volumes and tile each to match our single 7.5 Gpc, 6400 voxel simulation volume. We generate 50 independent realisations of the smaller simulations with 200, 400 and 600 voxels per side length, corresponding to lengths of $\sim 234.4,468.8$ and 703.1 Mpc. We tile the same individual simulation to the large volume footprint, thus we obtain 50 independent large volume tiled simulations. For each, we generate the full 21-cm light-cone and only consider a line-of-sight depth of 30 MHz, consistent with large volume simulation.

In Figure 8 we compare the 21-cm PS obtained from our 7.5 Gpc simulation (black curve) to the 21-cm PS obtained from an equivalent tiled 7.5 Gpc simulation using the 200 (red), 400 (blue) and 600 (purple) voxel smaller simulations. In each, the vertical dashed line corresponds to the largest scale (fundamental frequency) sampled by the smaller simulations ($k_{f}=2\pi/L_{\rm small}$). The shaded region corresponds to the 68th percentile scatter obtained from our 50 independent realisations. As expected, the 21-cm power drops away on scales below the fundamental frequency. This highlights that information on modes not originally simulated by the smaller volume simulations cannot be artificially created to match the full large volume simulation. Thus, a tiled simulation volume has a limited regime of validity restricted by the scale of the smaller simulations.

On scales just above the fundamental frequency, we find the 21-cm power in the tiled box underestimates the power from the full simulation. In some instances this can be as much as a factor of $\sim 2-3$. Although the physical scale of these spherically averaged Fourier modes are sampled by the original smaller volume simulation, on these scales we are still missing the contribution from Fourier modes which contain at least one component mode longer than the small volume simulation (e.g. $k_{x},k_{y}$ or $k_{z}<k_{f}$). These individual modes contain no power, but contribute to the total number of modes within the bin, resulting in an overall reduction to the averaged power in that bin. Thus, the 21-cm power from the tiled box will underestimate the full box power on all scales until the fraction of Fourier modes with individual components larger than the smaller box is sufficiently small. Therefore, care must be taken to adequately produce sufficiently sized smaller volumes prior to tiling to minimise the loss of power on the largest scales.

In addition to underestimating the 21-cm power following the tiling of smaller periodic simulations, it is also important to note that we will underestimate the cosmic variance uncertainty in our tiled simulation volume. Since in the act of tiling we produce repeated copies of the spatial information from our smaller simulations we cannot artificially increase the variance within a sampled Fourier bin. That is, the measured sample variance is fixed to the variance (i.e. Poisson fluctuations) measured from the smaller volume simulation. However, since we have scaled up our total volume, following from Equation 11 the sample variance decreases as $\propto 1/V$. Thus our estimated cosmic variance from our tiled, large-volume simulation must be scaled by the increased volume. In effect, the tiled simulation underestimates the cosmic variance error by $\propto\sqrt{N_{\perp}^{2}\times N_{\parallel}}$, where $N_{\perp}$ is the number of copies required to achieve the desired transverse scale and $N_{\parallel}$ is the number of copies required along the line-of-sight.

In Figure 9 we highlight this scaling by showing the ratio of the cosmic variance uncertainty from the full simulated volume compared to that obtained from tiling the smaller volume simulations to achieve the same volume. In each, the tiled simulation cosmic variance underestimates the true cosmic variance, as indicated by the ratio sitting well above unity. The dashed horizontal lines correspond to the expected scaling ($\propto\sqrt{N_{\perp}^{2}\times N_{\parallel}}$) by which we underestimate the true cosmic variance in each tiled simulation volume. That is, the cosmic variance in the tiled simulation must be decreased by this amount to match the true cosmic variance from the full simulation volume. For larger $k$ (i.e. $k>0.5$ Mpc^-1) the impact of the non-Gaussian (Trispectrum) contribution to the cosmic variance becomes increasingly important, modifying this simple scaling relation.

Importantly, here we have only considered the power extracted from the full simulation volume. However, if one also considers the instrumental response of the beam (i.e. beam multiplied by the tiled signal), then additional artefacts in the estimated power could be expected due to mode mixing and/or power leaking in from the side-lobes. We leave an investigation into these effects to future work.