Confusion noise due to clustered extragalactic point sources. Application of logarithmic cumulants for parameter estimation

TL;DR

This paper develops a method using logarithmic cumulants to estimate parameters of clustered extragalactic sources affecting confusion noise, improving accuracy over previous techniques by accounting for source correlations.

Contribution

It extends Zolotarev's theorem to compute cumulants of the logarithm of intensity, enabling better parameter estimation for clustered source populations.

Findings

Cumulants effectively determine number count parameters and correlation scale.

Method outperforms previous techniques by an order of magnitude in accuracy.

Estimates remain robust even with uncertain correlation scale.

Abstract

The calculation of the characteristic function of the signal fluctuations due to clustered astrophysical sources is performed in this paper. For the typical case of power-law differential number counts and two-point angular correlation function, we present an extension of Zolotarev's theorem that allows us to compute the cumulants of the logarithm of the absolute value of the intensity. As a test, simulations based on recent observations of radio galaxies are then carried out, showing that these cumulants can be very useful for determining the fundamental parameters defining the number counts and the correlation. If the angular correlation scale of the observed source population is known, the method presented here is able to obtain estimators of the amplitude and slope of the power-law number counts with mean absolute errors that are one order of magnitude better than previous techniques, that did not take into account the correlation. Even if the scale of correlation is not well known, the method is able to estimate it and still performs much better than if the effect of correlations is not considered.

Tables

Table 1. Table 1: Mean absolute percentage error (MAPE) for the estimation of the α 𝛼 \alpha and k 𝑘 k parameters as a function of α 𝛼 \alpha for one of the two kinds of astrophysical correlation described in the main text (Blake & Wall, 2002 ) . As a comparison, in the fourth and fifth columns the same quantities are computed using the estimators for the case of no spatial correlation (Herranz et al., 2004 ) .

$\alpha$	MAPEα	MAPEk	MAPEα,nc	MAPEk,nc
0.6	0.12	1	2.2	15
0.7	0.2	1.2	6	37
0.8	0.4	4	14	65
0.9	0.6	10	32	91

Equations

$$\varphi (\omega )=\exp \left[-i\mu \omega -\gamma |\omega {|}^{\alpha }\left(1+i\beta sgn(\omega )\tan (\alpha \pi /2)\right)\right],$$

$${\Omega }_e=\int h(\theta ,\varphi {)}^{\eta -1}d\Omega ,$$

$$<\rho ({\mathbf{x}}_1)\rho ({\mathbf{x}}_2)>={\mu }^2(1+\omega ({\mathbf{x}}_1,{\mathbf{x}}_2))$$

$$\omega (\theta )={\left(\frac{\theta }{{\theta }_0}\right)}^{-\delta }$$

$$\varphi (\omega )=\exp \left[\sum _{n=1}^{\infty }\frac{1}{n!}{\int }_{{\Omega }_b}\dots {\int }_{{\Omega }_b}\prod _{j=1}^n{\int }_{S_j}(e^{-i\omega S_{j}h({\mathbf{x}}_{\mathbf{j}})}-1)R(S_j)d{S}_j{w}_n({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{\mathbf{n}})\mathbf{d}{\mathbf{x}}_1\dots \mathbf{d}{\mathbf{x}}_{\mathbf{n}}\right],$$

$$\varphi (\omega )=\exp \left[\sum _{n=1}^{\infty }\frac{1}{n!}{\int }_{{\Omega }_b}\dots {\int }_{{\Omega }_b}\prod _{j=1}^n{h}^{\eta -1}({\mathbf{x}}_{\mathbf{j}})w_n({\mathbf{x}}_1,\dots ,{\mathbf{x}}_{\mathbf{n}})\mathbf{d}{\mathbf{x}}_1\dots \mathbf{d}{\mathbf{x}}_{\mathbf{n}}\prod _{j=1}^n{\int }_{s_j}(e^{-i\omega s_j}-1)k{s}_j^{-\eta }d{s}_j\right].$$

$${\int }_s(e^{-i\omega s}-1)k{s}^{-\eta }ds=k|\omega {|}^{\eta -1}\Gamma (1-\eta )\left(\sin \pi \eta /2-i\mathrm{sgn}(\omega )\cos \pi \eta /2\right),$$

$${\Omega }_n={\int }_{{\Omega }_b}\dots {\int }_{{\Omega }_b}\prod _{j=1}^n{h}^{\eta -1}({\mathbf{x}}_{\mathbf{j}})w_n({\mathbf{x}}_1,...,{\mathbf{x}}_{\mathbf{n}})\mathbf{d}{\mathbf{x}}_1\dots \mathbf{d}{\mathbf{x}}_{\mathbf{n}}.$$

$$\varphi (\omega )=\exp \sum _{n=1}^{\infty }-{\gamma }_n|\omega {|}^{n\alpha }(1+i\mathrm{sgn}(\omega )\tan (n\pi \alpha /2)),$$

$${\gamma }_n=-\frac{k^n{\Gamma }^n(-\alpha ){\Omega }_n\cos (n\pi \alpha /2)}{n!}.$$

$$h(\theta )=e^{-4\log 2(\theta /{\theta }_b{)}^2},$$

$${\Omega }_1=\frac{\pi {\theta }_b^2}{4\log 2\alpha }.$$

$$w_2({\mathbf{x}}_1,{\mathbf{x}}_2)={\left(\frac{{\theta }_{12}}{{\theta }_0}\right)}^{-\delta },$$

$${\Omega }_2=2\pi {\theta }_0^{\delta }{\left(\frac{{\theta }_b^2}{4\log 2\alpha }\right)}^{(4-\delta )/2}\Gamma \left(\frac{4-\delta }{2}\right)G(\delta )$$

$$G(\delta )={\int }_0^{\pi /2}{\int }_0^{\pi /2}\frac{\sin xdxdy}{(1+\sin x-2\sin x{\sin }^{2}y{)}^{\delta /2}}.$$

$$G(\delta )=\frac{\pi 2^{-\delta /2}}{(2-\delta )}.$$

$${\Omega }_2=\frac{{\pi }^2}{4}{\theta }_0^{\delta }{\left(\frac{{\theta }_b^2}{2\log 2\alpha }\right)}^{(4-\delta )/2}\Gamma \left(1-\frac{\delta }{2}\right).$$

$$E\left((\log |I|{)}^k\right)=\frac{d^k}{d{p}^k}{\left(E\left(|I{|}^p\right)\right)}_{p=0}.$$

$$E\left(|I{|}^p\right)=\frac{p}{\alpha \Gamma (1-p)\cos (p\pi /2)}\mathrm{Real}\left[{\int }_0^{\infty }\frac{1-e^{-au-b{u}^2}}{u^{p/\alpha +1}}du\right]$$

$$E(|I{|}^p)=\frac{\Gamma (1-p/\alpha )}{\Gamma (1-p)\cos (p\pi /2)}\mathrm{Real}\left[b^{p/(2\alpha )}H(p/\alpha ,z)\right],$$

$$E\left(|I{|}^p\right)=\frac{p}{\alpha \Gamma (1-p)\cos (p\pi /2)}\mathrm{Real}\left[M\left(\alpha ,p\right)\right],$$

$$M\left(\alpha ,p\right)=b^{p/(2\alpha )}\left[z\Gamma \left(\frac{1}{2}-\frac{p}{2\alpha }\right){}_1{F}_1\left(\frac{1}{2}-\frac{p}{2\alpha };\frac{3}{2};z^2\right)-\frac{1}{2}\Gamma \left(-\frac{p}{2\alpha }\right){}_1{F}_1\left(-\frac{p}{2\alpha };\frac{1}{2};z^2\right)\right]$$

$$E\left(|I{|}^p\right)=\frac{\Gamma (1-p/\alpha )}{\Gamma (1-p)\cos (p\pi /2)}\mathrm{Real}\left[a^{p/\alpha }\right],$$

$$E(\log |I|)=\left(\frac{1}{\alpha }-1\right)\gamma +\frac{1}{2\alpha }\log \frac{k^2{\Gamma }^2(-\alpha ){\Omega }_2}{2}+\frac{1}{\alpha }\mathrm{Real}\left[\frac{dH}{du}(0,z)\right]$$

$$\mathrm{var}(\log |I|)$$

$$z=\frac{{\Omega }_1}{\sqrt{2{\Omega }_2}}i={\left(\frac{{\theta }_b^2}{2\log 2\alpha {\theta }_0^2}\right)}^{\delta /4}\frac{1}{\sqrt{2\Gamma (1-\delta /2)}}i.$$

$$E(\log |I|)=\left(\frac{1}{\alpha }-1\right)\gamma +\frac{1}{\alpha }\log (k|\Gamma (-\alpha )|{\Omega }_1)+\frac{1}{\alpha }\left\{\mathrm{Real}\left[\frac{dH}{du}(0,z)\right]-\log (2|z|)\right\}.$$

$$r=\frac{k^2{\Gamma }^2(-\alpha ){\Omega }_2}{2}$$

$$\log E\left(|I{|}^p\right)$$

$${\kappa }_1=\psi (1)\left(1-\frac{1}{\alpha }\right)+\frac{1}{\alpha }\left(\frac{\log r}{2}+H_1^{\prime }(0,z)\right)$$

Full text

Extragalactic point sources, Statistical astronomy

††journal: PASP

1 Introduction

One of the main goals of modern Astronomy is the study of the formation and evolution of galaxies, as they constitute an extremely useful probe for Cosmology –the study of the origin, structure and ultimate fate of the whole Universe– and for extragalactic astronomy. Modern telescopes are able to reach epochs at which the Universe was ten percent of its present age and to observe billions of galaxies; these numbers make it imperative for astronomers to study galaxies statistically. In absence of distance information, extragalactic point sources (i.e., galaxies seen as point-like objects inside the observational beam, usually assumed as a circular Gaussian distribution whose angular dispersion determines its Full Width Half Maximum, FWHM) appear as luminous spots in the particular sky area actually surveyed. Their observed flux density depends on many factors including the distance to that particular galaxy, its evolutionary stage, the way part or all of it has been obscured by matter along the line of sight, etc. As a result of these many effects, the distribution of the number of galaxies as a function of their observed flux densities, $S(\nu)$, (differential number counts,dN/dS) is a very useful cosmological probe. In many cases of interest, and at least for some ranges of flux densities, the differential number counts can be well approximated by a simple power law $dN/dS=R(S)=k\,S^{-\eta}$ where the amplitude $k$ and the slope $\eta$ contain a wealth of information about the physical and statistical properties of the underlying galaxy populations, and hence of the whole Universe.

Differential number counts of extragalactic sources are usually measured (or estimated) by different analysis techniques, depending on the sensitivity of the detection instrument and on the observational band of interest: direct resolved detection, Bayesian detection, stacking analysis or probability of deflection analysis, $P(D)$, (see, e.g. De Zotti et al., 2010; Dole et al., 2004; Berta et al., 2011; Planck Collaboration et al., 2011, 2013; Harrison et al., 2016, etc.). Although the amplitude and slope of the number counts are usually known (or at last predictable according to current models of galaxy formation and evolution) for a wide range of wavelengths and at high flux densities, there are still many remaining gaps in our current knowledge, in particular at very faint flux densities, that are much more difficult to reach. Thus, astronomers who try to determine these two numbers across different regions of the observable electromagnetic spectrum can encounter great difficulties in doing that.

To make things more interesting, due to the fundamental law of diffraction and also to inevitable instrumental limitations, telescopes cannot observe the sky with infinite image definition. Any point of the sky observed through a telescope is blurred by the point spread function of such telescope. With a given beam, the final signal in a sky image is a mixing of bright and faint sources in which the brightest ones are still individually detectable above a given detection limit, $S_{lim}$. Sources at flux densities below this limit generate the so called confusion noise (see, e.g., Scheuer, 1957).

The probability distribution of the measured intensity, $f(I)$, where I is the signal response at a point produced by the sum of source fluxes convolved with the beam 111$I(x)=\Sigma\,S_{n}\,h(x-x_{n})$ with $h$ the beam profile and $x_{n}$ the position of the source with flux $S_{n}$, is also called $P(D)$ in the astrophysical jargon, and gives the total contribution of sources above and below the detection limit. However, P(D) is specially useful to estimate and partially reconstruct the number counts at very faint flux densities, below the detection limit of the instrument, where they cannot be directly measured.This situation is very common in astronomical images and it has been studied first at radio and X-ray frequencies.

The $P(D)$ distribution given by unresolved sources is strongly non-Gaussian and shows a long, positive tail. A wealth of information is encoded in the indistinct part of the distribution far from the tail. Hence the interest in Astronomy for the study and statistical characterization of this confusion noise.

The confusion noise phenomenon has been widely studied in the astronomical literature since the pioneering work of Scheuer (1957). That work was subsequently expanded by Hewish (1961). Condon (1974) studied the signal observed with a beam for differential number counts that follow a power law and are uncorrelated, $dN/dS=R(S)=k\,S^{-\eta}$. He wrote an analytical expression, formula (8) of his paper, for the characteristic function $\phi(\omega)$ of the observed intensity. However, the associated probability density function (pdf) $f(I)$, defined as the inverse Fourier transform of $\phi(\omega)$, cannot be calculated analytically. Many years later, Herranz et al. (2004) showed that Condon’s characteristic function was a particular case of an alpha-stable distribution $S_{\alpha}(\beta,\gamma,\mu)$, first introduced in Lévy (1925):

[TABLE]

For a general alpha-stable distribution: $\gamma>0,\,-1\leq\beta\leq 1$ and $0<\alpha\leq 2$ . In our case, $\alpha=\eta-1$, $\beta$ is equal to one222That is, the confusion noise generated by a power law-distributed population of compact sources corresponds to the particular case of an alpha-stable distribution of maximum positive skewness parameter. For $\alpha=1$ the alpha-stable expression is different from (1) and $\mu$ and $\gamma$ can be obtained in a simple way from $k$, $\eta$ and the effective beam size of the experiment

[TABLE]

where $h(\theta,\phi)$ is the beam profile. See formulas (10) and (11) in Herranz et al. (2004). Taking advantage of the theory of alpha-stable distributions (Lévy (1925); Zolotarev (1986); Shao & Nikias (1993); Kuruoğlu (2001); Samoradnitsky (2017)), Herranz et al. also defined suitable logarithmic moment estimators based on Zolotarev’s theorem (this theorem allowed the calculation of the expectation of $|I|^{p}$, $E\left(|I|^{p}\right)$). Then, they applied these estimators to the determination of $k$ and $\eta$ for power law number counts.

However, the presence of clustered –that is, spatially correlated– sources, so common in Astronomy, is not taken into account in any of the cited works. Barcons (1992) wrote a pioneering work in which he found a general expression for the characteristic function of the observed intensity, considering $n$-point correlations. In a similar way, by using random field theory, Takeuchi & Ishii (2004) wrote the Laplace transform of the same pdf. The outline of the paper is as follows. In Section 2 we briefly present and discuss the differential number counts of astrophysical (extragalactic) sources and their correlation functions. Then, inspired by the two last above-quoted papers, we manage to obtain a general $\phi(\omega)$ for the case of power law counts and general correlations in Section 3. In order to obtain precise results, finally we only consider two-point correlations and prove that, in this case, the observed intensity can be written as the sum of two independent random variables that follow alpha-stable distributions with parameters $\alpha$ and $2\alpha$. In Section 4, we generalize Zolotarev’s theorem for the clustering case and calculate the corresponding logarithmic moments and cumulants. The proof of this generalization is provided in Appendix A. In Section 5, we explore the symmetric case, defined by subtracting intensities at different points, distant enough to maintain statistical independence. This symmetrization enables us to simplify the formulas. In Section 6, we illustrate with simulations the usefulness of the logarithmic cumulants to compute the parameters defining the counts and the correlation. Finally, we draw our main conclusions in section 7.

2 Statistics of the distribution of extragalactic point sources

2.1 Differential number counts

The study of differential $dN/dS(\nu)$, and of integral, $N(>S(\nu))$, number counts of extragalactic point sources (i.e., galaxies seen as point-like objects in the sky given that their projected angular dimension is less than the beam size of the telescope/radio-telescope) is a field of great historical relevance in astrophysics. In fact, these statistics can be used to efficiently constrain the evolution with cosmic time of the observed galaxy populations at each selected frequency bands. This is easily performed by comparing the slope, $\eta$, of the observed number counts with the one foreseen for the static Euclidean universe, $\eta=2.5$, filled with galaxies uniformly distributed in the surveyed volume, which constitutes the benchmark comparison model. Any deviation in the observed number counts from the predictions of the static Euclidean universe implies evolution with cosmic time of the underlying galaxy population (see, e.g., Longair (1966) for a general and basic discussion on this fundamental comparison and on other related topics; see also Blain & Longair (1993) for an interesting application of the number counts distribution and of related statistics to the submillimeter wavebands of the electromagnetic spectrum).

The rapid development of radio antennas after the Second World War allowed that number counts of radio selected –at cm to mm wavelengths– extragalactic sources were the first ones to be measured: the very first results can be dated back to Mills (1952). Soon after, number counts of radio sources achieved instant notoriety when Ryle (1955) and Ryle & Scheuer (1955) announced that the slope of the observed counts at very bright flux densities ($S(\nu)\geq 1$ Jy) was steeper –i.e., $\eta>2.5$– than the one corresponding to a static Euclidean universe, that implied extragalactic radio sources must be evolving with cosmic time in space density or luminosity. Since then, many more surveys have been carried out to measure the source counts at various radio frequencies, both extended to the whole sky or to a limited sky area, which confirmed the importance of cosmological evolution (see, e.g., De Zotti et al. (2010) for a recent and comprehensive review on the subject).

More recent works at mm/sub-mm wavelengths or in the far-infrared domain of the electromagnetic spectrum (see, e.g., Scott et al., 2012; Hatsukade et al., 2013a, b; Carniani et al., 2015; Vernstrom et al., 2014; Whittam et al., 2017a, b) show that the differential number counts of faint extragalactic point sources can be well fitted usually by sub-euclidean count slopes $\eta\leq 2.5$ down to sub-mJy or even microJansky fluxes. Investigating the number counts at these faint fluxes is important for understanding the evolution of sources at low luminosities and/or high redshifts $(z>2)$. There is an open debate about what might be happening at flux densities fainter than the current limits ($S(\nu)<1$-10 $\mu$Jy) of the source number counts in the radio bands (see the works of Condon et al. (2012) and Vernstrom et al. (2014) at 3 GHz). New source populations, related to the epoch of galaxy formation, could be detected and their early cosmological evolution could be studied by the analysis of sub-$\mu$Jy number counts. Recent works in this field will be used in Section 6 to simulate fluxes in order to test our new technique for the determination of $k$ and $\eta$.

2.2 n-point and 2-point angular correlation functions

The angular distribution of sources on the celestial sphere is represented mathematically by a random field, $\rho(\mathbf{x})$, with $\rho(\mathbf{x})$ the number of sources per solid angle unit at a given point. All the information about the field is given by its mean $\mu=<\rho(\mathbf{x}>=\frac{N}{4\pi}$, with N the total number of sources (we assume that the field is homogeneous), the two-point correlation function $\omega(\mathbf{x_{1}},\mathbf{x_{2}})$, related to the two-point moment by

[TABLE]

and the n-point correlations $\omega_{n}(\mathbf{x_{1}},.........,\mathbf{x_{n}})$, that can be expressed in terms of higher-order cumulants of the field, (see Barcons, 1992), and are zero unless all positions at which they are calculated are correlated.

Few data have been so far collected on the correlations among galaxies of order higher than two (see, e.g., Meiksin et al., 1992, for the IRAS 1.2 Jy sample). This is obviously due to the lack of sufficient data (point sources) to calculate /estimate three- or n-point correlation functions than could be reliable, although in a limited angular scale interval. Historically, the determination of the three-point correlation function (3PCF) of galaxies was pioneered by Peebles & Groth (1975) using the Lick and Zwicky angular catalogs of galaxies. After this first outcome, among other results we can single out Verde et al. (2002) that computed the bispectrum –i.e. the Fourier counterpart of the 3PCF – using galaxies in the 2-Degree Field Galaxy Redshift Survey, and concluded that the non-linear bias was consistent with zero. More recent results still seem to find consistency with the usual hierarchical clustering model –the process by which larger structures are formed through the continuous merging of smaller structures– but it will be necessary to rely on future very large surveys (e.g., Euclid satellite, Large Synoptic Sky Telescope (LSST), ALMA, JWST, etc.) for obtaining new and more conclusive data on the 3PCF (see, e.g., Slepian & Eisenstein, 2015).

In view of the relative paucity of current data on the 3PCF and on higher order correlations, although we will write a general formula for the $n$-point correlation case, we will basically work with the two-point correlation function, that is, we neglect all higher-order contributions. The two-point correlation only depends on the angular distance $\theta$, –i.e. the underlying field is assumed to be homogeneous and isotropic–, and has been usually approximated by a power-law with a characteristic correlation angular distance $\theta_{0}$

[TABLE]

after the first evidence of it was found by Totsuji & Kihara (1969) and subsequently confirmed, with scaling relations demonstrated, by Groth & Peebles (1977). For instance, Loan et al. (1997) found $\omega(\theta)=0.01\theta^{-0.8}$, for sources at 5 GHz and Blake & Wall (2002) obtained $\omega(\theta)=0.001\,\theta^{-0.8}$, that is $\theta_{0}=0.64^{\prime\prime}$, at 1.4 GHz. More recently, Vernstrom et al. (2014) approximated the correlation of sources at 3 GHz with a similar power-law $\omega(\theta)=(\theta/\theta_{0})^{-0.8}$ but with $\theta_{0}=0.06^{\prime\prime}$ .

These results at radio frequencies are also confirmed in other frequency bands of the electromagnetic spectrum. For example, recent measurements of the angular correlation function of sub-millimeter galaxies identified in four of the five fields observed in the Herschel Astrophysical Terahertz Large Area Survey works (H-ATLAS, Eales et al., 2010) show that a power-law approximation of $\omega(\theta)$ with $\delta=0.8$ is still a very good approximation up to the highest redshift interval in which a signal is detected (Amvrosiadis et al., 2019). In all these cases, $\theta_{0}$, the characteristic correlation angle, is the only parameter that changes whereas $\delta=0.8$ remains fixed, which is related to the usual slope of the 3-dimensional spatial correlation, $\gamma=1.8$ (see, e.g., Davis & Peebles, 1977). Therefore, although our general results do not depend on the specific value of $\delta$, we will use $\delta=0.8$ in our simulations and different values of $\theta_{0}$, chosen from relevant observations.

3 Characteristic function of the observed signal distribution

As commented before, we try to calculate the confusion noise produced by filtered point sources when we include the effect of clustered sources. A general approach to this problem was developed in Takeuchi & Ishii (2004) by using the mathematical techniques of point field theory. They were able to write a general expression, formula (46) of their paper, for the Laplace transform of the probability density function (pdf) of the intensity fluctuations, $f(I)$, being $I$ the signal response observed with a certain beam.

We use a similar formula to that of Takeuchi & Ishii (2004), but for the Fourier transform (characteristic function), since we intend to obtain a generalization for clustered sources of the results found in Herranz et al. (2004) for the case of unclustered sources. Therefore, our formula can be written

[TABLE]

where n is the number of points (positions of sources $\mathbf{x_{1}},...,\mathbf{x_{n}}$ ) involved in the n integrals over the beams, taking into account the n-point correlation; and j, that goes from 1 to n, indicates the flux and position of each individual source. We multiply the integrals over the fluxes, with $R(S_{j})$ the differential source number counts with a flux density $S_{j}$, $h(\mathbf{x_{j}})$ the beam pattern of the instrument, $\Omega_{b}$ the beam area and $w_{n}(\mathbf{x_{1}},...,\mathbf{x_{n}})$ the $n$-point angular correlation, that encapsulates the statistical information about the spatial distribution of the sources, (see Barcons, 1992). To summarize, we calculate the exponential of a sum of integrals that include, in principle, the contribution of all n-point correlations. This sum is the characteristic function of the intensity , $I$.

From (5) we can find a simpler expression if we assume that the fluxes follow a power-law distribution $R(S)=kS^{-\eta}$. When we introduce this power-law formula in (5) and perform the change of variable $S_{j}h(\mathbf{x_{j}})=s_{j}$ we have

[TABLE]

The integral in $s_{j}$ can be readily calculated

[TABLE]

where $\Gamma$ is the gamma function. The integral of the real part converges for $1<\eta<3$ and the integral of the imaginary part for $0<\eta<2$. Therefore, both integrals are well defined for $1<\eta<2$. In the case without correlations, if we subtract independent intensities, we obtain a symmetrized random variable and are able to remove the imaginary part. Then, the characteristic function is real and if we define $\alpha=\eta-1$, it is an alpha-stable distribution within the standard range $0<\alpha<2$. However, in the general clustering case we have to resign ourselves to the convergence for $0<\alpha<1$.

Now, we define

[TABLE]

$\Omega_{n}$ is an integral incorporating the information about the beam pattern and the correlation. Putting together all these integrals, we arrive at the following formula for the Fourier transform

[TABLE]

with

[TABLE]

If we consider the special case $n=1$ (no correlations) we recover, after some simple algebra, formula (7) of Herranz et al. (2004). Handling the formula above in a general case is very complicated. In order to simplify matters we will only consider the two-point correlation, i.e. we neglect higher-order correlations from now on (see Subsection 2.2), $\gamma_{n}=0$ for $n>2$. All the formulas derived in this paper are valid for generic beam profiles and generic two-point correlation functions, provided that they are written in terms of $\Omega_{1}$ and $\Omega_{2}$. However, for the sake of simplicity, we will assume that the beam profile is Gaussian, (we will also assume a small beam and therefore work on the plane)

[TABLE]

with $\theta_{b}$ the full width half maximum (FWHM) of the beam. With this assumption the calculation of $\Omega_{1}$ is straightforward

[TABLE]

To compute $\Omega_{2}$ we need to specify a formula for the two-point angular correlation. We will use a typical power law, see formula (67) in Takeuchi & Ishii (2004) and Section 2

[TABLE]

where $\theta_{12}$ is the angular distance between the points and $\theta_{0}$ the characteristic correlation angular distance. By substituting this formula in $\Omega_{2}$, we obtain

[TABLE]

[TABLE]

This last integral converges for $0<\delta<2$ and can be worked out analytically

[TABLE]

We further simplify the formula for $\Omega_{2}$ and find

[TABLE]

Thus, if we knew the parameters $\alpha,k,\theta_{0},\delta$ and $\theta_{b}$, we could obtain $\phi(\omega)$ and, applying the inverse Fourier tranform, the pdf $f(I)$. $\phi(\omega)$, as given in formula (9) with $\gamma_{n}=0$ for $n>2$, is the characteristic function of the sum of two independent random variables with alpha-stable distibutions $S_{\alpha}(1,\gamma_{1},0)$ and $S_{2\alpha}(1,\gamma_{2},0)$ respectively. For the characteristic function to be invertible $\gamma_{2}$ must be positive, i.e. $1/2<\alpha<3/2$. $\gamma_{1}$ is always positive, since $0<\alpha<1$ in our calculations.

4 Logarithmic estimators

Our goal is to determine the unknown parameters $\alpha$, $k$, $\theta_{0}$, and $\delta$ from simulations based on real data and, in order to do that, we will follow the scheme presented in Herranz et al. (2004), that is, we will make use of logarithmic estimators. The logarithmic moments of our distribution can be found by using

[TABLE]

The right-hand side of this formula can be obtained in the case of an alpha-stable distribution by means of Zolotarev’s Theorem Zolotarev (1986); Kuruoğlu (2001), that allows us to calculate $E\left(|I|^{p}\right)$. We generalize Zolotarev’s theorem in Appendix A for the sum of two independent random variables with alpha-stable distributions. For our particular case, $\alpha_{2}=2\alpha_{1}$, we must compute the following integral

[TABLE]

with $a=\gamma_{1}(1+i\tan(\alpha\pi/2))$, $b=\gamma_{2}(1+i\tan(\alpha\pi))$, the integral converges if $\gamma_{1}$ and $\gamma_{2}$ are positive ($1/2<\alpha<1$) and for $0<p/\alpha<1$. In Zolotarev’s theorem, $b=0$. After performing this integral (see Appendix A), we find an extension of Zolotarev’s theorem

[TABLE]

where $H(p/\alpha,z)$ is the Hermite function Lebedev & Silverman (1972) and $z=a/2\sqrt{b}$. We can also express this formula in the following way

[TABLE]

where

[TABLE]

with ${}_{1}F_{1}$ the confluent hypergeometric function or Kummer function Lebedev & Silverman (1972). However, we prefer to use the Hermite function since the notation is more compact and the calculations simpler. If we perform the limit of (20) as $|z|$ tends to $\infty$ (|b| goes to zero), we find

[TABLE]

and we recover the formula for the case without correlations, Zolotarev’s theorem, Zolotarev (1986); Kuruoğlu (2001). Our next step is the differentiation of (20) to obtain the expectation and variance of $\log(|I|)$:

[TABLE]

[TABLE]

In these formulas $\gamma=0.5772...$ is the Euler constant, $\theta=\arctan(\tan(\alpha\pi))$, $u=p/\alpha$ and $z$ can be written as

[TABLE]

Equation (24) can be expanded by substituting $\Omega_{2}$ in terms of $\Omega_{1}$ and $z$

[TABLE]

Both (25) and (27) can be computed as $|z|$ approaches $\infty$, leading us to the formulas corresponding to the distribution of uncorrelated sources (formulas (21) and (22) in Herranz et al. (2004)). If we are interested in calculating the cumulants of $\log{|I|}$ in general, we can use the expression for $\log{E\left(|I|^{p}\right)}$ in terms of the Hermite functions and take the successive derivatives at $p=0$. If we define

[TABLE]

and separate the Hermite function in its real and imaginary part, $H_{1}+iH_{2}$, we can write

[TABLE]

Now, we can differentiate at $p=0$ and obtain the corresponding cumulants

[TABLE]

This is the mean, the same as (24). The following cumulant, the variance, is given by

[TABLE]

where $\psi^{(n-1)}(x)$ is the polygamma function, defined as

[TABLE]

The third-order cumulant is

[TABLE]

The expression for the fourth-order cumulant is too long to be written here, but quite straightforward to obtain. Besides, we can use the series expansion of the Hermite functions Lebedev & Silverman (1972) for calculating the derivatives of $H_{1}$ and $H_{2}$.

[TABLE]

For instance

[TABLE]

[TABLE]

[TABLE]

Higher order derivatives can be computed by iterative differentiation, giving rise to cumbersome formulas.

5 Symmetric case

For the sake of simplicity, we will now consider symmetrized sequences defined by

[TABLE]

for $k=1,\ldots,N/2$ and $N$ the number of independent data. It can be readily proved that in this case $a=2\gamma_{1}$ and $b=2\gamma_{2}$, so that the imaginary part of formula (9) vanishes. As before, we assume $\gamma_{n}=0$ for $n>2$, i.e. we only take into account two-point correlations. Please note that, in order to avoid correlations in the case of real data, we must consider data from distant enough pixels before applying (38). We can obtain the expressions for the logarithmic cumulants by using Hermite functions. Now, since $a$ and $b$ are real, so are

[TABLE]

and the Hermite function. This fact greatly simplifies the formulas, making these difference maps very useful for the determination of the parameters. If we define

[TABLE]

then $\Omega_{2}=(\Omega_{1}/\sigma)^{2}$ and $z=\sigma\cos{(\alpha\pi/2)}/\sqrt{-\cos{(\alpha\pi)}}$, and now we have

[TABLE]

with

[TABLE]

Bear in mind that since $1/2<\alpha<1$, $b$ is positive. Now, we can write the cumulants in a very simple way by taking succesive derivatives of (41) at $p=0$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (44) and (45) we could obtain $\alpha$ and $\sigma$ ($\theta_{0}$ and $\delta$ cannot be separated, and we must fix one, typically $\delta$, to determine the other). The variable $k$ can be isolated in (43), because it only appears in $\kappa_{1}$ inside $b$.

It would be possible to write more cumulants, by following the general pattern given by the polygamma functions, the $n$-th derivative of $\tan(p\pi/2)$ at $p=0$ and taking into account that the part involving the derivatives of $H(p/\alpha,z)$ follows the same pattern as the cumulants, expressed in terms of the derivatives of $H$ instead of the raw moments. The derivatives of $H$ can be calculated from the corresponding series (34). For instance

[TABLE]

[TABLE]

6 Simulations and Results

From now on, we will try to determine the parameters $k,\alpha,\sigma$ ( from $\sigma$ we could obtain $\theta_{0}$ assuming that we know $\delta$ ) by using our formulas for the cumulants and comparing their values with those obtained from simulations based on real data. The pdf of the symmetrized intensity $f(I^{s})$ can be calculated as the inverse Fourier transform of $\phi(\omega)$

[TABLE]

Please note that this formula is valid for the symmetrized intensity $I^{s}$, but not in the general case, in which we have to deal with the imaginary part. If we integrate this expression, we can calculate the cdf

[TABLE]

The cumulative $g(I^{s})$ can be computed numerically and used for generating simulations. However, we can also use the fact that our pdf is that of the sum of two symmetric alpha-stable variables with distributions $S_{\alpha}(0,a,0)$ and $S_{2\alpha}(0,b,0)$ respectively. Therefore, we simulate the corresponding alpha-stable variables Chambers et al. (1976); Weron (1996) and add them up. We simulate the first variable as

[TABLE]

and the second as

[TABLE]

where $t$ follows a uniform distribution between $-\pi/2$ and $\pi/2$ and w an exponential distribution with mean 1. We generate these data $x=x_{1}+x_{2}$, changing the values of $k$, $\alpha$, $\theta_{0}$ and $\delta$ and try to obtain $k$, $\alpha$ and $\sigma$ by using the logarithmic cumulants (43), (44), (45). The only caveat is that our results are only valid for $1/2<\alpha<1$ (however, see Condon et al. (2012); Vernstrom et al. (2014); Whittam et al. (2017a, b)). In the following, we will simulate symmetrized intensities $I_{s}$ by using the sum of two alpha-stable random numbers, $I_{s}=x_{1}+x_{2}$ , generated according to the counts and correlation parameters chosen from recent observations. We assume that the FWHM of the beam , $\theta_{b}$, is also known.

Among recent results on the analysis of the $P(D)$ – the confusion probability distribution, $P(D)$, of the deflections, $D$ (in Jy/beam), produced by point sources at flux densities below the flux detection limit of the survey – a partially new approach has been presented by Vernstrom et al. (2014), who analyzed a new sample of very faint radio sources. By this approach, “with a more robust model, and a comprehensive error analysis”, Vernstrom et al. (2014) were able to estimate the $\mu$Jy and sub–$\mu$Jy source counts by using new deep wide-band 3-GHz data in the Lockman Hole from the Karl G. Jansky Very Large Array (VLA). These authors found that the differential source counts between $0.05\mu$Jy and $0.2\mu$Jy present a slope of $\eta=1.79$. At the same time, they demonstrated that the source number counts can be constrained down to $\sim$50 nJy, a factor of $\sim$20 below the rms confusion noise.

In the Table 4 of Vernstrom et al. (2014) the differential number counts follow a power law with $k=1000$ and $\alpha=0.79$ at the lowest density fluxes probed by their observations. The correlation parameters are $\theta_{0}=0.06^{\prime\prime}$ and $\delta=0.8$; finally the FWHM of the experiment is $8^{\prime\prime}$, (Vernstrom et al., 2014). As a first application of our method, we use these numbers in a sort of toy model by applying the logarithmic cumulants in two different ways:

We carry out 100 simulations with the above parameters. Each simulation consists of $10^{6}$ symmetrized intensities, generated as the sum of alpha-stable random numbers as explained before. The number is large enough to have representative samples, though we have checked that numbers greater than $10^{4}$ could equally work. In principle, due to the simulation technique, we do not impose any cut in the simulated fluxes although the minimum fluxes are in the nanoJansky range.

For each simulation, we calculate, from our simulated data, the expectation and variance of $\log|I_{s}|$ and, assuming that we know the correlation parameters $\delta=0.8$ and $\theta_{0}=0.06^{\prime\prime}$, we obtain $\alpha$ from (44) and then $k$ from (43). When we compare the real values of $k$ and $\alpha$ with those obtained with our method we see that the mean of the absolute value of the relative error, taking into account our 100 simulations, is less than $0.1\%$ for $\alpha$ ($0.7\%$ for $k$) . If we calculate the expectation and variance of $\log|I_{s}|$ from our simulations and use instead formulas (24) and (25) of Herranz et al. (2004) to determine $\alpha$ and $k$,i. e. we assume no correlations, the error in $\alpha$ is still low, $3\%$, though much greater than with our formulas, but the error in $k$ increases clearly, $22\%$. If we simulate with other values of $\alpha$, $0.5<\alpha<1$,maintaining the value of the other parameters,the errors are very similar.

In this first case (knowledge of the correlation), we have chosen a second example, increasing the correlation distance $\theta_{0}=0.64^{\prime\prime}$ as in Blake & Wall (2002) and keeping the other parameters. Now the errors in $\alpha$ and $k$ are $0.4\%$ and $4\%$ respectively, for $\alpha=0.79$, $k=1000$, whereas they are $14\%$ and $65\%$ if we do not consider the correlations and apply the formulas in Herranz et al. (2004) to estimate the parameters.

We also calculate the relative errors for different values of $\alpha$ with the Blake & Wall (2002) correlation parameters. Table 1 shows that the errors are much higher if we do not take into account the real correlation and use the formulas in Herranz et al. (2004).

Finally, in Figures 1 and 2 we plot $\kappa_{1}$ and $\kappa_{2}$ obtained from (43) and (44), for the Vernstrom et al. (2014) and Blake & Wall (2002) cases, with values of $\alpha$ ranging from $0.5$ to $1$ and using the cited values for the other parameters. We can see that these cumulants differ more from the no correlation case the greater $\alpha$ and the correlation distance are. For $\alpha=0.5$, we recover the no correlation result, since $\gamma_{2}=0$. These figures are plotted just by substituting the parameters in the formulas, without the use of simulations.

In these examples, we have performed the simulations and recovered $\alpha$ and $k$ assuming a perfect knowledge of the correlations, which will not be, in general, the real case. However, we want to emphasize with our toy model that we could make large errors in the determination of these parameters if we do not consider the correlations at all, as can be seen in Table 1.

Now, we will assume, more realistically, that we know perfectly $\delta$ but we do not know $\theta_{0}$. As in case 1, we carry out 100 simulations with the Vernstrom et al. (2014) parameters, each consisting of $10^{6}$ symmetrized intensities and now, we calculate the cumulants for each simulation, solve (45) and (44) to obtain $\theta_{0}$ and $\alpha$, assuming $\delta=0.8$, and, as before, apply (43) to compute $k$. The mean absolute percentage errors in $\alpha$ and $k$ are $0.3\%$ and $2\%$ respectively, but the error in $\theta_{0}$ is greater, $13\%$.

For the Blake & Wall (2002) correlation, we simulate with $\delta=0.8$ and $\theta_{0}=0.64^{\prime\prime}$. We proceed as before, calculating the logarithmic cumulants from the simulations and solving (45), (44) and (43) to determine $\alpha$, $k$ and $\theta_{0}$, assuming $\delta=0.8$. The errors are now $1.6\%$, $10\%$ and $43\%$ for $\alpha$, $k$ and $\theta_{0}$.

These examples show that, in general, a certain knowledge of the correlation, even if we only know its slope, enables us to obtain a much better determination of the two parameters $k$ and $\alpha$, though in some cases the determination of the correlation characteristic angle is not very precise.

For the application of this method to real data, we should symmetrize the data, calculate the logarithmic cumulants from the symmetrized data and solve the system consisting of (43), (44) and (45), assuming a certain value of $\delta$.

An important issue is the influence of the instrumental noise in our results. It could be incorporated in the calculation of the characteristic function, formula (9), just multiplying by the characteristic function of the noise, e.g. in the case of Gaussian noise $\phi_{g}(\omega)=e^{-\sigma_{g}^{2}\omega^{2}/2}$. However, the calculation of the logarithmic cumulants in this case would be quite complicated. Logarithmic cumulants and minimum distance estimators (Ilow & Hatzinakos, 1998) were used in Herranz et al. (2004) to study the influence of Gaussian noise in their results. Their main conclusion was that if the noise dispersion $\sigma_{g}$ is low enough, the logarithmic cumulants, calculated without taking the noise into account, can be safely used to determine the parameters defining the counts. Otherwise, minimum distance estimators, that try to fit the characteristic function, perform better. In our paper, we intend to present the first method that allows the computation of the logarithmic cumulants in the clustering case and their use to estimate the counts and correlation parameters. We leave a careful study of the determination of the parameters in the clustering plus noise case for a further paper.

7 Conclusions

We study the characteristic function –the Fourier transform of the probability density function (pdf)– of the intensity produced by clustered extragalactic point sources filtered with a typical observational beam. We follow previous works by Barcons (1992) and Takeuchi & Ishii (2004) and find a new general expression for this characteristic function (9), (10), in the case of a power-law distribution for the differential number counts. 2. 2.

We analyze carefully the two-point correlation case $\gamma_{n}=0$ for $n>2$ in (9). Then, the observed intensity can be written as the sum of two independent random variables that follow alpha-stable distributions $S_{\alpha}(\gamma_{1},1,0)$, $S_{2\alpha}(\gamma_{2},1,0)$. We generalize Zolotarev’s theorem (20) and prove this generalization in Appendix A. Finally, we write formulas for the cumulants of the logarithm of the absolute value of the intensity (30), (31), (33). 3. 3.

We pay special attention to the symmetric case, where we subtract intensities at different points. In this situation, the general logarithmic cumulants happen to be quite simple (43)-(46). These formulas lead us to a general method for determining the parameters defining the differential number counts and the correlation. 4. 4.

Finally, we simulate intensity data with power law differential number counts and typical two-point correlation functions and show that the logarithmic cumulants enable us to calculate the parameters $k$, $\alpha$ and $\theta_{0}$ with small errors, assuming that we know $\delta$. We use the counts and correlations of Vernstrom et al. (2014) and Blake & Wall (2002), giving the relative errors in the text and in Table 1. In general, taking into account the correlation in the use of the logarithmic cumulants allows a better determination of the basic parameters. 5. 5.

The method can be applied to real data in a straightforward way. The data must be symmetrized, then the logarithmic cumulants can be calculated from the data and the parameters defining the counts $k$, $\alpha$ and the correlation angle $\theta_{0}$ found by solving (43), (44) and (45).

As a final remark, we can say that we have defined a new method to determine essential parameters related to the distribution of clustered astrophysical sources based on the calculation of logarithmic cumulants. The important issue of how to apply or extend the method to take into account the observational noise deserves a future detailed analysis.

Acknowledgements

We thank the Spanish MINECO for financial support under projects AYA2015-64508-P and AYA2015-65887-P. D. Herranz also thanks funding from the European Union’s Horizon 2020 research and innovation programme (COMPET-05-2015) under grant agreement number 687312 (RADIOFOREGROUNDS). Finally, we thank an anonymous referee for his/her very useful suggestions.

Appendix A AN EXTENSION OF ZOLOTAREV’S THEOREM.

An extension of Zolotarev’s theorem

If $X_{1}$ and $X_{2}$ are random independent variables that follow alpha-stable distributions $S_{\alpha_{1}}(\beta_{1},\gamma_{1},0)$, $S_{\alpha_{2}}(\beta_{2},\gamma_{2},0)$, respectively, with $k=\alpha_{2}/\alpha_{1}>1$ and both $\alpha_{1}$, $\alpha_{2}\neq 1$, then if $X=X_{1}+X_{2}$

[TABLE]

with $0<p<\alpha_{1}$, $a=\gamma_{1}(1+i\beta_{1}\tan{(\alpha_{1}\pi/2)})$, $b=\gamma_{2}(1+i\beta_{2}\tan{(\alpha_{2}\pi/2)})$, $z=a/(2b^{1/k})$ and

[TABLE]

For $k=2$, $G_{k}(\nu,z)=H(\nu,z)$, the Hermite function and we recover (20). If $p=1$ we must substitute $\pi/2$ for $\Gamma(1-p)\cos(p\pi/2)$ in (A1). Finally, (A1) is also valid if $\alpha_{1}=1\,(\alpha_{2}=1)$ provided $\beta_{1}=0$ $(\beta_{2}=0)$.

Proof.

First we calculate

[TABLE]

with the change$|X|t=u$, we obtain

[TABLE]

This result is valid for $0<p<2,\,p\neq 1$, however for $p=1$ we must substitute $\pi/2$ for $\Gamma(1-p)\,\cos(p\pi/2)$ in the calculations. We assume $p\neq 1$ from now on. Then

[TABLE]

Taking into account that $E(\cos(Xt))=E(\mathrm{Real}(e^{-iXt}))$ we rewrite (A5) as

[TABLE]

Since $E(e^{-iXt})$ is the characteristic function of $X$ and this is the product of the characteristic functions of $X_{1}$ and $X_{2}$, we obtain

[TABLE]

With the change $u=t^{\alpha_{1}}$, we have

[TABLE]

The integral converges for $0<p<\alpha_{1}$. In the case $k=2$, this is equation (19) of our paper. Finally, we make another substitution, $x=b^{1/k}u$. The integral, whose real part we have to calculate, can then be written as

[TABLE]

with $z=a\,b^{-1/k}/2$. After this last change, we have to add another integral in the complex plane, but its value is zero, since we integrate along an arc at infinite distance from the origin and the integrand approaches zero in this case. Our last integral can be readily worked out by expanding $e^{-2zx}$ as a power series. Then we obtain

[TABLE]

with

[TABLE]

and, integrating by parts

[TABLE]

Puting all the terms together

[TABLE]

multiplying by $\displaystyle\frac{p}{\Gamma(1-p)\cos{(p\pi/2)}\alpha_{1}}$, and taking the real part, we reach the final result

[TABLE]

For $k=2$, we recover the Hermite function (34), taking into account that

[TABLE]

If we only have an alpha-stable distribution $X=X_{1}$, then $b=0$ in (A8), this integral can be readily calculated and we find (23), i.e. Zolotarev’s theorem . As a final note, the theorem can be extended to general alpha-stable distributions $S_{\alpha_{1}}(\gamma_{1},\beta_{1},\mu_{1})$, $S_{\alpha_{2}}(\gamma_{2},\beta_{2},\mu_{2})$ if we write $E(|X-\mu_{1}-\mu_{2}|^{p})$ instead of $E(|X|^{p})$ in (A1).