Gravitational Wave Searches with Pulsar Timing Arrays. I: Cancellation of Clock and Ephemeris Noises
Massimo Tinto





n
p
k

apsrev

 

mtinto@ucsd.edu
University of California San Diego,

  Center for Astrophysics and Space Sciences, 

  9500 Gilman Dr, La Jolla, CA 92093, 

  U.S.A.  


  We propose a data processing technique to cancel monopole and dipole
  noise sources (such as clock and ephemeris noises respectively) in
  pulsar timing array searches for gravitational radiation. These
  noises are the dominant sources of correlated timing fluctuations in
  the lower-part (≈ 10^-9 - 10^-8 Hz) of the
  gravitational wave band accessible by pulsar timing
  experiments. After deriving the expressions that reconstruct these
  noises from the timing data, we estimate the gravitational wave
  sensitivity of our proposed processing technique to single-source
  signals to be at least one order of magnitude higher than that
  achievable by directly processing the timing data from an equal-size
  array. Since arrays can generate pairs of clock and ephemeris-free
  timing combinations that are no longer affected by correlated
  noises, we implement with them the cross-correlation statistic to
  search for an isotropic stochastic gravitational wave background. We
  find the resulting optimal signal-to-noise ratio to be more than one
  order of magnitude larger than that obtainable by correlating pairs
  of timing data from arrays of equal size.


04.80.Nn, 95.55.Ym, 07.60.Ly
Gravitational Wave Searches with Pulsar Timing Arrays.
 I:
  Cancellation of Clock and Ephemeris Noises
    Massimo Tinto
    December 30, 2023
=======================================================================================================





§ INTRODUCTION



The first LIGO <cit.> detection of a gravitational wave (GW)
signal from two medium mass Binary Black Hole (BBH) merger
<cit.>, followed by the observation of mergers from further
BBHs <cit.> and one binary neutron star <cit.>, mark
the beginning of GW astronomy. Because of seismic noise and limited
arm length issues, the lower part (below 10 Hz) of the GW spectrum
will only be accessible by space-based detectors such as the European
LISA mission <cit.>, and pulsar timing arrays.

Pulsar timing GW experiments entail timing highly stable millisecond
pulsars in our own Galaxy. These experiments have been performed for
decades now, and they aim to detect gravitational radiation in the nHz
frequency band complementary to those accessible by ground and future
space-based GW detectors <cit.>. The basic principle
underlining the pulsar timing technique is the same as that of other
GW detector designs <cit.>: to monitor the
frequency variations of a coherent electromagnetic signal exchanged by
two or more “point particles” separated in space. As a pulsar
continuously emits a series of radio pulses that are received at
Earth, a GW passing across the pulsar-Earth link introduces
fluctuations in the time-of-arrival (TOA) of the received
electromagnetic pulses. By comparing the pulses TOAs against those
predicted by a model, it is in principle possible to detect the
effects induced by any time-variable gravitational fields present,
such as the transverse-traceless metric curvature of a passing plane
GW train <cit.>. 

The frequency band in which the pulsar timing technique is most
sensitive to ranges from about 10^-9 to 10^-7 Hz, with the
lower-limit essentially determined by the overall duration of the
experiment (10^-9  Hz≃ 1/ 30  years), and the
upper limit identified by the signal-to-noise ratio (SNR) of the
received radio pulses. To attempt observations of GWs in this way, it
is thus necessary to control, monitor and minimize the effects of
other sources of timing fluctuations, and, in the data analysis, to
use optimal algorithms based on the different characteristics of the
pulsar timing response to GWs (the signal) and to other sources of
timing fluctuations (the noise).

A quantitative analysis of the noise sources affecting millisecond
pulsar timing searches for gravitational radiation
<cit.> has shown that, in the
region (10^-9 - 10^-7) Hz of the accessible frequency band, they
are due to:


  * finiteness of the SNR in the raw observations, resulting in what
  is usually referred to as thermal noise at the receiver;

  * uncertainties in solar system ephemeris, which are used to
  correct TOAs at the Earth to the barycenter of the solar
  system;

  * variation of the index of refraction in the interstellar and
  interplanetary plasma; 

  * intrinsic rotational stability of the pulsar, and

  * instability of a combination of the local clock and of the
  International Time Standard against which pulsars are timed, and
  noise in time transfer if the clock is not located at the
  observatory site.

Although the timing fluctuations induced by some of these noises can
be in principle either reduced or calibrated out, the fundamental
noise-limiting sensitivity of pulsar timing experiments is imposed by
the timing-fluctuations inherent to the pulsar, the reference clocks
that control the TOA measurements, and the noise affecting the
ephemeris used for referring to the Solar System Barycenter (SSB) the
timing measurements performed at the observatory. The magnitudes of
these noises can be comparable to or larger <cit.>
than a GW stochastic background possibly present in the timing
data. For instance, clocks such as the Linear-Ion-Trapped-Standard
(LITS) (presently in-the-field state-of-the art atomic clock with
long-time-scale timing stability) <cit.> would result in a
sinusoidal strain sensitivity of pulsar timing searches for GWs to a
level of about 10^-15 after coherently integrating the data for a
period of 10 years <cit.>. Since the characteristic wave
amplitude associated with a super-massive black-hole-binaries
background is predicted to be of comparable magnitude at the frequency
3 × 10^-9 Hz
<cit.>, it is clear
that a single telescope will not be able to unambiguously detect such
a GW signal.

A method for statistically enhancing the SNR of pulsar timing
experiments to an isotropic background of GW was first proposed by
Hellings and Downs <cit.> and improved by Jenet et al.
<cit.>. This technique relies on cross-correlating pairs of
TOA residuals data taken with an array of highly-stable millisecond
pulsars.  Since the GW background is common to all timing residuals
while the noises affecting them may be uncorrelated, the
cross-correlation technique should enhance the strength of the GW
signal over that of the noises. In particular, it was shown that the
correlation of pulsar timing data as a function of their enclosed
angle has a characteristic signature that should enhance the
likelihood of detection by correlating timing data from a sufficiently
large ensemble of millisecond pulsars.

Although the correlation of the noise induced by the interstellar
medium and the intrinsic timing noise of each individual pulsar can be
disregarded since pulsars are generally widely separated on the sky,
some of the other above-mentioned noises are common to the array data
and result in non-zero correlations
<cit.>. This in turn may
prevent us from detecting the angular dependence of the Hellings-Downs
curve induced by a GW background. Among the various noise sources, the
instability of the International Time Standard against which pulsars
are timed and the noise associated with the SSB ephemeris are the most
important ones, showing strong correlations in the lower part of the
accessible frequency band <cit.>. Although data
processing methods have been discussed in the literature to mitigate
the effects of correlated noises
<cit.> in pulsar timing arrays
searches for a GW stochastic background, here we propose a method to
cancel them. 

This article is organized as follows. In Sec. <ref> we present
the mathematical formulation of the problem after deriving the
transfer functions of the clock- and ephemeris noises in the TOA
residuals from an array of pulsars and noticing that they are
different from that of a GW signal. From these considerations we then
show that an array of 5 pulsars allows us to: (i) identify a linear
combination of its TOA residuals that simultaneously cancels both
clock- and ephemeris noises; (ii) optimally reconstruct these
noises. In Sec. <ref> we then estimate the GW sensitivity of
clock- and ephemeris-free combinations to single-source GW
signals. Since single-source searches are limited by clock and
ephemeris noises in the lower part of the accessible frequency band,
we estimate the sensitivity enhancement of our method over
single-pulsar searches to be at least one order of magnitude in this
part of the band. In Sec. <ref> we then turn our attention to
searches for an isotropic stochastic GW background implemented by
cross-correlating pairs of clock- and ephemeris-free combinations
generated by arrays of 10 or more pulsars. As the remaining noises
in some pairs of combinations are now uncorrelated
<cit.>, we derive the expressions of
the variance-covariance matrix associated with the correlation
statistic built with pairs of clock- and ephemeris-free data
combinations. We estimate that an array of 10 pulsars, with its six
clock- and ephemeris-free combinations, can generate three pairs of
clock- and ephemeris-free combinations whose noises are uncorrelated
and achieve a sensitivity to an isotropic stochastic GW background
that is more than one order of magnitude better than that achievable
by cross-correlating TOA residuals from an equal size array.  In
Sec. <ref> we finally present our conclusions and considerations
about the technique we propose and the advantages it offers to the nHz
GW search efforts.



§ CLOCK AND EPHEMERIS NOISE-FREE COMBINATIONS



Let us consider an array of M pulsars timed by M observatories
[The reason for using M telescopes simultaneously tracking
  M pulsars (where simultaneity is of course relative to the period
  of the waves searched for) rather than, for instance, N telescopes
  tracking M pulsars (with N < M) in a “switching mode” (i.e.,
  alternating between a subset of pulsars) is that the
  M-antennas-M-pulsars scenario will clearly show how our
  processing technique works.]. The SSB ephemeris, which are estimated
at the Jet Propulsion Laboratory (JPL) through analysis of tracking
data from interplanetary spacecraft <cit.>, are used
to correct TOAs at the observatory to the SSB. The expression relating
the TOA at the SSB, t̂^(i)_SSB ,  i = 1, ... M to the TOA
at one of the M observatories, t^(i)_obs ,  i = 1, ... M,
can be written in the following form <cit.>

    t̂^(i)_SSB = t^(i)_obs - 
    r̂^(i) (t) ·^(i)/c + ν^(i) (t)  ,

where the Römer delay has been shown explicitly, the symbol ·
represents the operation of scalar product between two vectors, c is
the speed of light, and ν^(i) (t) includes all other
contributors to the difference between the Earth and SSB TOAs
<cit.>. In Eq. (<ref>),
r̂^(i) (t) and ^(i) denote the
position of the radio telescope i and sky location of pulsar i
w.r.t. the SSB respectively, and the symbol  on both
observables emphasizes that they are affected by errors.  By rewriting
them as sums of their “true” values,
(r^(i) , ^(i)), and their errors,
(e^(i) (t), Δ^(i)), Eq.(<ref>) can be
rewritten in the following form (in which now the speed of light has
been taken to be equal to 1)

    t̂^(i)_SSB≃ t^(i)_SSB - 
    e^(i) (t) ·^(i)
    - r^(i) (t) ·Δ^(i)  + ν^(i) (t)  ,

where a term quadratic in the errors has been disregarded. Note that
the error vector e^(i) (t) can be decomposed into the sum
of two error terms, e (t) + η^(i) (t), with the
first term representing the error of the center of the Earth relative
to the SSB and the second the error of the position of the observatory
i relative to the center of the Earth. Since the position of
observatory i with respect to the center of the Earth is known with
accuracy and precision that are orders of magnitude better than those
associated with the position of the center of the Earth relative to
the SSB <cit.>, in what follows our focus will be on the
SSB ephemeris used to convert TOAs from the Earth' s center to the
SSB. Under this assumption, Eq. (<ref>) can be written in
the following form

    t̂^(i)_SSB≃ t^(i)_SSB - e (t) ·^(i)
    - r^(i) (t) ·Δ^(i) + ν^(i) (t)  .

Note that the angular error associated with the sky location of the
pulsar w.r.t. SSB, Δ^(i), can be as large as a few
tens of mas in both RA and DEC. This inaccuracy results in a
sinusoidal timing error of period one year and amplitude a few
microseconds <cit.>. Although the magnitude of such an
error is much larger than the other errors affecting the TOA
residuals, its well-defined frequency allows us to remove it
from the timing data and should not be regarded as a limiting
factor.

Let now R^(i) (t)   , i = 1, ...M be the TOA residual measured at
time t by the M radio telescopes timing M pulsars
[Although TOA residuals from an array are generally sampled
  unevenly and at different times, by applying Fractional-Delay
  Filtering (FDF) <cit.> to the timing data it is possible
  to reconstruct, with an exquisitely high accuracy, data points from
  the surrounding samples. As an example application of its use, FDF
  is integral part of the data processing technique used by LISA to
  digitally suppress (more than seven orders of magnitude) the laser
  noise by properly time-shifting and linearly combining the
  heterodyne measurements.]. From the noise-considerations made
earlier, the M-residuals can be described by the following
expressions

    R^(i) (t)  =  H^(i) (t)  +   C (t) 
     -  e (t) ·^(i) +  ζ^(i) (t)   ,   i =  1, ..., M  ,

where the first term on the right-hand-side represents the
contribution from a possibly present GW signal
<cit.>, C (t) is a monopole random process
associated with noises affecting all timing residuals at time t
(with the clock being probably the dominant one),
e (t) ·^(i) is the ephemeris noise and
ζ^(i) (t) corresponds to the timing fluctuations due to all
other noise sources affecting the timing residual i. In what follows
we will assume the random processes ζ^(i) (t) to be of
zero-mean, and regard them as being uncorrelated to each other
<cit.>. 

Although the number of random processes we want to cancel are 4 in
total, namely (C (t), e (t)), the minimum number of TOA
residuals needed is actually 5. This is because the exact removal of
the vector random process e (t) alone requires four timing
data. Since three directions to three pulsars are in general linearly
independent, it is possible to identify a linear combination of four
TOA residuals in which the resulting Römer terms add up to zero.

Before proceeding with the identification of the linear combinations
that simultaneously cancel the clock and ephemeris noises, we first
introduce an orthonormal basis (a⃗_1, a⃗_2, a⃗_3)
centered on the SSB, and denote with
n̂^(i)_j  ,  i=1,… M  ,  j = 1, 2, 3 the components
of the unit vectors associated with the directions to the pulsars in
this coordinate system.

To identify the linear combinations that simultaneously cancel the
clock and ephemeris noises, we rewrite Eq.(<ref>) in the
following matrix form

    [ R^(1); R^(2);     ⋮; R^(M);       ]
    =
    [    1 ^(1);    1 ^(2);    ⋮;    1 ^(M) ][    C; -  e ]
    +
    [ H^(1); H^(2);     ⋮; H^(M) ]
    +
    [   ζ^(1);   ζ^(2);       ⋮; ζ^(M) . ]

Since four or more unit vectors associated with the directions to the
pulsars are linearly dependent, we will assume three of them, say
(^(1), ^(2), ^(3)), to be
linearly independent and use them as a new basis. This means that the
remaining M - 3 unit vectors can be written as linear combinations
of them as follows

    ^(i) = ∑_j=1^3 α^(i)_j ^(j)  ,    i = 4,… M  ,

where the (M-3) × 3 matrix elements α^(i)_j are 
equal to

    α^(i)_j = ∑_k=1^3 n̂^(i)_k (𝒩^-1)^k_j   , i =
    4,… M   , j = 1, 2, 3  ,

and the 3 × 3 matrix 𝒩 is given in terms of the
components of the three vectors (^(1),
^(2), ^(3)) by the following expression

    𝒩
    =
    [    n̂^(1)_1    n̂^(1)_2    n̂^(1)_3;    n̂^(2)_1    n̂^(2)_2    n̂^(2)_3;    n̂^(3)_1    n̂^(3)_2 n̂^(3)_3  . ]

After substituting Eq. (<ref>) into the matrix
multiplying the vector (C, - e)^T in Eq.(<ref>), the
problem becomes one of finding the generators of the Kernel
<cit.> of this matrix. This means finding the vectors
λ⃗≡ (λ_1, λ_2, …, λ_M) that, once
applied to the left of both sides of Eq.(<ref>), satisfy the
following homogeneous linear system

    (λ_1, λ_2, …, λ_M)
    [                      1                   ^(1);                      1                   ^(2);                      1                   ^(3);                      ⋮;                      1   ∑_j=1^3 α^(i)_j ^(j);                      ⋮;                      1 ∑_j=1^3 α^(M)_j ^(j) . ]
    = 0

The above equation translates in a corresponding homogeneous linear
system after noticing that a linear combination of three linearly
independent unit vectors is equal to zero iff the coefficients
of the combination are identically null. After some algebra we obtain
the following homogeneous linear system of 4 equations in M
unknowns

    [        1        1        1        …        …        1;        1        0        0  α^(4)_1  α^(5)_1 …α^(M)_1;        0        1        0  α^(4)_2  α^(5)_2 …α^(M)_2;        0        0        1  α^(4)_3  α^(5)_3 …α^(M)_3;          ][ λ_1; λ_2;   ⋮; λ_M ]
    = 0   .

To find the generators of the null-space of the above 4 × M
matrix we have relied on the software Mathematica
<cit.>. We have verified, for instance, that when M = 4
the Kernel is empty and the Image space <cit.> has dimensionality
4. This means that with 4 pulsars we can estimate the vector
(C , - e)^T. To this end, if we treat again the three unit
vectors (^(1), ^(2), ^(3)) as
basis, Eq.(<ref>) can be rewritten in the following form

    [ R^(1); R^(2); R^(3); R^(4) ]
    =
    [       1       1       0       0;       1       0       1       0;       1       0       0       1;       1 α^(4)_1 α^(4)_2 α^(4)_3 ][         C; -  e·^(1); -  e·^(2); -  e·^(3) ]
    +
    [ H^(1); H^(2); H^(3); H^(4) ]
    +
    [ ζ^(1); ζ^(2); ζ^(3); ζ^(4);       ] .

The inverse of the matrix multiplying the vector
(C, -  e·^(1) , -  e·^(2) , -  e·^(3))^T
in Eq. (<ref>) can easily be derived and its expression is
equal to

    A ≡1/α^(4)_1 + α^(4)_2 + α^(4)_3 -
      1 [                   α^(4)_1                   α^(4)_2                   α^(4)_3                       - 1;      α^(4)_2 + α^(4)_3 -1                 - α^(4)_2                 - α^(4)_3                         1;                 - α^(4)_1     α^(4)_1 + α^(4)_3 - 1                 - α^(4)_3                         1;                 - α^(4)_1                 - α^(4)_2 α^(4)_1 +
        α^(4)_2 - 1                         1 ] .

By left-applying the matrix A to both sides of Eq.(<ref>),
we obtain the following estimate for the clock and ephemeris noises

    A  
    [ R^(1); R^(2); R^(3); R^(4) ]
    =
    [         C; -  e·^(1); -  e·^(2); -  e·^(3) ]
    +
    A 
    [ H^(1); H^(2); H^(3); H^(4) ]
    + A 
    [ ζ^(1); ζ^(2); ζ^(3); ζ^(4) ] .

Finally, the expressions for the components of the ephemeris noise
w.r.t. the orthonormal basis (a⃗_1, a⃗_2, a⃗_3)
can be obtained by applying the matrix 𝒩^-1 to the left
of the reconstructed components of the ephemeris given by
Eq.(<ref>).

The configurations with M > 4 have a Kernel of dimension equal to
M-4 and an Image that is 4-dimensional. This means that with
M > 4 pulsars it is possible to reconstruct the clock and ephemeris
noises in four different and independent ways. Although the noises
affecting these reconstructions may be correlated, there should exist
an optimal combination of the reconstructed clock and
ephemeris noises that improves upon each individual reconstruction. We
will analyze the problem of optimally reconstructing the clock and
ephemeris noises from linearly combining pulsar timing data in a
future publication as the focus of this article is on searches for GWs
with clock- and ephemeris-free combinations of timing data.

In the specific case M = 5, the components of the single generator
of the Kernel, λ⃗, are equal to the following expressions
(up to a scaling constant chosen to be equal to 1)

    λ_1     =     (α^(5)_2 + α^(5)_3 - 1) α^(4)_1 + (1 - α^(4)_2 -
                    α^(4)_3) α^(5)_1  ,
    
    λ_2     =     (α^(5)_3 + α^(5)_1 - 1) α^(4)_2 + (1 - α^(4)_3 -
                    α^(4)_1) α^(5)_2  ,
    
    λ_3     =     (α^(5)_1 + α^(5)_2 -1) α^(4)_3 + (1 - α^(4)_1 -
                    α^(4)_2) α^(5)_3  ,
    
    λ_4     =     (1 - α^(5)_1 - α^(5)_2 - α^(5)_3)  ,
    
    λ_5     =     (α^(4)_1 + α^(4)_2 + α^(4)_3 - 1)  .

To exemplify the efficacy of the time-series
I(t) ≡∑_i=1^5 λ_i R^(i) (t) in canceling clock
and ephemeris noises with five pulsars, we have numerically
synthesized five timing residuals R^(r) (t)  ,  r=1,...5 by
generating clock and ephemeris noises through the use of a Gaussian
random number generator. Both noises have been assumed to have
zero-mean and a root-mean-squared (r.m.s.) error equal to 100 ns and
100/√(3) ns for each of the three components of the ephemeris
noise respectively. In addition, we have added a GW signal
characterized by two sinusoidal polarization components with strain
amplitudes each equal to 6.0 × 10^-16 and frequency equal to
3.3 × 10^-9 Hz. The wave's propagation direction and the
directions to the five pulsars were randomly selected, while the
contribution from all other noises (denoted
ζ^(r)   , r=1,...5 in the above equations) affecting the
measured timing residuals were not included to visually emphasize the
clock and ephemeris noise cancellation in I(t). Figure <ref>
shows the five residuals (insert (a) through (e)) simulated over a
period of 50 years and sampled every two weeks. The black-colored
lines represent the response of each timing residual to the above GW
signal; insert (f) shows the effectiveness of the noise-canceling
algorithm by synthesizing the TOA residual combination I(t).


Although our analysis treats the pulsars' sky locations as constants,
it is clear that the linear combination I(t) works with pulsars that
may move across the sky. I(t) should in fact be regarded as an
example application of the more general data processing technique
called Time-Delay Interferometry (TDI)
<cit.>. By properly time-shifting and linearly
combining data measured by a network of GW detectors, TDI provides the
mathematical framework for deriving new time-series that are
unaffected by correlated noises while retaining sensitivity to GWs.



§ SENSITIVITY TO SINGLE-SOURCE SIGNALS
 


To quantify the advantages brought by our data processing technique,
we will first derive the expression of the sensitivity to individual
GW signals when M=5, and compare it against that of a single pulsar.
To this end, let us first rewrite the expression of I(t) in the
following form

    I(t) ≡ H(t) + N(t) = ∑_r=1^5 λ_r H^(r) + ∑_r=1^5 λ_r ζ^(r) (t)  ,

where H(t) corresponds to the first term on the right-hand-side (the
GW signal in I), and N(t) to the second (the noise in I). To
estimate the sensitivity of I to individual GW signals, we first
derive the expression of the GW signal power averaged over sources
randomly distributed on the sky and polarization states

    ⟨ |H (f)|^2 ⟩ =
      ∑_i,j=1^5 λ_i λ_j 
      ⟨H^(i) (f) H^*(j) (f) ⟩ = ∑_i,j=1^5 λ_i λ_j Θ^ij h^2(f)  ,
 
where the angle-brackets denote the averaging operation over the
celestial sphere and wave's polarization states, 
represents the operation of Fourier transform and ^* that of complex
conjugation. In Eq. (<ref>) Θ^ij is the Hellings-Downs
<cit.> correlation function of the angle enclosed by the
directions to the two pulsars (i,j), and
h^2 (f) ≡⟨ |H^(i) (f)|^2 ⟩ , i=1,...5,
is the averaged GW power in each timing residual. This can be written
as h^2(f) = ρ^2(f)  h^2_0(f), where h_0(f) is the Fourier
transform of the wave amplitude and ρ(f) is the resulting wave's
r.m.s. transfer function to the pulsar response
<cit.>. 

Since the noises ζ^(r) , r=1,...5 are uncorrelated and each
can be characterized by its own one-sided power spectral density,
P_ζ^(r)(f), the expression of the one-sided power spectral
density of the noise N, P_N (f), is equal to

    P_N (f) = ∑_r=1^5 λ^2_r P_ζ^(r)(f)  .

The GW sensitivity of the I(t) combination, γ (f), defined as
the ratio between the square-root of its noise spectrum, P_N(f), and
the r.m.s. transfer function of its GW response,
√(∑_i,j=1^5 λ_i λ_j Θ^ijρ^2 ), is
then equal to <cit.>

    γ (f) = √(∑_r=1^5 λ^2_r P_ζ^(r)(f)/∑_i,j=1^5 λ_i λ_j Θ^ijρ^2(f)) .


To get some insights about γ (f), let us consider the case of
timing residual noises ζ^(r) being characterized by the same
spectrum, i.e.
P_ζ (f) ≡ P_ζ^(r)(f)    , r=1,...5. Under this
assumption Eq.(<ref>) assumes the following form

    γ (f) = 
    √(∑_r=1^5 λ^2_r/∑_i,j=1^5 λ_i λ_j Θ^ij)√(P_ζ (f))/ρ(f) .

If we denote with (Θ_min , Θ_max) the minimum and
maximum values of the Hellings-Downs curve <cit.> when
i ≠ j, and use the identity
∑_r=1^5 λ^2_r = - ∑_i ≠ j = 1^5 λ_i
λ_j
(which follows from the first equation fulfilled by λ⃗ in
Eq.(<ref>)), from Eq. (<ref>) it is then
possible to derive the following inequality

    γ_max > γ (f)  > γ_min  ,  γ_(min,max)≡1/√(1 - Θ_(min,max)) √(P_ζ
        (f))/ρ(f) .

If we now multiply and divide the right-hand-side of
Eq.(<ref>) by the square-root of the one-sided power
spectral density of the noise of each timing residual, P_R(f), we
derive the following upper-limit for the function γ(f)

    γ (f)  < γ_max = 1/√(1 - Θ_max)(√(P_ζ(f)/P_R(f))) (√(P_R(f))/ρ(f))  .

Since Θ_min≃ -0.15 and Θ_max≃ 0.5
<cit.>, the factor that determines the sensitivity gain of the
clock- and ephemeris-free combination I(t) over that of a single
pulsar is the ratio √(P_ζ(f)/P_R(f)). This function
of the Fourier frequency can be significantly smaller than one in the
lower part of the band where clock and ephemeris noises dominate the
noise budget. Based on the noise-model discussed in <cit.>, in
Fig. <ref> we plot the estimated GW sensitivity of single-pulsar
experiments together with the sensitivity bounds achievable by the
noise-canceling combination I(t) given in Eq. (<ref>).
In this figure we have assumed (i) use of multiple-frequency
measurements to adequately calibrate timing fluctuations from
intergalactic and interplanetary plasma, and (ii) disregarded the
timing fluctuations due to the pulsars.  Stability analysis of known
millisecond pulsars <cit.> have shown that there
exist some displaying frequency stabilities superior to those of the
most stable operational clocks in the (10^-9 - 10^-8) Hz
frequency band.


As a final remark, arrays with M ≥ 5 pulsars are characterized by
Kernels of dimensionality equal to M-4 (as they have M-4
generators) and four-dimensional Images. When M > 5 in particular,
the sensitivity of single-source searches can be coherently improved
over that with five pulsars by diagonalizing the correlation matrix
associated with the resulting M-4 combinations that are clock- and
ephemeris- free. The anticipated sensitivity enhancement is somewhat
larger than √(M-4), and we refer the reader to
<cit.> for details.



§ SENSITIVITY TO AN ISOTROPIC STOCHASTIC GW BACKGROUND



The technique presented in this article for canceling correlated
noises (such as clock and ephemeris) affecting data from pulsar timing
arrays becomes particularly effective when searching for a stochastic
background of gravitational radiation with an array of at least 10
pulsars. This is because arrays with 10 or more pulsars can generate
pairs of clock- and ephemeris-free combinations whose noises are
expected to be uncorrelated by not sharing data from the same
pulsars. Since the kernel of clock- and ephemeris- free combinations
generated by M pulsars has dimensionality equal to M-4, we need to
determine the number of pairs of combinations that can be constructed
with the M-4 generators of the kernel and whose noises are
uncorrelated.

To better understand the problem, let us first consider the case of
10 pulsars. The associated kernel space is defined by 6
generators, i.e. any clock- and ephemeris-free combination can be
written as a linear combination of them. If we label the pulsars as
(1, 2, … 10), we can choose, for instance, the generators
constructed by combining the following six set of 5 pulsars
(1,2,3,4,5),(6,7,8,9,10),(2,3,4,5,6),(7,8,9,10,1),
(3,4,5,6,7),(8,9,10,1,2),
which were obtained by “circularly-right-shifting” to the right the
(1, 2, … 10) indices.  With these generators we can only form
the following three pairs of combinations whose noises are
uncorrelated, [(1,2,3,4,5),(6,7,8,9,10)] ;
[(2,3,4,5,6),(7,8,9,10,1)] ; [(3,4,5,6,7),(8,9,10,1,2)], and
implement with them the correlation statistic to search for a GW
stochastic background.

Let us consider one of the above three pairs and denote with
(^(a)I(t) , ^(a)I̅(t)) the corresponding two clock-
and ephemeris-free combinations [Latin indices
  (a, b, c,… h) will be used to label pairs of clock- and
  ephemeris-free data combinations.]. Here ^(a)I̅(t)
contains TOA residuals that are not entering in the
^(a)I(t) combination and, from Eq. (<ref>),
both combinations can be written in the following forms

    ^(a)I(t)    ≡   ^(a)H (t) + ^(a)ζ (t) = ∑_j=1^5 ^(a)λ_j ^(a)H^(j) + ∑_j=1^5
                        ^(a)λ_j ^(a)ζ^(j) (t)   ,
    
    ^(a)I̅(t)    ≡   ^(a)H̅ (t) + ^(a)ζ̅ (t) = ∑_j=1^5
                           ^(a)λ̅_j ^(a)H̅^(j) + ∑_j=1^5
                           ^(a)λ̅_j ^(a)ζ̅^(j) (t)   .

Since the noises in these two TOA residual combinations are
uncorrelated, we can implement with them the cross-correlation
statistic <cit.> to search for a
stochastic GW background. To keep as close as possible to the
literature on the implementation of the correlation statistic for GW
searches, we will formulate it in the Fourier domain [
  Although TOA residuals are not sampled at even rates, use of
  fractional-delay filtering <cit.> allows us to resample
  the data by “interpolating” the needed samples at an exquisite level
  of accuracy.]. In what follows we will assume a GW background that is
an isotropic, unpolarized, stationary, and Gaussian random process
with zero mean. Such a background can be characterized by a one-sided
power spectral density, P_h (|f|), defined through the following
expressions <cit.>

    ⟨^(a)H^(i) (f)  ^(a)H^*(j) (f')
      ⟩   ≡   1/2δ (f - f') ^(a)Θ^ij P_h (|f|)  ,
    
    ⟨^(a)H̅^(i) (f)  ^(a)H̅^*(j) (f')
      ⟩   ≡   1/2δ (f - f') ^(a)Θ̅^ij P_h (|f|)  ,
    
    ⟨^(a)H^(i) (f)  ^(a)H̅^*(j) (f')
      ⟩   ≡   1/2δ (f - f') ^(a)Γ^ij P_h (|f|)  ,

where angle-brackets, ⟨⟩, now denote the operation of
ensemble averaging of the random process and averaging over the
celestial sphere and polarization states, and the three coefficients
(^(a)Θ^ij, ^(a)Θ̅^ij, ^(a)Γ^ij)
are the Hellings-Downs functions <cit.> associated with the
correlation of pairs of timing residuals in ^(a)I(t),
^(a)I̅ (t), and from both combinations respectively
[Here ^(a)Γ^ii is not equal to 1 as
  it corresponds to the correlation of two different pulsars]. Note
that the spectrum P_h (|f|) appearing in Eq. (<ref>)
depends on the spectrum of the GW background, Ω_gw (f), as
well as on the r.m.s. transfer function, ρ (f), of the GW
background to the TOA residual response in the following way
<cit.>

    P_h (|f|) = ρ^2 (f)  3 H_0^2/32 π^3 |f|^-3Ω_gw(f)  ,

where H_0 is the Hubble constant (today). The typical functional
form for the GW spectrum is a power law, i.e.
Ω_gw (f) ≡Ω_α f^α.  The spectral index
α characterizes the shape of the GW spectrum and is the unknown
to be determined by the filtering procedure described below.

We will further assume the random processes ^(a)ζ^(i),
^(a)ζ̅^(j) to be stationary, Gaussian distributed
with zero-mean, and characterized by the one-sided power spectral
densities,
(P_^(a)ζ^(j)(|f|)  ,  P_^(a)ζ̅^(j)(|f|)), defined as follows

    ⟨^(a)ζ^(i) (f)  ^(a)ζ^*(j) (f')
      ⟩   ≡   1/2δ^ij δ (f - f')  P_^(a)ζ^(i) (|f|) 
    
    ⟨^(a)ζ̅^(i) (f)  ^(a)ζ̅^*(j) (f')
      ⟩   ≡   1/2δ^ij δ (f - f')                      P_^(a)ζ̅^(i) (|f|) 
    
    ⟨^(a)ζ^(i) (f)  ^(a)ζ̅^*(j) (f')
      ⟩   ≡    0  ,

where δ^ij is the Kronecker symbol and δ (f - f') is
the Dirac's delta function. 

The sample cross-correlation statistic is defined as

    ^(a)𝒮≡∫_-T/2^T/2 dt  ∫_-∞^∞ dt'
     ^(a)I(t) ^(a)I̅(t')
    ^(a) Q(t-t')  ,

where ^(a) Q(t-t') is the optimal filter that is non-zero over the
time interval (-T/2, T/2) (a property we have already taken into account
in Eq.(<ref>) by extending the t' integration over the entire real
axis), and reflects the stationarity of both the stochastic background
and the instrumental noises by depending on time differences
<cit.>.

From the statistical properties of both the GW stochastic background
and the noises, and by virtue of the central-limit theorem, it follows
that ^(a)𝒮 is also a Gaussian random process
<cit.> that can be fully characterized by estimating its
mean, ^(a)μ_𝒮, and variance,
^(a)σ_𝒮^2.  Although their derivations are long
(see appendix (<ref>) for details), they are straightforward
and result in the following expressions

    ^(a)μ_𝒮≡T/2∫^+∞_-∞^(a)Λ(|f|)  ^(a)Q (f)  df 
    = T/2∫^+∞_-∞∑_r,s=1^5 ^(a)λ_r
      ^(a)λ̅_s ^(a)Γ^rs P_h (|f|)  ^(a)Q (f)
     df  ,


    ^(a)σ_𝒮^2    ≡   T/4∫^+∞_-∞Δ(|f|)  |^(a)Q (f)|^2  df
    
        =    T/4∫^+∞_-∞[ 
    ∑_r,s,p,q=1^5 [ ^(a)λ_r ^(a)λ_s ^(a)λ̅_p ^(a)λ̅_q  ^(a)Θ^rs ^(a)Θ̅^pq 
    + ^(a)λ_r ^(a)λ̅_s ^(a)λ_p ^(a)λ̅_q 
     ^(a)Γ^rs ^(a)Γ^pq ]  P^2_h (|f|) .
    
        +    . ∑_r,s,q=1^5 [ ^(a)λ_r ^(a)λ_s
          ^(a)λ̅^2_q  ^(a)Θ^rs P_^(a)ζ̅^(q)(|f|) +
    ^(a)λ̅_r ^(a)λ̅_s ^(a)λ^2_q
          ^(a)Θ̅^rs P_^(a)ζ^(q)(|f|) ]  P_h (|f|) .
    
        +    . ∑_r,s=1^5 ^(a)λ_r^2  ^(a)λ̅^2_s
              P_^(a)ζ^(r)(|f|)  P_^(a)ζ̅^(s)(|f|) ]  |^(a)Q (f)|^2  df  .

By simple inspection of Eq. (<ref>), and following
<cit.>, it is convenient to define the following
operation of inner product between two arbitrary complex functions,
say A, and B

    (A,B) ≡∫^+∞_-∞ A(f) B^* (f) Δ(|f|)  df  .

With this newly defined inner product, ^(a)μ_𝒮 and
^(a)σ_𝒮 can be
rewritten in the following forms

    ^(a)μ_𝒮 = T/2 (^(a)Q, ^(a)Λ/^(a)Δ)    ;     ^(a)σ^2_𝒮 = T/4 (^(a)Q (f) , ^(a)Q (f) )  ,

while the squared signal-to-noise ratio of ^(a)𝒮,
^(a) SNR^2, is equal to

    ^(a) SNR^2 ≡^(a)μ_𝒮^2/^(a)σ_𝒮^2 = T   (^(a)Q, ^(a)Λ/^(a)Δ)^2/(^(a)Q , ^(a)Q ) .

From a simple geometrical interpretation of the inner product defined
through Eq. (<ref>), from Eq. (<ref>) we conclude that the
maximum of the SNR of ^(a)𝒮 is achieved by choosing
the filter function ^(a)Q to be equal to

    ^(a)Q (f) = ^(a)ξ ^(a)Λ/^(a)Δ ,

where ^(a)ξ is an arbitrary real number and, from the
definition of the functions ^(a)Λ and ^(a)Δ (Eqs. <ref>,
<ref>), it also follows that ^(a)Q (f) is
real. We can take advantage of the arbitrariness of ^(a)ξ
by making all mean values of the cross-correlations equal to the
following constant μ_𝒮≡^(a)μ_𝒮 = Ω_α T  ,  a = 1, 2, 3
<cit.>. As it will become clearer later on in this
section, this choice simplifies the derivation of the coherent SNR
achievable by properly combining the cross-correlations of pairs of
clock- and ephemeris-free combinations.

The expression for the filter ^(a)Q given by Eq. (<ref>)
implies the following maximum SNR achievable by cross-correlating the
pair (^(a)I(t) , ^(a)I̅(t))

    ^(a) SNR_^(a)I≡√(T) √(∫^+∞_-∞^(a)Λ^2(|f|)/^(a)Δ(|f|) df ) .

In the limit in which the noise spectra of the TOA residuals are
larger than that of the GW background, the integrand of
Eq.(<ref>) becomes equal to

    ^(a)Λ^2(|f|)/^(a)Δ(|f|)≃( ∑_r,s=1^5 ^(a)λ_r ^(a)λ̅_s Γ^rs)^2 P^2_h (|f|)/( ∑_r=1^5 ^(a)λ^2_r P_^(a)ζ^(r)(|f|) ) ( ∑_s=1^5 ^(a)λ̅^2_s P_^(a)ζ̅^(s)(|f|)
     )  .

Note that the above expression of the optimal SNR depends on the
pulsars' relative sky locations through the Hellings and Downs
correlation function and the vectors,
^(a)λ⃗ , ^(a)λ̅⃗̅
identifying the pulsars' clock- and ephemeris-free combinations. To
quantify the angular dependence of the optimal SNR, we can assume the
timing residual noises
(^(a)ζ^(r), ^(a)ζ̅^(r) ,  r = 1,
… 5)
to be characterized by the same spectrum, P_ζ (|f|). Under this
assumption the optimal SNR becomes equal to

    ^(a) SNR_^(a)I≡√(T)  | ∑_r,s=1^5 ^(a)λ_r ^(a)λ̅_s
      Γ^rs | /√(∑_r,s=1^5 ^(a)λ^2_r
        ^(a)λ̅^2_s)√(∫^+∞_-∞P^2_h (|f|)/P_ζ^2 (|f|) df ) ,

where the dependence of the SNR on the pulsars' relative sky locations
can now be factored out of the integral. To quantify the magnitude of
the angular function appearing in Eq. (<ref>) we have randomly
generated 10^7 sets of 10 pulsars' sky locations and compute it
for each set. We found the above angular function to assume values
within the interval (0, 1) and to have an r.m.s. value equal to
about 1/5.

To compare the effectiveness of our cross-correlation statistic
against that based on a pair of TOA residuals,
(R^(i) (t) , R^(j) (t)), we provide below the expression for the
optimal SNR, SNR^ij_R, associated with the cross-correlation
statistic of two timing residuals [The expression of the
  optimal SNR of the cross-correlation statistic of a pair of timing
  residuals follows from their definitions (Eq. (<ref>)),
  and by performing a calculation similar to that presented in
  Appendix <ref>.]

    SNR^ij_R = √(T)√(∫^+∞_-∞ df   
    [ Γ^ij P_h(|f|) + P_C(|f|) + P_e^ij(|f|) ]^2
    /Π (|f|)) ,

where:

    Π (|f|)    ≡    2P^2_C(|f|) + P_e^ii (|f|) P_e^jj (|f|) +
      P^2_e^ij(|f|) + [P_e^ii (|f|)  + P_e^jj (|f|) ] P_ζ(|f|)  
    
        +     [P_e^ii (|f|)  + P_e^jj (|f|)  + 2
      P_e^ij (|f|) + 2 P_ζ (|f|) ] P_C(|f|) + P^2_ζ(|f|)  .

Eqs. (<ref>, <ref>) reflect the assumption on the
spectra of the ζ-noises to be equal to each other and larger
than the GW background, and where we have also denoted with
(P_C(|f|) , P_e^ij(|f|)) the spectra of the clock and ephemeris
noises respectively. Eqs. (<ref>, <ref>) allow us to
quantify two aspects of the degradation in the likelihood of detection
due to correlated-noises in pulsar timing data
<cit.>. First, the cross-correlation
statistic of pairs of TOA residuals from different pulsars is affected
by clock and ephemeris noises through their contribution to the mean
value of the cross-correlation (i.e. the numerator of the integrand in
Eq. (<ref>)). Since these noises are characterized by
relatively large spectral components in the same part of the band
associated with the presence of a stochastic GW background, they
result in an increased false-alarm probability . Second, clock
and ephemeris noises contribute to the overall noise variance of the
TOA residuals and therefore reduce (by a factor larger than 10
<cit.>) the optimal cross-correlation SNR
(Eq. (<ref>)) over that associated with pairs of clock and
ephemeris-free combinations (Eq. (<ref>)).

Following <cit.>, we now provide the expression for
the SNR achievable by optimally combining the three pairs of clock-
and ephemeris-free combinations that can be synthesized with the
timing data from an array of 10 pulsars. This requires the
calculation of the inverse of the variance-covariance matrix of the
cross-correlation statistic,
^(ab) C ≡⟨^(a)𝒮^(b)𝒮⟩ - ⟨^(a)𝒮⟩ ⟨^(b)𝒮⟩.
The discussion on how to derive ^(ab) C is presented in the
appendix, and its expression can be written in the following form

    ^(ab) C     =    T/4∫^+∞_-∞[ ∑_r,s,p,q=1^5 [ ^(a)λ_r ^(b)λ_s
          ^(a)λ̅_p ^(b)λ̅_q  
          ^(ab)Θ^rs^(ab)Θ̅^pq + 
    ^(a)λ_r ^(b)λ̅_s ^(a)λ̅_p ^(b)λ_q
    ^(ab)Γ^rs^(ab)Γ̅^pq ]  P^2_h (|f|)
    .
    
        +    . ∑_r,s,p,q=1^5  [ ^(a)λ_r ^(b)λ_s
          ^(a)λ̅_p ^(b)λ̅_q ^(ab)Θ^rsη^(ab)_pq ^(ab) P_ζ̅_p(|f|) +
    ^(a)λ_r ^(b)λ_s
          ^(a)λ̅_p ^(b)λ̅_q ^(ab)Θ̅^pqη^(ab)_rs ^(ab) P_ζ_r(|f|) 
    .
    
        +    . ^(a)λ_r  ^(b)λ̅_s ^(a)λ̅_p
          ^(b)λ_q 
    ^(ab)Γ^rsη^(ab)_pq ^(ab) P_ζ_p (|f|) +
    ^(a)λ_r  ^(b)λ̅_s ^(a)λ̅_p
          ^(b)λ_q 
    ^(ab)Γ^pqη^(ab)_rs ^(ab) P_ζ̅_r(|f|) ]  P_h (|f|) 
    .
    
        +    . ∑_r,s,p,q=1^5 ^(a)λ_r  ^(b)λ̅_s ^(a)λ̅_p
          ^(b)λ_q η^(ab)_rsη^(ab)_pq ^(ab)
          P_ζ̅_r(|f|)  ^(ab) P_ζ_p(|f|) ]  ^(a)Q (f) ^(b)Q^* (f)   df  .

In the above equation the multi-indices symbol η^(ab)_pq is
either equal to 0 or 1 depending on the particular outcome of the
correlation of two ζ noises entering in the combinations' pair
(ab), and it reduces to δ_pq when a=b.  The expression
for the optimal SNR achievable by combining the cross-correlations of
our three pairs of clock- and ephemeris-free combinations is equal to
<cit.>

    SNR_opt^2 = (^(1)μ_𝒮, ^(2)μ_𝒮, ^(3)μ_𝒮)
     C^-1 
    (^(1)μ_𝒮, ^(2)μ_𝒮, ^(3)μ_𝒮)^𝒯  = (Ω_α T)^2  ∑_a,b = 1^3^(ab) (C^-1)  ,

where the mean values ^(i)μ_𝒮 ,  i=1, 2, 3 have
been normalized to the same constant value Ω_α T.

The analysis for a set of 10 pulsars presented above can be
generalized to an arbitrary array of size M. This is done by first
selecting the M-4 generators of the kernel, and then identifying the
set of all pairs of generators that do not have pulsars in common,
i.e. pairs of clock- and ephemeris-free combinations whose noises are
uncorrelated. As an example of how to identify such set of generators'
pairs, let us consider an array with 45 pulsars, equal in number to
the array recently analyzed by the NANOGRAV consortium
<cit.>.  Let us first consider the numerical vector
(1, 2, … 45), which can be used to label the 45 pulsars. If we
use the “right-circular-shifting” procedure described earlier for
selecting a set of generators of the 10-pulsar kernel, we obtain the
following 41 generators

    (     1,… 5) , (6, … 10) … (41, … 45) 
    
        (     45, … 4) , (5, … 9) … (40, … 44) 
    
        (    44, … 3) , (4, … 8) … (39, … 43) 
    
        (     43, … 2) , (3, … 7) … (38, … 42) 
    
        (     42, … 1) , (2, … 6) , (7, … 11) , (12, … 16) , (17, … 21)

This set of 41 clock- and ephemeris-free combinations of timing
residuals results in a total of 672 pairs that do not share data
from the same pulsars, and can therefore be used to implement the
cross-correlation statistic[The total number of 672 pairs
  of clock- and ephemeris-free combinations generated by an array of
  45 pulsars was calculated numerically by using the program Mathematica <cit.>]. Since this number is comparable
to the number of pairs (990) of timing residuals given by an array
45 pulsars (and upon which the usual cross-correlation statistic is
built), it follows that the optimal SNRs associated with both
cross-correlation statistics will scale roughly by the same amount
over their respective single-pair SNRs <cit.>.



§ CONCLUSIONS



The data processing technique presented in this article allows us to
cancel clock and ephemeris noises affecting nHz GW pulsar timing
experiments. This is done by properly constructing linear combinations
of TOA residuals generated by arrays of millisecond pulsars.

We have found that searches for single-source GW signals will benefit
from this technique when implemented with arrays of at least 5
millisecond pulsars. The estimated sensitivity enhancement over that
from individual pulsar experiments is of at least one order of
magnitude in the lower-part of the accessible frequency band. This is
the frequency region where clock and ephemeris noises are leading
noise sources equally effecting the array's timing measurements.

Searches for an isotropic stochastic GW background can also be
performed with clock- and ephemeris-free combinations from an array of
10 or more pulsars. This is done by implementing the
cross-correlation statistic with pairs of clock- and ephemeris-free
combinations that do not share timing residuals from the same pulsars
to prevent noise correlations. With an array of 10 pulsars we have
found the associated cross-correlation statistic to be characterized
by an optimal SNR that is more than an order of magnitude larger than
the optimal SNR achievable by cross-correlating pairs of timing
residuals from an equal-size array.

As a final note, we have shown that clock and ephemeris noises can be
reconstructed with the timing data from an array of 4 or more
pulsars. In a future article we will estimate the accuracies by which
these observable can be reconstructed as functions of the number of
pulsars, their relative sky locations, and the magnitudes of the
remaining noises affecting the timing measurements.



§ ACKNOWLEDGMENTS


I would like to thank Dr. Stephen Taylor for the many useful
conversations during the development of this work, Dr. Lee Lindblom
for several comments during the early development phase of the idea
presented in this article, and Dr. Frank B. Estabrook and John
W. Armstrong for their constant encouragement.



§ DERIVATION OF THE MEAN AND VARIANCE-COVARIANCE OF 𝒮



In what follows we derive the expressions of the mean and the
variance-covariance matrix of the cross-correlation function
^(a) S (Eqs. (<ref>, <ref>, <ref>) given in
Section (<ref>)) for an array of 10 pulsars. We will assume
the noises ζ_r to be Gaussian-distributed with zero-mean, and to
have one-sided power-spectral densities P_ζ_r(|f|). We will
also take the GW stochastic background to be a Gaussian random process
of zero-mean, unpolarized, and stationary.  As a consequence of these
assumptions such a background is uniquely characterized by a one-sided
power spectral density, P_h (|f|) as given by Eqs.  (<ref>,
<ref>).

From the expressions of the clock- and ephemeris-free data
combinations, ^(a) I(t), their complements,
^(a)I̅ (t), (i.e. those combinations that do not include
timing data from the pulsars entering in ^(a) I(t)), and the
definition of their cross-correlation statistic,
^(a)𝒮, we have

    ^(a)𝒮≡∫_-T/2^T/2 dt ∫_-∞^∞ dt'  ^(a)I(t) ^(a)I̅(t')
      ^(a) Q(t-t')  ,

where ^(a) Q(t-t') is the optimal filter. Since this is non-zero
over the interval (-T/2, T/2), we have extended the t' integration
over the entire real axis. Equivalently, Eq. (<ref>) can be
rewritten in the Fourier domain as follows

    ^(a)𝒮 = ∫_-∞^∞ df  ∫_-∞^∞ df'
     δ_T (f - f')  ^(a)I^* (f) ^(a)I̅(f')
    ^(a)Q (f')  ,

where δ_T (f) ≡ T sinc(π f T) is the finite-time
approximation of the Dirac's delta function.

Since the noises in the combination ^(a) I(t) do not enter in
^(a)I̅ (t), and because a GW stochastic background is
uncorrelated with the measurement noises, we infer that the ensemble
average of the cross-correlation statistic, Eq. (<ref>), contains
contribution only from the GW stochastic background in the following
form

    ^(a)μ_𝒮≡⟨^(a)𝒮⟩    =    ∫_-∞^∞ df  ∫_-∞^∞ df'
     δ_T (f - f')  ⟨^(a)I^* (f) ^(a)I̅(f') ⟩^(a)Q (f') 
    
        =    ∫_-∞^∞ df  ∫_-∞^∞ df'
     δ_T (f - f')  ⟨^(a)H^* (f) ^(a)H̅(f') ⟩^(a)Q (f')  ,

where ^(a)H^* (f), ^(a)H̅(f')
are the contributions of the GW stochastic background to the
clock- and ephemeris-free combinations ^(a)I^* (f),
^(a)I̅(f') respectively. Since
^(a)H^* (f) = ∑_r=1^5 ^(a)λ_r
^(a)H^*r (f), and 
^(a)H̅ (f') = ∑_r=1^5 ^(a)λ̅_r
^(a)H̅^r (f), we have

    ⟨^(a)H^* (f) ^(a)H̅(f') ⟩    =    ∑_r,s=1^5 ^(a)λ_r ^(a)λ̅_s ⟨^(a)H^*r (f) ^(a)H̅^s (f') ⟩
        =    1/2δ(f - f')
    ∑_r,s=1^5 ^(a)λ_r ^(a)λ̅_s ^(a)Γ^rs ^(a)P_h (|f|)  ,

where ^(a)Γ^rs is the Hellings and Downs correlation function
<cit.>. After substituting Eq. (<ref>) into
Eq. (<ref>) and exercising the Dirac's delta function, we finally
obtain Eq. (<ref>).

From the definition of the variance-covariance matrix of the cross-correlation
statistic,

    ^(ab) C ≡⟨^(a)𝒮^(b)𝒮⟩ - ⟨^(a)𝒮⟩⟨^(b)𝒮⟩ ,

we have

    ^(ab) C      =    ∫_-∞^∞ df  ∫_-∞^∞ df'  ∫_-∞^∞ dk  ∫_-∞^∞ dk'
     δ_T (f - f')  δ_T (k - k') 
    
       ×    [ ⟨^(a)I^* (f) ^(a)I̅(f')
     ^(b)I (k) ^(b)I̅^*(k') ⟩ - 
    ⟨^(a)I^* (f) ^(a)I̅(f') ⟩⟨^(b)I (k) ^(b)I̅^*(k') ⟩]
    
       ×    ^(a)Q (f') ^(b)Q^* (k')  .

Since the random processes associated with the noises and the GW
background are Gaussian, we infer that the linear combinations
^(a)I (f), ^(a)I̅(f'),
^(b)I (f), ^(b)I̅(f') are
also Gaussian. This implies that the term in the integrand containing
four I-combinations, inside the ensemble average operator, can be
written in the following form <cit.>

    ⟨^(a)I^* (f) 
    ^(a)I̅(f')
    ^(b)I (k) 
    ^(b)I̅^*(k')
    ⟩    =    ⟨^(a)I^* (f) 
    ^(a)I̅(f') 
    ⟩⟨^(b)I (k) 
    ^(b)I̅^*(k') 
    ⟩
        +    ⟨^(a)I^* (f) 
    ^(b)I (k) 
    ⟩⟨^(a)I̅(f') 
    ^(b)I̅^*(k') 
    ⟩
        +    ⟨^(a)I^* (f) 
    ^(b)I̅^*(k') 
    ⟩⟨^(a)I̅(f') 
    ^(b)I (k) 
    ⟩

After substituting Eq. (<ref>) into Eq. (<ref>)
we get

    ^(ab) C     =    ∫_-∞^∞ df  ∫_-∞^∞ df'  ∫_-∞^∞ dk  ∫_-∞^∞ dk'
     δ_T (f - f')  δ_T (k - k') 
    
       ×    [ 
    ⟨^(a)I^* (f) 
    ^(b)I (k) 
    ⟩⟨^(a)I̅(f') 
    ^(b)I̅^*(k') 
    ⟩
    +
    ⟨^(a)I^* (f) 
    ^(b)I̅^*(k') 
    ⟩⟨^(a)I̅(f') 
    ^(b)I (k) 
    ⟩]
    
       ×    ^(a)Q (f') ^(b)Q^* (k')  .

By replacing the Is in terms of their GW signal and noise terms, and
applying the ensemble average operation (see Eqs. (<ref>,
<ref>) on the resulting expressions entering in
Eq. (<ref>), Eq. (<ref>) can be obtained after some
straightforward algebra.

