Proving the short-wavelength approximation in Pulsar Timing Array gravitational-wave background searches

Gravitational waves (GWs) are ripples in the fabric of space-time, originating from some of the most violent events in the Universe, including the mergers of supermassive black holes. High frequency GWs from the merger of stellar-mass black holes were first detected by the Laser Interferometer Gravitational-wave Observatory (LIGO) in September 2015 [1], hailing the dawn of gravitational-wave astronomy. However, LIGO can only detect high frequency GWs, in the 100–1000Hz range. Similarly to electromagnetic radiation, different GW detectors are needed to probe different GW frequencies. Currently there are plans to launch a space-based GW detector in 2034—the Laser Interferometer Space Antenna (LISA) [2]—which will probe the millihertz GW frequency regime, thought to be populated primarily by merging supermassive black holes (SMBHs) in the 10⁵–10⁶ M_⊙ range. At the very low-frequency end of the GW spectrum, one expects to find nanohertz GWs from very massive inspiraling SMBHs, in the ${10}^{8}\mbox{--}{10}^{9}\,{M}_{\odot }$ range. These can be detected by timing millisecond pulsars, called a Pulsar Timing Array (PTA) [3–6]. Millisecond pulsars are excellent clocks, and delays or advances in their arrival times—inducing a timing residual—could signal the presence of GWs. PTA experiments are very active, and have been taking data for over a decade [7–9]. With a PTA, one can detect not only GWs from inspiraling SMBH binaries (SMBHBs), see e.g. [10, 11], but the GW background (GWB) from the cosmic merger history of SMBHBs [12–14]. This GWB is expected to be detected in the next few years [15, 16], with the details depending on the underlying astrophysics of the SMBH mergers. More details on PTAs can be found in recent review articles, e.g. [17–20], and an outline GW astrophysics covering nanohertz to kilohertz frequencies can be found in [21].

Indeed, a rigorous exploration and examination of the tools which will be used to make the first detection of a GWB is crucial. An isotropic GWB will induce characteristic correlations in the pulsar timing residuals. By cross-correlating these residuals, one expects to see a characteristic correlation called the Hellings and Downs curve [5]. Deviations from an isotropic GWB can be induced by nearby and/or particularly loud SMBHBs, inducing anisotropy in the GWB. Anisotropic GWBs will induce different correlations patterns, and have been explored by [22–26].

Here we prove analytically, and for the first time, that the Hellings and Downs curve can be extracted from the cross-correlated pulsar residuals, without making assumptions that the pulsars are all at the same distance L from the Earth. Part of this proof is a consequence of the application of the Riemann-Lebesgue lemma and the Lebesgue Dominated Convergence Theorem—well-known in the mathematics community, but less well-known in the field of GWs. We emphasize that no previous work has been able to do this analytically, though computer-aided integration has been used to verify one's intuition numerically.

The International PTA (IPTA) published combined data on 49 millisecond pulsars in their first data release [14]. These millisecond pulsars are the most stable natural astrophysical clocks known [27], and are regularly monitored by 8 radio telescopes: 5 in Europe [8], 2 in North America [7] and one in Australia [9]. PTAs take advantage of the precise arrival times of millisecond pulsars to enable GW detection.

The GWB is described in terms of its characteristic strain, h_c(f), with amplitude A at a reference frequency of 1/yr (e.g. [28]):

$$\begin{eqnarray}&&{h}_{c}(f)=A{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{-2/3}.\end{eqnarray} \tag{ 2.1 }$$

The current upper limits on A are difficult to compare, since it was recently discovered that errors in planetary masses and positions (called the solar system ephemeris model) can directly affect the limit on A [12, 29], and in some cases mimic a GWB signal.

While the current upper limit on A from NANOGrav can take this into account, and limit A < 1.35 × 10⁻¹⁵, other PTAs have not yet published updates to their limits. Projections the characteristic strain accessible with future IPTA and Square Kilometer Array (SKA) [30–32] detectors are shown in figure 1.

Zoom In Zoom Out Reset image size

The observed residuals due to the presence of a GWB with characteristic strain h_c(f) is described by the cross-power spectral density of pulsar 1 and pulsar 2 by

$$\begin{eqnarray}&&{S}_{\mathrm{1,2}}(f)=\displaystyle \frac{{{\rm{\Gamma }}}_{\mathrm{1,2}}({{fL}}_{1},{{fL}}_{2},\zeta ){h}_{c}^{2}(f)}{12{\pi }^{2}{f}^{3}},\end{eqnarray} \tag{ 2.2 }$$

see e.g. [34], where Γ_1,2 is the so-called overlap reduction function, which describes the GWB-induced correlation signature in the pulsar residuals. This is a function of the frequency of the GWB, the distance to the pulsars L_1,2, and the angular separation of the pulsars, ζ. PTA geometry is explored in detail in figure 2. For an isotropic GWB, this is called the Hellings and Downs curve [5], and for anisotropic GWBs see [22–25].

Zoom In Zoom Out Reset image size

In analogy with [22, 24], we present an overview of how one arrives to the Hellings and Downs curve. A source of GWs in direction $-\hat{{\rm{\Omega }}}$, see figure 2, generates a metric perturbation ${h}_{{ij}}(t,\hat{{\rm{\Omega }}})$, which we describe as a plane wave:

$$\begin{eqnarray}&&{h}_{{ij}}(t,\vec{x})=\displaystyle \sum _{A}{\int }_{-\infty }^{\infty }{df}{\int }_{{S}^{2}}d\hat{{\rm{\Omega }}}{h}_{A}(f,\hat{{\rm{\Omega }}}){e}^{i2\pi f(t-\hat{{\rm{\Omega }}}\cdot \vec{x})}{e}_{{ij}}^{A}(\hat{{\rm{\Omega }}}).\end{eqnarray} \tag{ 3.1 }$$

This can be decomposed over two polarization tensors ${e}_{{ij}}^{A}(\hat{{\rm{\Omega }}})$, and two independent polarization amplitudes ${h}_{A}(t,\hat{{\rm{\Omega }}})$ [35, 36]:

$$\begin{eqnarray}&&{h}_{{ij}}(t,\hat{{\rm{\Omega }}})={e}_{{ij}}^{+}(\hat{{\rm{\Omega }}}){h}_{+}(t,\hat{{\rm{\Omega }}})+{e}_{{ij}}^{\times }(\hat{{\rm{\Omega }}}){h}_{\times }(t,\hat{{\rm{\Omega }}}).\end{eqnarray} \tag{ 3.2 }$$

We note that General Relativity predicts only two independent polarizations, plus +, and cross ×, while other theories predict additional polarizations, such as breathing modes [37–39]. He we restrict ourselves to the well-known tensor polarizations, A = +, ×.

The ${e}_{{ij}}^{A}(\hat{{\rm{\Omega }}})$ polarization tensors are uniquely defined by specifying $\hat{m}$ and $\hat{n}$—the GW principal axes, illustrated in figure 2:

$$\begin{eqnarray}&&{e}_{{ij}}^{+}(\hat{{\rm{\Omega }}})={\hat{m}}_{i}{\hat{m}}_{j}-{\hat{n}}_{i}{\hat{n}}_{j},\hspace{1cm}{e}_{{ij}}^{\times }(\hat{{\rm{\Omega }}})={\hat{m}}_{i}{\hat{n}}_{j}+{\hat{n}}_{i}{\hat{m}}_{j}.\end{eqnarray} \tag{ 3.3 }$$

For a stationary, Gaussian, and unpolarized GWB, the polarization amplitudes satisfy (see e.g. [40]):

$$\begin{eqnarray}&&\langle {h}_{A}^{* }(f,\hat{{\rm{\Omega }}}){h}_{A^{\prime} }(f^{\prime} ,\hat{{\rm{\Omega }}}^{\prime} )\rangle ={\delta }^{2}(\hat{{\rm{\Omega }}},\hat{{\rm{\Omega }}}^{\prime} ){\delta }_{{AA}^{\prime} }\delta (f-f^{\prime} )H(f).\end{eqnarray} \tag{ 3.4 }$$

The metric perturbation will change the proper distance between the Earth and the pulsars, inducing an advance or delay in the pulsar pulse's arrival time at the Earth. Consider for example a millisecond pulsar with frequency ν₀ whose location in the sky is described by $\hat{p}$, at a distance L from the Earth. The metric perturbation affects the frequency of the radio pulses, ν, received at the radio telescope. This frequency shift is given by

$$\begin{eqnarray}&&z(t,\hat{{\rm{\Omega }}})\equiv \displaystyle \frac{\nu (t)-{\nu }_{0}}{{\nu }_{0}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\hat{p}}^{i}{\hat{p}}^{j}}{1+\hat{{\rm{\Omega }}}\cdot \hat{p}}{\rm{\Delta }}{h}_{{ij}}(t,\hat{{\rm{\Omega }}}),\end{eqnarray} \tag{ 3.5 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}{h}_{{ij}}(t,\hat{{\rm{\Omega }}})\equiv {h}_{{ij}}({t}_{e},\hat{{\rm{\Omega }}})-{h}_{{ij}}({t}_{p},\hat{{\rm{\Omega }}})\end{eqnarray} \tag{ 3.6 }$$

is the difference between the GW-induced metric perturbation at the Earth ${h}_{{ij}}({t}_{e},\hat{{\rm{\Omega }}})$, the Earth term, with coordinates $({t}_{e},{\vec{x}}_{e})$, and at the pulsar ${h}_{{ij}}({t}_{p},\hat{{\rm{\Omega }}})$, the pulsar term, with coordinates $({t}_{p},{\vec{x}}_{p})$:

$$\begin{eqnarray}&&{t}_{p}={t}_{e}-L,\quad {\vec{x}}_{p}=L\hat{p},\quad {\vec{x}}_{e}=0.\end{eqnarray} \tag{ 3.7 }$$

The indices 'e' and 'p' refer to the Earth and the pulsar, however, it is standard write t_e = t, see e.g. [22, 24, 41, 42].

We can now write (3.6), using (3.1) and (3.2) as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{h}_{{ij}}(t,\hat{{\rm{\Omega }}}) & = & \displaystyle \sum _{A}{\displaystyle \int }_{-\infty }^{\infty }{{dfe}}_{{ij}}^{A}(\hat{{\rm{\Omega }}})\ {h}_{A}(f,\hat{{\rm{\Omega }}})\\ & & \times \,{e}^{i2\pi {ft}}[1-{e}^{-i2\pi {fL}(1+\hat{{\rm{\Omega }}}\cdot \hat{p})}].\end{array}\end{eqnarray} \tag{ 3.8 }$$

The fractional frequency shift, z(t), produced by a stochastic GWB is simply given by integrating equation (3.5) over all directions. Using (3.1) and (3.8), we obtain:

$$\begin{eqnarray}&&z(t)=\displaystyle \int d\hat{{\rm{\Omega }}}\,z(t,\hat{{\rm{\Omega }}})\end{eqnarray} \tag{ 3.9 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & \displaystyle \sum _{A}{\displaystyle \int }_{-\infty }^{\infty }{df}{\displaystyle \int }_{{S}^{2}}d\hat{{\rm{\Omega }}}{F}^{A}(\hat{{\rm{\Omega }}}){h}_{A}(f,\hat{{\rm{\Omega }}}){e}^{i2\pi {ft}}\\ & & \times [1-{e}^{-i2\pi {fL}(1+\hat{{\rm{\Omega }}}\cdot \hat{p})}],\end{array}\end{eqnarray} \tag{ 3.10 }$$

where ${F}^{A}(\hat{{\rm{\Omega }}})$ are the antenna beam patterns for each polarization A, which we write as

$$\begin{eqnarray}&&{F}^{A}(\hat{{\rm{\Omega }}})=\left[\displaystyle \frac{1}{2}\displaystyle \frac{{\hat{p}}^{i}{\hat{p}}^{j}}{1+\hat{{\rm{\Omega }}}\cdot \hat{p}}\ {e}_{{ij}}^{A}(\hat{{\rm{\Omega }}})\right].\end{eqnarray} \tag{ 3.11 }$$

Searches for the GWB rely on looking for correlations induced by GWs in the timing residuals of pulsar pairs. Indeed, the observed quantity in PTA experiments is the timing residual $r(t)$, which is simply the integral of equation (3.9) in time:

$$\begin{eqnarray}&&r(t)={\int }^{t}{dt}^{\prime} z(t^{\prime} ).\end{eqnarray} \tag{ 3.12 }$$

The expected value of the correlation between a residual from pulsar 1 at time t_j, with that from a different pulsar, say pulsar 2 at time t_k, depends on terms of the form:

$$\begin{eqnarray}&&\langle {r}_{1}^{* }({t}_{j}){r}_{2}({t}_{k})\rangle =\left\langle {\int }^{{t}_{j}}{dt}^{\prime} {\int }^{{t}_{k}}{dt}^{\prime\prime} {z}_{1}^{* }(t^{\prime} ){z}_{2}(t^{\prime\prime} )\right\rangle ,\end{eqnarray} \tag{ 3.13 }$$

$$\begin{eqnarray}&&\langle {r}_{1}^{* }({t}_{j}){r}_{2}({t}_{k})\rangle ={\int }^{{t}_{j}}{dt}^{\prime} {\int }^{{t}_{k}}{dt}^{\prime\prime} {\int }_{-\infty }^{+\infty }{{dfe}}^{-i2\pi f(t^{\prime} -t^{\prime\prime} )}H(f){\rm{\Gamma }}({{fL}}_{1},{{fL}}_{2},\zeta ),\end{eqnarray} \tag{ 3.14 }$$

where H(f) contains the information of the spectrum of radiation. In analogy with [22, 23, 36], we define the quantity above that depends on the angular separation of the pulsars, ζ, their distances from the Earth, L₁, L₂, and the GW frequency f, as the overlap reduction function

$$\begin{eqnarray}&&{\rm{\Gamma }}({{fL}}_{1},{{fL}}_{2},\zeta )\equiv \int d\hat{{\rm{\Omega }}}\,\kappa (f,\hat{{\rm{\Omega }}})\left[\displaystyle \sum _{A}{F}_{1}^{A}(\hat{{\rm{\Omega }}}){F}_{2}^{A}(\hat{{\rm{\Omega }}})\right],\end{eqnarray} \tag{ 3.15 }$$

where

$$\begin{eqnarray}&&\kappa ({{fL}}_{\mathrm{1,2}},\hat{{\rm{\Omega }}})\equiv [1-{e}^{i2\pi {{fL}}_{1}(1+\hat{{\rm{\Omega }}}\cdot {\hat{p}}_{1})}][1-{e}^{-i2\pi {{fL}}_{2}(1+\hat{{\rm{\Omega }}}\cdot {\hat{p}}_{2})}].\end{eqnarray} \tag{ 3.16 }$$

In order to write a closed-form, analytic solution to (3.15), we choose a reference frame where one pulsar is placed along the z-axis and the other in the x-z plane as seen in figure 2. Specifically, we write

$$\begin{eqnarray}&&{\hat{p}}_{1}=(0,0,1),\end{eqnarray} \tag{ 3.17a }$$

$$\begin{eqnarray}&&{\hat{p}}_{2}=(\sin \zeta ,0,\cos \zeta ),\end{eqnarray} \tag{ 3.17b }$$

$$\begin{eqnarray}&&\hat{{\rm{\Omega }}}=(\sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \theta ),\end{eqnarray} \tag{ 3.17c }$$

$$\begin{eqnarray}&&\hat{m}=(\sin \phi ,-\cos \phi ,0),\end{eqnarray} \tag{ 3.17d }$$

$$\begin{eqnarray}&&\hat{n}=(\cos \theta \cos \phi ,\cos \theta \sin \phi ,-\sin \theta ).\end{eqnarray} \tag{ 3.17e }$$

We remind the reader that ${\hat{p}}_{1}$ and ${\hat{p}}_{2}$ are the unit vectors pointing to pulsars 1 and 2, respectively, $\hat{{\rm{\Omega }}}$ is the direction of GW propagation and $\hat{m}$ and $\hat{n}$ are the GW principal axes, see figure 2. Note that in this reference frame ${F}_{a}^{\times }=0$ by (3.11), making it a convenient choice.

For an isotropic GWB one is free to choose whichever coordinate system is most convenient, as was done here. However, one must be more careful when considering reference frames which are used to describe pulsar locations in an anisotropic GWB, as was done by [22].

We choose the coordinate system defined in equation (3.17), and apply it to equation (3.15). The result is equations (4.1) and (4.2).

Claim. Let ${L}_{1},{L}_{2},f$, be real positive constants. Then, for each $\zeta \in (0,\pi ]$, as ${{fL}}_{1}\to \infty$ and ${{fL}}_{2}\to \infty$, we have

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\pi }d\theta {\displaystyle \int }_{0}^{2\pi }d\phi \,{K}_{2}(\zeta ,\phi ,\theta )[1-{e}^{2\pi {\rm{i}}{{fL}}_{1}(1+\cos \theta )}]\\ \ \ \ \times \,[1-{e}^{-2\pi {\rm{i}}{{fL}}_{2}(1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi )}]\,\longrightarrow \,{\displaystyle \int }_{0}^{\pi }d\theta {\displaystyle \int }_{0}^{2\pi }d\phi {K}_{2}(\zeta ,\phi ,\theta )\end{array}\end{eqnarray} \tag{ 4.1 }$$

except when $\zeta =0$ and ${L}_{1}={L}_{2}$, a case covered in [24]. Here,

$$\begin{eqnarray}&&{K}_{2}(\zeta ,\phi ,\theta )=\displaystyle \frac{\sin \theta (1-\cos \theta )\ ({\sin }^{2}\phi {\sin }^{2}\zeta -{(\sin \zeta \cos \theta \cos \phi -\sin \theta \cos \zeta )}^{2})}{1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi }.\end{eqnarray} \tag{ 4.2 }$$

Note that the above integrals are now written in terms of the coordinate system constructed in equation (3.17), and illustrated in figure 1, which was applied to (3.15) and (3.16).

Until now, one was only able to show this result by picking some values of pulsar distance 1, L₁, and pulsar distance 2, L₂ and solve (4.1) numerically assuming some GW frequency f. In the literature, e.g. [22, 24], the authors invoke the reader's physical intuition to support the numerical result—that if the exponents in (4.1) are large, ${fL}\gg 1$, these oscillatory pieces rapidly converge to zero. This is often referred to as the 'short wavelength approximation', and has been used without proof, which we will now provide.

4.1. Proof of claim

To prove this result, we estimate each of the four integrals (4.4)–(4.8) below which make up (4.1) separately. We apply the Lebesgue Dominated Convergence theorem (see appendix), Fubini's Theorem, and the two-dimensional Divergence theorem to get the required limiting value (4.1). A key result used in the proofs which follow is a variant of the Riemann-Lebesgue Lemma in harmonic analysis (see [43], p. 277 and [44], p.2): let a, b be finite (though this is not necessary). Then for a Lebesgue integrable function f (a comprehensive definition and examples of this are given in appendix) over [a, b],

$$\begin{eqnarray}&&{\int }_{a}^{b}{dt}\ {e}^{{itx}}f(t)\to 0\hspace{0.5cm}\mathrm{as}\hspace{0.5cm}x\to \infty .\end{eqnarray} \tag{ 4.3 }$$

The aforementioned Dominated Convergence Theorem, equation (A.1), basically gives us conditions under which we can interchange the operation of taking the limit of an integral with the integral of the limit.

First, we show that K₂(ζ, ϕ, θ) can be made continuous—and so absolutely integrable over its domain of definition—for all values of θ ∈ [0, π], ζ ∈ (0, π], and ϕ ∈ [0, 2π].

We use the identity

$$\begin{eqnarray*}&&\begin{array}{l}\quad 1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi \\ =1+\cos (\theta +\zeta )+\sin \theta \,\sin \zeta (1+\cos \phi ),\end{array}\end{eqnarray*}$$

to show that the denominator of (4.2), i.e.,

$$\begin{eqnarray*}&&1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi \geqslant 0\end{eqnarray*}$$

for all θ ∈ [0, π], ζ ∈ (0, π], and ϕ ∈ [0, 2π]. It follows that the only singularities of K₂ must occur when the denominator vanishes, and this occurs precisely when $1+\cos (\theta +\zeta )=0$ and $\sin \theta \,\sin \zeta (1+\cos \phi )=0$, since both these quantities are necessarily non-negative. This, in turn, implies that for given ζ, θ = π − ζ and ϕ = π or ζ = 0, ϕ any, or ζ = π, ϕ any. Each of these cases is handled by limiting arguments.

For example, we note that

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{\zeta \to {0}^{+}}{K}_{2}(\zeta ,\phi ,\pi -\zeta )=0,\quad \mathop{\mathrm{lim}}\limits_{\zeta \to {\pi }^{-}}{K}_{2}(\zeta ,\phi ,\pi -\zeta )=0,\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{\phi \to {\pi }^{-}}{K}_{2}(\zeta ,\phi ,\pi -\zeta )=\displaystyle \frac{2{\sin }^{3}\zeta }{1-\cos \zeta }.\end{eqnarray*}$$

The previous equation gives a zero limit as ζ → 0⁺ and is otherwise finite. It follows from this that K₂ can be defined to be a continuous function for any given value of ζ ∈ (0, π] and all values of θ ∈ [0, π] and ϕ ∈ [0, 2π]. Thus K₂ is Lebesgue integrable over the region [0, 2π] × [0, π]. Next,

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\pi }d\theta {\displaystyle \int }_{0}^{2\pi }d\phi \,{K}_{2}(\zeta ,\phi ,\theta )[1-{e}^{2\pi {{ifL}}_{1}(1+\cos \theta )}]\\ \ \ \times \,[1-{e}^{-2\pi {{ifL}}_{2}(1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi )}]=\end{array}\end{eqnarray} \tag{ 4.4 }$$

$$\begin{eqnarray}&&-\,{\int }_{0}^{\pi }d\theta {\int }_{0}^{2\pi }d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{2\pi {\rm{i}}{{fL}}_{1}(1+\cos \theta )}\end{eqnarray} \tag{ 4.5 }$$

$$\begin{eqnarray}&&-\,{\int }_{0}^{\pi }d\theta {\int }_{0}^{2\pi }d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{-2\pi {\rm{i}}{{fL}}_{2}(1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi )}\ d\phi \,d\theta \end{eqnarray} \tag{ 4.6 }$$

$$\begin{eqnarray}&&+\,{\int }_{0}^{\pi }d\theta {\int }_{0}^{2\pi }d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{2\pi {\rm{i}}({{fL}}_{1}(1+\cos \theta )-{{fL}}_{2}(1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi ))}\ \end{eqnarray} \tag{ 4.7 }$$

$$\begin{eqnarray}&&+\,{\int }_{0}^{\pi }d\theta {\int }_{0}^{2\pi }d\phi {K}_{2}(\zeta ,\phi ,\theta )\end{eqnarray} \tag{ 4.8 }$$

the last of which is identical to the required integral, (4.16). Note that each of the previous four integrals is necessarily finite since the region of integration is finite and K₂ is absolutely integrable over it.

Write $\lambda := {{fL}}_{1},\mu := {{fL}}_{2}$. Now we treat each of the previous three integrals (4.5)–(4.7) separately, and fix ζ ∈ (0, π].

4.2. Equation (4.5) tends to zero

Here we show that the first of the three equations with the exponential pulsar terms, (4.5), tends to zero. Using the above notation, and using the fact that K₂ is absolutely integrable over ${ \mathcal R }$, Fubini's theorem on the interchange of iterated integrals yields the equality,

$$\begin{eqnarray*}\begin{array}{rcl}{I}_{1}(\lambda ) & := & {\displaystyle \int }_{0}^{\pi }d\theta {\displaystyle \int }_{0}^{2\pi }d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{2\pi {\rm{i}}{{fL}}_{1}(1+\cos \theta )}\\ & = & {\displaystyle \int }_{0}^{2\pi }d\phi \left\{{\displaystyle \int }_{0}^{\pi }d\theta \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{2\pi {\rm{i}}\lambda (1+\cos \theta )}\right\},\end{array}\end{eqnarray*}$$

which, after the change of variable $u=1+\cos \theta$, gives us,

$$\begin{eqnarray*}&&{I}_{1}(\lambda )={\int }_{0}^{2\pi }d\phi \left\{{\int }_{0}^{2}{du}\,{K}_{2}^{* }(\zeta ,\phi ,u)\,{e}^{2\pi {\rm{i}}\lambda u}\right\},\end{eqnarray*}$$

where ${K}_{2}^{* }(\zeta ,\phi ,u)={K}_{2}(\zeta ,\phi ,\theta )/\sin \theta$ in the new variables is still absolutely integrable. Next, since ${K}_{2}^{* }$ is absolutely integrable over its domain, the ordinary two-dimensional version of the Riemann-Lebesgue Lemma, equation (4.3), implies that

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{\lambda \to \infty }{\int }_{0}^{2}{du}\,{K}_{2}^{* }(\zeta ,\phi ,u)\,{e}^{2\pi {\rm{i}}\lambda u}=0.\end{eqnarray*}$$

Since the previous integral is itself ${ \mathcal O }(| | {K}_{2}| {| }_{\infty })$, the Lebesgue Dominated Convergence theorem (A.1) can be used to interchange the order of the limit and the integral. We find:

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{\lambda \to \infty }{I}_{1}(\lambda )={\int }_{0}^{2\pi }d\phi \left\{\mathop{\mathrm{lim}}\limits_{\lambda \to \infty }{\int }_{0}^{2}{du}\,{K}_{2}^{* }(\zeta ,\phi ,u)\,{e}^{2\pi {\rm{i}}\lambda u}\right\}=0.\end{eqnarray*}$$

Thus, (4.5) tends to zero as $\lambda \to \infty$.

4.3. Equation (4.6) tends to zero

Preamble. In order to extend the previous idea to more general exponents, we apply integration by parts to double integrals via the Divergence Theorem. In order to prove either (4.6) or (4.7) it suffices that we obtain the decay estimates ${ \mathcal O }(1/\mu )$ or ${ \mathcal O }(1/\lambda )$ as $\mu \to \infty$ or $\lambda \to \infty$. What follows is the general idea which we then apply to the various cases. We need to estimate limits of the form

$$\begin{eqnarray}&&I(\omega ):= {\iint }_{{ \mathcal D }}\,d\theta \,d\phi \,f(\theta ,\phi ){e}^{i\omega g(\theta ,\phi )}={\iint }_{{ \mathcal D }}{dA}\,{{fe}}^{{\rm{i}}\omega g},\end{eqnarray} \tag{ 4.9 }$$

as $\omega \to \infty$. (Note that we used Fubini's theorem to justify the interchange of the order of integration in the iterated integral (4.6)). Here ${ \mathcal D }$ along with its boundary (or perimeter), ${ \mathcal C }$, are completely contained in ${ \mathcal R }$ and are chosen so that ${\rm{\nabla }}g(\theta ,\phi )\ne 0$ on and inside ${ \mathcal D }\cup { \mathcal C }$ (which necessarily has no points in common with ${ \mathcal R }$). By construction, the gradient of $g$, ${\rm{\nabla }}g$ does not vanish on ${ \mathcal D }\cup { \mathcal C }$ and therefore the quantity

$$\begin{eqnarray}&&{\bf{u}}=\displaystyle \frac{{\rm{\nabla }}g}{| {\rm{\nabla }}g{| }^{2}}\,f\end{eqnarray} \tag{ 4.10 }$$

is well-defined on ${ \mathcal D }\cup { \mathcal C }$.

We need to estimate the integral in (4.9) for large ω. First, observe that (suppressing the variables for clarity of exposition)

$$\begin{eqnarray*}&&{\rm{\nabla }}\cdot ({\bf{u}}\,{e}^{{\rm{i}}\omega g})=({\rm{\nabla }}{e}^{{\rm{i}}\omega g})\cdot {\bf{u}}\,+{e}^{{\rm{i}}\omega g}\,({\rm{\nabla }}\cdot {\bf{u}}),\end{eqnarray*}$$

where we assume, in addition, that f is sufficiently smooth so that ${\rm{\nabla }}\cdot {\bf{u}}$ is defined. Since ${\rm{\nabla }}g\cdot {\bf{u}}=f$ we have ${\rm{\nabla }}({e}^{{\rm{i}}\omega g})\cdot {\bf{u}}={\rm{i}}\omega {{fe}}^{{\rm{i}}\omega g},$ which when inserted into the previous display and integrated over ${ \mathcal D }$ yields,

$$\begin{eqnarray*}\begin{array}{rcl}{\iint }_{{ \mathcal D }}{dA}\,{\rm{\nabla }}\cdot ({\bf{u}}\,{e}^{{\rm{i}}\omega g})\ & = & {\iint }_{{ \mathcal D }}{dA}\,({\rm{\nabla }}\cdot {\bf{u}}){e}^{{\rm{i}}\omega g}+{\rm{i}}\omega {\iint }_{{ \mathcal D }}{dA}\,{{fe}}^{{\rm{i}}\omega g}\\ & = & {\iint }_{{ \mathcal D }}{dA}\,({\rm{\nabla }}\cdot {\bf{u}}){e}^{{\rm{i}}\omega g}+{\rm{i}}\omega I(\omega ).\end{array}\end{eqnarray*}$$

An application of the divergence theorem to the integral on the left gives us,

$$\begin{eqnarray*}&&{\iint }_{{ \mathcal D }}{dA}\,{\rm{\nabla }}\cdot ({\bf{u}}\,{e}^{{\rm{i}}\omega g})={\int }_{{ \mathcal C }}d\sigma \,({\bf{u}}\cdot {\bf{n}})\,{e}^{{\rm{i}}\omega g},\end{eqnarray*}$$

where ${\bf{n}}$ is the unit normal to ${ \mathcal C }$, itself oriented in the positive direction, and σ is arc length. Combining the two previous displays we get,

$$\begin{eqnarray}&&I(\omega )=-\displaystyle \frac{{\rm{i}}}{\omega }{\int }_{{ \mathcal C }}d\sigma \,({\bf{u}}\cdot {\bf{n}}){e}^{{\rm{i}}\omega g}+\displaystyle \frac{{\rm{i}}}{\omega }{\iint }_{{ \mathcal D }}{dA}\,({\rm{\nabla }}\cdot {\bf{u}}){e}^{{\rm{i}}\omega g}.\end{eqnarray} \tag{ 4.11 }$$

Once we know that both integrands are absolutely integrable over ${ \mathcal C }$ and ${ \mathcal D }$ respectively, we get $I(\omega )={ \mathcal O }(1/\omega )$ or $I(\omega )\to 0$ as $\omega \to \infty$, over ${ \mathcal D }\cup { \mathcal C }$. The results (4.6) or (4.7) are obtained by a careful limiting analysis of the case where ${ \mathcal D }\cup { \mathcal C }$ approaches ${ \mathcal R }$ which then gives us the desired decay estimate.

Proof, Case 1.  $\zeta \in (0,\pi )$. Set $g(\theta ,\phi )=-(1+\cos \theta \,\cos \zeta +\sin \theta \sin \zeta \cos \phi )$, $\omega =2\pi \,{{fL}}_{2}$ in (4.9). Then ${\rm{\nabla }}g(\theta ,\phi )=(\sin \theta \sin \zeta \cos \phi ,\sin \theta \cos \zeta -\cos \theta \sin \zeta \cos \phi ),$ so that ${\rm{\nabla }}g(\theta ,\phi )={\bf{0}}$ if and only if $\sin \theta \sin \phi =0$ and $\sin \theta \cos \zeta -\cos \theta \sin \zeta \cos \phi =0$. For $\zeta \in (0,\pi )$ this yields the eight (8) critical (or stationary) points

$$\begin{eqnarray}\begin{array}{rcl}(\theta ,\phi ) & = & \left(0,\displaystyle \frac{\pi }{2}\right),\left(0,\displaystyle \frac{3\pi }{2}\right),\left(\pi ,\displaystyle \frac{\pi }{2}\right),\left(\pi ,\displaystyle \frac{3\pi }{2}\right),\\ & & (\zeta ,0),(\pi -\zeta ,0),(\zeta ,2\pi ),(\pi -\zeta ,2\pi ),\end{array}\end{eqnarray} \tag{ 4.12 }$$

all of which are located on the perimeter of ${ \mathcal R }$. Since we want ${ \mathcal D }$ to be critical-point-free, for given $\varepsilon \gt 0$, choose

$$\begin{eqnarray*}&&{ \mathcal D }=\{(\theta ,\phi ):\varepsilon \lt \theta \lt \pi -\varepsilon ,\varepsilon \lt \phi \lt 2\pi -\varepsilon \},\end{eqnarray*}$$

and its perimeter,

$$\begin{eqnarray*}&&{ \mathcal C }=\{(\theta ,\phi ):\phi =\varepsilon ,\phi =2\pi -\varepsilon ,\theta =\varepsilon ,\theta =\pi -\varepsilon \}.\end{eqnarray*}$$

Then, by construction, ${\rm{\nabla }}g\ne 0$ in ${ \mathcal D }$ as well as on its perimeter, ${ \mathcal C }.$ Defining ${\bf{u}}$ as in (4.10) we then obtain (4.11) for suitably smooth functions $f,g$, i.e., $I(\omega )\to 0$ on ${ \mathcal D }\cup { \mathcal C }$, for every $\varepsilon \gt 0$. Now set $f={K}_{2}$ and note that both integrals in (4.11) are finite on their respective region of integration. Thus, for every $\varepsilon \gt 0$,

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{{{fL}}_{2}\to \infty }{\iint }_{{ \mathcal D }}\ d\theta \,d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}=0,\end{eqnarray*}$$

i.e., so taking the limit as ε approaches zero, we must have

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{\varepsilon \to 0}\mathop{\mathrm{lim}}\limits_{{{fL}}_{2}\to \infty }{\iint }_{{ \mathcal D }}\ d\theta \,d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}=0.\end{eqnarray} \tag{ 4.13 }$$

All that remains to be shown is that the interchange of the limits in the next expression is justified, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\mathrm{lim}}\limits_{\varepsilon \to 0}\mathop{\mathrm{lim}}\limits_{{{fL}}_{2}\to \infty }{\iint }_{{ \mathcal D }}\ d\theta \,d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}\\ \quad =\,\mathop{\mathrm{lim}}\limits_{{{fL}}_{2}\to \infty }\mathop{\mathrm{lim}}\limits_{\varepsilon \to 0}{\iint }_{{ \mathcal D }}\ d\theta \,d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )},\end{array}\end{eqnarray} \tag{ 4.14 }$$

as the right hand side of (4.14) is necessarily equal to (4.6) and so must vanish as well by (4.13), which is what we set out to prove. To this end, we note that, by continuity of the integrals,

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\mathrm{lim}}\limits_{\varepsilon \to 0}{\iint }_{{ \mathcal D }}d\theta \,d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}\\ \quad ={\iint }_{{ \mathcal R }}d\theta \,d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}\end{array}\end{eqnarray} \tag{ 4.15 }$$

and that, in fact, the convergence is uniform in ω, for $\omega \in [0,\infty )$.

Indeed, observe that

$$\begin{eqnarray*}&&\begin{array}{l}\left|{\iint }_{{ \mathcal R }}\ d\theta d\phi {K}_{2}(\zeta ,\phi ,\theta ){e}^{{\rm{i}}\omega g(\theta ,\phi )}-{\iint }_{{ \mathcal D }}\ d\theta d\phi {K}_{2}(\zeta ,\phi ,\theta )\,{e}^{{\rm{i}}\omega g(\theta ,\phi )}\right|\\ \quad =\,\left|{\iint }_{{ \mathcal R }\setminus { \mathcal D }}\ d\theta d\phi {K}_{2}(\zeta ,\phi ,\theta ){e}^{{\rm{i}}\omega g(\theta ,\phi )}\right|\\ \quad \leqslant {\iint }_{{ \mathcal R }\setminus { \mathcal D }}| \ d\theta \,d\phi \,{K}_{2}(\zeta ,\phi ,\theta )| \end{array}\end{eqnarray*}$$

and this last integral may be made arbitrarily small, independently of ω, if ε is sufficiently restricted. Hence the convergence in (4.15) is uniform in ω (actually for any $\omega \in {\bf{R}}$ but we only require this for ω on the half axis, $[0,\infty )$).

We are now in a position to apply a fundamental theorem on the interchange of such limits (see citepf, p. 395) to validate the equality in (4.14) and complete the proof in the case where $\zeta \in (0,\pi ).$

Case 2.  $\zeta =\pi$. In this case $g(\theta ,\phi )=-1-\cos \theta$ is independent of ϕ and the resulting double integral can be handled in a similar way as (4.5), the only difference being the presence of a negative sign in the exponent. This, however, causes no difficulty with the argument in that section and so we omit the details.

4.4. Equation (4.7) tends to zero

We use the same basic technique as in the proof of (4.6). The proof of the limiting result for (4.7) can be obtained by reducing it to the case of (4.6) just proved. For example, (4.7) may be rewritten in the form

$$\begin{eqnarray*}&&\begin{array}{l}{\displaystyle \int }_{0}^{\pi }d\theta {\displaystyle \int }_{0}^{2\pi }d\phi \{{K}_{2}(\zeta ,\phi ,\theta )\,{e}^{2\pi {\rm{i}}\lambda (1+\cos \theta )}\}\\ \quad \times \,{e}^{-2\pi {\rm{i}}\mu (1+\cos \theta \cos \zeta +\sin \theta \sin \zeta \cos \phi )}\,\end{array}\end{eqnarray*}$$

where now it is ${K}_{2}(\zeta ,\phi ,\theta )\,{e}^{t}2\pi {\rm{i}}\lambda \,(1+\cos \theta )$ that is absolutely integrable over ${ \mathcal R }$, since K₂ is and the exponential term has modulus equal to one. So, it follows from the methods above leading to (4.6) approaching zero as $\mu \to \infty$ that (4.7) also tends to zero as $\mu \to \infty$. Similarly, interchanging the μ and λ terms in the preceding integral we obtain that (4.7) tends to zero as $\lambda \to \infty$ as well.

4.5. Final result: the Hellings and Downs curve

We have shown that for an isotropic GWB, and for ${{fL}}_{i}\to \infty$, that the pulsar terms tend to zero. We can now write down the final form of the overlap reduction function: the 'Hellings and Downs' curve [5]:

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\pi }d\theta \,{\displaystyle \int }_{0}^{2\pi }d\phi \,{K}_{2}(\zeta ,\phi ,\theta )\\ \quad =\,\displaystyle \frac{\sqrt{\pi }}{2}\left\{1+\displaystyle \frac{\cos \zeta }{3}+4(1-\cos \zeta )\mathrm{ln}\left(\displaystyle \frac{\sin \zeta }{2}\right)\right\}(1+{\delta }_{\mathrm{1,2}}),\end{array}\end{eqnarray} \tag{ 4.16 }$$

for $\zeta \in (0,\pi ]$ by [5, 22, 41]. Several comments should be made about equation (4.16) regarding a choice of normalization for the Hellings and Downs curve, the failure of the short-wavelength approximation, and the subsequent approximation of the pulsar term by a delta function for the autocorrelation.

When evaluating the autocorrelation, one can easily see that the value of the overlap reduction function is $4\sqrt{\pi }/3$. Indeed, it is a choice of normalization to set the autocorrelation equal to one when $\zeta =0$, which requires that equation (4.16) be multiplied by $3/(4\sqrt{\pi })$.

Next, one will note the $(1+{\delta }_{\mathrm{1,2}})$ term: this takes into account the failure of the short-wavelength approximation when evaluating the autocorrelation term. We approximate this is a delta function, however, in appendix C of [25], equation C4, it is shown that the exact solution is

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{auto}}(f,{L}_{1}={L}_{2})=\displaystyle \frac{1}{{(2\pi f)}^{2}}\left\{\displaystyle \frac{8\pi }{3}-\displaystyle \frac{1}{\pi {({fL})}^{2}}[1-{j}_{0}(4\pi {fL})]\right\},\end{eqnarray} \tag{ 4.17 }$$

where we adopt natural units (c = 1). Here ${j}_{0}(x)=\sin x/x$ is a spherical Bessel function of the first kind. Since the oscillatory piece is suppressed by a factor of $1/{({fL})}^{2}$, approximating the pulsar term by a multiplicative factor of 2 for the $\zeta =0$ autocorrelation case is justified.

Mingarelli and Sidery 2014 [24] also explored analytic expressions for the autocorrelation, and found that the $\kappa ({{fL}}_{\mathrm{1,2}},\hat{{\rm{\Omega }}})$ term, equation (3.16), can be well approximated as $2-2\cos M$, where $M=2\pi {fL}(1+\cos \theta )$. They showed numerically that the $2\cos M$ term contributed very little for large values of fL.

We note that for small values of ζ and fL, there is an intermediate regime where one requires more than $(1+{\delta }_{\mathrm{1,2}})$ to approximate $\kappa ({{fL}}_{\mathrm{1,2}},\hat{{\rm{\Omega }}})$. This case was explored in detail by Mingarelli and Sidery 2014 [24], who found that there are strong additional correlations from the pulsar term when the pulsars are separated by less than one gravitational wavelength. The authors also gave first and second order corrections for these cases.

We have shown analytically that the Hellings and Downs curve approaches the Earth-term only solution, even when the pulsars are arbitrarily distant from the Earth, and not themselves at the same distance L from the Earth. Of course, the case when ${{fL}}_{1}={{fL}}_{2}:= {fL}$ is easily recovered, since in this case, $\lambda =\mu$ and all terms (4.5)–(4.7) approach zero as the common value of this parameter fL approaches infinity.

The proofs indicate that the asymptotic estimate, equation (4.1), holds for sufficiently smooth kernels that are absolutely integrable over the region ${ \mathcal R }$, and not just kernels of the form ${K}_{2}(\zeta ,\phi ,\theta )$ as considered here.

The astrophysical interpretation of this result is that if one monitors any galactic millisecond pulsar, and cross-correlates it with a pulsar in e.g. the Large Magellanic Cloud, the Hellings and Downs curve would still be correct correlation function to use, under the assumption that the GWB is isotropic. Anisotropic GWBs can be handled similarly, but care is required when evaluating the new kernel.

To summarize, we have shown that for pulsars at distances L₁ and L₂ from the Earth, that the pulsar terms tend to zero as the ${{fL}}_{i}\to \infty$. The asymptotic estimate (4.1) is false if ${{fL}}_{2}$ is fixed as there is no reason for the integral (4.6) to tend to zero as ${{fL}}_{1}\to \infty$, since it is independent of ${{fL}}_{1}$. While this result is consistent with the previous intuition developed in the field of nanohertz GW astronomy, and indeed verified numerically for a few values of fL in [24], it has never before been proven analytically or generally for any ${{fL}}_{i}$. This result is an important validation of a fundamental result in the field, and lends credibility and rigor to current GWB searches as we enter detection era in nanohertz GW astronomy.

The authors thank Yacine Ali-Haïmoud for useful comments and a thorough reading of this manuscript. Figure 1 was generated in part with the online tool 'gwplotter'[32]. The Flatiron Institute is supported by the Simons Foundation.

Here we give a brief overview of the Lebesgue integral on the real line, ${\bf{R}}$, see ([45]), Chapter 10. Readers familiar with this concept may proceed immediately to the Main Results in section 3 and their proofs.

One of the great advantages of the Lebesgue integral over the classical (Riemann) integral is in the handling of limiting processes such as limits of integrals of sequences of functions which, in the classical case, usually requires uniform convergence of the sequence in question—this is not the case for the Lebesgue integral.

Fundamental to this now-standard integral in mathematics is the notion of Lebesgue measure. One intuitive notion of measure assigns the value b − a for the length of an interval $[a,b]$ or, more generally, $(b-a)(d-c)$, for the area of a plane rectangle ${ \mathcal R }=[a,b]\times [c,d]=\{(x,y):x\in [a,b]\ \ \mathrm{and}\ \ y\in [c,d]\}$ formed by the Cartesian product of the two intervals $[a,b]$ and $[c,d]$. The notion of measure allows one to extend this property (length, area, volume, etc) to Lebesgue measurable sets. Aside traditional examples such as rectangles, unions of disjoint intervals, etc one can consider the set ${\bf{Q}}$ of rational numbers on the real line, or its complement, the set of irrational numbers there. Both the latter two sets are Lebesgue measurable, the former has Lebesgue measure zero while the latter has positive Lebesgue measure.

More precisely, a set $E\subset {\bf{R}}$ is said to have Lebesgue measure zero if for any given $\varepsilon \gt 0$ there is a sequence of intervals $[{a}_{n},{b}_{n}]$, n = 1, 2, ... such that E is contained within the union of all these intervals where, in addition, ${\sum }_{n=1}^{\infty }({b}_{n}-{a}_{n})\lt \varepsilon$. For example, the set of all rational numbers ${\bf{Q}}\subset {\bf{R}}$ has measure zero.

One of the connections between the Lebesgue theory and the ordinary Riemann theory of the integral is the following result: If a function f is bounded on an interval $[a,b]$ then f is Riemann integrable over $[a,b]$ if and only if the set of its discontinuities is a set of measure zero. By its very definition, the Lebesgue integrability of f forces that the absolute value of the function, $| f|$, in question be Lebesgue integrable, which is not so in the case of the Riemann integral. Still, the advantage in using Lebesgue integrals is huge in that we can extend the class of functions and sets over which we are integrating and still get meaningful results.

Now we offer an extremely brief introduction to the Lebesgue integral, with a few key definitions and results. Briefly, with the Lebesgue integral we subdivide the range of a given function f, look for those parts where horizontal lines intersect the graph of f and then drop rectangles onto the domain of f. (Recall that we subdivide the domain of f in the case of the Riemann integral).

More generally, the measure of a set $E\subset {\bf{R}}$ is the greatest lower bound (see e.g. [45], p. 11) of the set of all numbers of the form ${\sum }_{n=1}^{\infty }({b}_{n}-{a}_{n})$ where the union of all the intervals $[{a}_{n},{b}_{n}]$ contains E. A set E is said to be Lebesgue measurable if it has finite measure. A function f is said to be measurable if the special set $F=\{x:f(x)\gt a\}$ is measurable for every real number a. These are precisely the functions that we can integrate in the Lebesgue sense, i.e., the Lebesgue integrable functions over E. First, the integral is defined for simple functions (i.e., those functions whose range is a finite set of points). Then, using the fact that for a given measurable function $f(x)\geqslant 0$ there is a monotonically increasing sequence of simple functions ${s}_{n}(x)$ that converges to $f(x)$ and whose integrals, ${\int }_{E}{s}_{n}(x){dx}$, are bounded by a constant C (that depends on f), we define the Lebesgue integral of f over E by

$$\begin{eqnarray*}&&{\int }_{E}f(x){dx}=\mathop{\mathrm{lim}}\limits_{n\to \infty }{\int }_{E}{s}_{n}(x){dx},\end{eqnarray*}$$

Finally, when f is a general measurable function (i.e., not necessarily non-negative) we define its integral using its decomposition into positive and negative parts, that is, we know that $f(x)={f}^{+}(x)-{f}^{-}(x)$ where ${f}^{\pm }(x)=\max \ \pm f(x),0\}$. Thus, the Lebesgue integral of a general measurable function f over E is by definition,

$$\begin{eqnarray*}&&{\int }_{E}f(x){dx}={\int }_{E}{f}^{+}(x){dx}-{\int }_{E}{f}^{-}(x){dx},\end{eqnarray*}$$

It is then an easy matter to see that

$$\begin{eqnarray*}&&{\int }_{E}| f(x)| \,{dx}={\int }_{E}{f}^{+}(x){dx}+{\int }_{E}{f}^{-}(x){dx}.\end{eqnarray*}$$

The Lebesgue Dominated Convergence theorem states that if ${f}_{n}(x)$ is a sequence of measurable functions with ${f}_{n}(x)\to f(x)$ almost everywhere on E. (This means that ${f}_{n}(x)\to f(x)$ at all points x except for a set of measure zero.) In addition, let $g\geqslant 0$ be Lebesgue integrable on E with $| {f}_{n}(x)| \leqslant g(x)$ then

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{n\to \infty }{\int }_{E}{f}_{n}(x){dx}={\int }_{E}f(x){dx}.\end{eqnarray} \tag{ A.1 }$$

This theorem is the main one being used in this paper in order to guarantee convergence of the various integrals. A similar theory and similar results hold in 2 or more dimensions.

The space ${L}^{1}[0,\pi ]$ is by definition the space of all complex valued integrable functions f such that

$$\begin{eqnarray*}&&{\int }_{0}^{\pi }| f(x)| \,{dx}\lt \infty .\end{eqnarray*}$$

In addition, the space ${L}^{\infty }({{\bf{R}}}^{+})$ is the space of all such functions f for which there exists a constant, C, depending on f, such that $| f(x)| \lt C$ almost everywhere on ${{\bf{R}}}^{+}$.