MNRAS 506, 4121–4130 (2021) https://doi.org/10.1093/mnras/stab1941 Advance Access publication 2021 July 8 Future radio continuum cosmology clustering surveys Jacobo Asorey 1,2 and David Parkinson 2,3‹ 1Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas (CIEMAT), Av. Complutense, 40, E-28040 Madrid, Spain 2Korea Astronomy and Space Science Institute, Yuseong-gu, Daedeok-daero 776, Daejeon 34055, Korea 3University of Science and Technology, Daejeon 34113, Korea Accepted 2021 July 2. Received 2021 July 2; in original form 2021 February 10 ABSTRACT The use of continuum emission radio galaxies as cosmological tracers of the large-scale structure will soon move into a new phase. Upcoming surveys from the Australian Square Kilometre Array Pathfinder (ASKAP), MeerKAT, and the Square Kilometre Array project (SKA) will survey the entire available sky down to an ∼ 100μJy flux limit, increasing the number of detected extra-galactic radio sources by several orders of magnitude. External data and machine learning algorithms will also enable some low-resolution radial selection (photometric redshift binning) of the sample, increasing the cosmological utility of the sample observed. In this paper, we discuss the flux limit required to detect enough galaxies to decrease the shot-noise term in the error to be 10 per cent of the total. We show how future surveys of this type will be limited by available technology. The confusion generated by the intrinsic sizes of galaxies may have the consequence that surveys of this type eventually reach a hard flux limit of ∼100 nJy, as is predicted by the current modelling of AGN sizes by simulations such as the Tiered Radio Extragalactic Continuum Simulation (T-RECS). Finally, when considering the multitracer approach, where galaxies are split by type to measure some bias ratio, we find that there are not enough AGN present to achieve a reasonable level of shot noise for this kind of measurement. Key words: large-scale structure of Universe – radio continuum: galaxies. 1 IN T RO D U C T I O N The distribution of matter in the Universe, originally seeded as tiny fluctuations during the initial expansion, provide an important cosmological probe of both the physics of the early time, and the subsequent expansion and evolution. The late-time1 matter distribution is often measured through the statistics of luminous matter, and this is most commonly in the form of galaxies. The use of clustering statistics of continuum radio sources as a cosmological probe (Blake & Wall 2002; Jarvis et al. 2015; Chen & Schwarz 2016) has lagged somewhat behind that of optical galaxies. There are two reasons for this: the relative difficulty of obtaining a equivalently large sample (given the available technology), and the comparative lack of spectral features at those wavelengths, making localization in the radial direction much harder. To close this gap, the next generation of radio observatories will need to have a high sensitivity and capable of both precise angular localization and fast coverage of a wide-area. Fortunately, with the development of the SKA pathfinders such as ASKAP (the Australian Square Kilometre Array Pathfinder; Hotan et al. 2021), MWA (Murchison Widefield Array; Tingay et al. 2013), and MeerKAT, this will be possible. The next few years should see a rapid growth in the number of galactic radio sources useful for cosmology as wide-area surveys reach an � E-mail: davidparkinson@kasi.re.kr 1Here ‘late-time’ commonly refers to any survey or probe sampling the Uni- verse after the cosmic microwave background was emitted at recombination. In this regard, even an AGN survey extending up to redshift z ∼ 5 can be considered a late-time cosmology experiment. rms noise of ∼ 10μJy beam−1 (Norris et al. 2011; Hurley-Walker et al. 2016; Santos et al. 2017). However, while this door of opportunity may be opening now, it may not be possible to keep it open into the distant future. In attempting to reach sub μJy flux limits for the next generation of instruments, all-sky radio continuum surveys beyond the current SKA design plan may find the accessible area of the sky to be reduced. In the SKA technical memo by Condon (2009), the author discussed the effect of confusion noise and natural confusion providing a hard bound on the achievable flux limit of the survey, given the radio telescope being used. Even if this very faint limit can be reached for small, pinpoint surveys, extending these to the whole sky will be limited by the necessity of masking bright radio sources, and the presence of diffuse galactic foreground emission increasing the difficulty of detecting faint sources. In Hopkins et al. (2000), the authors made a detailed analysis of what technology would be required to detect radio continuum galaxies in the sub-μJy regime. In this paper, we build on those ideas and subsequent work to explore the optimal manner of sampling both wide areas and reaching deep flux limits that will be required for the cosmology use of radio-continuum surveys. Much of this work is based on the simulation of the radio sky from the Tiered Radio Extragalactic Continuum Simulation (T-RECS; Bonaldi et al. 2019). Though this is the current ‘state of the art’ simulation, it is still based on observational data. This means that the extrap- olations that we make to as yet unexplored regions of parameter space (e.g. flux densities � 1μJy) may not be reflected in future observations. In Section 2, we show how the flux limit for radio contin- uum surveys influences the measurement of the angular clustering C© 2021 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 http://orcid.org/0000-0002-6211-499X http://orcid.org/0000-0002-7464-2351 mailto:davidparkinson@kasi.re.kr 4122 J. Asorey and D. Parkinson statistics, and how the separation of the sample into redshift bins and by population can reduce the precision of these clustering measurements. In Section 3, we show how achieving the maximal area that minimizes sample variance can be technically challenging due to the sources of noise and confusion present in the sky. In Section 4, we discuss the designs for current and future wide-area continuum surveys, and the instrumentation and data processing requirements to reach the proposed cosmology limits. 2 C L U STERIN G SURVEYS The design of a future radio continuum survey will depend on previous surveys and the science goals that we have prioritized. In the area of cosmology the science topics of the mysterious acceleration of the Universe (through some possible dark energy or modified gravity), and the initial conditions of the early Universe (generated by cosmological inflation or some other mechanism) will both return optimal results from all-sky extra-galactic surveys. Previous work has identified that testing inflation through non-Gaussianity will require an all-sky survey, with a number density high enough for shot-noise effects to be subdominant, which can be subdivided into at least five redshift bins (Raccanelli et al. 2012, 2015; Bernal et al. 2019). We set this as the standard required to measure the non-Gaussian component on the initial curvature fluctuation to the same or better accuracy than Planck. In this section, we discuss how the measured clustering power spectrum and its errors change as the survey achieves deeper flux limits, and the sample is subdivided into redshift bins. 2.1 Clustering statistics The two-point distribution of radio galaxy positions in angular space can be decomposed into a series of spherical harmonics. For a given direction on the sky θ , if the continuous density field in that direction σ (θ) is Gaussian and randomly distributed, then it can be decomposed into its multiple moment using spherical harmonics Y�m, such that the amplitudes a�m are given by a�m = ∫ dθY ∗ �mσ (θ) . (1) The amplitudes of the spherical harmonic contributions for each of the multipole moments � and m can be averaged in one of the two spherical directions (assuming istoropy) giving the angular correlation power spectrum C�, such that C� = 〈a�ma∗ �m〉 . (2) Assuming that the galaxies trace some underlying matter distribution, which can be represented as a power spectrum of fluctuations, then the prediction for the angular power spectrum from theory is given by C� = 4π ∫ dk k �2(k)[Wg � (k)]2 , (3) where k is the wavenumber, �2(k) is the logarithmic matter power spectrum, and W g � (k) is the radio galaxy window function. This form assumes a single redshift bin with the galaxies distributed across this bin, giving the window function as (e.g. Giannantonio et al. 2008; Raccanelli et al. 2008) W�(k) = ∫ dN (χ ) dχ b(z)D(z)j�[kχ ]dχ . (4) Here dN/dχ is distribution of sources per steradian with comoving radial distance χ (z) within the redshift bin (brighter than some survey magnitude or flux limit), b(χ ) is the bias factor relating tracer overdensity to matter overdensity, D(χ ) is the growth factor of density perturbations, and j�(x) is the spherical Bessel function. The power spectrum of density fluctuations �(k) and the cosmo- logical distances χ (z) are sensitive to the values of the cosmological parameters and the cosmological model. With good knowledge of the bias and number distribution of radio galaxies, accurate measurements of the radio galaxy angular power spectrum can be used to constrain the cosmological parameters, and test the standard concordance cosmology �CDM. But in order to do so, we will need very accurate measurement of this power spectrum. We assume that the field describing the inhomogeneities is Gaussian-distributed in the amplitudes, and so we can describe the covariance of the observed C� as (Asorey et al. 2012) Cov(C�) = 2(C� + 1/n̄)2 N (�) , (5) where N(�) = (2� + 1)fsky is the number of modes sampled for a given �, for fsky as the fraction of the sky being surveyed, and n̄ is the angular number density of sources (in units of number count per steradian) given by n̄ = Ngal,bin/��. Therefore, the variance on an angular power spectrum measure- ment, in the completely Gaussian case, is given by σ 2 C� = 2 fsky(2� + 1) [ C� + 1 n̄ ]2 . (6) The two sources of error can therefore be broken down into sample variance, which is limited by the number of independent modes on the sky that the power can be measured in, and shot noise, which is limited by the number density of sources. Finally, the signal in the angular power spectrum can also be increased by measuring tracers with a greater galaxy bias, which will therefore have a larger angular clustering amplitude. For a target sample of radio continuum galaxies, detected through their synchrotron emission, the number of sources will increase as the flux-limit decreases, and the redshift distribution will change. As future radio surveys reach a greater depth (smaller flux limit), equation (6) implies that the number density of sources will reach a sufficient number density such that the shot noise will be subominant to the sample variance error. Using a model of the total number and redshift distribution of sources, we can compute the corresponding flux limit that gives a fixed fraction of the shot noise contribution to the total errors. The predicted normalized number distribution of galaxies, multi- plied by the linear bias of that tracer, is the kernel that we integrate over to find the window function in equation (4), and is shown in Fig. 1. In this figure, we see that the overall amplitude changes as the number of less-biased SFGs come to dominate over the more highly biased AGN. But the localization changes as well, as we detect more high-redshift SFGs at lower fluxes. To extract cosmological-relevant data from the angular power spectrum, we need to be able to both accurately measure and accurately model it, which gives a limit to range of multipoles we can utilize, as modelling the non-linear power-spectrum accurately becomes more difficult as k increases. Following the approach in Bernal et al. (2019), we define a minimum and maximum multipole number. For the largest scales, �min is limited by the fraction of sky surveyed, �min = π /(2fsky). We assume a 75 per cent total sky coverage, though address the question of sky fraction as a function of flux limit in Section 3. For the small-scale limit, we need to consider the distance, χ (z), to the particular redshift bin in consideration. We MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 Future continuum cosmology clustering surveys 4123 Figure 1. The evolution of the product of the bias and the normalized redshift distribution of total radio continuum populations, AGNs and SFGs for different survey flux cuts. Here the number distribution is drawn from the T-RECS simulation (Bonaldi et al. 2019), and the bias is estimated using the Colossus code (Diemer 2018) with the relationship between halo mass and type the same as adopted by the S simulations (Wilman et al. 2008). This quantity is averaged over when computing the theoretical value of the angular power spectrum, C�, and the change in amplitude and localization of n(z)b(z) relates directly to the amplitude and shape of the C� spectrum. define the relevant scale r� as r� = ∫ χ dN (χ ) dχ dχ, (7) and then we define �max = kmaxr�, where kmax = 0.1 Mpc−1 in order to consider only the linear regimes. In order to realistically simulate the predicted C� and shot noise, we need to model the galaxy bias b(z) and the redshift distribution N(z) of the samples we will observe, for a given flux cut. We use a combination of information from two sets of simulations. We define the redshift distributions of a given population with the T-RECS simulation (Bonaldi et al. 2019). However, for the galaxy bias we use the Colossus (Diemer 2018) suite to estimate the bias from the halo model, using the Tinker model (Tinker et al. 2008). The values for the halo masses assigned to each type of radio galaxy taken by the S simulations (Wilman et al. 2008), which follows the approach of Raccanelli et al. (2012) and Bernal et al. (2019), associating the names in that work with the names given by T-RECS (Bonaldi et al. 2019). In Table 1, we show the different names given to different classifications of radio galaxies for different works. As the relevant quantity is the combination of the bias b(z)n(z), we need to compute a combined bias b(z) for the ensemble. We combine the different populations in the following way: b(z) = ∑type α bα(z)Nα(z) Nall(z) , (8) where α = {AGNRQQ, AGNFR I, AGNFR II, SFGSFG, SFGSB} following the types used in the S3 simulation. We derive this form of the total bias of the ensemble in Appendix A. As seen in equation (4), the combination of the bias and the normalized redshift distribution b(z)nnorm(z) is the main component of the angular power spectrum window function. We show in Fig. 1 the kernel of this window function for different flux cuts. For bright flux limits, where the redshift distribution is dominated by AGN, there is a inherent noisiness to the redshift distribution of this species, caused by the small number of galaxies and the finite size of the simulation. However, this distribution is firstly transformed through equation (4) to get the window function, and the square of the window function is convolved with the three-dimensional matter power spectrum to give the angular power spectrum through equation (3), and both operations will smooth over any initial oscillations in n(z). We therefore expect this noisiness to have almost no effect on the predicted power spectrum at bright flux limits, and no consequence for our conclusions. In Fig. 2, we show how the angular power spectra changes with flux limit Scut. We find that as the flux limit decreases and more sources are detected, the power increases until a maximum at around Scut ∼ 10−5Jy, owing to the increased radial localization of the tracers, as the number of low-redshift AGN and SFG rises, as well as a higher overall bias. For flux limits smaller than this, the amplitude decreases again slightly, as the number of higher redshift tracers increases, and the power becomes more radially smeared. In Fig. 2, we also show how the error on the angular power spectrum changes with flux limit, for the sample variance alone, and the sample variance and shot noise combined. We find that the shot noise contribution to the error becomes subominant to the sample variance at a flux limit of roughly 30μJy, and the two lines converge. Note that since the shot noise component is inversely proportional to number density, rather than raw number, this change is independent of the total area surveyed. If all the galaxies that can be detected are used in a single ‘redshift bin’, then integrating down to a flux limit fainter than 30μJy will increase the number of galaxies, but will not decrease the fractional error on the measured C� power spectrum. 2.2 Redshift binning Subdividing the sample of galaxies by redshift will allow the angular power spectrum to be measured in more than one redshift bin, and so increase the cosmological utility (e.g. Camera et al. (2012). However, as continuum sources are normally detected at the ≈1-GHz frequency through their synchrotron emission, the featureless power-law nature of their flux spectrum makes for a greater challenge than for optical photometric redshift binning. As such, future surveys may need to make use of data from other surveys, to either cross-identify with an optical/NIR source with redshift information, or make use of statistical clustering redshifts (Kovetz, Raccanelli & Rahman 2017), MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 4124 J. Asorey and D. Parkinson Table 1. A comparison given to the names of different types of radio galaxies, based on the way that they are classified. SKADS name Spectral classification T-RECS name (Wilman et al. 2008) (Hale et al. 2018) (Bonaldi et al. 2019) Star-forming galaxy (SFG) SFG SFG Starburst (SB) SB RQQ (radio-quiet quasar) Moderate- to low-luminosity AGN (MLAGN) BL Lac FRI (Fanaroff–Riley type I) Flat-spectrum radio quasar (FSRQ) FRII (Fanaroff–Riley type II) High- to moderate-luminosity AGN (HLAGN) Steep-spectrum AGN (SS-AGN) Figure 2. Predicted angular power spectra C� (top panel) and error on the angular power spectrum (bottom panel) of the clustering of radio continuum galaxies of a single redshift bin, for different values of the flux limit scut. We see the increase and decrease in overall power amplitude as the flux limit changes, as well as a relatively small change in shape, owing to the number distribution and bias evolution of the tracers. For the error, the sample variance contribution in orange, and the combined sample variance and shot noise effect combined in blue. The sample variance error increases on the right hand side as the amplitude also increases with decreasing scut,. again using known redshifts of sources over the same area of sky. In either case though, much more information is needed to subivide the sample by redshift will be available at low redshift than high. In this paper, rather than assuming a particular set of redshift bins, we assume that the subdivision will be roughly equal in log (1 + z), with a maximum of z = 5. In Fig. 3, we show the redshift ranges spanned by the different bins for each redshift binning case considered in this paper. As the galaxies being used in the sample are localized into smaller bins, the minimum scale that can accurately be modelled also changes, and will be different for each bin. We compute the maximum multiple number �max using equation (7) for each bin. We show the values of �max for the extreme case of eight redshift bins in Table 2. When the sample of galaxies can be subdivided, either by redshift or population, or both, the required number of galaxies increases. Figure 3. The ranges in redshift of the different redshift binning configura- tion considered in this paper. We assumed that there would be a greater amount of multispectrum, morphological, and clustering information available at lower redshifts than higher, which would allow more precise subivisions in redshift. In this paper we assume this would manifest as for equal binning in log (1 + z), as shown here. Throughout this paper, we assume top-hat redshift bins. Table 2. Scales used in the analysis for the configuration with eight redshift bins that correspond to a comoving scale of k = 0.1 Mpc−1. Redshift bin Mean distance (Mpc) �max 0–0.25 504 50 0.25–0.57 1546 155 0.57–0.96 2649 265 0.96–1.45 3756 376 1.45–2.06 4822 482 2.06–2.83 5815 581 2.83–3.8 6725 673 3.8–5 7548 755 We assume a number of redshift bins, with the size of each scaled by the logarithm of 1 + z and spaced to cover the range 0 < z < 5, as shown in Fig. 3. In the left-hand panel of Fig. 4, we show how the flux limit for the total population required to achieve 10 per cent shot noise increases with number of bins. We find that there is a ‘plateau’ between nbins = 5 and 8 for a flux limit of 0.5μJy, such that subdividing the sample into more bins does not significantly impact on the shot noise contribution. 2.3 Multitracer sampling One way to remove or reduce sample variance consists of using two differently biased tracers of the same matter field (McDonald & Seljak 2009; Seljak 2009; Asorey, Crocce & Gaztañaga 2014), and measuring the ratio of the density fluctuation of the two different samples, over the same volume. If the two tracers are following the same underlying matter field (i.e. there is no stochasticity), it follows MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 Future continuum cosmology clustering surveys 4125 Figure 4. Fraction of the signal noise that is shot noise for different redshift configurations (described by the number of redshift bins nbins) and flux limits, for all galaxies (left-hand panel) and AGN only (right-hand panel). The fraction is shown by the colour bar, with dark blue being mostly shot noise in the error budget, and bright yellow being mostly sample variance. The red line shows the flux limit that is needed for each redshift case for the shot noise to be less than 10 per cent of the total signal noise. There is no red line on the right-hand panel of the figure, because in no cases do we achieve a less than 10 per cent shot-noise contribution to the error budget when using AGN only. This assumes an all-sky survey, fsky = 1. that this ratio will be independent of the underlying fluctuation. Then the error on the power spectrum measurement of this ratio will be independent of the modes of the underlying fluctuation, effectively cancelling the sample variance, allowing for better than cosmic variance uncertainty on large-scale physical effects that are not dependent on the underlying power spectra, such as the non-linear bias generated by primordial non-Gaussianity (e.g. Ferramacho et al. (2014), Bernal et al. (2019), Gomes et al. (2019)). For this technique to be effective, a large number density of tracers with different biases must be surveyed over the entire survey volume. In the radio continuum it is possible to separate the sample into AGN and star-forming galaxies, though doing so will reduce the number density, in comparison to the combined cohort. To achieve the same shot-noise accuracy as we defined for Section 2.2, it will be necessary to increase the numbers of these individual populations by pushing down to even smaller flux limits. In the right-hand panel of Fig. 4, we show the flux limit required in order to avoid shot noise saturation of the error for a given set for redshift bins, when only considering AGN galaxies. We only consider AGN here because their number density is much lower than that of SFGs below 10 mJy, as shown in Fig. 1, and so this population will be the limiting factor. It is shown that we can only achieve the same shot noise error by reaching a flux limit of scut = 0.1 μJy, and in this case, only for a single continuous redshift bin. This is similar to the result found in Bernal et al. (2019), where the best constraint on the non-Gaussianity parameter fNL from the EMU survey would come from the analysis with an SFG and AGN multitracer approach but only a single redshift bin. Here we state more strongly that no future continuum surveys will be able to subdivide in both redshift and tracer and see any improvement in cosmological constraints, even those achieving a fainter flux limit than EMU. 3 O N-SKY RADIO SURVEY LIMITS If the design goal of large-area angular clustering surveys is to minimize the shot-noise contribution to the error that comes from a low number density of tracers, it must also be to minimize the sample variance by maximizing the sky area covered. However, even if there is enough observing time to integrate down to the required flux limit uniformly over the entire available sky, the resulting sample of tracers may not have the same flux limit and sky fraction as was planned. This is because there are on-sky limitations to the effectiveness of the survey that are generated by bright sources, faint diffuse foreground emission, and confusion between tracers being resolved in the same beam. In this section, we discuss each of them in turn. 3.1 Foreground maps The signal in a given beam would be the combination of the source temperature and the system (instrument) noise: TA = Tsys + Tsource (9) where the system temperature is a combination of the background, which is mostly given by foreground emission (Tsky), the atmospheric emission (Tatm), ground signal scattered from the feed and support structure (Tscat), and electronics noise from the receiver (Trecv), i.e. Tsys = Tsky + Tatm + Tscat + Trecv . (10) The expected flux limit of a radio continuum survey is derived from the root mean square noise on the flux intensity map. For a given integration on the sky of time τ , the rms noise in the flux measurement (Srms) from the temperature of the system is given by the radiometer equation, Srms = 2kBTsys AeffNdish √ 2Bτ , (11) where kB is the Boltzmann constant, B is the bandwidth, Aeff is the effective collecting area, and Ndish is the number of dishes. Note that this noise is assumed to be Gaussian and generated by the internal temperature fluctuations of the radio telescope and some homogeneous external radio field (background sky temperature). Sources are thus identified as being fluctuations of an amplitude significantly larger than this flux rms, and normally associated with a cut at (for example) 10σ . While random fluctuations of this size MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 4126 J. Asorey and D. Parkinson Figure 5. Fraction of the sky that will be inaccessible as a function of integration time (τ ) for different flux limits and different system temperatures, assuming a 10σ detection threshold. We assume a background sky modelled by the GSM (de Oliveira-Costa et al. 2008; Zheng et al. 2017), a frequency range of 800 MHz to 1.4 GHz, a brightness efficiency of ηB = 1.0, and a beam size of 60-arcsec FWHM. can happen, the probability of such a fluctuation (10−23) in a sample of 1010–1012 pixels is small enough to avoid consideration. The diffuse sky brightness can overwhelm the emission from a localized extragalactic source, reducing the signal-to-noise ratio (S/N) of the detection. The closer to the galactic equator that observations are made, the brighter the sky will be, and longer integrations will be needed to achieve a uniform deepness. Reversing this argument, we can estimate the amount of sky that will be inaccessible for a given flux limit and integration time. We use the Global Sky Model (GSM; de Oliveira-Costa et al. 2008; Zheng et al. 2017) to model the sky temperature across the sky, and the radiometer equation (equation 11) to model the S/N for detection. In Fig. 5, we show how the sky fraction varies with exposure time and flux limit, for two possible system temperatures. To reach the required depth across the sky, such that the number of galaxies gives < 10 per cent shot noise for at least two redshift bins, we would need integrations of at least 2 h for a 50-K detector. Greater depths can be reach with longer integrations or more sensitive equipment. This argument is applied to the entire sky, not to any particular observing patch. There are, of course quieter regions of the sky that can be observed that have relatively low sky temperature (e.g. the Lockman Hole), and so can be used for deeper observations with shorter integration time. In cosmology, since the sample variance will decrease as the survey area increases, the aim will always be to maximize the areal coverage. For a fixed total observing time, there will therefore be a trade-off between longer integrations and a wider area. 3.2 Confusion limits Even if the integration time is increased indefinitely, there are natural limits to the depth at which radio continuum surveys can operate. The homogeneous radio sky approximation breaks down when considering the existence of faint extra-galactic radio sources, which generate ‘confusion noise’ (Condon et al. 2012). The number of sources above some flux limit is given by N (> S0) ∫ ∞ S0 n(S)dS . (12) The solid angle for a particular beam is �beam = πθ2 4 ln(2) , (13) Figure 6. Flux limit of different surveys generated by the confusion of multiple detectable sources in the same beam, as a function of angular resolution. Here we assume a frequency of 1 GHz. We see that the EMU survey, which has a confusion limit of ∼ 30μJy (5σ ) is operating close to the design limit given by the angular resolution, as it lacks the very long 150-km baselines of SKA-MID1. where θ is the beam FWHM. If the angular resolution of the telescope is too small, the number of sources per beam will be too large, and the images of these faint sources will overlap, causing the images generated to merge. So even if they are bright enough to be detected, the generated image will be too noisy for the sources to be distinguished. In this case the receiver may be sensitive enough to detect the flux from the fainter sources, but the source detection is limited by the ability to separate them on the sky. This limit is given by β = [N (> S0)�beam]−1 , (14) and is described in terms of some ‘number of sigma’ cut, where the sources are distinct enough on the sky to be non-overlapping. So for β = 25, this would a 5σ confusion limit for a flux cut of S0, assuming a slope of 2 for the differential source counts. We show how this translates into a flux limit for a given angular resolution, assuming the T-RECS simulated catalogue, in Fig. 6. At around 1 GHz, we see that the EMU survey, which has a confusion limit of ∼ 30μJy (5σ ) is operating close to the design limit given by the angular resolution of ASKAP. SKA-MID1 does a lot better here, as it has access to some very long baselines of 150 km. However, using these very long baselines and so generating very high-resolution images will increase the processing power required to run such a survey (though the discussion of processing power is beyond the current scope of this paper). An increase in the number of long-baselines that is needed to achieve some high-resolution imaging, and so lower the confusion limit, will also increase the computing power capacity needed to process the data. The ASKAP correlator is built to operate at 340 TFLOPS, and the raw data rate for ASKAP is approximately 100 Tbit s−1 (Hotan et al. 2021). However, the data rate would be expected to scale as (Bmax/D)4, where Bmax is the length of the maximum baseline, and D is the diameter of the dish. If an observatory was constructed with similar size dishes to ASKAP but the 150km baselines planed for SKA Mid, the increase in data rate would be by a factor of (150/6)4 ∼ 108, or approximately 100 exaFLOPS. This would be a more than significant upgrade in comparison to the current computational facilities. This confusion noise is generated even in the ‘scale-free’ case, where the clustering of these sources is not considered. There is also a ‘natural confusion’ limit, where the inherent angular size of the sources themselves is large enough for them to overlap on the sky, MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 Future continuum cosmology clustering surveys 4127 Figure 7. The median size of radio continuum galaxies, as a function of flux limit, for AGN and SFG, as taken from the T-RECS deep simulated catalogue (Bonaldi et al. 2019). We also add the sizes of galaxies (or the upper limits on their sizes) as measured by the Cotton et al. (2018) and Vernstrom et al. (2016) surveys at 3 GHz. We compare these sizes to the 5σ natural confusion limit and the best-fitting formula for the median source size taken from Windhorst (2003). Here we assume a frequency of 1.4 GHz. We find that at a flux limit of about 100 nJy the presence of a large number of ‘monster’ Steep-Spectrum AGN (SS-AGN) would place a limit on the possible depth of large-area surveys. Figure 8. The number of radio continuum galaxies as a function of flux limit, for AGN and SFG, as taken from the T-RECS deep simulated catalogue (Bonaldi et al. 2019). Here we assume a frequency of 1.4GHz. and increased angular resolution will not save you from noise-limited images. In this case, the limit is given by βn.c. = [N (> S0)�source]−1 . (15) If the source solid angle is given by �source ∼ πθ2 s /[4 ln(2)], where θ s is the median angular source size, then we can compute the value of θ s that corresponds to 5σ natural confusion limit (again assuming a slope of 2 for the differential source counts). If the average source radius is above this size, then the sources will be overlapping and confused. In Fig. 7, we show this source size limit as a function of flux limit, and compare with the source sizes present in the T-RECS catalogue (Bonaldi et al. 2019), for SFG and the different classifications of AGN present in the T-RECS simulated catalogue. We also compare these sizes to the best-fitting curve for the median sizes as a function of flux density taken from Windhorst (2003), which is given by θs = 2 ( S 1 mJy )0.3 , (16) where S is the flux and θ s is median source angular size, measured in arcseconds. We also add the sizes of galaxies (or upper limit on sizes) as measured by the Cotton et al. (2018) and Vernstrom et al. (2016) surveys at 3GHz, showing that sizes have been measured down to roughly 10μJy, but no deeper. The Windhost relation relation is a reasonable fit to the SFG sizes from T-RECS for almost all flux densities. However, this is not the case for AGN. The T-RECS simulation assumes a constant size-flux relation of the different AGN species. It therefore contains a number of ‘monster’ AGN (the largest species, Steep-Spectrum AGN), which have large sizes but low flux densities. Though these would be subdominant in number, they would be large enough to confuse any surveys operating at the nJy level. If this population is present in reality, this would be a hard limit beyond which it would not be possible to increase the number density by further integration. 3.3 Masking bright sources For a given faint flux limit, there is a corresponding maximum flux limit beyond which bright sources will cause problems for reaching the required faint limit. The relationship between the bright and faint limits is given by the dynamic range DR, such that DR = 10 log10 ( Smax Smin ) , (17) where the dynamic range is given in decibels (dB). MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 4128 J. Asorey and D. Parkinson Figure 9. Fraction of the sky that will be inaccessible as a function of dynamic range for size different flux limits (Sig), using the distribution of radio galaxy fluxes taken from the T-RECS simulated catalogue (Bonaldi et al. 2019). As the dynamic range decreases, the number of bright sources that needs to be masked increases, until the masks dominate the sky area. It is clear that to achieve an all-sky survey of radio continuum galaxies with a a sub-μJy flux density, a dynamic range of at least 70 dB (107 in flux ratio) would be required. Any sources in the field with a flux higher than this limit will generate imaging artefacts in the map, as the sidelobes of bright sources will be incompletely CLEANed (Högbom 1974; Chen & Schwarz 2016), and so obscure nearby fainter sources. Therefore these will need to masked before the angular correlation function can be accurately measured. In the NVSS catalogue (Condon et al. 1998) sources closer than about 35 arcmin to a strong source of flux density S and weaker than 10−3S (given the local dynamic range assumed to be 30) were excluded from the analysis (Blake & Wall 2002; Chen & Schwarz 2016). So as the dynamic range decreases the number of sources above the critical flux maximum Smax will increase, and, because of the need to mask more of these bright sources, a greater area of the sky will become unusable. Different analysis have made different assumptions about the amount of area that needs to be cut. In Chen & Schwarz (2016) they cut a circular region of radius 0.◦6 around the bright sources, which is an area of an area of roughly 1 deg2. But in Bengaly, Maartens & Santos (2018) this is widened to a 1◦ radius circle, which is an area of roughly 3 deg2. Even increasing the number or length of the baselines may not help reduce the size of the region to be removed, since, as shown in Fig. 7 and equation (16), the brighter sources that generate these artefacts also tend to be large. Fig. 9 shows the relationship between sky coverage and dynamic range for a number of different flux limits. Here we assume that we must remove about 1 deg2 around each bright source above the flux maximum Smax, as theCLEAN algorithm has not been significantly improved on for this action. 4 N E X T-G E N E R AT I O N C O N T I N U U M S U RV E Y S Table 3 summarizes the current and future radio continuum clustering surveys. These are the NRAO VLA Sky Survey (NVSS; Condon et al. 1998), the LOFAR Two-metre Sky Survey (LoTSS; Shimwell et al. 2017, 2019), the Evolutionary Map of the Universe survey (EMU; Norris et al. 2011), the MeerKAT Large Area Synoptic Survey (MeerKLASS; Santos et al. 2017), and the Square Kilometer Array Mid medium and wide cosmology surveys (SKA; Jarvis et al. 2015; Square Kilometre Array Cosmology Science Working Group et al. 2020). Since the current construction plan for SKA does not contain a survey telescope, conducting a large-area continuum survey with SKA-MID will be more time-consuming, due to the smaller field of view. Therefore, the SKA Cosmology Red book (Square Kilometre Array Cosmology Science Working Group et al. 2020) baseline continuum survey is designed to cover only 5000 deg2, and the sample will be also used for radio weak lensing shear analysis. It may also be possible to use the SKA-MID Band 1 wide cosmology survey to detect continuum sources, though this wide survey is primarily designed as an H I Intensity Mapping survey. As noted in the Red book, further study is needed on the feasibility of conducting the two (continuum and H I intensity mapping) commensally. A true all-sky radio continuum survey with SKA, which can achieve a significant improvement on the EMU depth, may therefore need to wait until some future expansion of the SKA. In Table 3, we list such a successor survey and observatory, with a proposed flux limit of 1 μJy. This future survey would require an angular resolution of at least 1 arcsec, an integration time of at least two hours per field of view (assuming a 50-K system temperature), and a dynamic range of at least 60 dB. This would allow us to achieve sample- variance limited measurements of the clustering power spectrum in four redshift bins, assuming a single population (from Fig. 4), and reasonably small shot noise contribution for a multitracer analysis for a single bin (from the right-hand panel in Fig. 4). Beyond this, it would be possible to extend down further to a 100-nJy limit, which would still only need an angular resolution of 1 arcsec (to avoid confusion), a dynamic range of 70 dB, and a significant improvement in system temperature, as was previously described in Condon (2009). But, given the results from from Fig. 7, any fainter than this, and the sources would be confused, overlapping due to their natural size. This natural confusion limit, which is generated by a population of faint but large ‘monster’ AGN (if they are present), would not be resolvable by any improvement in angular resolution by increasing the length and number of long baselines. The next generation of radio continuum surveys, particularly very deep radio continuum observations from the SKA, would reveal or falsify such a population. 5 C O N C L U S I O N S In this paper we have studied the limitations of current and future radio continuum surveys, from the perspective of measuring the angular clustering statistics of radio galaxies to determine the large- scale structure of the Universe. We considered how sources of noise (both sky, and instrumental) can limit the possible flux limit, and so the number of available sky sources, and the region of sky available to survey. This limit on the number of sources then has consequences for how precisely the angular clustering statistics can be measured, especially in the cases where the sample is subivided by redshift and species. We summarize our conclusions as follows: (i) The shot noise term in the clustering error measurement can be reduced by going deeper, which adds more galaxies to the sample. Surveys such as EMU & MeerKLASS should reach the ∼ 100μJy flux limit needed for the shot noise to be subominant (< 10 per cent of the total) in a single redshift bin. (ii) For multiple redshift bins, the sample is subdivided, and so for a given flux the number density will be lower in each individual bin. The survey will need to go even deeper in flux to reduce the shot noise, and to achieve a < 10 per cent shot noise contribution in at least six redshift bins, a flux limit of roughly 500 nJy, corresponding to a 50-nJy rms, would be needed. (iii) For a survey such that only reaches an ∼ 100μJy limit, the multitracer approach will be not very effective, as the shot noise MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 Future continuum cosmology clustering surveys 4129 Table 3. Survey parameters of the current and future planned radio continuum surveys, based on exiting data sets and future plans. Survey Observatory Frequency (MHz) Flux limit fsky Optimal nbins NVSS VLA 1350–1450 2.5 mJy 0.65 – LoTSS LoFAR 120–168 500 μJy 0.65 – EMU ASKAP 800–1400 50 μJy 0.65 1 MeerKLASS MeerKAT 900–1670 25 μJy 0.1 1 SKA-Mid medium cosmology SKA-MID B2 950–1750 8 μJy 0.12 2 SKA-Mid wide cosmology SKA-MID B1 350–1050 22.8 μJy 0.45 1 SKA all-sky cosmology SKA Survey 800–1400 1 μJy 0.65 4 Notes. The appropriate references are found in the text. Here we define the optimal number of redshift bins based on those that give less than 10 PER CENT shot noise contribution to the error budget. for individual populations is too large (as was previously shown in Bernal et al. (2019)). (iv) Considering the less numerous AGN population as the limit- ing factor for a multitracer analysis, the shot noise is never smaller than 10 per cent of the sample variance, even for a single redshift bin. There are simply not enough AGN available in the sky (given current modelling of the population) to achieve the required accuracy. (v) To reach the required flux limits for a usable cosmological sample, instrument sensitivity (system temperature) is now less important than angular resolution, needed to minimize the confusion noise, and dynamic range, necessary to maximize the sky area. (vi) The requirements for a future large-scale survey with a flux limit of around 100nJy (needed for shot noise < 10 per cent in six redshift bins), the telescope needs to operate at 70dB with a angular resolution of 1”. So while 50K detectors, such as the Phased Array Feeds of ASKAP, would be sensitive enough, more detectors with longer baselines and a greater data processing capacity would be needed for such a future survey. Such would be the requirement for a future expansion of SKA with a survey telescope. (vii) The T-RECS (Bonaldi et al. 2019) simulation predicts a population of ‘monster’ AGN with median size greater than 1 arcsec hiding at depth of 100 nJy. It will not be possible to survey to much greater depths than this over a wider area because of these, as the natural confusion of overlapping sources will make surveys deeper than that impossible. If these exist, they provide a hard limit beyond which it would not be possible to increase the number density by longer observations. The next generation of large-area radio surveys, such as a future expansion of SKA with a survey telescope, would have the capabili- ties needed to measure the angular correlations of galaxies in multiple redshift bins with a small shot noise contribution to the error. Such measurements will provide a very useful data set to learn more about cosmology and fundamental physics. But for the generation of radio telescopes that come after SKA, and plan to measure the clustering of even fainter radio continuum galaxies over a large area of the sky, there will be challenges, and both in terms of instrumental limitations and the nature of the radio galaxy population. Such a survey may not produce any significant improvement in what will by then have already been measured. AC K N OW L E D G E M E N T S We would like to thank Stefano Camera, Ian Harrison, and Jeffrey Hodgson for helpful discussions when preparing this paper. We would like to thank the EMU cosmology working group for their continued support, including Ray Norris, Glen Rees, Song Chen, and Faisal Rahman. We also thank the referee for helpful comments that contributed to this paper. JA has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 776247 EWC. DP is supported by the project (‘Understanding Dark Universe Using Large Scale Structure of the Universe’), funded by the Ministry of Science. This research made use of ASTROPY,2 a community-developed core PYTHON package for astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), the HEALPIX and HEALPY package (Górski et al. 2005; Zonca et al. 2019), the NUMPY package (Oliphant 2006), the SCIPY package (Virtanen et al. 2019) and MATPLOTLIB package (Hunter 2007). DATA AVAILABILITY STATEMENT The data underlying this paper were accessed from the Tiered Radio Extragalactic Continuum Simulation (T-RECS), accessed from http://cdsarc.u-strasbg.fr/viz-bin/qcat?VII/282. The derived data generated in this research will be shared on reasonable request to the corresponding author. REFERENCES Asorey J., Crocce M., Gaztañaga E., Lewis A., 2012, MNRAS, 427, 1891 Asorey J., Crocce M., Gaztañaga E., 2014, MNRAS, 445, 2825 Astropy Collaboration et al., 2013, A&A, 558, A33 Bengaly C. A. P., Maartens R., Santos M. G., 2018, J. Cosmol. Astropart. Phys., 2018, 031 Bernal J. L., Raccanelli A., Kovetz E. D., Parkinson D., Norris R. P., Danforth G., Schmitt C., 2019, J. Cosmol. Astropart. Phys., 2019, 030 Blake C., Wall J., 2002, MNRAS, 329, L37 Bonaldi A., Bonato M., Galluzzi V., Harrison I., Massardi M., Kay S., De Zotti G., Brown M. L., 2019, MNRAS, 482, 2 Camera S., Santos M. G., Bacon D. J., Jarvis M. J., McAlpine K., Norris R. P., Raccanelli A., Röttgering H., 2012, MNRAS, 427, 2079 Chen S., Schwarz D. J., 2016, A&A, 591, A135 Condon J., 2009, SKA Memo 114. Available at: https://www.skatelescope.o rg/uploaded/12336 114 Memo Condon.pdf Condon J., Cotton W., Greisen E., Yin Q., Perley R., Taylor G., Broderick J., 1998, AJ, 115, 1693 Condon J. J. et al., 2012, ApJ, 758, 23 Cotton W. D. et al., 2018, ApJ, 856, 67 de Oliveira-Costa A., Tegmark M., Gaensler B. M., Jonas J., Landecker T. L., Reich P., 2008, MNRAS, 388, 247 Diemer B., 2018, ApJS, 239, 35 Ferramacho L. D., Santos M. G., Jarvis M. J., Camera S., 2014, MNRAS, 442, 2511 2http://www.astropy.org. MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 http://cdsarc.u-strasbg.fr/viz-bin/qcat?VII/282 http://dx.doi.org/10.1111/j.1365-2966.2012.21972.x http://dx.doi.org/10.1093/mnras/stu1955 http://dx.doi.org/10.1051/0004-6361/201322068 http://dx.doi.org/10.1088/1475-7516/2018/04/031 http://dx.doi.org/10.1088/1475-7516/2019/02/030 http://dx.doi.org/10.1046/j.1365-8711.2002.05163.x http://dx.doi.org/10.1093/mnras/sty2603 http://dx.doi.org/10.1111/j.1365-2966.2012.22073.x http://dx.doi.org/10.1051/0004-6361/201526956 https://www.skatelescope.org/uploaded/12336_114_Memo_Condon.pdf http://dx.doi.org/10.1086/300337 http://dx.doi.org/10.1088/0004-637X/758/1/23 http://dx.doi.org/10.3847/1538-4357/aaaec4 http://dx.doi.org/10.1111/j.1365-2966.2008.13376.x http://dx.doi.org/10.3847/1538-4365/aaee8c http://dx.doi.org/10.1093/mnras/stu1015 http://www.astropy.org 4130 J. Asorey and D. Parkinson Giannantonio T., Scranton R., Crittenden R. G., Nichol R. C., Boughn S. P., Myers A. D., Richards G. T., 2008, Phys. Rev. D, 77, 123520 Gomes Z., Camera S., Jarvis M. J., Hale C., Fonseca J., 2019, MNRAS, 492, 1513 Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759 Hale C. L., Jarvis M. J., Delvecchio I., Hatfield P. W., Novak M., Smolčić V., Zamorani G., 2018, MNRAS, 474, 4133 Högbom J. A., 1974, A&AS, 15, 417 Hopkins A., Windhorst R., Cram L., Ekers R., 2000, Exp. Astron., 10, 419 Hotan A. W. et al., 2021, Publ. Astron. Soc. Aust., 38, e009 Hunter J. D., 2007, Comput. Sci. Eng., 9, 90 Hurley-Walker N. et al., 2016, MNRAS, 464, 1146 Jarvis M., Bacon D., Blake C., Brown M., Lindsay S., Raccanelli A., Santos M., Schwarz D. J., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14) Kovetz E. D., Raccanelli A., Rahman M., 2017, MNRAS, 468, 3650 McDonald P., Seljak U., 2009, J. Cosmol. Astropart. Phys., 2009, 007 Norris R. P. et al., 2011, Publ. Astron. Soc. Aust., 28, 215 Oliphant T., 2006, Guide to NumPy. Trelgol Publishing, United States Price-Whelan A. M. et al., 2018, AJ, 156, 123 Raccanelli A., Bonaldi A., Negrello M., Matarrese S., Tormen G., de Zotti G., 2008, MNRAS, 386, 2161 Raccanelli A. et al., 2012, MNRAS, 424, 801 Raccanelli A. et al., 2015, J. Cosmol. Astropart. Phys., 2015, 042 Santos M. G. et al., 2017, preprint (arXiv:1709.06099) Seljak U., 2009, Phys. Rev. Lett., 102, 021302 Shimwell T. W. et al., 2017, A&A, 598, A104 Shimwell T. W. et al., 2019, A&A, 622, A1 Square Kilometre Array Cosmology Science Working Group et al., 2020, Publ. Astron. Soc. Aust., 37, e007 Tingay S. J. et al., 2013, Publ. Astron. Soc. Aust., 30, e007 Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709 Vernstrom T., Scott D., Wall J. V., Condon J. J., Cotton W. D., Kellermann K. I., Perley R. A., 2016, MNRAS, 462, 2934 Virtanen P. et al., 2019, Nat. Methods, 17, 261 Wilman R. J. et al., 2008, MNRAS, 388, 1335 Windhorst R. A., 2003, New Astron. Rev., 47, 357 Zheng H. et al., 2017, MNRAS, 464, 3486 Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, J. Open Source Softw., 4, 1298 APPENDI X A : BI AS O F A N ENSEMBLE The bias is an ansatz that relates the number density of galaxies or other tracers to that of matter , δg(x) = ng(x) n̄g − 1 = bδm(x) , (A1) where δg(x) is the normalized galaxy number density fluctuation at some position in space x, ng(x) is the galaxy number density at that position, δm(x) is the matter density fluctuation at that position, n̄g is the homogeneous number density, and b is the bias. If we have multiple galaxy species, each with its own bias value (i.e. galaxies that have formed differently and trace the underlying matter density in a different way), how can we estimate the total bias of the ensemble? First, we write down what the fluctuation will be in terms of the combined ensemble of tracers, and some total bias, nA g (x) + nB g (x) + . . . + nX g (x) n̄A g + n̄B g + . . . + n̄X g − 1 = btotalδm(x) , (A2) Now we assume a form of the total bias btotal, such that the total is the sum of the individual biases {bA, bB, . . . , bX} in the form btotal = bAn̄A g + bBn̄B g + . . . + bXn̄X g n̄A g + n̄B g + . . . + n̄X g . (A3) It can therefore easily be shown that this is the only form of the total bias that satisfies the definition of the linear bias given in equation (A1) for any and all of the tracers when considered individually. This paper has been typeset from a TEX/LATEX file prepared by the author. MNRAS 506, 4121–4130 (2021) D ow nloaded from https://academ ic.oup.com /m nras/article/506/3/4121/6317625 by Biblioteca U niversidad de Zaragoza user on 25 February 2026 http://dx.doi.org/ 10.1103/PhysRevD.77.123520 http://dx.doi.org/10.1093/mnras/stz3581 http://dx.doi.org/10.1086/427976 http://dx.doi.org/10.1093/mnras/stx2954 http://dx.doi.org/10.1017/pasa.2021.1 http://dx.doi.org/10.1109/MCSE.2007.55 http://dx.doi.org/10.1093/mnras/stw2337 http://dx.doi.org/10.1093/mnras/stx691 http://dx.doi.org/10.1088/1475-7516/2009/10/007 http://dx.doi.org/10.1071/AS11021 http://dx.doi.org/10.3847/1538-3881/aabc4f http://dx.doi.org/ 10.1111/j.1365-2966.2008.13189.x http://dx.doi.org/10.1111/j.1365-2966.2012.20634.x http://dx.doi.org/10.1088/1475-7516/2015/01/042 http://arxiv.org/abs/1709.06099 http://dx.doi.org/10.1103/PhysRevLett.102.021302 http://dx.doi.org/10.1051/0004-6361/201629313 http://dx.doi.org/10.1051/0004-6361/201833559 http://dx.doi.org/10.1017/pasa.2019.51 http://dx.doi.org/10.1017/pasa.2012.007 http://dx.doi.org/10.1086/591439 http://dx.doi.org/10.1093/mnras/stw1836 http://dx.doi.org/10.1111/j.1365-2966.2008.13486.x http://dx.doi.org/10.1016/S1387-6473(03)00045-9 http://dx.doi.org/10.1093/mnras/stw2525 http://dx.doi.org/10.21105/joss.01298