text:
# Cusped Mass Models Of Gravitational Lenses

## 1 Introduction

The families of density distributions used to model gravitational lenses are rooted in the lore of the early 1980s, where massive galaxies were believed to have finite core radii (e.g. Lauer 1985). Hence, surveys of lens models almost always use models related to the softened power-law surface density, $`\mathrm{\Sigma }(m^2+s^2)^{(\alpha 2)/2}`$ where $`m^2`$ is a two-dimensional ellipsoidal coordinate and $`s`$ is the core radius, as a template for exploring the effects of the radial mass distribution on lensing phenomenon (e.g. Kochanek 1995, Chen, Kochanek & Hewitt 1995, Grogin & Narayan 1996, Chae, Turnshek & Khersonsky 1998, Schmidt, Webster & Lewis 1998, Keeton & Kochanek 1998, Chae 1999, Keeton et al. 2000, Cohn et al. 2001). The asymptotic exponent $`\alpha `$ covers a range of historical models for galaxies from the singular isothermal ellipsoid (or SIE, $`\alpha =1`$ and $`s=0`$), to the modified Hubble profile ($`\alpha =0`$) and the Plummer model ($`\alpha =2`$). Viewed as projections of thin disks, these models include the Mestel (1963) disk (equivalent to the SIE), and the Kuzmin (1956) disk ($`\alpha =1`$). For $`\alpha =1`$ and $`1`$ there are simple analytic expressions for the potentials (see Keeton & Kochanek 1998 and references therein), and relatively fast approximation schemes have been developed for a broad range of exponents (Barkana 1998, Chae, Khersonsky & Turnshek 1998). Unfortunately, the only model from this broad class of potentials that can be regarded as astrophysical is the SIE.

Since the advent of HST and high resolution numerical simulations, our view of the central regions of galaxies has changed significantly, in that massive galaxies have central cusps rather than central core radii in both high resolution observations of luminosity profiles (e.g. Faber et al. 1997) and theoretical simulations (e.g. Navarro, Frenk & White 1997, Moore et al. 1998). Most simple, modern, analytic representations of galaxy mass distributions have density distributions of the form $`\rho r^\gamma (r+a)^{n+\gamma }`$ where $`\gamma `$ is the exponent of the central density cusp, $`n`$ is the exponent for the asymptotic slope of the density distribution, and $`a`$ is the break radius. Popular limits are the Hernquist (1990) model ($`\gamma =1`$, $`n=4`$), the Jaffe (1983) model ($`\gamma =2`$, $`n=4`$), the more general $`\eta `$-models (Dehnen 1993, Tremaine et al. 1994, $`n=4`$), the Navarro, Frenk & White (1997, NFW) models ($`\gamma =1`$, $`n=3`$) and the more general Zhao (1996) models. The classical de Vaucouleurs (1948) profile for early-type galaxies is similar to the Hernquist model, although its central “cusp” is better approximated by a slightly shallower exponent (Hernquist 1990).

Few of these more realistic models have been considered as lens models. Lehar et al. (2000) fit SIE and ellipsoidal de Vaucouleurs models as standard models which represent the two typical limits of central cusps. Only one analytic study by Evans & Wilkinson (1998) explored the general mathematics of models with different central cusps, focusing on circular models because there are no general analytic results for cusped ellipsoidal models. Bartelmann (1996) explored using the NFW density distribution to explain the radial arcs sometimes found in cluster lenses, and Wyithe, Turner & Spergel (2000) and Keeton & Madau (2001) explored the general statistical properties of a family of circular, cusped models for clusters. Cusps have qualitatively different lensing properties from finite cores, making it important to begin exploring their effects on gravitational lenses. In §2 we introduce a general class of ellipsoidal, cusped density distributions for use in gravitational lensing, and in §3 we apply the new models to the gravitational lenses APM 08279+5255 and B 1933+503. Our results are summarized in §4. The models are implemented and available as part of the lensmodel package of lens modeling codes (Keeton 2001).111The lensmodel package can be obtained from http://cfa-www.harvard.edu/castles

## 2 A New Class of Ellipsoidal Lens Model

Our models are defined from the spherical density distribution,

$$\rho =\frac{\rho _0}{r^\gamma }\frac{a^n}{(r^2+a^2)^{(n\gamma )/2}}$$
(1)
for cusp exponent $`\gamma `$, break radius $`a`$ and outer exponent $`n`$ (it is the same as a Zhao (1996) model with $`\alpha =1/2`$). This differs from the standard cusped models in using $`(r^2+a^2)^{(n\gamma )/2}`$ rather than $`(r+a)^{n\gamma }`$, leading to a smoother break in the density profile, more convenient analytic properties for the lensing calculations, and overlaps with several standard lens models. This form also has physical advantages in some cases. For example, the rotation curve of the Jaffe (1983) model, $`\rho r^2(r+a)^2`$, declines linearly for $`r<a`$ despite the $`1/r^2`$ central density profile, while our “pseudo-Jaffe” model with $`\rho r^2(r^2+a^2)^1`$ has a fairly flat rotation curve for $`r<a`$. Whenever the exponents of our density model match those of a standard cusped model we will refer to it as a “pseudo-” model. Hence for $`n=4`$, we have the pseudo-Jaffe model ($`\gamma =2`$), the pseudo-Hernquist model ($`\gamma =1`$), and the pseudo-$`\eta `$ models (any $`\gamma `$), and for $`n=3`$ and $`\gamma =1`$ we have the pseudo-NFW model. The model also matches the standard SIE model ($`\gamma =2`$, $`n=4`$ and $`ra`$) and can match all the softened power-law models (for $`\gamma =0`$ and viewing the break radius $`a`$ as a core radius).

Given the density distribution, the total mass of the model is

$$M_T=2\pi a^3\rho _0B(\frac{n3}{2},\frac{3\gamma }{2}),$$
(2)
where $`B(a,b)=\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)/\mathrm{\Gamma }(a+b)`$ is the Euler beta function. The mass is well-defined (i.e. finite) only for $`\gamma <3`$ and $`n>3`$. The surface mass density in units of the critical surface density for lensing is

$$\kappa (R)=\frac{\mathrm{\Sigma }(R)}{\mathrm{\Sigma }_{crit}}=\kappa _0B(\frac{n1}{2},\frac{1}{2})(1+\xi ^2)^{(1n)/2}{}_{2}{}^{}F_{1}^{}(\frac{n1}{2},\frac{\gamma }{2},\frac{n}{2};\frac{1}{1+\xi ^2})$$
(3)
where the normalization $`\kappa _0=\rho _0a/\mathrm{\Sigma }_{crit}`$, $`{}_{2}{}^{}F_{1}^{}`$ is the hypergeometric function and $`\xi =R/a`$ is the projected radius $`R`$ normalized by the break radius $`a`$. For a circular lens with $`n3`$, the deflection angle is

$$\alpha (R)=\frac{2\kappa _0a}{\xi }\left[B(\frac{n3}{2},\frac{3\gamma }{2})B(\frac{n3}{2},\frac{3}{2})(1+\xi ^2)^{(3n)/2}{}_{2}{}^{}F_{1}^{}(\frac{n3}{2},\frac{\gamma }{2},\frac{n}{2};\frac{1}{1+\xi ^2})\right]$$
(4)
For $`n=3`$ the deflection must be expressed in terms of generalized hypergeometric functions. The first derivative of $`\kappa `$, which is used in some of the integrals needed to compute the properties of ellipsoidal models, is

$`{\displaystyle \frac{d\kappa (R)}{dR}}`$ $`=`$ $`{\displaystyle \frac{\kappa _0}{a}}B({\displaystyle \frac{n+1}{2}},{\displaystyle \frac{1}{2}})\xi (1+\xi ^2)^{(n+3)/2}[n(1+\xi ^2){}_{2}{}^{}F_{1}^{}({\displaystyle \frac{n1}{2}},{\displaystyle \frac{\gamma }{2}},{\displaystyle \frac{n}{2}};{\displaystyle \frac{1}{1+\xi ^2}})+`$ (5)
$`+\gamma {}_{2}{}^{}F_{1}^{}({\displaystyle \frac{n+1}{2}},{\displaystyle \frac{\gamma +2}{2}},{\displaystyle \frac{n+2}{2}};{\displaystyle \frac{1}{1+\xi ^2}})]`$
In Figure 1 we plot the deflection angle, $`\alpha (R)`$, for several values of $`\gamma `$, a fixed outer exponent of $`n=4`$ and a fixed surface density scale $`M_T/\pi a^2\mathrm{\Sigma }_{crit}=1`$. The central deflection diverges for $`\gamma >2`$, is finite for $`\gamma =2`$ and passes through zero for $`\gamma <2`$.

The properties of the lens models in the core can be better understood by examining power series expansions of the deflection angle and the magnification as $`\xi =R/a0`$. If we define the constants $`A\kappa _0B[(n1)/2,(1\gamma )/2]`$ and $`B\kappa _0B[(\gamma 1)/2,1/2]`$, for $`\gamma 1`$ we find that

$$\kappa (R)=A+B\xi ^{1\gamma }+O(\xi ^2,\xi ^{3\gamma })$$
(6)
$$\alpha (R)=R\left[A+\frac{2B}{3\gamma }\xi ^{1\gamma }+O(\xi ^2,\xi ^{3\gamma })\right]$$
(7)
with logarithmic divergences when $`\gamma =1`$. For circular lenses, the inverse of the magnification, $`\mu ^1=12\kappa (R)\alpha ^2(R)/R^2+2\kappa (R)\alpha (R)/R`$, can be expressed in terms of $`\kappa (R)`$ and $`\alpha (R)`$ (e.g. Schneider, Ehlers and Falco 1992). Using eqns. (6) and (7), we find a power series expansion for the central magnification of

$$\mu ^1=(A1)(A1+2B\xi ^{1\gamma })+O(\xi ^{2(1\gamma )},\xi ^2,\xi ^{3\gamma }).$$
(8)
For $`\gamma 2`$ the cusp model produces only an even number of images; there is no central “odd” image. For $`\gamma <2`$ the model does produce a central, odd image, but for $`1<\gamma <2`$ the image can be significantly demagnified compared to the other images (e.g. Evans & Wilkinson 1998).

When we make the models ellipsoidal we must calculate the deflections and magnifications numerically. The projected radius $`R`$ is replaced by the ellipsoidal coordinate $`R^2\eta ^2=x^2+y^2/q^2`$ where $`q`$ is the projected axis ratio and the Cartesian coordinates $`x`$ and $`y`$ are rotated to match the principal axes of the lens. The elliptical symmetry reduces the computation of the lensing properties to one-dimensional integrals (see e.g. Schramm 1990). Modeling lenses using the numerical integrals is too slow for practical use, so we tabulate the potential, the deflections and the components of the magnification tensor as a function of scaled radius ($`x/a`$ and $`y/a`$) and axis ratio ($`q`$). The modeling is then done by interpolating over the tabulated models. We have implemented the ellipsoidal cusped models as part of the lensmodel software package (Keeton 2001).

## 3 Application of the New Model to APM 08279+5255 and B 1933+503

We now apply these models to two gravitational lenses. The first lens, APM 08279+5255 (Irwin et al. 1998), is the first clear example of a gravitational lens with a bright central image (see Ibata et al. 1999, Egami et al. 2000). In this system we can explore the role of a shallow central cusp, rather than a finite core radius, in producing a visible central image. The second lens, B 1933+507 (Sykes et al. 1998), is a 10 image gravitational lens which is well-suited to determining the radial mass distribution of the lens galaxy. A detailed analysis by Cohn et al. (2001) found that models with cusps slightly shallower than isothermal produced the best fits, but these models were not globally consistent due to their infinite masses and perpetually rising rotation curves.

### 3.1 APM 08279+5255

When Irwin et al. (1998) discovered the $`z_s=3.87`$ BAL quasar APM 08279+5255 they also noted that the image was extended and likely to be a compact double (confirmed by Ledoux et al. 1998). HST images in the infrared with NICMOS (Ibata et al. 1999) and improved ground-based images (Egami et al. 2000) demonstrated that the source was actually a triple. There are two bright point sources separated by $`0\stackrel{}{\mathrm{.}}38`$ with a flux ratio of 0.8, and between the two lies a relatively bright third image with a flux ratio of 0.2 relative to the brightest image. A lens galaxy has yet to be detected, although this is not surprising given the extraordinary brightness of the quasar images (about 15 R mag!). The image configuration is easily modeled as an isothermal ellipsoid with a finite core radius (Ibata et al. 1999, Egami et al. 2000). However, given our expectation that the lens galaxy should not have a core radius, we will investigate whether the system is equally well modeled by a shallow cusp.

While it is likely that the system is composed of three lensed images given the fluxes and the similarity of the infrared colors (see Ibata et al. 1999), we can use optical flux ratios for further confirmation. The CfA/Arizona Space Telescope Lens Survey (CASTLES) obtained I-band (F814W) and V-band (F505W) images of APM 08279+5255 on 1999.11.09 with the Hubble Space Telescope (HST). Each observation consisted of 20 dithered 40 sec exposures designed to avoid saturating the extraordinarily bright quasars. The data were reduced and analyzed following the standard procedures of Lehar et al. (2000) for the CASTLES observations. Table 1 presents the resulting astrometry and fluxes for the system where we also include the photometry obtained from Ibata et al. (1999) with the NICMOS infrared images. We now see that not only are the infrared colors of the images similar, but the spectrum remains similar over the full range from the V-band to K-band, which makes it virtually impossible for the three images to be anything other than three images of the BAL quasar. We use these results to constrain our lens models with a more conservative error of 0$`\stackrel{}{\mathrm{.}}`$01 for the quasar positions.

We first fit the system using a softened isothermal ellipsoid. Like Ibata et al. (1999) and Egami et al. (2000) we find that the model produces a perfect fit to the data given a core radius of $`a=0\stackrel{}{\mathrm{.}}16_{0\stackrel{}{\mathrm{.}}05}^{+0\stackrel{}{\mathrm{.}}10}`$, which would be $`0.7_{0.2}^{+0.4}(\mathrm{\hspace{0.17em}0.9}_{0.3}^{+0.6})h^1`$ kpc for a lens redshift of $`z_l=0.5`$ ($`1.0`$) in an $`\mathrm{\Omega }_0=0.3`$ flat cosmological model. Since the bulges of massive galaxies usually show cusps down to smaller scales (see e.g. Faber et al. 1997), this solution appears to be inconsistent with the known properties of galaxies. The best model reaches $`\chi ^2=0`$, as we would expect given a system with $`6`$ constraints and $`6`$ parameters leaving $`N_{dof}=0`$ degrees of freedom. We also fit a model where the lens was an exponential disk without a bulge or dark matter halo. Exponential disks have little density contrast between their centers and the Einstein ring, which makes it relatively easy to fit the system. We show the parameters of these two models in Table 2. Lens models with a de Vaucouleurs profile were not able to fit the data for any effective radius. The central density of the $`R^{1/4}`$ law roughly corresponds to a central density distribution with $`\rho r^{3/4}`$ as $`r0`$ (see Hernquist 1990), and a central exponent $`\gamma =3/4`$ is too steep to reproduce the central bright image (see below).

Next we fit a sequence of our cusped models to the system (see Figures 2 & 3). The models are underconstrained when we allow the cusp exponent and the break radius both to vary ($`N_{dof}=1`$ given 6 constraints and 7 parameters). Nevertheless, we could only find successful models for cusps outside the standard range of $`1\gamma 2`$ found for stellar luminosity profiles and dark matter simulations. We found that the cusp exponent has to satisfy $`\gamma <0.4`$ to achieve $`\chi ^2<1`$, and the models can produce perfect fits only for $`\gamma <0.2`$. These cusps are so shallow that we are almost at the limit of reproducing the softened power law models ($`\gamma =0`$). While the cusp models with a small break radius have total magnifications comparable to the softened isothermal models ($`\mu _T100`$), there is a second solution branch with a very large break radius where the total magnification slowly diverges. For example a model with $`\gamma =0.2`$ and $`a1\mathrm{}`$ (see in Fig. 3) provides a perfect fit combined with a total magnification of $`\mu _T1800`$. Without any correction for lensing, the implied luminosity of the system is $`5\times 10^{15}L_{}`$ ($`\mathrm{\Omega }=1`$, Egami et al. 2000), and in the previous models and most of our new models the lens magnification of $`\mu _T100`$ reduces this to a more physically plausible $`5\times 10^{13}L_{}`$. On the branch of solutions with $`\mu _T1800`$, however, the intrinsic luminosity of the system drops to $`3\times 10^{12}L_{}`$ and the source ceases to be remarkable in its intrinsic properties. We can find models with steeper central cusps ($`1<\gamma <2`$) by adding an external shear, but these models are very underconstrained (9 parameters for 6 constraints) and require very high lens ellipticities. They are probably unphysical, and we do not present any examples.

Table 2 summarizes a representative set of models. Since the models are underconstrained and do not include an external shear (which we know is required in a complete model, see Keeton, Kochanek & Seljak 1997), we do not include error estimates for the parameters. It is impossible to distinguish between the models using the current constraints, in large part because the actual deflection profiles of the models are very similar (see Fig. 4). The position of the lens galaxy is predicted to vary by only a few milli-arcseconds between the models, and it will be extraordinarily difficult to detect the lens galaxy at all given the enormous flux from the lensed images. With the exception of the high magnification models, the time delays vary by only 20% and a precise knowledge of the Hubble constant would be required to distinguish between the models. The delay ratios, which are independent of the Hubble constant, vary by only 1% between the models. The only prospect for significantly improving the constraints in this system is to search for the lensed images of the quasar host galaxy (see Kochanek, Keeton & McLeod 2001), although the brightness of the quasar images makes this a challenging observation.

### 3.2 B 1933+503

B 1933+507 (Sykes et al. 1998, Nair 1998, Marlow et al. 1999, Browne et al. 1999, Biggs et al. 2000, Norbury et al. 2000) has 10 lensed images of a $`z_s=2.62`$ radio source generated by a $`z_l=0.76`$ lens galaxy. The source has three radio components, two of which have 4 lensed images and one of which has two lensed images (see Fig. 5). The large number of constraints provided by the 10 images makes this lens an ideal laboratory for studying the radial mass distribution of the lens galaxy. Extensive models by Cohn et al. (2001) found that the best mass distribution they considered was a power law model with a singular core and an exponent slightly shallower than isothermal ($`s=0`$ and $`\alpha 1.1`$). The $`\alpha >1`$ power-law models with zero core radius favored by the Cohn et al. (2001) analysis correspond to our cusped models with $`\gamma =3\alpha `$ in the limit of infinite break radius, $`a\mathrm{}`$. The Cohn et al. (2001) models are unphysical because they have infinite mass and a rising rotation curve. Our new cusp models allow us to examine these solutions using the same central cusp exponents combined with an outer break radius to produce a physical model. We used the same constraints as Cohn et al. (2001) to study the system.

Figure 6 shows the goodness of fit as a function of the break radius for a representative range ($`1.1\gamma 2.1`$) of cusp exponents. Cusps shallower than $`\gamma =1.1`$ are ruled out because they produce observable central images. The best fit models for a range of cusp exponents are presented in Table 3, and Figure 7 shows goodness of fit and the permitted range of break radii as a function of the cusp exponent. The best fit model has $`\gamma =1.85`$ and a large break radius (which we fixed to $`a=10\mathrm{}`$) with a total fit statistic $`\chi ^2=10.8`$ for 17 degrees of freedom given that the cusped model has 9 parameters and the data supplies 26 constraints (see Cohn et al. 2001). Figure 5 compares the data and the model images and shows the locations of the source components relative to the model caustics. The largest contribution to the $`\chi ^2`$ (4.3) is due to the position of image 1a (see Fig. 5). This best fit model is nearly identical to that of Cohn et al. (2001) who considered models with the same central density exponent but an infinite break radius.

The cusp exponent is limited to the range $`1.6<\gamma <2.0`$, marginally excluding the isothermal exponent $`\gamma =2`$ and strongly excluding the popular $`\gamma =1`$ cusps of the NFW and Hernquist models. Shallower cusps, $`\gamma <1.8`$, require finite break radii that become steadily smaller as the cusp becomes shallower. In other words, to maintain a large enough surface density contrast between the center of the lens and the Einstein ring to avoid creating visible central images, the shallower cusps require global profiles that rapidly break onto the steeper outer profile. As pointed out by Rusin & Ma (2000), shallow cusps of infinite extent ($`\gamma <1.8`$ with $`a\mathrm{}`$) can be ruled out because they will produce visible central images in too many lensed systems to be consistent with the data. Here, we can rule out the shallower cusps for any break radius, rather than the limiting case of a break radius which is large compared to the Einstein ring of the lens. Steeper central cusps, $`\gamma >2`$, intrinsically lack central images (see §2), but are also inconsistent with the constraints.

Figure 8 shows the monopole deflection of lens models which are consistent with the data. Despite the enormous freedom in the models, only models with nearly flat rotation curves are consistent with the data. The dashed and solid lines are used to distinguish between the models with $`\gamma <1.8`$ where the best lens model has a finite break radius and $`\gamma 1.8`$ where the best model is reached for $`a\mathrm{}`$. This discontinuity produces two different trends in the monopole deflection as well as a small shift in the crossing point.

The location of the outer break radius does affect the time delays between the images, where for a given $`\gamma `$ we find shorter delays as we increase the break radius. Figure 9 shows the longest time delay between the four images of the compact radio-core (it corresponds to images 1 and 6, see Figure 5) as a function of the $`\chi ^2`$. The points correspond to the best lens models for different values of $`\gamma `$. The dotted lines show the variation of the time delay as a function of the $`\chi ^2`$ when we vary the break radius for each value of $`\gamma `$. The dashed line shows the limit $`\mathrm{\Delta }\chi ^21`$. Given the generality of the cusped model we can consider the values inside this region as a very robust prediction of the time delay. This figure also shows more clearly the consequences of the jump in the optimal break radii from finite to infinite at $`\gamma 1.8`$.

## 4 Summary

Gravitational lenses are used both to estimate the mass distribution of galaxies (e.g. Kochanek et al. 1989, Kochanek 1991, Kochanek & Narayan 1992, Wambsganss & Paczynski 1994, Kochanek 1995, Chen et al. 1995, Wallington, Kochanek & Narayan 1996, Ellithorpe, Kochanek, & Hewitt 1996, Chae et al. 1998, Schmidt et al. 1998, Kochanek et al. 2000, Cohn et al. 2001) and the Hubble constant from measurements of the time delays between the lensed images (e.g. Grogin & Narayan 1996, Kundic et al. 1997, Schechter et al. 1997, Keeton & Kochanek 1997, Impey et al. 1998, Chae 1999, Barkana et al. 1999, Koopmans & Fassnacht 1999, Fassnacht et al. 1999, Keeton et al. 2000, Kochanek, Keeton & McLeod 2001, Oscoz et al. 2001). Almost all these analyses use parametric models for the mass distribution of the lens. Parametric models have three advantages: first, all standard parametric models represent physically self-consistent, dynamically stable equilibrium systems; second, parametric models are relatively easy to interpret physically; and third, they are computationally efficient. The disadvantage of parametric models, particularly for estimating the systematic uncertainties in the mass distribution or the Hubble constant, is that they may overly restrict the true freedom in the mass distributions of real astronomical systems.

Our new family of ellipsoidal models for gravitational lenses is a significant advance for two reasons. First, the new models are far more realistic models of galaxies than the standard lens models given our current understanding of the central stellar density profiles of galaxies (e.g. Faber et al. 1997) and simulations of dark matter halos (Navarro et al. 1997, Moore et al. 1998), where the central density follows a power-law cusp with $`\rho r^\gamma `$ and $`1<\gamma <2`$. Our model, which is a subcase of the more general Zhao (1996) models, includes analogues of many currently popular analytic density distributions such as the Jaffe (1983) model, the Hernquist (1990) model, the $`\eta `$ model (Dehnen 1993, Tremain et al. 1994) and the NFW model (Navarro et al. 1997). Second, independent of the models being physically more realistic than the standard lens models, they provide a genuinely different class of density distributions from the standard lens models. Since the freedom in the radial mass distribution appears to be the dominant source of systematic uncertainties in determining the Hubble constant from gravitational lenses (see e.g Impey et al. 1998, Fassnacht et al. 1999, Keeton et al. 2000), the new models can be used as an additional probe of the uncertainties. The cusped models are implemented and publically available as part of the lensmodel (Keeton 2001) package for modeling gravitational lenses.

Non-parametric models of the mass distribution avoid the restrictions of parametric models. In particular, the elegant linear programming method introduced by Saha & Williams (1997) assumes only that the surface mass density of the lens is non-negative. This method will find a broader range of systematic uncertainties in the mass distribution or the Hubble constant than any parametric model or set of parametric models. Unfortunately, the models also lose the first advantage of the parametric models, that they correspond to physical, dynamically stable systems. Physical systems must have non-negative density distributions and phase space distribution functions, and many non-negative surface densities will correspond to unphysical models with regions of negative density or phase space density (see Cohn et al. 2001). Models which are physical may still be improbable because they are dynamically unstable, with lifetimes of only a few dynamical (crossing) times.222It is possible to detail these additional conditions on the surface density $`\mathrm{\Sigma }(R)`$ for spherical systems (see Binney & Tremaine 1987). The mass density, $`\rho =\pi ^1_r^{\mathrm{}}𝑑R(d\mathrm{\Sigma }/dR)(R^2r^2)^{1/2}`$, must be positive for all $`r`$, which is trivially satisfied if the surface density distribution is a monotonically declining function, $`d\mathrm{\Sigma }/dR<0`$. For an isotropic distribution function, which depends only on the binding energy $`ϵ`$, the distribution function corresponding to a given density distribution is positive if $`_0^ϵ(d\rho /d\psi )(ϵ\psi )^{1/2}𝑑\psi `$ is an increasing function of $`ϵ=\psi v^2/2`$. This is trivially satisfied if $`d\rho /d\psi >0`$, which is true for all the parametric models. The system is dynamically stable if the derivative of the distribution function $`df/dϵ<0`$. Each step in this process imposes a smoothness criterion (through the constraints on the derivatives of the previous variable) which goes well beyond simply requiring a positive surface density. Viewed in terms of the surface density $`\mathrm{\Sigma }`$, a positive density is a restriction on $`d\mathrm{\Sigma }/dR`$, a positive distribution function is a restriction on $`d^2\mathrm{\Sigma }/dR^2`$, and stability is a restriction on $`d^3\mathrm{\Sigma }/dR^3`$! Unfortunately, there are no analytic expressions even for spherical systems with distribution functions depending on energy and angular momentum, let alone systems with more complicated symmetries. Although galaxies are more complicated than simple spherical systems, the problems we illustrate using the spherical systems apply to all galaxies and mean that a positive surface density is only a first, minimal criterion for a physical mass model.  These additional criteria mean that the current non-parametric methods greatly over-estimate the freedom in the mass distribution and frequently produce unphysical models. These issues make the two approaches to modeling lenses complementary. The parametric models are always physical but can underestimate the freedom in the mass distribution, while the non-parametric models are frequently unphysical and overestimate the freedom in the mass distribution.

We explored these new parameterized mass distributions by using them to model two lenses, APM 08279+5255 and B 1933+503. APM 08279+5255 (Irwin et al. 1998, Ibata et al. 1999, Egami et al. 2000) is the only lens which clearly possesses three images of the background source. We find further evidence of this from the CASTLES V-band and I-band HST images of the system, which show that the flux ratios of the images change little from the K-band in the infrared to the V-band in the optical. Although the ability of the cusped models to model a three image gravitational lens like APM 08279+5255 without a finite core radius was one of the motivations for exploring the cusped models, we actually found that APM 08279+5255 cannot be modeled by cusps in the favored range $`1<\gamma <2`$. Only models with $`\gamma <0.4`$, where the cusps become very similar to finite core radii ($`\gamma =0`$), could successfully explain the properties of this lens. This suggests that the lens galaxy in APM 08279+5255 is a spiral galaxy, where rotation curves are sometimes consistent with finite central densities rather than steep central density cusps. Alternatively, this could be an example of a “disk” geometry lens (Keeton & Kochanek 1998) where we see three images produced by a nearly edge on disk galaxy. In a disk lens geometry, the lens galaxy should be offset perpendicular to the A–B image separation. The small image separation ($`0\stackrel{}{\mathrm{.}}38`$) also argues for a spiral lens (see Fukugita & Turner 1991, Kochanek 1996), but the similar colors of the three images suggests that there is little dust along the three lines of sight through the lens. Unfortunately, the extreme brightness of the lensed images makes it nearly impossible to find the far fainter lens galaxy, let alone measure its position and morphology.

The second lens we considered, B 1933+503 (Sykes et al. 1998), already had an extensive model survey suggesting that the best models would have steep central cusps (Cohn et al. 2001). Here we show using physically realistic cusped models rather than simple power law density distributions that the system requires a steep central cusp, $`\rho r^{1.8\pm 0.2}`$, with a break radius located well outside the Einstein ring of the lens. Shallow cusps, particularly the $`\rho r^1`$ cusps of the Hernquist and NFW models, are inconsistent with the data for any break radius. Lenses with shallow central cusps and break radii larger than the Einstein radius of the lens generically predict visible central images, and are ruled out by the absence of such images in well-defined lens samples (Rusin & Ma 2000). Shallow cusps combined with a break radius smaller than the Einstein ring can demagnify the central image, but the monopole deflection curve (the rotation curve) becomes centrally peaked and is ruled out by the need for a flat deflection profile to fit the positions of the observed images. A massive central black hole can also demagnify central images (Mao, Witt & Koopmans 2000), thereby weakening the limits on shallow cusps (or a central core radius) based on the absence of central images. In the case of B 1933+507, the restrictions on shallow cusps are only partially due to the limits on the presence of central images, and when we eliminate the flux limits on the presence of a central image the restriction on the cusp exponent is only reduced to $`\gamma >1.5`$ from $`\gamma >1.6`$. The generality of the new developed cusped models yields a very robust prediction of the time delay for this system. In particular the lower limit obtained in the time delay computation for any cusped exponent $`\gamma `$ and any break radius (see Figure 9) would produce a strong constraint in the upper value of the Hubble constant provided an observational measurement of the time delay.

Acknowledgments: We thank the CASTLES collaboration (E.E. Falco, C.D. Impey, C.R. Keeton, C.S. Kochanek, J. Lehar, B.A. McLeod, J.A. Muñoz & H.-W. Rix) for permission to include their APM 08279+5255 data. We also thank Joanne Cohn for providing us with her constraints for B 1933+507 so that we could exactly match her models of B 1933+503. CSK is supported by the Smithsonian Institution and NASA grant NAG5-9265.


format: text

text:
# Space missions to detect the cosmic gravitational-wave background

## 1 Introduction

Before embarking on a quest to discover the Cosmic Gravitational-wave Background (CGB), it is worth reflecting on how much has been learned from its electromagnetic analog, the Cosmic Microwave Background (CMB). The recent Boomerang and Maxima experiments have furnished detailed pictures of the universe some 300,000 years after the Big Bang. These measurements of the CMB anisotropies have been used to infer that the visible universe is, to a good approximation, spatially flat. The earlier COBE-DMR measurements of the CMB anisotropy fixed the scale for the density perturbations that formed the large scale structure we see today, and showed that the density perturbations have an almost scale free spectrum. With all this excitement over the anisotropy measurements it is easy to forget how much the CMB taught us before the anisotropies were detected. The initial observation that the CMB is highly isotropic is one of the major observational triumphs of the Big Bang theory. Equally important was the discovery by COBE-FIRAS that the CMB has an exquisite black body spectrum with a temperature of 2.728 Kelvin.

The lesson that we take from the CMB is that while it would be wonderful to have a COBE style map of the early universe in gravitational waves, a great deal can be learned by detecting the isotropic component. The main focus of this work is on the fixing the amplitude of the CGB at some frequency. Some attention will also be given to the prospects of measuring the the energy spectrum. Detecting anisotropies in the CGB is a considerably harder problem that we address elsewhere.

There are two major obstacles that stand in the way of detecting the CGB. The first is the extreme weakness of the signal, and the second is the competing stochastic background produced by astrophysical sources. Figure 1 shows various predictions for the CGB energy spectra in several early universe scenarios including standard inflation, a particular M-theory model and quintessence based inflation. Also shown is the sensitivity curve for the proposed Laser Interferometer Space Antenna (LISA), and a compilation of possible extragalactic astrophysical foregrounds taken from the work of Schneider et al.. The spectrum is expressed in terms of $`\mathrm{\Omega }_{\mathrm{gw}}(f)`$, the energy density in gravitational waves (in units of the critical density) per logarithmic frequency interval, multiplied by $`h_0^2`$, where $`h_0`$ is the Hubble constant in units of 100 km s-1 Mpc-1. The predictions for the CGB power spectra are fairly optimistic, and some of them come close to saturating existing indirect bounds on $`\mathrm{\Omega }_{\mathrm{gw}}(f)`$ – see the review by Maggiore for details.

With the exception of the M-theory motivated model, we see that the major obstacle to detecting the CGB is not the sensitivity of our detectors, but competition from astrophysical foregrounds. The astrophysical foregrounds are produced by the combination of a great number of weak sources that add together to form a significant stochastic signal. Based on current estimates, it is thought that the contribution from astrophysical sources falls off fairly rapidly below 1 $`\mu `$Hz. For that reason, we will focus our attention on detecting the CGB in the sub-micro Hertz regime. To put this number into perspective, an interferometer built using LISA technology would need spacecraft separated by one-tenth of a light year to reach a peak sensitivity at 1 $`\mu `$Hz! A more practical approach is to use two smaller interferometers that lack the intrinsic sensitivity to detect the CGB on their own, but together are able to reach the necessary sensitivity. The idea is to cross correlate the outputs from the two interferometers and integrate over some long observation time $`T`$. Since the noise in the two interferometers is uncorrelated while the signal is correlated, the signal to noise ratio will steadily improve as the observation time is increased.

In principle, it is possible to detect an arbitrarily weak signal by observing for an arbitrarily long time. In practice, the prospects are not so rosy as the sensitivity of the detector pair only improves on the sensitivity of a single detector as $`(T\mathrm{\Delta }f)^{1/2}`$, where $`\mathrm{\Delta }f`$ is the frequency bandpass. The frequency bandpass is typically taken to equal the central observing frequency. Suppose we try and look for the CGB at a frequency of $`f\mathrm{\Delta }f=10^3`$ Hz using the two LIGO detectors. Over an observation time of $`T=1`$ year, the sensitivity is improved by a factor of $`2\times 10^5`$ by cross correlating the Washington and Louisiana detectors. For a pair of LISA detectors with a bandpass of $`\mathrm{\Delta }f10^2`$ Hz, cross correlating for a year would result in a 500 fold improvement in sensitivity. However, for detectors operating in the $`\mu `$Hz range, cross correlating two detectors for one year only improves on the sensitivity of a single detector by a factor of $`6`$. So why use two interferometers when one will do? The reason is simple: with a single interferometer it is not possible to tell the difference between the CGB and instrument noise since both are, to a good approximation, stationary Gaussian random processes. Moreover, it is impossible to shield a detector from gravitational waves in order to establish the noise floor. With two or more independent interferometers the cross correlated signal can be used to seperate the stochastic bacground from instrument noise. The LISA system has three arms, and the signal from these arms can be used to form three different Michelson interferometers. However, these interferometers share spacecraft, and thus common sources of noise, so nothing is gained by cross correlating the outputs. However, Tinto et al. have recently suggested another way of combining the signals, known as a Sagnac system, that creates an interferometer that is fairly insensitive to gravitational waves. Since the response of this “bad” interferometer is noise dominated, its output can be used to estimate the noise floor for the standard interferometer configuration. Thus, a single LISA type observatory might be able to discriminate between the CGB and instrument noise.

We will show that cross correlating two LISA interferometers results in a signal to noise ratio of

$$\mathrm{SNR}=1.4\left(\frac{T}{\mathrm{year}}\right)^{1/2}\left(\frac{\mathrm{\Omega }_{\mathrm{gw}}h_0^2}{10^{14}}\right)$$
(1)
for a scale invariant stochastic gravitational wave background. This represents a five hundred fold improvement over the sensitivity of a single LISA interferometer. If it were not for the astrophysical foregrounds that are thought to dominate the CGB signal in the LISA band, a pair of LISA detectors would stand a good chance of detecting the CGB. In an effort to avoid being swamped by astrophysical sources, we focus our efforts on detecting the CGB below 1 $`\mu `$Hz. We show that a pair of interferometers can reach a signal to noise of

$$\mathrm{SNR}=3.1\left(\frac{T}{\mathrm{year}}\right)^{1/2}\left(\frac{\mathrm{\Omega }_{\mathrm{gw}}h_0^2}{10^{14}}\right)\left(\frac{L}{\sqrt{3}\mathrm{AU}}\right)^2\left(\frac{3\times 10^{16}\mathrm{ms}^2}{\delta a}\right)^2\left(\frac{f}{\mu \mathrm{Hz}}\right)^{3/2}$$
(2)
for a scale invariant CGB spectrum. The above result is scaled against a possible LISA follow-on mission (LISA II) that calls for two identical interferometers, comprising six equally spaced spacecraft that form two equilateral triangles overlayed in a star pattern. The constellations follow circular orbits about the Sun in a common plane at a radius of 1 AU, so that each interferometer arm is $`L=\sqrt{3}`$ AU in length. The acceleration noise, $`\delta a`$, is scaled against a value that improves on the LISA specifications by one order of magnitude. The LISA II design could detect a stochastic background at 99% confidence with one year of observations if $`\mathrm{\Omega }_{\mathrm{gw}}h_0^2=10^{14}`$. Ideally we would like to reach a sensitivity of at least $`\mathrm{\Omega }_{\mathrm{gw}}h_0^210^{20}`$ in order to detect inflationary spectra with a mild negative tilt. It is difficult to achieve this with spacecraft orbiting at 1 AU as thermal fluctuations due to solar heating of the spacecraft make it hard to reduce the acceleration noise below $`\delta a10^{16}\mathrm{ms}^2`$. The best strategy is to increase the size of the orbit, $`R`$, as this increases the size of the interferometer, $`L=\sqrt{3}R`$, and decreases the thermal noise by $`R^2`$.

## 2 Outline

We begin by deriving the response of a pair of cross correlated space-based laser interferometers to a stochastic background of gravitational waves. Similar calculations have been done for ground based detectors, but these omit the high frequency transfer functions that play an important role in space based systems. Our calculations closely follow those of Allen and Romano, and reduce to theirs in the low frequency limit111Here high and low frequencies are defined relative to the transfer frequency of the detector $`f_{}=c/(2\pi L)`$, where $`L`$ is the length of one interferometer arm.. Many of the calculational details that we omit for brevity can be found in their paper. Having established the general formalism for cross correlating space based interferometers we use our results to identify some promising configurations for a GrAvitational Background Interferometer (GABI) mission to detect the CGB.

## 3 Detector Response

We attack the problem of cross correlating two space based interferometers in stages. We begin by reviewing how the Doppler tracking of a pair of spacecraft can be used to detect gravitational waves. Using this result we derive the response of a two arm Michelson interferometer and express the result in a convenient coordinate-free form. The response of a single interferometer is then used to find the sensitivity that can be achieved by cross correlating two space based interferometers.

### 3.1 Single arm Doppler tracking

Following the treatment of Hellings, we consider two free spacecraft, one at the origin of a coordinate system and the other a distance $`L`$ away at an angle $`\theta `$ from the $`z`$ axis. Both spacecraft are at rest in this coordinate system. The spacecraft at the origin sends out a series of photons, while a weak plane gravitational wave is passing through space in the $`+z`$ direction. To leading order, the spacetime metric in the transverse-traceless gauge is:

$$ds^2=c^2dt^2+dz^2+(1+h\mathrm{cos}2\psi )dx^2+(1h\mathrm{cos}2\psi )dy^2h\mathrm{sin}2\psi dxdy,$$
(3)
where $`h(tz/c)`$ is the wave amplitude and $`\psi `$ is the angle between the principal polarization vector and the $`x`$ axis. If we choose our coordinates such that the second spacecraft is in the $`xz`$ plane, the path of the photons can be parameterized by

$$x=\rho \mathrm{sin}\theta z=\rho \mathrm{cos}\theta .$$
(4)
The photon path in the perturbed spacetime is given by $`ds^2=0`$, or

$$cdt\pm (1+\frac{1}{2}h\mathrm{cos}2\psi \mathrm{sin}^2\theta )d\rho .$$
(5)
The round trip journey from spacecraft 1 to spacecraft 2 and back again is given by

$`\mathrm{}(t_2t_0)={\displaystyle _{t_0}^{t_2}}c𝑑t`$ $`=`$ $`2L+{\displaystyle \frac{1}{2}}\mathrm{cos}2\psi \mathrm{sin}^2\theta ({\displaystyle _0^L}h[t_0+{\displaystyle \frac{\rho }{c}}(1\mathrm{cos}\theta )]d\rho `$ (6)
$`+{\displaystyle _0^L}h[t_1+{\displaystyle \frac{L}{c}}{\displaystyle \frac{\rho }{c}}(1+\mathrm{cos}\theta )]d\rho ).`$
Here $`t_0`$ is the time of emission, $`t_1`$ is the time of reception at spacecraft 2 and $`t_2`$ is the time of reception back at the first spacecraft. Using $`h(q)=h_0\mathrm{exp}(i\omega q)`$, the varying portion of the round-trip distance is

$`\delta \mathrm{}(t_2)=Lh(t_2)\mathrm{cos}2\psi \mathrm{sin}^2\theta {\displaystyle \frac{1}{2}}`$ $`(\mathrm{sinc}\left({\displaystyle \frac{f}{2f_{}}}(1\mathrm{cos}\theta )\right)\mathrm{exp}(i{\displaystyle \frac{f}{2f_{}}}(3+\mathrm{cos}\theta ))`$ (7)
$`+\mathrm{sinc}\left({\displaystyle \frac{f}{2f_{}}}(1+\mathrm{cos}\theta )\right)\mathrm{exp}(i{\displaystyle \frac{f}{2f_{}}}(1+\mathrm{cos}\theta ))),`$
where $`f=\omega /(2\pi )`$ is the frequency of the gravitational wave and $`f_{}=c/(2\pi L)`$ is the transfer frequency. To zeroth order in $`h`$, the optical path length is $`\mathrm{}=2L`$. The general, coordinate independent version of this expression is given by

$$\frac{\delta \mathrm{}}{\mathrm{}}=\frac{1}{2}(𝐮𝐮):𝐡(\widehat{\mathrm{\Omega }})𝒯(𝐮\widehat{\mathrm{\Omega }},f),$$
(8)
where $`𝐮`$ is a unit vector pointing from the first to the second spacecraft and $`\widehat{\mathrm{\Omega }}`$ is a unit vector in the direction the gravitational wave is propagating. The colon denotes the double contraction $`𝐚:𝐛=a_{ij}b^{ij}`$. The transfer function $`𝒯`$ is given by

$`𝒯(𝐮\widehat{\mathrm{\Omega }},f)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\mathrm{sinc}\left({\displaystyle \frac{f}{2f_{}}}(1𝐮\widehat{\mathrm{\Omega }})\right)\mathrm{exp}(i{\displaystyle \frac{f}{2f_{}}}(3+𝐮\widehat{\mathrm{\Omega }}))`$ (9)
$`+`$ $`\mathrm{sinc}\left({\displaystyle \frac{f}{2f_{}}}(1+𝐮\widehat{\mathrm{\Omega }})\right)\mathrm{exp}(i{\displaystyle \frac{f}{2f_{}}}(1+𝐮\widehat{\mathrm{\Omega }}))].`$
This expression for $`𝒯`$ agrees with the one derived by Schilling. The gravitational wave is described by the tensor

$$𝐡(t,𝐱)=h^+(\omega t\omega \widehat{\mathrm{\Omega }}𝐱)ϵ^+(\widehat{\mathrm{\Omega }},\psi )+h^\times (\omega t\omega \widehat{\mathrm{\Omega }}𝐱)ϵ^\times (\widehat{\mathrm{\Omega }},\psi ),$$
(10)
with polarization tensors

$`ϵ^+(\widehat{\mathrm{\Omega }},\psi )=𝐞^+(\widehat{\mathrm{\Omega }})\mathrm{cos}2\psi 𝐞^\times (\widehat{\mathrm{\Omega }})\mathrm{sin}2\psi `$
$`ϵ^\times (\widehat{\mathrm{\Omega }},\psi )=𝐞^+(\widehat{\mathrm{\Omega }})\mathrm{sin}2\psi +𝐞^\times (\widehat{\mathrm{\Omega }})\mathrm{cos}2\psi .`$ (11)
The basis tensors can be written as

$`𝐞^+(\widehat{\mathrm{\Omega }})=\widehat{m}\widehat{m}\widehat{n}\widehat{n}`$
$`𝐞^\times (\widehat{\mathrm{\Omega }})=\widehat{m}\widehat{n}+\widehat{n}\widehat{m}`$ (12)
where $`\widehat{m}`$, $`\widehat{n}`$ and $`\widehat{\mathrm{\Omega }}`$ are an orthonormal set of unit vectors. The expression in (7) can be recovered from (8) by setting

$`h^+=h,h^\times =0,\widehat{\mathrm{\Omega }}=\widehat{z},𝐮=\widehat{x}\mathrm{sin}\theta +\widehat{z}\mathrm{cos}\theta ,`$
$`𝐞^+=\widehat{x}\widehat{x}\widehat{y}\widehat{y},\mathrm{and}𝐞^\times =\widehat{x}\widehat{y}+\widehat{y}\widehat{x}.`$ (13)
The transfer function $`𝒯`$ approaches unity for $`ff_{}`$, and falls of as $`1/f`$ for $`ff_{}`$ by virtue of the sinc function. For a ground based detector such as LIGO the transfer function can be ignored since instrument noise keeps the operational range below LIGO’s transfer frequency of $`f_{}10^4`$ Hz (not to mention that LIGO is a Fabry-Perot interferometer so our results do not really apply). For a space based detector such as LISA, the instrument noise does not rise appreciably at high frequencies, and it is the transfer function that limits the high frequency response. The sinc function is responsible for the wiggly rise in the LISA sensitivity curve at high frequency that can be seen in Figure 4.

### 3.2 Interferometer response

We can form an interferometer by introducing a third spacecraft at a distance $`L`$ from the corner spacecraft and differencing the outputs of the two arms222More complicated differencing schemes have to be applied if the arms of the interferometer have unequal lengths in order to cancel laser phase noise.. The measured strain is given by

$`s(\widehat{\mathrm{\Omega }},f,𝐱,t)`$ $`=`$ $`{\displaystyle \frac{\delta \mathrm{}_𝐮(t)\delta \mathrm{}_𝐯(t)}{\mathrm{}}}`$ (14)
$`=`$ $`𝐃(\widehat{\mathrm{\Omega }},f):𝐡(\widehat{\mathrm{\Omega }},f,𝐱,t),`$
where

$$𝐃(\widehat{\mathrm{\Omega }},f)=\frac{1}{2}\left((𝐮𝐮)𝒯(𝐮\widehat{\mathrm{\Omega }},f)(𝐯𝐯)𝒯(𝐯\widehat{\mathrm{\Omega }},f)\right)$$
(15)
is the detector response tensor and $`𝐮`$ and $`𝐯`$ are unit vectors in the direction of each interferometer arm, directed out from the vertex of the interferometer. The above expression gives the response of the interferometer to a plane wave of frequency $`f`$ propagating in the $`\widehat{\mathrm{\Omega }}`$ direction. A general gravitational wave background can be expanded in the terms of plane waves:

$`h_{ij}(t,𝐱)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle 𝑑\widehat{\mathrm{\Omega }}\stackrel{~}{h}_{ij}(\widehat{\mathrm{\Omega }},f,𝐱,t)}`$ (16)
$`=`$ $`{\displaystyle \underset{A}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle 𝑑\widehat{\mathrm{\Omega }}\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})e^{2\pi ift}e^{2\pi if\widehat{\mathrm{\Omega }}𝐱/c}e_{ij}^A(\widehat{\mathrm{\Omega }})}.`$
Here $`𝑑\widehat{\mathrm{\Omega }}=_0^{2\pi }𝑑\varphi _0^\pi \mathrm{sin}\theta d\theta `$ denotes an all sky integral. The Fourier amplitudes obey $`\stackrel{~}{h}_A(f)=\stackrel{~}{h}_A^{}(f)`$ since the waves have real amplitudes. In the final expression we have chosen $`𝐞^+`$ and $`𝐞^\times `$ as basis tensors for the decomposition of the two independent polarizations. The response of the detector (which we are free to locate at $`𝐱=\mathrm{𝟎}`$) to a superposition of plane waves is then

$`s(t)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle 𝑑\widehat{\mathrm{\Omega }}s(\widehat{\mathrm{\Omega }},f,\mathrm{𝟎},t)}`$ (17)
$`=`$ $`{\displaystyle \underset{A}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle 𝑑\widehat{\mathrm{\Omega }}\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})e^{2\pi ift}𝐃(\widehat{\mathrm{\Omega }},f)}:𝐞^A(\widehat{\mathrm{\Omega }}).`$
From this we infer that the Fourier transform of $`s(t)`$ is given by

$$\stackrel{~}{s}(f)=\underset{A}{}𝑑\widehat{\mathrm{\Omega }}\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})𝐃(\widehat{\mathrm{\Omega }},f):𝐞^A(\widehat{\mathrm{\Omega }}).$$
(18)
In the low frequency limit the transfer function $`𝒯`$ approaches unity and $`𝐃(\widehat{\mathrm{\Omega }},f)`$ only depends on the geometry of the detector:

$$𝐃=\frac{1}{2}\left(𝐮𝐮𝐯𝐯\right).$$
(19)
To a first approximation, a stochastic gravitational wave background can be taken to be isotropic, stationary and unpolarized. It is fully specified by the ensemble averages:

$$\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})=0,\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})\stackrel{~}{h}_A^{}(f^{},\widehat{\mathrm{\Omega }}^{})=0,$$
(20)
and

$$\stackrel{~}{h}_A^{}(f,\widehat{\mathrm{\Omega }})\stackrel{~}{h}_A^{}(f^{},\widehat{\mathrm{\Omega }}^{})=\frac{1}{2}\delta (ff^{})\frac{\delta ^2(\widehat{\mathrm{\Omega }},\widehat{\mathrm{\Omega }}^{})}{4\pi }\delta _{AA^{}}S_h(f).$$
(21)
Here $`S_h(f)`$ is the spectral density of the stochastic background and the normalization is chosen such that

$`\stackrel{~}{h}_A^{}(f)\stackrel{~}{h}_A^{}(f^{})`$ $``$ $`{\displaystyle 𝑑\widehat{\mathrm{\Omega }}𝑑\widehat{\mathrm{\Omega }}^{}\stackrel{~}{h}_A^{}(f,\widehat{\mathrm{\Omega }})\stackrel{~}{h}_A^{}(f^{},\widehat{\mathrm{\Omega }}^{})}`$ (22)
$`=`$ $`{\displaystyle \frac{1}{2}}\delta (ff^{})\delta _{AA^{}}S_h(f).`$
The spectral density has dimension Hz-1 and satisfies $`S_h(f)=S_h(f)`$. It is related to $`\mathrm{\Omega }_{\mathrm{gw}}(f)`$ by

$$\mathrm{\Omega }_{\mathrm{gw}}(f)=\frac{4\pi ^2}{3H_0^2}f^3S_h(f).$$
(23)
The time averaged response333The limit $`\tau \mathrm{}`$ is approximated in practice by observing for a period much longer than the period of the gravitational wave. of the interferometer,

$$s(t)_\tau =\underset{\tau \mathrm{}}{lim}\frac{1}{\tau }_{\tau /2}^{\tau /2}s(t)𝑑t,$$
(24)
when applied to a stochastic background, is equivalent to the ensemble average

$`s(t)_\tau `$ $`=`$ $`s(t)`$ (25)
$`=`$ $`{\displaystyle \underset{A}{}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle 𝑑\widehat{\mathrm{\Omega }}\stackrel{~}{h}_A(f,\widehat{\mathrm{\Omega }})e^{2\pi ift}𝐃(\widehat{\mathrm{\Omega }},f)}:𝐞^A(\widehat{\mathrm{\Omega }})`$
$`=`$ $`0.`$
Since the expectation value of $`s(t)`$ vanishes we need to consider higher moments such as $`s^2(t)`$. Using equations (14) and (16) we find

$$s^2(t)=_0^{\mathrm{}}𝑑fS_h(f)(f),$$
(26)
where the transfer function $`(f)`$ is given by

$$(f)=\frac{d\widehat{\mathrm{\Omega }}}{4\pi }\underset{A}{}F^A(\widehat{\mathrm{\Omega }},f)F^A(\widehat{\mathrm{\Omega }},f)^{},$$
(27)
and $`F^A(\widehat{\mathrm{\Omega }},f)=𝐃(\widehat{\mathrm{\Omega }},f):𝐞^A(\widehat{\mathrm{\Omega }})`$ is the detector response function. In the low frequency limit, $`ff_{}`$, it is easy to show that $`(f)=2/5\mathrm{sin}^2\beta `$, where $`\mathrm{cos}\beta =𝐮𝐯`$ is the angle between the interferometer arms.

The response of the interferometer can be expressed in terms of the strain spectral density, $`\stackrel{~}{h}_s(f)`$, which has units of Hz-1/2 and is defined by

$$s^2(t)=_0^{\mathrm{}}𝑑f\stackrel{~}{h}_s^2(f)=_0^{\mathrm{}}𝑑fS_h(f)(f)$$
(28)
We see that the strain spectral density in the interferometer is related to the spectral density of the source by

$$\stackrel{~}{h}_s(f)=\sqrt{S_h(f)(f)}.$$
(29)
The total output of the interferometer is given by the sum of the signal and the noise:

$$S(t)=s(t)+n(t).$$
(30)
Assuming the noise is Gaussian, it can be fully characterized by the expectation values

$$\stackrel{~}{n}(f)=0,\mathrm{and}\stackrel{~}{n}^{}(f)\stackrel{~}{n}(f^{})=\frac{1}{2}\delta (ff^{})S_n(f),$$
(31)
where $`S_n(f)`$ is the noise spectral density. The total noise power in the interferometer is thus

$$n^2(t)=_0^{\mathrm{}}𝑑fS_n(f)=_0^{\mathrm{}}𝑑f\stackrel{~}{h}_n^2(f)$$
(32)
where $`\stackrel{~}{h}_n(f)`$ is the the strain spectral density due to the noise. Comparing equations (28) and (32), we define the signal to noise ratio at frequency $`f`$ by

$$\mathrm{SNR}(f)=\frac{\stackrel{~}{h}_s^2(f)}{\stackrel{~}{h}_n^2(f)}=\frac{S_h(f)(f)}{S_n(f)}.$$
(33)
Sensitivity curves for space-based interferometers typically display the effective strain noise

$$\stackrel{~}{h}_{\mathrm{eff}}(f)=\sqrt{\frac{S_n(f)}{(f)}}.$$
(34)
At high frequencies it is the decay of the transfer function $`(f)`$ leads to a rise in the effective noise floor. The actual noise power, $`S_n(f)`$, does not rise significantly at high frequencies for space based systems.

### 3.3 The LISA interferometer

As a concrete example of the general formalism described above we derive the sensitivity curve for the LISA interferometer using the design specifications quoted in the LISA Pre-Phase A Report. Using the coordinate system shown in Figure 2, the unit vectors along each arm are given by

$`𝐮={\displaystyle \frac{1}{2}}\widehat{x}+{\displaystyle \frac{\sqrt{3}}{2}}\widehat{y}`$
$`𝐯={\displaystyle \frac{1}{2}}\widehat{x}+{\displaystyle \frac{\sqrt{3}}{2}}\widehat{y},`$ (35)
and the gravitational wave is described by

$`\stackrel{~}{\mathrm{\Omega }}=\mathrm{cos}\varphi \mathrm{sin}\theta \widehat{x}+\mathrm{sin}\varphi \mathrm{sin}\theta \widehat{y}+\mathrm{cos}\theta \widehat{z}`$ (36)
$`\widehat{m}=\mathrm{sin}\varphi \widehat{x}\mathrm{cos}\varphi \widehat{y}`$ (37)
$`\widehat{n}=\mathrm{cos}\varphi \mathrm{cos}\theta \widehat{x}+\mathrm{sin}\varphi \mathrm{cos}\theta \widehat{y}\mathrm{sin}\theta \widehat{z}.`$ (38)
The angle between the interferometer arms is $`\mathrm{arccos}(𝐮𝐯)=\beta =\pi /3`$. The various ingredients we need to calculate $`(f)`$ are

$`𝐮\widehat{\mathrm{\Omega }}=\mathrm{sin}(\varphi +\pi /6)\mathrm{sin}\theta `$ (39)
$`𝐯\widehat{\mathrm{\Omega }}=\mathrm{sin}(\varphi \pi /6)\mathrm{sin}\theta `$ (40)
and

$$(𝐮𝐮):𝐞^+=\frac{1}{4}\mathrm{sin}^2\theta +\frac{1}{2}\mathrm{cos}^2\theta \mathrm{cos}^2\varphi \frac{\sqrt{3}}{4}\mathrm{sin}2\varphi (1+\mathrm{cos}^2\theta ).$$
(41)
$$(𝐯𝐯):𝐞^+=\frac{1}{4}\mathrm{sin}^2\theta +\frac{1}{2}\mathrm{cos}^2\theta \mathrm{cos}^2\varphi +\frac{\sqrt{3}}{4}\mathrm{sin}2\varphi (1+\mathrm{cos}^2\theta ).$$
(42)
$$(𝐮𝐮):𝐞^\times =\mathrm{cos}\theta \mathrm{sin}(2\varphi +\pi /3).$$
(43)
$$(𝐯𝐯):𝐞^\times =\mathrm{cos}\theta \mathrm{sin}(2\varphi \pi /3).$$
(44)
Before proceeding to find $`(f)`$ it is instructive to evaluate the detector response functions, $`F^A(\widehat{\mathrm{\Omega }},f)`$, at zero frequency:

$`F^+(\widehat{\mathrm{\Omega }})`$ $`=`$ $`{\displaystyle \frac{1}{2}}(𝐮𝐮𝐯𝐯):𝐞^+`$ (45)
$`=`$ $`{\displaystyle \frac{\sqrt{3}}{4}}\mathrm{sin}2\varphi (1+\mathrm{cos}^2\theta ),`$
and

$`F^\times (\widehat{\mathrm{\Omega }})`$ $`=`$ $`{\displaystyle \frac{1}{2}}(𝐮𝐮𝐯𝐯):𝐞^\times `$ (46)
$`=`$ $`{\displaystyle \frac{\sqrt{3}}{2}}\mathrm{cos}2\varphi \mathrm{cos}\theta ,`$
The magnitudes of these response functions are shown in Figure 3. They can be thought of as the polarization dependent antenna patterns for the LISA detector appropriate to a stochastic background of gravitational waves444These shapes are not coordinate independent. They depend on how we choose $`𝐮`$ and $`𝐯`$ in relation to $`𝐞^+`$ and $`𝐞^\times `$.

Inserting equations (39) through (44) into equation (27) gives the transfer function $`(f)`$ in terms of an integral over the angles $`\theta `$ and $`\varphi `$. Our explicit expression for LISA’s transfer function agrees with the one found by Larson, Hiscock and Hellings using an alternative approach. We were unable to perform the angular integral to arrive at a general form for $`(f)`$, but for low frequencies the integrand can be expanded in a series to give

$`(f)={\displaystyle \frac{2}{5}}\mathrm{sin}^2\beta (1{\displaystyle \frac{37\mathrm{cos}^2(\beta /2)10\mathrm{cos}^4(\beta /2)1}{84\mathrm{cos}^2(\beta /2)}}\left({\displaystyle \frac{f}{f_{}}}\right)^2`$
$`+{\displaystyle \frac{20\mathrm{cos}^6(\beta /2)+163\mathrm{cos}^2(\beta /2)80\mathrm{cos}^4(\beta /2)6}{2268\mathrm{cos}^2(\beta /2)}}\left({\displaystyle \frac{f}{f_{}}}\right)^4\mathrm{}).`$ (47)
At frequencies above $`f\frac{3}{2}f_{}`$, the sinc function takes over and the transfer function falls of as $`1/f^2`$. Setting $`\beta =\pi /3`$, we find that a good approximation for the transfer function is given by

$$(f)=\{\begin{array}{cc}\frac{3}{10}\left(1\frac{169}{504}\left(\frac{f}{f_{}}\right)^2+\frac{425}{9072}\left(\frac{f}{f_{}}\right)^4\frac{165073}{47900160}\left(\frac{f}{f_{}}\right)^6\right),\hfill & f<\frac{3}{2}f_{}\hfill \\ \frac{16783143}{126156800}\left(\frac{3f_{}}{2f}\right)^2,\hfill & f\frac{3}{2}f_{}\hfill \end{array}.$$
(48)
The coefficient in front of the high frequency term is chosen so that the transfer function is continuous at $`f=\frac{3}{2}f_{}`$.

Now that we have calculated the transfer function, our next task is to estimate the the detector noise power, $`S_n(f)`$. There are many noise contributions discussed in the LISA Pre-Phase A Report. The dominant ones are thought to be acceleration noise from the inertial sensors and position noise due to laser shot noise. A position noise of $`\stackrel{~}{\delta x}=2\times 10^{11}`$ m Hz-1/2 is quoted for each LISA spacecraft. There are two such contributions per arm, giving a total of 4 contributions for the interferometer. Since the contributions are uncorrelated they add in quadrature to give a total position noise of $`2\stackrel{~}{\delta x}`$. Dividing this by the optical path length of $`2L`$ and squaring gives the position noise power:

$$S_n^{\mathrm{pos}}(f)=\left(\frac{\stackrel{~}{\delta x}}{L}\right)^2.$$
(49)
An acceleration noise of $`\stackrel{~}{\delta a}=3\times 10^{15}`$ ms-2 Hz-1/2 is quoted for each inertial sensor. This noise acts coherently on the incoming and outgoing signal, for a combined acceleration noise of $`2\stackrel{~}{\delta a}`$ per spacecraft. There are four such contributions in the interferometer that add in quadrature for a total acceleration noise of $`4\stackrel{~}{\delta a}`$. Dividing this by the square of the angular frequency of the gravitational wave yields the effective position noise due to spurious accelerations. Dividing the effective position noise by the optical path length and squaring gives the acceleration noise power:

$$S_n^{\mathrm{accl}}(f)=\left(\frac{2\stackrel{~}{\delta a}}{(2\pi f)^2L}\right)^2.$$
(50)
Adding together the acceleration and position noise power gives the total noise power. Using the nominal values for the LISA mission we obtain

$$S_n(f)=\left(9.24\times 10^{40}\left(\frac{\mathrm{mHz}}{f}\right)^4+1.6\times 10^{41}\right)\mathrm{Hz}^1.$$
(51)
Inserting the above expression for $`S_n(f)`$, along with the approximate expression for the transfer function (48), into equation (34) yields a useful analytic approximation for the effective strain noise, $`\stackrel{~}{h}_{\mathrm{eff}}(f)`$, in the LISA interferometer. The analytic approximation is compared to the full numerical result in Figure 4.

LISA reaches a peak sensitivity in the frequency range $`3\times 10^310^2`$ Hz. Using equations (33) and (34) we have

$$S_h(f)=\stackrel{~}{h}_{\mathrm{eff}}^2(f)\mathrm{SNR}(f).$$
(52)
For a SNR of 2, this translates into a sensitivity of

$$\stackrel{~}{h}_{\mathrm{eff}}(f)=6.2\times 10^{21}\left(\frac{\mathrm{\Omega }_{\mathrm{gw}}(f)h_0^2}{10^{13}}\right)^{1/2}\left(\frac{f}{\mathrm{mHz}}\right)^{3/2}\mathrm{Hz}^{1/2}.$$
(53)
Thus, a stochastic background with $`\mathrm{\Omega }_{\mathrm{gw}}(f)h_0^2>7\times 10^{12}`$ should dominate LISA’s instrument noise for frequencies near 3 mHz. The difficulty would come in deciding if the interferometer response was due to instrument noise or a stochastic background, as both are Gaussian random processes. One way to be sure is to fly two LISA interferometers and cross correlate their outputs.

## 4 Cross correlating two interferometers

Having demonstrated that our formalism recovers the standard results for a single interferometer, we now investigate how the sensitivity can be improved by cross correlating two interferometers. The most general form for the cross correlation of two detectors is given by

$`C`$ $`=`$ $`{\displaystyle _{T/2}^{T/2}}𝑑t{\displaystyle _{T/2}^{T/2}}𝑑t^{}S_1(t)S_2(t^{})Q(tt^{})`$ (54)
$``$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f^{}\delta _T(ff^{})\stackrel{~}{S}_1^{}(f)\stackrel{~}{S}_2(f^{})\stackrel{~}{Q}(f^{}).`$
Here $`Q(tt^{})`$ is a filter function and the $`S_i`$’s are the strain amplitudes that we read out from the $`i^{\mathrm{th}}`$ detector:

$$S_i(t)=h_i^{\mathrm{astro}}(t)+s_i(t)+n_i(t).$$
(55)
The contributions are from resolvable astrophysical sources $`h_i^{\mathrm{astro}}`$, the stochastic background $`s_i`$, and intrinsic detector noise $`n_i`$. The function $`\delta _T`$ that appears in the Fourier space version of the correlation function is the “finite time delta function”

$$\delta _T=_{T/2}^{T/2}𝑑te^{2\pi ift}=\frac{\mathrm{sin}(\pi fT)}{\pi f}.$$
(56)
It obeys

$$\delta _T(0)=T\mathrm{and}\delta (x)=\underset{T\mathrm{}}{lim}\delta _T(x).$$
(57)
After removing resolvable astrophysical sources, the ensemble average of $`C`$ is given by

$$C=s_1,s_2+s_1,n_2+n_1,s_2+n_1,n_2=s_1,s_2.$$
(58)
Notice that terms involving the noise vanish as the noise is uncorrelated with the signal, and the noise in different detectors is also uncorrelated. Thus,

$$C=s_1,s_2=_{\mathrm{}}^{\mathrm{}}𝑑f_{\mathrm{}}^{\mathrm{}}𝑑f^{}\delta _T(ff^{})\stackrel{~}{s}_1^{}(f)\stackrel{~}{s}_2(f^{})\stackrel{~}{Q}(f^{}).$$
(59)
Our next task is to evaluate the quantity $`\stackrel{~}{s}_1^{}(f)\stackrel{~}{s}_2(f^{})`$ that appears in equation (59). Taking the Fourier transform of the expression in (14) and performing the ensemble average we find

$$\stackrel{~}{s}_1^{}(f)\stackrel{~}{s}_2(f^{})=\frac{1}{2}\delta (ff^{})S_h(f)\gamma (f)\frac{2}{5}\mathrm{sin}^2\beta ,$$
(60)
where $`\gamma (f)`$ is the overlap reduction function

$$\gamma (f)=\frac{5}{2\mathrm{sin}^2\beta }\frac{d\widehat{\mathrm{\Omega }}}{4\pi }\left(F_{1}^{+}{}_{}{}^{}(\widehat{\mathrm{\Omega }},f)F_2^+(\widehat{\mathrm{\Omega }},f)+F_{1}^{\times }{}_{}{}^{}(\widehat{\mathrm{\Omega }},f)F_2^\times (\widehat{\mathrm{\Omega }},f)\right)e^{2\pi if\widehat{\mathrm{\Omega }}(𝐱_1𝐱_2)/c}.$$
(61)
The normalization is chosen so that $`\gamma (0)=1`$. The function $`\gamma (f)`$ is Hermitian in the sense that $`\gamma _{12}(f)=\gamma _{21}^{}(f)`$. The inverse Fourier transform of $`\gamma (f)`$ is a real function since $`\gamma ^{}(f)=\gamma (f)`$. The term “overlap reduction function” refers to the fact that $`\gamma (f)`$ takes into account the misalignment and separation of the interferometers. In addition, $`\gamma (f)`$ contains the transfer functions for each detector, as can be seen by setting $`F_1=F_2`$ in (61) and comparing with (27):

$$\gamma _{1=2}(f)=\frac{5}{2\mathrm{sin}^2\beta }(f).$$
(62)
In other words, $`\gamma (f)`$ is an overlap reduction function and detector transfer function all rolled into one. Inserting equation (60) into equation (59) yields the expectation value of $`C`$:

$`C`$ $`=`$ $`{\displaystyle \frac{T}{5}}\mathrm{sin}^2\beta {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑fS_h(f)\gamma (f)\stackrel{~}{Q}(f)`$ (63)
$`=`$ $`{\displaystyle \frac{2T}{5}}\mathrm{sin}^2\beta {\displaystyle _0^{\mathrm{}}}𝑑fS_h(f)\mathrm{}[\gamma (f)\stackrel{~}{Q}(f)].`$
In the final line we have switched from the mathematically convenient range of integration $`f(\mathrm{},\mathrm{})`$ to the physically relevant range $`f[0,\mathrm{})`$.

The noise in a measurement of $`C`$ is given by $`N=CC`$, and the signal to noise ratio (squared) is given by

$$\mathrm{SNR}^2\frac{C^2}{N^2}=\frac{C^2}{C^2C^2}.$$
(64)
Our task is to find the filter function $`\stackrel{~}{Q}(f)`$ that maximizes this signal to noise ratio. A lengthy but straightforward calculation yields

$$N^2=\frac{T}{4}_{\mathrm{}}^{\mathrm{}}|\stackrel{~}{Q}(f)|^2M(f)𝑑f.$$
(65)
where

$`M(f)`$ $`=`$ $`S_{n1}(f)S_{n2}(f)+S_{n1}(f)S_h(f)_2(f)+S_{n2}(f)S_h(f)_1(f)`$ (66)
$`+S_h^2(f)\left(|\gamma (f)|^2+{\displaystyle \frac{_1(f)_2(f)}{(2/5\mathrm{sin}^2\beta )^2}}\right).`$
Here $`S_{n,i}`$ and $`_i`$ denote the spectral noise and the transfer function for the $`i^{\mathrm{th}}`$ interferometer. The square of the signal to noise ratio can be written as

$$\mathrm{SNR}^2=\frac{8T}{25}\mathrm{sin}^4\beta \frac{\{\stackrel{~}{P},\stackrel{~}{Q}\}^2}{\{\stackrel{~}{Q},\stackrel{~}{Q}\}},$$
(67)
where $`\{A,B\}`$ denotes the inner product

$$\{A,B\}=_{\mathrm{}}^{\mathrm{}}𝑑fA^{}(f)B(f)M(f),$$
(68)
and

$$\stackrel{~}{P}(f)=\frac{S_h(f)\gamma ^{}(f)}{M(f)}.$$
(69)
The signal to noise ratio is maximized by choosing the optimal filter “parallel” to $`\stackrel{~}{P}`$. Since the normalization of $`\stackrel{~}{Q}`$ drops out, we set $`\stackrel{~}{Q}(f)=\stackrel{~}{P}(f)`$. Using this filter, the optimal signal to noise ratio for the cross correlated interferometers is given by

$$\mathrm{SNR}^2=\frac{8T}{25}\mathrm{sin}^4\beta _0^{\mathrm{}}𝑑f\frac{|\gamma (f)|^2S_h^2(f)}{M(f)},$$
(70)
or, equivalently, as

$$\mathrm{SNR}^2=\frac{9H_0^4\mathrm{sin}^4\beta }{50\pi ^4}T_0^{\mathrm{}}𝑑f\frac{|\gamma (f)|^2\mathrm{\Omega }_{\mathrm{gw}}^2(f)}{f^6M(f)}.$$
(71)
In the limit that the noise power in each interferometer is very much larger than the signal power we find

$$\mathrm{SNR}^2=\frac{9H_0^4\mathrm{sin}^4\beta }{50\pi ^4}T_0^{\mathrm{}}𝑑f\frac{|\gamma (f)|^2\mathrm{\Omega }_{\mathrm{gw}}^2(f)}{f^6S_{n1}(f)S_{n2}(f)},$$
(72)
which recovers the expressions quoted by Flanagan and Allen when we set $`\beta =\pi /2`$. The factor of $`\mathrm{sin}^4\beta `$ was missed in these papers, as the implicit assumption that $`\beta =\pi /2`$ crept into supposedly general expressions. The main new ingredient in our expression are the transfer functions $`𝒯_i(f)`$ that reside in the overlap reduction function $`\gamma (f)`$. The transfer functions prevent us from performing the integral over the two-sphere in equation (61) in closed form except at low frequency. However, it is a simple matter to numerically evaluate $`\gamma (f)`$ for a given detector configuration.

If we are trying to detect a weak stochastic background with a noisy detector, then equation (72) is the appropriate expression to use for the SNR. However, there may be some range of frequencies where the signal power dominates the noise power in each detector. The contribution to the SNR from a clean frequency window of this sort is given by

$$\mathrm{SNR}^2(f,\mathrm{\Delta }f)=T\frac{8\mathrm{sin}^4\beta }{25}_{f\mathrm{\Delta }f/2}^{f+\mathrm{\Delta }f/2}𝑑f^{}\frac{|\gamma (f^{})|^2}{\left(|\gamma (f^{})|^2+_1(f^{})_2(f^{})/(2/5\mathrm{sin}^2\beta )^2\right)}.$$
(73)
The integrand is approximately equal to $`1/2`$ for all frequencies, so that

$$\mathrm{SNR}(f,\mathrm{\Delta }f)\frac{2}{5}\mathrm{sin}^2\beta \sqrt{T\mathrm{\Delta }f}.$$
(74)
This expression provides a useful lower bound for the SNR that can be achieved for detectors with a clean frequency window of width $`\mathrm{\Delta }f`$ centered at some frequency $`f`$. For example, the spectral density due to white dwarf binaries in our galaxy is expected to exceed LISA’s noise spectral density for frequencies in the range $`10^43\times 10^3`$ Hz. To detect this background at 90% confidence requires a signal to noise ratio of SNR=1.65, which can be achieved with less than three hours of integration time $`T`$.

### 4.1 Cross correlating two LISA interferometers

One possible modification to the current LISA proposal would be to fly six spacecraft instead of three, and thereby form two independent interferometers. The cost of doing this is considerably less than twice the cost of the current proposal, and it has the advantage of providing additional redundancy to the mission. How would we best use a pair of LISA interferometers? If our main concern is getting better positional information on bright astrophysical sources, then we would fly the two interferometers far apart, eg. with one leading and the other trailing the Earth. However, if we want to maximize the cross correlation then we need the interferometers to be coincident and coaligned. However, a configuration of this type is likely to share correlated noise sources, which would defeat the purpose of cross correlating the interferometers. A better choice is to use a configuration that is coaligned but not coincident. This can be done by placing six spacecraft at the corners of a regular hexagon as shown in Figure 5. Notice that interferometers have parallel arms and that the corner spacecraft are separated by the diameter of the circle.

Using the same coordinate system that we used earlier for a single interferometer, the unit vectors $`𝐮_i`$ and $`𝐯_i`$ along each interferometer arm are:

$`𝐮_1=𝐮_2={\displaystyle \frac{1}{2}}\widehat{x}+{\displaystyle \frac{\sqrt{3}}{2}}\widehat{y}`$
$`𝐯_1=𝐯_2={\displaystyle \frac{1}{2}}\widehat{x}+{\displaystyle \frac{\sqrt{3}}{2}}\widehat{y},`$ (75)
where the index $`i=1,2`$ labels the two interferometers. The displacement of the corner spacecraft is given by

$$𝐱_1𝐱_2=2R\widehat{y}$$
(76)
where $`R`$ is the radius of the orbit. The angle between the interferometer arms is $`\beta =\pi /3`$, and the length of each arm is $`L=\sqrt{3}R`$. The light crossing time between the two interferometers, $`2R/c`$, is almost equal to the light crossing time along each interferometer arm, $`\sqrt{3}R/c`$. Thus, the loss of sensitivity due to multiple wavelengths fitting between the interferometers occurs for frequencies near the transfer frequency $`f_{}`$ (the transfer frequency corresponds to wavelengths that fit inside the interferometer arms). For LISA the transfer frequency is $`f_{}=9.54\times 10^3`$ Hz.

The various ingredients we need to calculate $`\gamma (f)`$ follow from those given in equations (39) through (44). Putting everything together in (61) and working in the low frequency limit we find

$$\gamma (f)=1\frac{383}{504}\left(\frac{f}{f_{}}\right)^2+\frac{893}{3888}\left(\frac{f}{f_{}}\right)^4\frac{5414989}{143700480}\left(\frac{f}{f_{}}\right)^6+\mathrm{}$$
(77)
For frequencies above $`f_{}`$ the overlap reduction function decays as $`f^2`$. A numerically generated plot of $`\gamma (f)`$ is displayed in Figure 6. Scaling the signal to noise ratio (72) in units appropriate to a pair of LISA interferometers we have

$$\mathrm{SNR}^2=34.4\left(\frac{T}{\mathrm{year}}\right)_0^{\mathrm{}}\left(\frac{df}{\mathrm{mHz}}\right)|\gamma (f)|^2\left(\frac{\mathrm{\Omega }_{\mathrm{gw}}(f)h_0^2}{10^{15}}\right)^2\left(\frac{\mathrm{mHz}}{f}\right)^6\left(\frac{10^{41}\mathrm{Hz}^1}{S_n(f)}\right)^2.$$
(78)
Using equation (51) for $`S_n(f)`$ and assuming a scale invariant stochastic background yields a signal to noise of

$$\mathrm{SNR}=1.44\sqrt{\frac{T}{\mathrm{year}}}\left(\frac{\mathrm{\Omega }_{\mathrm{gw}}h_0^2}{10^{14}}\right).$$
(79)
This represents a 500 fold improvement on the sensitivity of a single LISA detector. Indeed, if it were not for the astrophysical foregrounds, a pair of LISA detectors would be well poised to detect a scale invariant gravitational wave background from inflation.

## 5 Missions to detect the CGB

We will work on the assumption that astrophysical foregrounds swamp the CGB for frequencies above a few $`\mu `$Hz, and design our missions accordingly. What we have in mind is a post-LISA mission based on the (by then) tried and tested LISA technology, but with some allowance for improvements in basic technologies such as the accelerometers. Since cross-correlating two interferometers at ultra-low frequencies does not buy us a major improvement in sensitivity, we need to start with a design that has excellent sensitivity at low frequencies. Basically this means building bigger interferometers with better accelerometers. To be concrete, we show the sensitivity curves for three generations of LISA missions in Figure 7. LISA I corresponds to the current LISA design with the spacecraft cart-wheeling about the Sun at 1 AU, separated by $`L=5\times 10^9`$ m. LISA II refers to a possible follow on mission with the three spacecraft evenly spaced around an orbit at 1 AU, so that the spacecraft are separated by $`L=\sqrt{3}`$ AU. The LISA II mission would use similar optics555The main difference is that we have to “lead our target” by a much larger amount for the LISA II mission. In other words, the angle, $`\mathrm{\Delta }\theta `$, between the received and transmitted laser beams is much larger for LISA II than for LISA I. A simple calculation yields $`\mathrm{\Delta }\theta =2v/c`$, where $`v`$ is the velocity around the circle shown in Figure 5. For LISA II $`v/c=(GM_{}/c^2/R)^{1/2}`$, while for LISA I $`v/c=2e(GM_{}/c^2/R)^{1/2}`$. Here $`e=0.01`$ is the eccentricity of the LISA I orbits and $`R=1`$ AU. This equates to lead of $`\mathrm{\Delta }\theta =4\times 10^6`$ radians for LISA I and $`\mathrm{\Delta }\theta =2\times 10^4`$ radians for LISA II. to LISA I (same laser power and telescope size), but allows for an order of magnitude improvement in accelerometer performance. LISA III is similar to LISA II, except that the constellation would orbit at 35 AU (between Neptune and Pluto) and the accelerometers would be improved by a further two orders of magnitude. The acceleration noise for LISA III would benefit from the three orders of magnitude reduction in solar radiation relative to LISAs I and II, but this would come at the cost of having to power the spacecraft using nuclear generators (RTGs).

For all three missions we see that the $`\mu `$Hz range lies well below the interferometer’s transfer frequency

$$f_{}=1.84\times 10^4\left(\frac{\sqrt{3}\mathrm{AU}}{L}\right)\mathrm{Hz}.$$
(80)
Indeed, it would take a mission orbiting at 180 AU ( somewhere inside the Kuiper belt) to achieve a transfer frequency of $`f_{}=10^6`$ Hz. In practical terms this allows us to work in the low frequency limit $`ff_{}`$ where it is easy to derive analytic expressions for the overlap reduction function $`\gamma (f)`$. Moreover, we need only consider acceleration noise when working below 1$`\mu `$Hz. All our calculations are based on the same hexagonal cross correlation shown in Figure 5 that we used for the LISA I cross correlation. Using equation (50) for the acceleration noise and ignoring the position noise yields

$$S_n(f)=3.42\times 10^{33}\left(\frac{\sqrt{3}\mathrm{AU}}{L}\right)^2\left(\frac{\mu \mathrm{Hz}}{f}\right)^4\left(\frac{\delta a}{3\times 10^{16}\mathrm{ms}^2}\right)^2\mathrm{Hz}^1.$$
(81)
In the limit that the noise power dominates the spectral density of the CGB we arrive at an estimate for the signal to noise ratio for frequencies below 1$`\mu `$Hz

$`\mathrm{SNR}^2`$ $`=`$ $`29.4\left({\displaystyle \frac{T}{\mathrm{year}}}\right)\left({\displaystyle \frac{L}{\sqrt{3}\mathrm{AU}}}\right)^4\left({\displaystyle \frac{3\times 10^{16}\mathrm{ms}^2}{\delta a}}\right)^4\times `$
$`{\displaystyle _0^f}|\gamma (f^{})|^2\left({\displaystyle \frac{\mathrm{\Omega }_{\mathrm{gw}}(f^{})h_0^2}{10^{14}}}\right)^2\left({\displaystyle \frac{f^{}}{\mu \mathrm{Hz}}}\right)^2\left({\displaystyle \frac{df^{}}{\mu \mathrm{Hz}}}\right).`$ (82)
The SNR is scaled against the LISA II specifications. Since $`f/f_{}1`$ it is a good approximation to set $`\gamma (f)=1`$ inside the integral. We model the CGB spectrum near 1$`\mu `$Hz by a simple power law:

$$\mathrm{\Omega }_{\mathrm{gw}}(f)=\mathrm{\Omega }_{\mu \mathrm{Hz}}\left(\frac{f}{\mu \mathrm{Hz}}\right)^\alpha ,$$
(83)
scaled relative to a reference value at 1$`\mu `$Hz. With these approximations we arrive at our final expression for the SNR for a pair of cross correlated interferometers in the ultra low frequency regime:

$`\mathrm{SNR}={\displaystyle \frac{3.13}{\sqrt{1+2\alpha /3}}}\left({\displaystyle \frac{T}{\mathrm{year}}}\right)^{1/2}\left({\displaystyle \frac{L}{\sqrt{3}\mathrm{AU}}}\right)^2\left({\displaystyle \frac{3\times 10^{16}\mathrm{ms}^2}{\delta a}}\right)^2\left({\displaystyle \frac{\mathrm{\Omega }_{\mu \mathrm{Hz}}h_0^2}{10^{14}}}\right)\left({\displaystyle \frac{f}{\mu \mathrm{Hz}}}\right)^{\alpha +3/2}`$
(84)
We can compare this to the SNR of a single interferometer in the ultra low frequency regime by combining equations (23), (33) and (81) to find

$$\mathrm{SNR}_s(f)=0.683\left(\frac{L}{\sqrt{3}\mathrm{AU}}\right)^2\left(\frac{3\times 10^{16}\mathrm{ms}^2}{\delta a}\right)^2\left(\frac{\mathrm{\Omega }_{\mu \mathrm{Hz}}h_0^2}{10^{14}}\right)\left(\frac{f}{\mu \mathrm{Hz}}\right)^{\alpha +1}.$$
(85)
We see that the individual and cross correlated sensitivities are related:

$$\mathrm{SNR}=\frac{0.82}{\sqrt{1+2\alpha /3}}(T\mathrm{\Delta }f)^{1/2}\mathrm{SNR}_s(f),\mathrm{for}f\mathrm{\Delta }f1\mu \mathrm{Hz}.$$
(86)
This supports our earlier assertion that cross correlating two interferometers improves on the sensitivity of a single interferometer by $`(T\mathrm{\Delta }f)^{1/2}`$ where $`\mathrm{\Delta }f`$ is approximately equal to the central observing frequency $`f`$.

Equation (5) tells us that a GABI detector built from two LISA II interferometers could detect the CGB at greater than 90% confidence with one year of data taking if $`\mathrm{\Omega }_{\mathrm{gw}}(1\mu \mathrm{Hz})h_0^2=10^{14}`$. The same equation also tells us what has to be done to achieve greater sensitivity. It is clear that increasing the duration of the mission is not the best answer as the SNR only improves as the square root of the observation time. In contrast, increasing the size of the interferometer or reducing the acceleration noise produces a quadratic increase in sensitivity. For example, a GABI detector built from a pair of LISA III interferometers could detect the CGB at 90% confidence for signals as small as $`\mathrm{\Omega }_{\mathrm{gw}}(1\mu \mathrm{Hz})h_0^2=4\times 10^{22}`$. In-fact, for signal strengths above this value the LISA III detectors are signal dominated and we need to use equation (74) instead of (5) to calculate the SNR.

We can also use equation (5) to determine how well a GABI mission can measure the CGB power spectrum. Suppose that we break up the frequency spectrum into bins of width $`\delta f`$. Nyquist’s theorem tells us that $`\delta f1/T`$, eg. for an observation time of one year the frequency resolution is $`3.17\times 10^8`$ Hz. Writing $`\delta f=n/T`$ where $`n1`$, and taking $`\delta ff`$, we find from (5) that in the frequency window $`(f\delta f/2,f+\delta f/2)`$ the SNR is

$$\mathrm{SNR}(f,\delta f)=0.965\sqrt{n}\left(\frac{L}{\sqrt{3}\mathrm{AU}}\right)^2\left(\frac{3\times 10^{16}\mathrm{ms}^2}{\delta a}\right)^2\left(\frac{\mathrm{\Omega }_{\mu \mathrm{Hz}}h_0^2}{10^{14}}\right)\left(\frac{f}{\mu \mathrm{Hz}}\right)^{\alpha +1}.$$
(87)
We see that the SNR in each frequency bin scales as $`f^{1+\alpha }`$. For a scale invariant spectrum ($`\alpha =0`$) this translates into poor performance at frequencies below 1 $`\mu `$Hz and limits the range over which we can measure the spectrum. For example, if the CGB has a scale invariant spectrum and an amplitude of $`\mathrm{\Omega }_{\mathrm{gw}}h_0^22\times 10^{13}`$, a pair of LISA II interferometers could measure the spectrum over the range $`10^710^6`$ Hz. For a pair of LISA III detectors the main limitation at low frequencies comes from Nyquist’s theorem. Since mission lifetimes are limited to tens of years, it will not be possible to measure the CGB spectrum much below $`10^8`$ Hz using space based interferometers. Indeed, it may be difficult to push much below $`10^7`$ Hz unless ways can be found to build detectors that are stable for many months. A pair of LISA III interferometers could measure the spectrum between $`10^810^6`$ Hz if $`\mathrm{\Omega }_{\mathrm{gw}}(10^8\mathrm{Hz})h_0^210^{18}`$. To measure the CGB spectrum below $`10^8`$ Hz requires a return, full circle, to the world of CMB physics. Detailed polarization measurements of the CMB can be used to infer the CGB power spectrum for frequencies in the range $`10^{18}10^{16}`$ Hz.

## 6 Discussion

The LISA follow-on missions we have described will be able to detect or place stringent bounds on the CGB amplitude and spectrum between $`10^8`$ and $`10^6`$ Hz. But how can we be sure that it is the CGB we have detected and not some unresolved astrophysical foreground? The answer can be found in the statistical character of the competing signals. Most early universe theories predict that the CGB is truly stochastic. In contrast, the astrophysical signal is only approximately stochastic, in a sense that can be made precise by appealing to the central limit theorem. The line in Figure 1 indicating the predicted amplitude of various astrophysical foregrounds corresponds to what is known as the confusion limit. It marks the amplitude at which we can expect to find, on average, one source per frequency bin666The confusion limit goes down when the frequency resolution goes up. Most plots of the confusion limit assume a one year observation period so the bins are $`3.17\times 10^8`$ Hz in width.. It is not until we reach signal strengths considerably smaller than the confusion limit that the astrophysical signal starts to look stochastic. And therein lies the answer to our question. A “stochastic” signal due to astrophysical sources will always have bright outliers that sit just above the amplitude of the stochastic signal. The number and distribution of the outliers can be predicted on statistical grounds. If we do not see the outliers, then we can safely conclude that we have detected the CGB.

## Acknowledgements

We thank Peter Bender and Raffaella Schneider for sharing their thoughts on possible astrophysical sources of gravitational waves.

## References