Impact of Instrument Responses on the Detectability of Onepoint Statistics from Redshifted 21 cm Observations
Abstract
We study the impact of instrumental systematics on onepoint statistics (variance, skewness, and kurtosis) of redshifted 21 cm intensity fluctuation observations from the Epoch of Reionization. We simulate realistic 21 cm observations based on the Murchison Widefield Array (MWA) Phase I reionization experiment, using the array’s point spread function (PSF) and antenna beam patterns, fullsky 21 cm models, and the FHD imaging pipeline. We measure the observed redshift evolution of pixel probability density functions (PDF) and onepoint statistics from the simulated maps, comparing them to the measurements derived from simpler simulations that represent the instrument PSFs with Gaussian kernels. We find that both methods yield onepoint statistics with similar trends with greater than 80% correlation for all statistics. We perform additional simulations based on the Hydrogen Epoch of Reionization Array (HERA), using Gaussian kernels as the instrument PSFs, and study the effect of frequency binning on the statistics. We find that PSF smoothing and sampling variance from measuring the statistics over limited field of view dilute intrinsic features and add fluctuations to the statistics but, at the same time, reveal new detectable features. Particularly, observed kurtosis will increase when a few extremely high or low temperature regions are present in the maps. Frequency binning reduces the thermal uncertainty but can also blur regions along the frequency dimension, resulting in kurtosis peaks that only appear in statistics derived from maps of certain frequency bins. We further find that the kurtosis peaks will reach their maxima when the angular resolution of the PSFs match the size scale of the extreme regions that produce the peaks. The HERA array should be capable of charting the evolution of the observed skewness and kurtosis of the 21 cm fluctuations with high sensitivity while the MWA Phase I will likely only be capable of detecting the peak in variance.
1. Introduction
Considerable effort is underway to constrain the Epoch of Reionization (EoR), the era when radiation from the first stars and galaxies transformed gas in the intergalctic medium (IGM) from neutral to ionized. Observations of the Lyman forest in highredshift quasars have set a limit to the end of reionization of (Fan et al., 2006), and measurements of the cosmic microwave background (CMB) optical depth by the Planck experiment have indicated that reionization is still progressing at (Planck Collaboration et al., 2016).
The 21 cm emission from the hyperfine transition in the ground state of neutral hydrogen is arguably the most direct probe to detect the EoR (Sunyaev & Zeldovich, 1972; Scott & Rees, 1990; Madau et al., 1997; Tozzi et al., 2000; Iliev et al., 2002). Fullsky observations of 21 cm spectra, redshifted to meterwave, will produce tomographic maps of neutral hydrogen throughout the reionization era and beyond, allowing the study of the evolution of this structure and its implication for the underlying ionizing sources.
Many telescopes have been built or are being built to conduct experiments hoping to map this signal; these include the MWA (Murchison Widefield Array; Tingay et al. (2013); Bowman et al. (2013)), LOFAR (Low Frequency Array), PAPER (Donald C. Backer Precision Array for Probing the Epoch of Reionization; Parsons et al. (2010)), HERA (Hydrogen Epoch of Reionization Array; DeBoer et al. (2016)), and the SKA (Square Kilometer Array; Mellema et al. (2013)).
These telescopes utilize many compact antennas to yield wide fields of view and point spread functions (PSF) with approximately arcminute angular resolutions. These characteristics are ideal for EoR tomographic mapping but also give rise to widefield beam chromaticity and strong sidelobe interferences that complicate mitigation of bright astrophysical foregrounds.
Thus, much attention has been focused on the statistical analysis of redshifted 21 cm EoR observations, in particular with power spectrum measurements. Preliminary results from current observations have recently been released (Beardsley et al., 2016; Trott et al., 2016; Dillon et al., 2015; Ali et al., 2015; Pober et al., 2015; Jacobs et al., 2015), including robust characterization of the foreground contamination and instrumental systematics in the power spectrum (Thyagarajan et al., 2015a, b; Jacobs et al., 2013). Lessons from the first generation of experiments have led to the design of HERA, which will be a highlypacked, redundantbaseline array that is optimized for the power spectrum analysis, while preserving imaging performance on subdegree scales.
As reionization progresses, ionized regions (H II) will form around groups of sources with highenergy UV radiation, causing the distribution of 21 cm intensity field to deviate from Gaussian, and the power spectrum alone will be insufficient to fully describe the signal (Lidz et al., 2007). This limitation has led to studies of several other statistics.
One promising alternative is the onepoint probability distribution function (PDF) and higherorder moments (variance and skewness) of the 21 cm brightness temperature fluctuations. Simulations have suggested strong evolution of the PDF throughout the reionization redshifts (Furlanetto et al., 2004), and the nonGaussianity in the distribution can be captured by measuring the skewness of the PDF (Wyithe & Morales, 2007; Harker et al., 2009). Ichikawa et al. (2010) has also developed a maximumlikelihood method to directly measure the shape of the PDF. Gluscevic & Barkana (2010) further used such a method to differentiate the PDF from six different reionization models. Recently, Watkinson & Pritchard (2014, 2015) and Watkinson et al. (2015) studied the sensitivities of the variance and skewness on different reionization scenarios, including the global and local reionization models, sinks of reionization radiation, and Lyman coupling and Xray heating. They suggest that signatures in the statistics could be detected by upcoming instruments. Although not discussed here, other efforts to constrain nonGaussianity in the 21 cm signal have explore morecomplicated statistics, including the twopoint difference probability density function (Barkana & Loeb (2008)) and the 21 cm bispectrum (e.g. Saiyad Ali et al. (2006), Pillepich et al. (2007), Yoshiura et al. (2015), and Shimabukuro et al. (2016) among others).
Yet characterization of foregrounds and instrumental systematics on the onepoint statistics have not been well studied. The complicated PSFs of the reionization arrays are often ignored from the analysis or approximated by convolving the model images with ideal Gaussian PSF. Noise from foreground residuals and other systematics are usually added to the simulation as Gaussian random noise or analytically derived from theory. Real observations also pose other limitations such as the limited field of view, the effect of frequency bandwidth, and the evolution of the signal along the frequency dimension. These issues must be robustly studied before predictions can be made for the upcoming data.
In this paper, we study some of the above issues in detail by performing realistic simulations that closely mimic actual EoR observations, combining fullsky 21 cm models, readily available radio array simulators, and other software packages from the data reduction pipelines of the MWA and HERA. We focus our analysis and discussion on the detectability of the variance, skewness and kurtosis, taking into account the impact of instrument PSF, field of view, and frequency bandwidth. We describe our simulation method for the MWA in Section 2. We present our results and discuss the effect of the instrument PSF and field of view on the statistics in Section 3. In Section 4, we build HERA simulations based on results from previous sections and discuss the impact of frequency bandwidth on the statistics, as well as the potential of using kurtosis to detect ionized regions in 21 cm maps.
2. Simulations
We construct our MWA simulation pipeline based on existing 21 cm models, a telescope simulator and the actual MWA imaging pipeline. We describe each component in the following subsections. All simulations are noiseless to simplify the interpretation of the results although mathematically derived uncertainties are included in the analysis. We also ignore foreground contamination, postponing it to a future work. Although we only discuss the MWA simulation here, our simulation pipeline can be adapted to different instruments.
2.1. FullSky 21 cm Models
Neglecting foregrounds, the EoR sky consists of 21 cm brightness temperature intensity fluctuations (), which are characterized by the 21 cm spin temperature (), the CMB temperature (), the fractional overdensity of baryons , and the gradient of the proper velocity along the line of sight (),
(1) 
The current stateoftheart 21 cm simulators such as the 21cmFAST use semianalytic methods to produce threedimensional cubes of 21 cm brightness temperature fluctuation at different redshifts (Mesinger et al., 2010). These cubes are represented in rectangular comoving coordinates and can be up to a couple of (Gpc/h) in size, roughly equivalent to deg volume at relevant redshifts.
In contrast, EoR observations produce treedimensional data where two of the dimensions map the spatial dimensions of the sky in sine projected coordinates and one dimension conflates the telescopes’ lineofsight distances with time and redshift/frequency. The MWA observes the EoR from 138.915 MHz to 195.235 MHz with 80 kHz frequency channel resolution, equivalent to 704 maps that span . Thus, the observed data form a single ”lightcone” cube, where slices across the frequency dimension evolve with redshift.
In order to match existing 21 cm models to an instrumental observation, we transform the fourdimensional (three spatial and one redshift) outputs of theoretical 21 cm simulations into a threedimensional (two angular and one frequency) 21 cm cube. We use a tileandgrid method that maps the simulation cubes from different redshifts to fullsky maps in HEALPix^{1}^{1}1http://healpix.jpl.nasa.gov coordinates with NSIDE=4096, capturing the full spatial structure and evolution of the signal in frequencyspace.
We start with an original suite of “raw” 21 cm cubes from the semianalytic simulations of Malloy & Lidz (2013). The cubes span the redshift range , with steps, representing the universe from to ionized respectively. The simulation volume is in a pixel box with periodic boundary. We linearly interpolate the simulated 21 cm cubes across redshift to produce new cubes that match the redshifts observed by the MWA at each of its frequency channels, yielding 704 cubes. For each of the 704 redshifts, we tile the interpolated cube with itself to construct an arbitrarily large simulated volume for that redshift. Then, we draw an observable sky at that redshift as a sphere of radius equal to the comoving distance of that redshift from a fixed origin inside the volume and linearly interpolate the four nearest neighboring pixels from the cube to the corresponding HEALPix pixel location on the sphere. The HEALPix pixel area is times smaller than the resolution of the simulations to avoid Nyquist sampling. This process is repeated for every observed redshift to produce a suit of fullsky maps that accurately represent the 21 cm lightcone model. The HEALPix maps are later reprojected into halfsky sine projected maps with pixel area arcminute and passed to our observation simulator.
2.2. Instrument Model
The main instrumental effects of a radio interferometer result from the array’s PSF response which is both directional and spectral dependent. The size and shape of the PSF depend on the response of individual antenna elements (primary beam), the antenna layout, and the frequency and pointing of the observations.
We model the PSF effects of the MWA Phase I (Tingay et al., 2013), which consists of 128 phasedarray antenna tiles arranged over a region with a diameter of 3 km. Each tile is a fourbyfour grid of dualpolarization dipoles. Nearly half of the antenna tiles are contained in a compact core with a diameter of 100 meters. EoR observations only use this compact core to improve EoR power spectrum sensitivity, foreground subtraction and calibration (Beardsley et al., 2013; Bowman et al., 2009). This yields naturallyweighted PSF resolution degree across the MWA EoR band. The primary beam extends across the whole sky, but the usable main lobe is about 30 degrees diameter centered around the pointing center.
The MWA acquires EoR data in a driftandshift mode. The telescope is pointed at a specific coordinate in a radio quiet region of the sky (for example, h, ) and observes as the sky drifts throughs the primary beam. Then the telescope is repointed, and the process repeats. We adopt a simpler singlepoint observation where the telescope is always pointed at the zenith for our simulation. We model the limiting case of a single snapshot image. This results in the worstcase PSF for the MWA with no rotation synthesis to improve UV coverage.
We use the MIT Array Performance Simulator (MAPS^{2}^{2}2http://www.haystack.mit.edu/ast/arrays/maps/) to simulate MWA visibilities. MAPS generates beam responses from userconfigurable antenna configurations, array configurations and observing parameters, and performs convolution with sky models. We configure MAPS to use the MWA 128tile array layout, using a short crossdipole on an infinite ground plane as an antenna model. MAPS also applies wprojection to correct for widefield effects. We run MAPS with our fullsky models for each observing frequency to take into account spectral dependency of the PSF. Data is integrated for two seconds.
The visibilities are then regridded and Fourier transformed by FFT into “dirty” maps of brightness temperature fluctuation with Fast Holographic Deconvolution (FHD) (Sullivan et al., 2012), one of the two primary calibration and imaging pipelines of the MWA (Jacobs et al., 2016). We assume that the derived dirty maps from FHD are the final data products that are used for EoR analysis. This is consistent with data from real observations in an absence of foregrounds. In principle, a deconvolution through a matrix inversion or an iterative algorithm such as CLEAN can be utilized to remove the PSF structure and recover the âtrueâ intensity map. However, the sky structures in maps from reionization arrays are too complicated and not sufficiently sparse to utilize such an approach. For current EoR experiments, foregroundsubtracted dirty images are generally where derived statistics are calculated.
2.3. Gaussian Approximation of PSFs
The use of idealized Gaussian models to represent telescope PSFs is common in previous studies. Gaussian kernels with FWHM equivalent to the angular resolution of the telescopes’ PSFs are assumed as the telescope PSFs and convolved with sky models to produce “observed” intensity maps. This method neglects sidelobe responses and directional dependency of the PSFs. Figure 1 shows these differences for the MWA. The grayed solid line shows a cross section of the PSF of the full MWA Phase I array that includes all antennas. The dotted line shows the PSF of the compact core that only includes antennas within the 100 meter diameter, which is the PSF used in EoR anlysis and in this work. The solidline with squares shows a Gaussian fit to the main lobe of the compact core PSF that would be used in a simple convolution method. To investigate the differences between the more realistic and idealized PSFs on our analysis, we carry out simulations with both the MWA Phase I Core PSF, using FHD and MAPS, and its fitted Gaussian PSF, using a simple convolution. We fit a Gaussian kernel to the main lobe of the MWA Phase I Core PSFs derived from the beam configuration in FHD. Then, we convolve our twodimensional Gaussian kernel to the sine projected 21 cm maps that we input to the MAPSFHD pipeline. This produces similar output maps to the pipeline but without the the effects of the full instrumental PSF sidelobes..
2.4. OnePoint Statistics
Our quantities of interest are the onepoint statistics of brightness temperature fluctuations in the dirty maps, in particular variance, skewness and kurtosis. For a map with data value and mean value , the pth order statistical moments are defined as,
(2) 
The variance (), skewness () and kurtosis^{3}^{3}3In statistics, the precise term of this definition is excess kurtosis, which subtracts 3 from the standard definition of kurtosis, , to yield zero for Gaussian distribution. Here, we simply use kurtosis to refer to excess kurtosis. () are standardizations of the 2nd, 3rd and 4th moments,
(3)  
(4)  
(5) 
Along with the mean, these three quantities describe simple deviations in the shape of the brightness temperature PDF relative to a standard Gaussian PDF. The variance measures the spread of the PDF. Both skewness and kurtosis describe the outliers, or the tails, of the PDF in different ways. Skewness measures asymmetry of the outliers, in which positive and negative skewness values indicate that the values of the outliers are greater and less than the mean value that would be the case for a Gaussian distribution. Kurtosis describes the heaviness, or density, of the outliers; a positive kurtosis indicates more outliers whereas a negative kurtosis indicate less outliers. Thus, a PDF with high kurtosis will look more peaked in comparison to a Gaussian PDF with the same variance. A perfect Gaussian PDF has zero skewness and kurtosis.
We measure the PDF, variance, skewness and kurtosis from both the FHD and Gaussian simulated dirty maps, as well as from the residual maps, defined as the difference between the FHD and Gaussian maps. In real measurements, signaltonoise ratio decreases substantially near the edge of the field due to instrument response. Therefore, we only use pixels within the FWHM of the primary beam to construct the PDFs and statistics. All PDFs are calculated using 60 linearly spaced bins from the minimum to maximum pixel values within each map.
It is common in radio astronomy to oversample the PSF resolution and produce dirty maps with higher angular resolution than the PSF’s angular resolution. The oversampled pixels contain no extra information as features smaller than the size of the PSF are filtered out from PSF smoothing. Mathematically, the oversampling of the PSF has no effect on the statistics measured in this work. Thus, we use all pixels in the dirty maps that fall within the FWHM to calculate onepoint statistics. In Section 4, we will calculate thermal noise uncertainties, and in those calculations, we take care to use the proper number of independent PSF samples in the maps.
As a reference, we plot the variance (dashed line), skewness (dotdash line) and kurtosis (dotted line) calculated from the raw 21 cm input model using all pixels, along with the ionized fraction of the cubes (solid line), as a function of frequency and redshift in Figure 2. The left yaxis shows corresponding statistical values, whereas the right yaxis shows the ionized fraction. Although not shown here, we find that statistics derived from the HEALPix and sine projected outputs of our gridandtile extrapolation closely match statistics from the raw model. The average fractional errors between statistics derived from the raw model and the HEALPix maps are less than 0.1% with Pearson Correlation coefficients (PCC) greater than 0.99 (99% correlated) for all statistics. Statistics derived from the sine projected maps diverge a bit more from the raw model with average fractional errors of 12%, 0.8% and 3.7% for variance, skewness and kurtosis but still with PCC greater than 0.99. Specifically, the variance and kurtosis decrease due to interpolation effects. We will discuss and compare features in the model statistics with statistics derived from our simulations in Section 3.1.
3. Effects of PSF
3.1. Observable Features in OnePoint Statistics
Figure 3 shows example PDFs derived from the FHD and Gaussian simulations at () along with the corresponding simulated maps. Two effects are readily apparent. First, the lack of zerospacing baselines in the MWA forces the observed sky brightness mean to equal zero. As a result, all PDFs derived from simulated observations are offset from the PDFs of the input theoretical 21 cm model. Second, the simulated observation PDFs are also more Gaussianlike than the input model. Since there is no thermal noise in our simulations, this effect arises solely from the PSF smoothing that reduces the amplitudes of the intrinsic 21 cm fluctuations and blurs the relatively sharp transitions from ionized to neutral IGM in the input model as seen in the simulated maps.
Figure 4 shows the variance, skewness and kurtosis derived from the FHD (solid line) and Gaussian (dashed line) simulations as a function of ionization fraction (or frequency). The statistics derived from the simulated maps exhibit fluctuations throughout the measured redshifts, as opposed to smoothly varying measurements in the models, due to the sample variance inherent in simulating an observation of a single field, rather than modeling a statistical ensemble.
The variance is seen to steadily increase as reionization progresses, reaching a maximum around ionized fraction , before falling off as reionization concludes. The intrinsic variance of the input model shows similar evolution, but reaches the highest point earlier, around (see Figure 2). We interpret the shift in the variance between input and output as a result of the PSF smoothing. The peak in the variance from the simulated observations does not occur until the ionized bubble size grows large enough to match the MWA Phase I Core’s angular resolution. This occurs later than the intrinsic peak in the input model. We can relate this effect to power spectrum measurements. The MWA Phase I Core only samples sizescales corresponding to low kmodes in Fourier space. Hence, it is when these modes reach their peak power that determines when the variance peaks in the MWA observations. The instrument is acting analogously to a matchedfilter to identify signal features of a particular size scale. In general, the redshift of the observed peak will be dependent on a telescope’s angular resolution, and hence, may occur earlier (at higher redshift) for larger telescopes with higher angular resolution since ionized bubbles are smaller when the intrinsic variance in the signal peaks . The PSF of the telescope also smooths out 21 cm fluctuations, lowering the signal amplitude. Thus, the observed variance peak has a lower amplitude than the intrinsic variance peak in the model.
The skewness derived from both simulated observations shows little evolution for most of reionization, remaining slightly negative around and mean values with standard deviation of and for FHD and Gaussian cases respectively. At the end of reionization, however, it becomes increasingly positive until the last redshift in the simulation. In contrast, the input model shows highly negative skewness early in the reionization due to the high filling factor of neutral hydrogen and rare regions of low emission, creating an offset positive peak in the PDF with a long tail extending toward zero emission. The model skewness becomes positive around when large ionized regions start to form and the PDF begins to become bimodal. The transition to positive skewness occurs later in the simulated observations, beginning at approximately the same time as the observed variance peaks, indicating that the skewness is diluted by the telescope’s PSF smoothing until the ionized regions grow larger than the angular resolution.
The kurtosis also shows little evolution throughout reionization, oscillating around and mean with standard deviation of and for FHD and Gaussian cases, but does exhibit a few occasional spikes that rise above the standard deviation. For example, at and . Upon examining the maps, we find that these spikes occur when the observed field contains a few large regions of extremely low or high brightness temperature (relative to the mean at that frequency). These outlier regions add more density to the tails of the PDFs, resulting in higher kurtosis.
Taking the nearzero values of the skewness and kurtosis together indicates that the observed distribution remains nearly Gaussian throughout most of reionization, in contrast to the evolving distribution intrinsic to the input model. Hence, for the MWA Phase I, the power spectrum is a sufficient measurement to fully describe the information that is retrievable by the telescope from the underlying 21 cm signal, except at the end of reionization.
3.2. Comparison of Gaussian and FHD PSFs
It is evident from Figures 3 and 4 that the PDFs and onepoint statistics resulting from the FHDsimulated maps and Gaussiansmoothed maps are qualitatively similar. Both PSFs produce similarly shaped PDFs. However, the FHDderived PDF is consistently broader than the Gaussianderived PDF and its variance is consistently higher by approximately . We interpret the increased variance in the FHD case to be due to the complicated sidelobe structure inherent in the full instrumental PSF. This structure results in a sidelobe confusion effect that scatters power throughout the map, raising the overall variance and diluting the signal’s skewness and kurtosis.
Despite the discrepency in variance between the two cases, the dependence on ionization fraction of each of the statistics is highly correlated between the two cases. We quantify the similarity by calculating the Pearson Correlation Coefficients (PCC) between the two cases for each onepoint statistic. We find coefficients of for variance, for skewness, and for kurtosis.
We attempt to further quantify the effect of sidelobe confusion by measuring statistics from the residual maps derived by subtracting FHD maps from Gaussian maps. The residual maps contain features from sidelobe confusion and other effects inherent in the full instrumental PSFs. Statistics derived from the residual maps are plotted in Figure 4 as dotted lines along with those from the FHD and Gaussian simulated maps. Variance of the residual maps closely follow variance from the Gaussian convolved maps as expected from the fact that the FHDderive variance is 56% higher than the Gaussianderive variance. The residual skewness and kurtosis remains noiselike through the redshifts, oscillating around mean and standard deviation for skewness and mean and standard deviation for kurtosis.
Based on these findings, we conclude that Gaussiansmoothing kernels are reasonable approximations to full instrument simulations for denselypacked antenna arrays such as the MWA Phase I Core. The skewness and kurtosis estimates derived from Gaussian smoothed maps are likely sufficient for many statistical investigations. For variance estimates, a simple transfer function of the form can be used to compensate for the factional difference between the estimates derived from realistic simulations and morecommonly modeled idealized cases. Inspection of Figure 4 suggests that this function may reduce to a multiplicative correction factor such that , with in our case. We have not explored different 21 cm input models, but given that most theoretical models yield similar features, we expect the correction factor to be generally independent of the detailed properties of the 21 cm signal for a given instrument.
4. Detectability
Now we explore the detectability of the onepoint statistics. We will find below that the MWA Phase I does not have sufficient thermal sensitivity to detect skewness and kurtosis. Hence, building on our findings above that idealized Gaussian beam simulations are sufficient for dense aperture arrays, we perform new simulations based on HERA using only idealized Gaussian beams to focus our investigation on the detectability of the onepoint statistics by a moresensitive array. We begin in Section 4.1 with an overview of the HERA telescope and our simulations of it. We will then review our thermal uncertainty estimates in Section 4.2 and discuss the effect of frequency bandwidth in Section 4.3 before presenting our results.
4.1. Hera
HERA is a secondgeneration radio interferometer optimized for redshifted 21 cm power spectrum detection. Presently under construction, HERA builds on the lessonslearned from the MWA and PAPER. It consists of zenithpointed, 14meter diameter dishes fed by dualpolarization dipoles that are closely packed into a hexagon shape. The telescope has a 10degree primary field of view and will operate between 100 and 200 MHz. Nineteen dishes have been constructed to date as part of a staged buildout that will culminate in 331 dishes. DeBoer et al. (2016) provides additional details about the design of HERA.
We model the HERA instrument response using its FWHM fieldofview and an idealized Gaussian smoothing kernel, corresponding to the angular resolution of the array, to approximate its PSF. We perform the Gaussian simulation with the final planned 331dish HERA array, as well as a 37dish HERA array that will be operable by 2017. The latter happens to match the angular resolution of the MWA Phase I Core, allowing for a relatively direct comparison of the two. Table 1 summarizes our simulation parameters for the MWA Phase I Core and both HERA configurations.
Simulation Parameters  MWA  HERA37  HERA331 

Phase I Core  
Number of antennas  50  37  331 
Effective area (m)  21.5  153.86  153.86 
Maximum baseline (m)  100  42  140 
Integration time (h)  1000  100  100 
Angular Resolution (deg)  11.6  11.6  0.30.5 
FoV diameter (deg)  2035  711  711 
4.2. Thermal Uncertainty
To estimate the impact of thermal noise on recovering onepoint statistics with the MWA Phase I Core and HERA, we calculate thermal uncertainties for the statistics derived from our noiseless simulations. Our derivation of the thermal uncertainties follows Watkinson & Pritchard (2014). The thermal uncertainty in interferometric measurements () is associated with the instrument and observing parameters (Furlanetto et al., 2006).
(6) 
Here, is the array filling factor, defined as the ratio of the total effective area () and the maximum baseline length () of the array. is the spectral channel size, and is the integration time of the observation. is the system temperature of the array, which is dominated by the Galactic synchrotron radiation in redshifted 21 cm observations below 200 MHz. This becomes the main contribution to the noise, giving . Using this assumption, Equation 6 can be rewritten to estimate the thermal uncertainty in 21 cm observations,
(7) 
Assuming that independent Gaussian random noise with zero mean and standard deviation is added to every independent “beam” pixel in the simulated maps, the estimator variance of each pth order statistical moment resulting from the noise can be calculated and later propagated to variance, skewness and kurtosis. The uncertainty for each statistic is then just the square root of the estimator variance. The noise propagation is described in detail in the Appendix, where we rederive the equations given by Watkinson & Pritchard (2014), along with an additional derivation for the kurtosis.
In summary, the estimator variance of the 2nd, 3rd and 4th moments can be described by,
(8)  
(9)  
(10) 
and the estimator variance of skewness and kurtosis are as follows,
(11)  
(12) 
Here, is the number of independent PSF samples in the data, which is limited by the angular resolution of the telescope. Because our maps are oversampled, we calculate the number of independent PSF per pixel from the ratio of the pixel and the beam area, , and multiply this factor to the number of pixels in our sample to obtain . and are the estimator covariance of the 2nd and 3rd moments, and the 3rd and 4th moments respectively. The estimator covariance can be derived in the same manner as the estimator variance of the moments resulting in,
(13)  
(14) 
Figure 5 shows statistics derived from the Gaussian simulation for all three telescopes along with corresponding estimated thermal uncertainty. We begin with this figure showing the thermal uncertainty for the 80 kHz (roughly native) spectral resolution of the arrays. In this case, the narrow spectral channels result in substantial thermal noise. We address frequency binning to reduce the uncertainty in the next section. The higher angular resolution of HERA331, compared to HERA37 and the MWA Phase I core, smooths out less signal, especially early in reionization when feature sizes are smaller, leaving more outliers in the field. This results in an overall higher observed variance that peaks at lower frequency, specifically at 185.475 MHz compared to at 188.835 MHz for HERA37. Similarly, skewness and kurtosis derived from the HERA331 simulation show more fluctuation, especially earlier in reionization, due to more outliers being left in the observed field. The overall trend of skewness and kurtosis, however, remain the same. In contrast, the larger PSFs of both HERA37 and the MWA Phase I Core smooth out more signal, resulting in overall lower variance and fewer spikes in skewness and kurtosis..
The thermal uncertainty rapidly decreases as the reionization progresses regardless of the statistics derived from the telescopes. This is because the synchrotron radiation, which dominates the system temperature, is brighter at lower frequencies. In terms of telescope sensitivity, it is clear from Figure 5 that the MWA Phase I Core will not have enough sensitivity to detect any statistical features with the current observational parameters using 80 kHz spectral resolution, where as HERA331 array will be able to detect most features in all statistics with high signaltonoise ratio (SNR) due to the much improved sensitivity of the array. Even at 80 kHz spectral resolution, we see that HERA37 will also be able to detect variance and place some decent limits on skewness and kurtosis near the end of reionization despite a comparable angular resolution to the MWA Phase I Core due to significantly larger collecting area. We show in Section 4.3 that frequency binning will significantly improve telescope sensitivity for all statistics, including enable the MWA Phase I Core to detect features in the variance with sufficient SNR.
4.3. Frequency Binning and Windowing
It is evident from Equations 7–10 that the thermal uncertainty on the onepoint statistics is sensitive to both the spectral channel size () of the observed maps and the number of samples () used to calculate a given estimate of a statistic. Both of these parameters can be tuned to achieve the bestpossible signaltonoise ratio (SNR) for a particular instrument and 21 cm signal model combination. We explore the benefits of varying spectral channel size through frequency binning (averaging) of raw maps, as well as the effects of increasing the number of samples per estimate by grouping neighboring maps into the calculation of a onepoint statistic without frequency averaging. We term the latter process frequency windowing.
From Equation 7, we see that increasing the spectral channel size, , reduces the thermal noise in a given observation. Both MWA and HERA have relatively narrow raw spectral channels ( kHz) yielding high thermal noise in raw maps. However, we can bin observed maps in neighboring frequency channels to yield larger effective and lower thermal noise per binned map. Frequency binning does have potentially detrimental effects. It comes at the expense of reducing the number of samples that can be used to chart the evolution of the signal, and it can also reduce the observed 21 cm signal amplitude if it is performed over a range larger than the coherence length of the signal. Unexpected amplitude gain or loss can also occur as we shall see in the results.
Frequency windowing uses maps created at the raw spectral channel width (80 kHz in our cases), but calculates statistics using multiple maps for each estimate. The motivation for this method is the inclusion of more samples into each statistical calculation. The noise per pixel is unchanged from the raw maps as is the underlying signal contribution to each pixel, but using more pixels per estimate lowers the overall uncertainty on the estimate as shown in Equation 8–10. The disadvantage is that frequency windowing has the potential to dilute evolution in the signal similar to frequency binning.
To investigate the tradeoff between the improvement in thermal uncertainty and the potential loss of signal of the two methods, we group our MWA and HERA Gaussian simulations into trial cases of different total spectral window size, spanning 130 MHz with 1 MHz increment between each case. Then, we measure the variance, skewness and kurtosis, with and without binning the simulation sets, using the same method as described in Section 2.4, and calculate the corresponding thermal uncertainty as above.
Figures 6, 7 and 8 show the variance, skewness and kurtosis and their uncertainty derived from the MWA (dotted), HERA37 (dashed) and HERA331 (solid) simulations with 2, 3, 4 and 8 MHz. The right and left columns show results from the frequency binning and frequency windowing trials respectively. In addition, the top left panel in all three figures show results from our raw spectral channel size (0.08 MHz window), same as those in Figure 5. We select these particular frequencies to show several trends that are apparent in our results.
Frequency binning produces similar effects to the PSF smoothing but along the frequency axis. Binning reduces observed statistical variance particularly at the larger bin sizes by smoothing out the signal over the frequency scale corresponding to the bin sizes. To illustrate this effect, a lightcone map showing the evolution of signal in the HERA331 simulated observations is shown in Figure 9. The typical sizes of the signal along the frequency axis grow as reionization progresses, reaching the maximum size of MHz near the end of reionization. Binning the signal with bin size greater than this size scale will significantly reduce the amplitude of signal since there are few features that are larger than this size. The observed variance declines significantly between 4 and 8 MHz bins.
Skewness and kurtosis do not suffer overall amplitude reduction from the binning due to both statistics being sensitive to the relative shape of distribution, not the amplitude of the signal. However, binning smooths out the extremely low and high regions in the maps to neighboring frequency channels within the bin. Thus, the binned skewness and kurtosis appear to have fewer occasional fluctuation throughout. Binning can also reduce or amplify the overall amplitude of these regions depending on their appearances in the neighboring frequency bins and the bin location. A large peak appears in the kurtosis near 145 MHz when the bin size is 3–4 MHz but disappears when the bin size is smaller or larger. Interestingly, this peak appears only in the HERA331 simulation where the angular resolution is 0.5 degree but not the HERA37 or MWA Phase I Core simulations where the angular resolution is 1.5 degree. This implies a correlation between angular resolution and the amplitude of the kurtosis peaks. We explore its origin in the next section.
Compared to frequency binning, frequency windowing simply adds more samples to the statistics. Windowing over a large scale does not directly reduce the signal amplitude but rather convolves the intrinsic variance across a range of redshifts. Thus, we see that the variance amplitudes remain at similar levels regardless of window size. Apart from retaining the variance amplitude, frequency windowing retains large numbers of pixels in each estimate, reducing the susceptibility to sample variance within each map, and allowing the statistics to evolve smoothly with redshift. A good example is the skewness in Figure 7 in which frequency binning results in larger variation in the skewness from one redshift to the next while frequency windowing produces skewness that smoothly evolves.
An important question to ask is whether we should use frequency binning or frequency windowing and at what bin or window size for each array. For the MWA Phase I Core, Figure 6 suggests that frequency binning is preferable since it tends to yield greater signal to noise. For HERA, an overall look at the signaltonoise ratio for all cases provides a more robust answer. Figures 10 and 11 show color charts that depict the SNR of skewness and kurtosis, derived from the HERA37 and HERA331 simulations, as a function of bin/window size. Each block in this diagram represent the bin/window for each cases, with windowing cases in the left column, and binning cases in the right column. The scale is normalized such that the color is white at SNR=1, and becomes bluer and redder at SNR and SNR respectively. Based on the chart, we find that statistics derived from the frequency binning approach have lower uncertainty than the frequencywindowed counterparts for the same bin/window size. This should be true for all other arrays, including the MWA Phase I Core. A mathematical explanation is that the variance for pth order moment is approximately proportional to the noise uncertainty to the power of p and the reciprocal of the number of beam, . Equation 7 implies that , leading to . Binning raw channels will reduce the moment uncertainty by while windowing raw channel will reduce the moment uncertainty by just . Thus, we suggest that spectral binning be use to detect features in the statistics. Figures 10 and 11 also suggest that a bin size of MHz be used as averaging over wider frequency range does not significantly improve the signaltonoise. This is in agreement with the previously discussed point that features reach the the maximum size of MHz along the frequency axis, and averaging over a range larger than this number would only reduce the signal.
4.4. Correlation in the Observed Kurtosis Peaks
Throughout this work, we see occasional spikes that rise above the fluctuations caused by sample variance in the kurtosis. We find that these spikes occur when the observed field contains a few large regions of extremely low or high brightness temperature (relative to the mean at that frequency) (see Section 3.1). These outlier regions add more density to the tails of the PDFs, resulting in higher kurtosis. Binning multiple maps of neighboring frequencies that contain these outlier regions will smooth them over the frequency bin, resulting in a broader peak in the kurtosis. The amplitude of the peak will be reduced or amplified, depending on whether the size of the outlier regions along the frequency axis matches the frequency bin size. Thus, we see kurtosis peaks at 145.115 MHz and 178.555 MHz in the 2 and 3 MHz bin cases, but not in other cases. These particular peaks, however, only appear in the kurtosis derived from HERA331 simulation, but not from other telescope simulations (see Section 4.3. We also see more peaks in the kurtosis derived from the unbinned maps of HERA331 simulation in comparison to unbinned maps from other telescopes (see 4.1). These imply a correlation between kurtosis peaks, the angular resolution of PSF and the average bin size.
To understand more about this correlation, we perform additional simulations for intermediate HERA array sizes, specifically with 61, 91, 127, 169, 217 and 271 antennas, corresponding to approximately 1.1, 0.9, 0.8, 0.7, 0.6, 0.5 degree angular resolutions at 150 MHz, to give us more coverage over a range of angular resolution. We then recalculate the kurtosis from maps of all additional configurations, binning and windowing the maps with 130 MHz bin size as before. Noise parameters are adjusted to match the PSF resolutions.
As expected, we see the correlation between the kurtosis peaks and the angular resolution of the telescopes. Figure 12 illustrates the relation. The top panel shows the kurtosis derived from the 3 MHz windowed and binned maps of HERA37, HERA91, HERA127, HERA217 and HERA331 simulations (left to right). Other array configurations are left off for simplicity. The bottom two panels show maps from the simulation at 145.115 MHz (2nd row) and 178.555 MHz (3rd row) from the respective telescopes. In addition, the PSF size of each telescope is drawn as a white ellipse at the bottom left corner of each map. It is clear from the diagram that the amplitude of peak at 145.115 MHz rises as the sizes of outlier regions in the map start to match the size of the PSF. Since the size scale of the signal at 145.555 MHz is small, the kurtosis peak reaches its maximum when observed with the matched PSF of HERA331. Similarly, the peak at 178.555 MHz reaches its maximum when observed with HERA91 where the larger PSF of the telescope matches the typical size scale at that frequency. Although not shown here, we see the same trend in the kurtosis derived from maps of other bin sizes. For example, the two peaks near 171 and 179 MHz in the kurtosis derived from 2 MHz binned maps (2nd row, right, in Figure 8) reach their maximum when observed with HERA217 and HERA91 respectively, rising and falling according to the sizes of the PSF. The peak near 145 MHz in the 4 MHz binned maps also rises as the PSF resolution increases, reaching its maximum when observed with HERA331 similar to the peak at 15.115 MHz in the 3 MHz bin case.
5. Conclusion
For this work, we have developed realistic simulations of EoR observations that produce dirty maps of 21 cm brightness temperature intensity fluctuations. Our simulation incorporates fullsky 21 cm models and realistic beam response. We build fullsky 21 cm models with a tileandgrid method that projects and interpolates small 21 cm simulation cubes to lightcone gridded HEALPix maps. We base our simulations on the MWA Phase I Core configuration. Simple Gaussian smoothing is also performed for a comparison with additional PSF and noise parameters from HERA37 and HERA331 array. We measure the PDF, variance, skewness and kurtosis from the simulated maps and study the effect of frequency binning and windowing on the statistics. We also study the correlation between the size of the array point spread function and the peaks in observed kurtosis by performing additional simulations based on different HERA configurations. Uncertainty from thermal noise is mathematically derived and propagated to the statistics. Errors from foreground contamination and other systematics are ignored in this work, postponing them to future works.
Base on our findings we conclude with the following.

Gaussian PSF kernels can reasonably approximate the full instrument simulations for arrays with denselypack antenna configuration such as the MWA Core. High correlation is maintained between statistics derived from simulations with Gaussian and full PSFs although sidelobe confusion in the full PSF can scale the measured statistics by as much as 4060. A transfer function could potentially be used to rescale the signal, but more studies with various array types and 21 cm models are necessary.

Apart from spatially smoothing the signal by the PSF shape, the instrument acts analogously to a match filter that emphasizes a spatial size scale that match the angular resolution of the PSF. The observed variance will have lower overall amplitude and reach the maximum peak at a different redshift depending on the angular resolution of the PSF. Skewness and kurtosis remains near zero for most of reionization due to the underlying distribution becoming more Gaussianlike from PSF smoothing, except near the end of reionization where the intrinsic signal is highly bimodal.

Measuring the statistics over a limited field of view introduces sampling variance, resulting in statistics that oscillate over a range of redshifts as opposed to the smoothly varying curves derived from the models. Coupling with the sizescale emphasis of the PSF, this can result in strong peaks that rise above variation from sampling variance as we have seen in the kurtosis measurements throughout this work.

Similarly, frequency binning smoothes and emphasizes certain size scales along the frequency axis apart from improving the signaltonoise of the measurement. Thus, observed variance will have lower overall amplitude if measured from frequency binned maps, and certain peaks can appear in the kurtosis derived from a specific bin size but not another. In contrast, frequency windowing (grouping neighboring frequency channel without averaging) preserves the variance of the signal and reduces sampling variance in the observed statistics due to higher statistical samples in the windowed data. Frequency windowing, however, is less effective in reducing thermal noise.

The origin of the peaks in the observed kurtosis can be physically explained by considering the distribution of regions with extremely low or high brightness temperature, relative to the mean temperature at the redshift. When these regions dominate the observing field, observed kurtosis will increase due to higher density being added to the tails of the PDF of the measured signal. Similarly, blurring of these regions due to mismatched bin size or PSF size, will reduce the amplitude of the peaks. Hence, we see peaks only at specific frequency binning and angular resolution combinations. The correlation between the kurtosis peaks, frequency bin size and angular resolution of the telescopes could be very beneficial in quantifying 21 cm signal. A similar study to this work performed over wider range of these variables can potentially result in an empirical relationship that quantifies the statistical correlation between 21 cm model and observed number of large kurtosis and skewness features. Although these regions may not be the exact representations of the ionized or overdense regions, due to the instrumental smoothing of the complex underlying ionization structure, they can be statistically related as demonstrated in Beardsley et al. (2015).

In terms of detectability, our study suggests that the MWA Phase I Core should be able to detect the variance peak with SNR, while HERA, even with the smaller 37 array, will be able to detect all studied statistics with high sensitivity throughout the reionization redshifts. Frequency binning significantly improves the sensitivity for all detections, and our study suggests that the bin size of MHz be used to optimize SNR.
Appendix A Uncertainty Propagation
The method of uncertainty propagation for onepoint statistics is first described in Watkinson & Pritchard (2014). Here, we summarize and expand on their work, deriving uncertainty propagation for the kurtosis in addition to variance and skewness.
To recap, for a map with pixel value , mean and pixels, the pth central moment of the map is defined as,
(A1) 
Then, the variance, skewness and kurtosis are standardizations of 2nd, 3rd and 4th central moments as follow.
variance:  (A2)  
skewness:  (A3)  
kurtosis:  (A4) 
Assuming that every pixel in a 21 cm brightness temperature fluctuation map represents an independent brightness temperature fluctuation “pixel” , we can simply substitute and to compute the “true” moments of the brightness temperature fluctuation.
Adding noise with standard deviation to the true signal, each pixel now consists of the true signal plus the noise, , and the noise will bias the moments. An unbiased estimator for the pth moments can be estimated by averaging the moment equations over noise realization. Assuming that the noise is Gaussian and independent on different pixels, the averaged noise terms can be rewritten as functions of standard deviation of Gaussian noise using the Gaussian moment identity,
(A5) 
where the angle bracket designates an average. In addition to Equation A5 and the identities given in the Table A1 in Watkinson & Pritchard (2014), Table 2 in this work provide more identities necessary for the derivation of the kurtosis uncertainty.
The uncertainty for each statistic that will arise from the added noise can be determined by calculating the estimator variance and covariance of the unbiased estimator of the moments,
(A6)  
(A7) 
Since skewness and kurtosis are functions of the 2nd and 3rd, and 2nd and 4th moments, their estimator variance can be calculated by propagating the estimator variance of the moments to respective statistics using Taylor expansion,
(A8) 
Here, is a function of two nonindependent variables and . In other word, and for skewness, and and for kurtosis.
Equations A9 to A14 summarize results from Watkinson & Pritchard (2014). Here, is assumed to be equal to in Equation 7 for all pixels.
(A9)  
(A10)  
(A11)  
(A12)  
(A13)  
(A14) 
We follow the procedure in Watkinson & Pritchard (2014) and are able to confirm their results. In addition, we derive the estimator variance for kurtosis as follows.
First the the 4th moment with added noise is constructed.
(A15) 
Then, we average over the noise.
(A16) 
This implies that the unbiased estimator of the 4th moment is,
(A17) 
Next we derive the estimator variance of the 4th moment. We substitute and to simplify the derivation.
(A18) 
Expanding this expression and moving the noise averaging brackets inside the summation gives,
(A19) 
Using the Gaussian noise identities reduces the expression to,
(A20) 
Doing the summation to perform index conversion via , substituting all terms with the pth moments and with , resubstituting , and cancelling out many terms will yield the estimator variance of the 4th moment,
(A21) 
The estimator covariance between 2nd and 4th moment can be found in a similar manner, resulting in,
(A22) 
The estimator variance of the moments can then be propagated to the kurtosis with Equation A8 to obtain,
(A23) 
References
 Ali et al. (2015) Ali, Z. S., Parsons, A. R., Zheng, H., et al. 2015, ApJ, 809, 61
 Barkana & Loeb (2008) Barkana, R., & Loeb, A. 2008, MNRAS, 384, 1069
 Beardsley et al. (2015) Beardsley, A. P., Morales, M. F., Lidz, A., Malloy, M., & Sutter, P. M. 2015, ApJ, 800, 128
 Beardsley et al. (2012) Beardsley, A. P., Hazelton, B. J., Morales, M. F., et al. 2012, MNRAS, 425, 1781
 Beardsley et al. (2013) —. 2013, MNRASL, 429, L5
 Beardsley et al. (2016) Beardsley, A. P., Hazelton, B. J., Sullivan, I. S., et al. 2016, arXiv, arXiv:1608.06281
 Bowman et al. (2009) Bowman, J. D., Morales, M. F., & Hewitt, J. N. 2009, ApJ, 695, 183
 Bowman et al. (2013) Bowman, J. D., Cairns, I., Kaplan, D. L., et al. 2013, PASA, 30, 31
 DeBoer et al. (2016) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2016, arXiv, arXiv:1606.07473
 Dillon et al. (2015) Dillon, J. S., Neben, A. R., Hewitt, J. N., et al. 2015, Phys. Rev. D, 91, 123011
 Fan et al. (2006) Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, AJ, 132, 117
 Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181
 Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 16
 Gluscevic & Barkana (2010) Gluscevic, V., & Barkana, R. 2010, MNRAS, 408, 2373
 Harker et al. (2009) Harker, G. J. A., Zaroubi, S., Thomas, R. M., et al. 2009, MNRAS, 393, 1449
 Ichikawa et al. (2010) Ichikawa, K., Barkana, R., Iliev, I. T., Mellema, G., & Shapiro, P. R. 2010, MNRAS, 406, 2521
 Iliev et al. (2002) Iliev, I. T., Shapiro, P. R., Ferrara, A., & Martel, H. 2002, ApJ, 572, L123
 Jacobs et al. (2013) Jacobs, D. C., Bowman, J., & Aguirre, J. E. 2013, ApJ, 769, 5
 Jacobs et al. (2015) Jacobs, D. C., Pober, J. C., Parsons, A. R., et al. 2015, ApJ, 801, 51
 Jacobs et al. (2016) Jacobs, D. C., Hazelton, B. J., Trott, C. M., et al. 2016, ApJ, 825, 114
 Lidz et al. (2007) Lidz, A., Zahn, O., McQuinn, M., et al. 2007, ApJ, 659, 865
 Madau et al. (1997) Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429
 Malloy & Lidz (2013) Malloy, M., & Lidz, A. 2013, ApJ, 767, 68
 Mellema et al. (2013) Mellema, G., Koopmans, L. V. E., Abdalla, F. A., et al. 2013, Exp Astron, 36, 235
 Mesinger et al. (2010) Mesinger, A., Furlanetto, S., & Cen, R. 2010, MNRAS, 411, 955
 Parsons et al. (2010) Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, AJ, 139, 1468
 Pillepich et al. (2007) Pillepich, A., Porciani, C., & Matarrese, S. 2007, ApJ, 662, 1
 Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
 Pober et al. (2015) Pober, J. C., Ali, Z. S., Parsons, A. R., et al. 2015, ApJ, 809, 62
 Saiyad Ali et al. (2006) Saiyad Ali, S. K., Bharadwaj, S., & Pandey, S. K. 2006, MNRAS, 366, 213
 Scott & Rees (1990) Scott, D., & Rees, M. J. 1990, MNRAS, 247, 510
 Shimabukuro et al. (2016) Shimabukuro, H., Yoshiura, S., Takahashi, K., Yokoyama, S., & Ichiki, K. 2016, MNRAS, 458, 3003
 Sullivan et al. (2012) Sullivan, I. S., Morales, M. F., Hazelton, B. J., et al. 2012, ApJ, 759, 17
 Sunyaev & Zeldovich (1972) Sunyaev, R. A., & Zeldovich, Y. B. 1972, A&A, 20, 189
 Thyagarajan et al. (2015a) Thyagarajan, N., Jacobs, D. C., Bowman, J. D., et al. 2015a, ApJL, 807, L28
 Thyagarajan et al. (2015b) —. 2015b, ApJ, 804, 14
 Tingay et al. (2013) Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, 7
 Tozzi et al. (2000) Tozzi, P., Madau, P., Meiksin, A., & Rees, M. J. 2000, ApJ, 528, 597
 Trott et al. (2016) Trott, C. M., Pindor, B., Procopio, P., et al. 2016, ApJ, 818, 139
 Watkinson et al. (2015) Watkinson, C. A., Mesinger, A., Pritchard, J. R., & Sobacchi, E. 2015, MNRAS, 449, 3202
 Watkinson & Pritchard (2014) Watkinson, C. A., & Pritchard, J. R. 2014, MNRAS, 443, 3090
 Watkinson & Pritchard (2015) —. 2015, MNRAS, 454, 1416
 Wyithe & Morales (2007) Wyithe, J. S. B., & Morales, M. F. 2007, MNRAS, 379, 1647
 Yoshiura et al. (2015) Yoshiura, S., Shimabukuro, H., Takahashi, K., et al. 2015, MNRAS, 451, 266