GRAVITATIONAL-WAVE EMISSION FROM COMPACT GALACTIC BINARIES

A subclass of DDWDs with periods of several minutes becomes semi-detached when angular-momentum loss due to gravitational radiation drives the two bodies sufficiently close that the lower-mass WD fills its Roche lobe and mass transfer begins, either via direct impact or via an accretion disk. The stability of mass transfer depends critically on the mass ratio of the two initial WDs (requiring typically q = M₂/M₁ ≲ 2/3) and on whether spin–orbit tidal coupling feeds angular momentum accreted by the primary star back to the orbit; the strength of this coupling is characterized by the synchronization timescale τ_s (Marsh et al. 2004).

Under optimistic assumptions, strong tidal coupling stabilizes mass transfer, favoring the formation of stable AM CVn systems and resulting in fewer mergers. By contrast, under pessimistic assumptions, there is no tidal coupling and many more DDWD systems merge. In the presence of a disk, the accreted material's angular momentum is mainly transferred back into the orbit by tidal forces at the outer edge of the disk. However, the disk's material at the inner edge still has sufficient angular momentum to spin up the accretor. Effects not considered here, such as magnetic braking (Farmer & Roelofs 2010), may also impact the formation of AM CVn systems and favor the merger of the WDs.

If stable mass transfer proceeds, then the donor WD loses mass and expands in size, while the system evolves to longer orbital periods. The minimum orbital period at which mass transfer can begin is ∼2 minutes for a 0.45 M_☉ and a 0.25 M_☉ binary. Maximum orbital periods for evolved systems in the Galaxy are estimated to be about 90 minutes; these correspond to very old systems whose orbits have slowly evolved over a Hubble time.

Surveys specifically targeted at low-mass WD (such as the ELM survey of WDs; Brown et al. 2010; Kilic et al. 2011a; Brown et al. 2011) have resulted in the discovery of numerous DDWD systems with q < 2/3, which are expected to get into contact within a Hubble time. Several of these newly discovered low-mass binaries will form stable mass-transferring AM CVn systems, even in the most pessimistic case of no tidal coupling (Brown et al. 2011).

In this channel, a WD accretes from an initially nondegenerate He star with mass 0.4–0.6 M_☉ before the He has been exhausted in the He-star core. He stars are produced from stars with masses greater than 2 M_☉ that lose their envelopes on the red giant branch, and typically have lifetimes comparable to their main-sequence progenitors, allowing for mass transfer before evolution to the WD stage. Mass transfer is stable when the mass ratio of the donor He star to the accretor WD is less than approximately 1.2 (Tutukov & Fedorova 1989; Ergma & Fedorova 1990; see also Nelemans et al. 2001a). Due to mass loss, He burning in the donor's core quenches and at masses ∼0.2 M_☉ the star becomes semi-degenerate. The minimum orbital period reached by these systems is ∼10 minutes.

It remains in question whether 0.1 M_☉ or 0.3 M_☉ of He accreted from an He star will destroy a CO WD by edge-lit detonations (henceforth denoted as ELDs; see Shen et al. 2010 and references therein). Nelemans et al. (2001a) distinguish between optimistic and pessimistic assumptions for this channel depending on whether double ELDs remove CO WD systems.

In analogy to traditional CVs, a so-called evolved CV comprises a main-sequence star which transfers mass to a WD or neutron star (Podsiadlowski et al. 2001). The difference is that mass transfer occurs as the donor star evolves off the hydrogen main sequence. However, it is thought that the evolved-CV channel contributes only marginally (<2%) to the AM CVn population (Nelemans et al. 2004; Roelofs et al. 2007b) for systems with periods less than 1500 s. Therefore, we do not consider this channel in this work, which focuses on the detectability of AM CVn systems by GW interferometers.

In our analysis, we considered several different assumptions on AM CVn formation, generated simulated binary catalogs and detector data sets for each, and estimated the resulting binary detection and subtraction prospects. Among the large number of cases that we investigated, here for conciseness we report on five scenarios that we believe are representative of the range of possibilities given the constraints of current observations. Two of the scenarios represent an upper bound (with maximally efficient formation channels) and conversely, a bound derived assuming inefficient formation channels on predictions from population-synthesis models (cases 1 and 2); two others attempt to reproduce SDSS-estimated AM CVn densities (cases 3 and 4 below); a fifth (case 5) explores a Galactic population with a larger scale height at times earlier than 7 Gyr ago. As discussed below, the modified Galactic-disk scenario (case 5) should be consistent with SDSS observations; however, this agreement is difficult to quantify as the SDSS-estimated densities given in Roelofs et al. (2007b) are derived assuming a standard disk model.

We generate DDWD and AM CVn catalogs with the SeBa population-synthesis code (Portegies Zwart & Verbunt 1996; Nelemans et al. 2001b, 2004; Toonen et al. 2011), with the following assumptions and parameters: (1) the time- and position-dependent star formation rate (SFR) based on Boissier & Prantzos 1999 (see Nelemans et al. 2004), (2) initial primary masses distributed with power-law index −2.5 (see Kroupa et al. 1993), (3) a flat initial mass-ratio distribution, (4) a flat distribution for the logarithm of the semimajor axis, up to a radius of 10⁶ R_☉, (5) eccentricities distributed as P(e)∝e², and (6) a 50% binary fraction in the initial population of main-sequence stars. We expect a factor of two difference in population-synthesis results obtained using different assumptions and initial parameters than those stated above (Portegies Zwart & Verbunt 1996; Nelemans et al. 2001b, 2004). Each system is described by seven parameters: the binary masses M₁ and M₂, the orbital period P_orb, the mass-transfer rate $\dot{M}$, the Galactic latitude l and longitude b, and the luminosity distance d to the solar system.

Following Nelemans et al. (2004), we distribute single and binary stars according to a Galactic model that comprises both disk and bulge components.⁶ The density of systems in the disk is

$$\begin{equation} \rho _\mathrm{disk} (R,Z) = \rho _{\rm BP} (R) \, \mathrm{sech}^2 \biggl (\frac{Z}{Z_{\mathrm{disk}}} \biggr) \, \mathrm{pc}^{-3}, \end{equation} \tag{ 3 }$$

where R is the cylindrical radius from the Galactic center, ρ_BP(R) results from integrating the time-dependent plane-projected SFR of Boissier & Prantzos (1999; see Equations (1) and (2) of Nelemans et al. 2004), Z is the perpendicular distance from the Galactic plane, and Z_disk is the disk's scale height. The density of systems in the spherical bulge is given by

$$\begin{equation} \rho _\mathrm{bulge} \, (r) \propto e^{-(r/r_\mathrm{bulge})^2} \,\mathrm{pc}^{-3}, \end{equation} \tag{ 4 }$$

where r is the spherical distance from the Galactic center, r_bulge is the characteristic radius of the bulge, and Equation (4) is normalized so that there are as many bulge systems as disk systems in the inner 3 kpc of the Galaxy, with a current total mass of 2.6 × 10⁶ M_☉.

For the first four scenarios, we set Z_disk = 300 pc and r_bulge = 0.5 kpc (Jurić et al. 2008; Nelemans et al. 2004; Binney & Tremaine 1987; Sofue et al. 2009); ρ_BP(R) is fit well by a density ∝exp (− R/R_disk) with R_disk ≈ 2.5 kpc for DDWD systems and R_disk ≈ 2.2–2.4 kpc for AM CVn systems. In the last scenario, we explore a modified model of the Galactic disk with a larger scale height for binaries that are produced at earlier times (>7 Gyr ago).

For all but one scenarios, we use the DDWD catalog analyzed in Nelemans et al. (2004), with local density 1.2 × 10⁻⁵ pc⁻³, and orbital periods as large as 5.5 hr (f_GW = 10⁻⁴ Hz). Binaries at longer periods are not typically detectable by GW instruments because detector noise is higher at low frequencies; they are therefore not considered further in this analysis. All the catalogs of AM CVn systems, by contrast, show the effects of our different assumptions on AM CVn formation:

• Case 1: upper bound (strong tidal coupling and inefficient ELD events). This catalog includes AM CVn systems that have formed via the DDWD channel operating with strong tidal coupling and the He-star–WD channel with few ELD events (see Nelemans et al. 2001b, 2004). Although the resulting AM CVn local density is higher (2.8 × 10⁻⁵ pc⁻³) than SDSS-estimated values, past investigations of Galactic WD binaries as LISA sources have all used variants of this scenario (see, e.g., Edlund et al. 2005; Timpano et al. 2006). Thus, this case provides a meaningful comparison between our study and previously published results, which may vary depending on their idealizations of "detection," "resolution," etc. (see Section 4.3).
• Case 2: inefficient formation channels (weak tidal coupling and efficient ELD events). This catalog includes AM CVn systems that have formed via the DDWD channel operating with minimal tidal coupling and the He-star–WD channel with many ELD events (see Nelemans et al. 2001b, 2004). The AM CVn local density is 6.1 × 10⁻⁶ pc⁻³, a factor of two greater than SDSS-estimated values (Roelofs et al. 2007a), and falls within the factor of two uncertainty of population-synthesis results.
• Case 3: mild tidal coupling and no ELD events. This catalog includes AM CVn systems that have formed via only the DDWD channel, where mild tidal coupling (τ_s = 2 years) between the WDs allows for stable mass transfer in up to 10% of systems compared to the optimistic channel of case 1. We assume that the He-star–WD channel is suppressed. (However, observations discussed in Roelofs et al. 2007a suggest the presence of hot donors, which imply either that the He-star–WD channel is active with core He burning having been quenched only recently, or that following the second common-envelope ejection, the initial DDWD orbital separation is so small that the donor WD cannot cool within the gravitational-radiation timescale.) The AM CVn local density is 1.1 × 10⁻⁶ pc⁻³.
• Case 4: fewer progenitor DDWDs. In this scenario, we explain the SDSS-estimated AM CVn densities by postulating a 1/5 reduction in the number of their DDWD progenitors. This reduction in DDWD systems is consistent with existing observations and with the factor of two uncertainty introduced in population-synthesis models by our limited understanding of binary evolution and Galactic structure. This catalog also includes AM CVn systems that have formed both via the DDWD channel with mild tidal coupling (τ_s = 1 year) and via the pessimistic He-star–WD channel (see Nelemans et al. 2001b, 2004); their densities are also reduced by 1/5 (across all periods). The DDWD and AM CVn local densities are 2.4 × 10⁻⁶ and 1.5 × 10⁻⁶ pc⁻³.
• Case 5: modified Galactic-disk model. In this scenario, we place older systems in a thicker disk. With the caveat of small-number statistics, AM CVn systems identified by emission lines in SDSS have inferred projected velocities >70 km s⁻¹, hinting that AM CVn systems with longer periods may originate from a disk with a thick component, or at least from a thin disk with a larger scale height. Indeed, such a distribution of AM CVns at birth is consistent with our current (if everchanging) understanding of the spatial and temporal evolution of thin and thick galaxy disks (Scannapieco et al. 2011; Jurić et al. 2008; Schönrich & Binney 2009a, 2009b; Siebert et al. 2011; Bovy et al. 2012). Observational constraints are unfortunately limited (see Napiwotzki 2009 for a discussion on the thin-disk, thick-disk, and halo populations of single WDs). As an illustrative example, we altered the disk component of Equation (3) to have Z_disk = 1250 pc for systems older than 7 Gyr, Z_disk = 300 pc for systems younger than 7 Gyr, and r_bulge = 0.5 kpc as for cases 1–4. This effectively places all systems older than 7 Gyr in a thick disk, which includes the great majority of longer period systems. In addition, we assumed that AM CVn systems only formed via the DDWD channel. The DDWD and AM CVn local densities are 1.2 × 10⁻⁵ and 3.8 × 10⁻⁶ pc⁻³.

For each scenario, Table 1 shows the numbers of AM CVn and DDWD systems in the catalogs, as well as their local density ρ₀ (i.e., the Galactic mid-plane density at 8.5 kpc from the Galactic center).

Compact Galactic binaries in the low-frequency band emit quasi-monochromatic GWs, with frequency drifts of (at most) several frequency bins for year-long observation. In this paper, we model compact-binary waveforms by including only the dominant quadrupolar emission (Peters & Mathews 1963), with instantaneous strain amplitude

$$\begin{equation} h = 10^{-21} \biggl (\frac{\mathcal {M}_z}{M_{\odot }} \biggr)^{5/3} \biggl (\frac{P_{\mathrm{orb}}}{\mathrm{hr}}\biggr)^{-2/3} \biggl (\frac{d}{\mathrm{kpc}} \biggr)^{-1}, \end{equation} \tag{ 5 }$$

where ${\mathcal M_z}= (1+z) \mu ^{3/5} M^{2/5}=(1+z){\mathcal M}$ is the redshifted chirp mass (with M the total mass, μ the reduced mass, $\mathcal M$ the chirp mass, and z the redshift), d is the luminosity distance to the source, and P_orb is the orbital period of the binary.

Figure 1 shows the number density as a function of h (the GW strain of Equation (5)), and f = 2/P_orb (the GW frequency), for the different DDWD and AM CVn catalogs of each of our scenarios. The colored dots indicate detected binaries (see Section 4.3): green shows DDWDs, red shows AM CVn systems (from the DDWD channel), and blue shows AM CVn systems (from the He-star–WD channel).

Zoom In Zoom Out Reset image size

Binaries with orbital periods shorter than approximately 20 minutes will appear as isolated signals in year-long data sets of space-based GW observatories, and will be detectable when the signals are sufficiently strong compared to instrument noise at the same frequency. By contrast, most binaries with longer orbital periods will not be detected individually, because they will be too dense in frequency space to tell apart. These will produce a noise-like confusion foreground that will affect the detection and parameter estimation of extra-Galactic GW sources such as massive black hole binaries. Nevertheless, the brightest binaries will be detected above the confusion foreground.

Looking at Figure 1, we see that DDWD systems have larger densities than AM CVn systems, and larger GW strains (because they have greater chirp masses). Thus, it is DDWD systems that are principally responsible for the shape and strength of the GW foreground from WD binaries (see Section 5.2). We note also that all AM CVn populations show an overall low-frequency cutoff at log₁₀f = −3.5, due to the age of the universe.

A simplified (but roughly correct) description of GW data analysis is that GW signals will be detectable when their signal-to-noise ratio (S/N) is greater than a detection threshold S/N_thr chosen to yield a minimal number of false alarms. S/N is the ratio of signal strength to rms instrument noise (plus confusion noise, if present) at the same frequency. For LISA-like interferometric detectors, the yearly orbit and rotation of the spacecraft constellation imprint the signals with frequency and amplitude modulations that encode the sky position of the source. Detector response is further affected in frequency-dependent fashion by Time-Delay Interferometry (TDI), the technique used to subtract the otherwise overwhelming laser frequency noise (thus, S/N would be evaluated for the TDI signal with respect to TDI noise). See Appendix B for an overview of LISA-like detectors, their response to GWs, and the theory of probabilistic signal detection in noise.

To estimate the number of Galactic-binary systems that we will be able to detect individually with our reference detectors, as well as the confusion noise from the residual unsubtracted systems, we performed simulations of detection and subtraction using a variant of the method introduced by (Timpano et al. 2006; see Section 5.1). To wit:

• 1.  
  For each of our scenarios, we compute the frequency-domain TDI response for each system in the corresponding catalog, using an optimized variant (Vallisneri 2011) of the Mock LISA Data Challenge fast-binary code (Cornish & Littenberg 2007); we then superpose the responses to produce a full year-long data set.
• 2.  
  We run through the catalog and evaluate the S/N of each system (see Appendix B) with respect to a background of instrumental noise plus confusion noise from all nearby systems, estimated as a running median of their power spectral density. We flag as detected the systems with $\mathrm{{\rm S/N}} > {\rm{S/N}}_{\rm thr}$.
• 3.  
  We subtract the signals detected in step 2 from the data set. We then go back to step 2 and iterate, using the resulting lower estimate of confusion noise.
• 4.  
  The process converges after a few iterations, with the number of detected systems decreasing better than exponentially. We then evaluate and fit the power spectral density from the residual confusion foreground.

This process represents a rather optimistic model of source detection and subtraction that neglects two important effects: source detection will be confused by the presence of nearby signals beyond what is predicted by the increased noise level, and subtraction will be degraded by the imperfect estimation of parameters in noise (see, e.g., Crowder & Cornish 2004). Nevertheless, we can adjust the level of optimism by varying our assumptions regarding the duration of observation T_obs, the detection threshold S/N_thr, and the number N_obs of TDI observables that are available (two, in effect, for the standard LISA configuration with five or six operating links; one for a two-arm configuration).

It turns out that our results for the number of detected sources and for the residual confusion foregrounds are fit remarkably well by power laws of a single effective-S/N parameter ρ_eff = S/N_thr(T_obs/yr)^−1/2(N_obs)^−1/2. As a baseline, we take a standard assumption, consistent with prior literature, of one year, one observable (the standard unequal-arm Michelson X), and S/N_thr = 5. A conservative assumption would correspond to taking one year, one observable, and S/N_th = 10, yielding ρ_eff = 10; an optimistic assumption to five years, two independent observables, and S/N_th = 5, yielding ρ_eff = 1.58. Results are not sensitive to details such as the choice of running-median window, and whether signals are subtracted immediately or after each full catalog run-through.

Tables 2 (for the 5 Mkm detector) and 3 (for the 1 Mkm detector) show the number of recovered systems of each type in our five scenarios. The numbers quoted correspond to simulation results for one year, one interferometric observable, and S/N_thr = 5; the power laws give approximate fits as functions of ρ_eff = S/N_thr(T_obs/yr)^−1/2(N_obs)^−1/2.

Consider first the number of detectable sources as seen by the longer, standard LISA-like detector. In all cases except the "fewer progenitors" scenario 4, the number of detected DDWDs remains comparable to prior optimistic projections (∼12, 000–14, 000). Even in case 4, where the DDWD systems in the catalog are reduced by a factor of five, detections drop only by half, thanks to a weaker confusion foreground. By contrast, the number of detected AM CVn varies considerably depending on their density and formation scenarios (∼15–2000). In cases 2 and 3, more DDWDs are detectable at relatively high frequencies because of the reduced confusion noise from fewer AM CVn.

Case 1 was chosen to be directly comparable to the standard population catalogs of earlier works (Nelemans et al. 2001a, 2001b, 2001c, 2004), but the claimed number of detections for AM CVn systems is significantly lower than those earlier estimates, derived with a simple resolvability criterion of a single source per frequency bin with instrument-noise S/N > 1 (Nelemans et al. 2001b, 2004). This is unsurprising, given that in our analysis we model detection by identifying and subtracting bright sources iteratively, and that we require a higher S/N, defined with respect to both instrument and confusion noise. Furthermore, in this paper, we use proper interferometric observables instead of "raw" GW strain.

Downscaling the detector's armlength to 1 Mkm approximately halves the number of DDWD systems detected for ρ_eff = 5, except for case 4, where the decrease is ∼3. The number of AM CVn system detections is likewise reduced by factors ∼3–4. Note however that the ρ_eff scalings of Tables 2 and 3 are slightly different. The ESA eLISA/NGO study (Amaro-Seoane et al. 2012) finds numbers compatible with ours, once we take into account our noise assumptions for the 1 Mkm detector, which are somewhat more favorable than NGO's.

In Appendix A, we apply a familiar analytical treatment of the compact-binary GW foreground to a broad range of Galactic models examined in the literature, and we conclude that the foreground depends strongly on the Galactic-binary merger rate and on their characteristic chirp masses, but only weakly on other parameters. Namely,

$$\begin{eqnarray} S_h(f) \, & \simeq & \, 1.9 \times 10^{-44} (f/\mathrm{Hz})^{-7/3} \, \mathrm{Hz}^{-1}\nonumber \\ &&\times \biggl (\frac{\mathcal {D}_{\rm char}}{6.4\ {\rm kpc}}\biggr)^{\!-2} \! \biggl (\frac{\mathcal {R}_{\rm gal}}{0.015\ {\rm yr}^{-1}}\biggr) \biggl (\frac{\mathcal {M}_{z, \mathrm{char}}}{0.35\, M_\odot }\biggr)^{5/3}, \quad \end{eqnarray} \tag{ 6 }$$

where $\mathcal {D}_\mathrm{char}$ is the characteristic distance from Earth to Galactic binaries (a weak function of the model), $\mathcal {R}_\mathrm{gal}$ is the binary merger rate in the Galaxy, and ${\cal M}_{z, {\rm char}}$ is the characteristic chirp mass $\langle {\cal M}_z^{5/3}\rangle ^{-5/3}$. The f^−7/3 slope describes the power emitted at GW frequency f by a single binary evolving due to GWs, and it applies to a stationary ensemble where new binaries are born with a constant rate (equal to $\mathcal {R}_\mathrm{gal}$) at frequencies below the detector bandwidth. Equation (6) is normalized to values representative of our assumptions; as shown in the left panel of Figure 3, it is in good agreement with the results of our simulations for cases 1–3 and 5, while the coefficient is reduced by a factor of five for case 4.

Zoom In Zoom Out Reset image size

The residual confusion foreground consists of the unresolvable DDWD and AM CVn systems that are left once relatively brighter or isolated sources have been subtracted iteratively from the original population catalog. The residual foregrounds found in our simulations are very similar among cases 1–3 and 5, while they are a factor of five weaker for case 4. We note that the overall number of DDWD systems in cases 1–3 and 5 differs by a factor of five, thereby confirming earlier work and our predictions that the DDWD population provides the dominant component to the strength of the confusion foregrounds (both unsubtracted and subtracted). As shown in the middle and right panels of Figure 3, for both the 5 Mkm and 1 Mkm detectors, the total equivalent strain noise (instrument noise plus confusion foreground) is fit well by the expression

$$\begin{equation} S^\mathrm{total}_h(f) = S_h^\mathrm{inst} + S^\mathrm{conf}_h \times \mathrm{tanh}^\alpha \bigl ({\textstyle \frac{1}{2}} \beta \, {d}N/df \bigr), \end{equation} \tag{ 7 }$$

where α ≃ 1, β ≃ 1 yr⁻¹,

$$\begin{equation} S_h^\mathrm{conf}(f) = 1.4 \times 10^{-45} (f / \mathrm{Hz} )^{-8/3} \, \mathrm{Hz}^{-1}, \end{equation} \tag{ 8 }$$

and where

$$\begin{equation} \frac{dN}{df} = 5 \times 10^{-3} (f / \mathrm{Hz} )^{-11/3} \, \mathrm{Hz}^{-1} \end{equation} \tag{ 9 }$$

is the frequency-domain density of binary systems, which modulates the transition from the confused to the resolved regime. Equation (7) is formulated in analogy to the standard expression Equation (C1) (see Barack & Cutler 2004a, 2004b), with a different functional form that fits our simulations better. In Appendix C, we give more precise α and β fits as functions of ρ_eff, and we provide expressions for the total TDI X noise and for the seasonal modulation of the GW foreground due to detector motion (Giampieri & Polnarev 1997; Seto 2004; Edlund et al. 2005).

More than two dozen papers have made estimates of the GW signal from individually detectable Galactic compact binaries and from the continuum of unresolvable sources at low frequencies (variously called "confusion background" or "Galactic foreground") for a LISA-like GW detector with 5 Mkm arms. An earlier informative review of the literature is given by Ruiter et al. (2010). We note that for a 1 Mkm detector, the foreground decreases substantially, to a level comparable to instrument noise (see Figure 3 in Section 5.2).

Table 4 compares a subset (large, but possibly not complete) of relevant published estimates, grouped by population model. For each entry, we show whether the estimates include resolved sources (RSs), the continuum background (CB), or both; we show also which classes of binaries are included. In Table 5, we provide relevant scale lengths for the Galactic disk and bulge distributions, where available; these are indicative of the range of assumed Galactic models, but they do not provide a complete picture since different density laws also appear. We do not give quantitative information on possible halo components, which are also frequently used. The major differences between the various estimates arise from two sets of assumptions:

• 1.  
  Assumptions about the relative numbers of systems in terms of total mass, mass ratio, and separation. These arise from different treatments of compact-binary progenitor evolution and of the stability of mass-transferring systems. Approaches include both analytic estimates and Monte Carlo population-synthesis calculations.
• 2.  
  Assumptions about Galactic structure and population space densities. These can affect estimates in two ways: first, the magnitude of GW signals will depend on the scale size assumed for the Galactic population; second, the "calibration" of the local number density of sources using observations of a particular source type will depend sensitively on the assumptions about Galactic structure.

Table 5. Estimates of GW from Populations of Compact Binaries—Spatial Distributions

Reference	Model	Zdisk	Rdisk	$\mathcal {D}_{\rm char}$ (Disk)	rbulge
		(kpc)	(kpc)	(kpc)	(kpc)
This work, Nelemans et al. (2004)	See Section 4.1	0.3	2.5a	6.4b	0.5
This work (thick disk)		1.25	2.5a	6.8b	0.5
Nelemans et al. (2001b)	$e^{-R/R_\mathrm{disk}}\,\mathrm{sech}^2(Z/Z_\mathrm{disk})$	0.2	2.5	5.3	...
Edlund et al. (2005)
Timpano et al. (2006)	See Hils & Bender (2000) and Nelemans et al. (2004)
Seto (2002)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.09	3.5	4.8	...
Bender & Hils (1997)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.09	3.5	4.8	...
Hils et al. (1990)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.09	3.5	4.8	...
Hils & Bender (2000)	ns	ns	ns	...	ns
Hils (1998)	ns	ns	ns	...	ns
Yu & Jeffery (2010)	$e^{-R/R_\mathrm{disk}}\,\mathrm{sech}^2(Z/Z_\mathrm{disk})$	0.352	2.5	5.6	0.5
	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	1.158	2.5	6.5	...
Liu et al. (2010); Liu (2009)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.09	2.5	5.0	...
Farmer & Phinney (2003)	...	...	...	...	...
Webbink & Han (1998)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.09	4.0	4.7	...
Ruiter et al. (2010)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.2	2.5	5.4	0.5
Ruiter et al. (2009)	(halo)
Benacquista et al. (2004)	$e^{-R/R_\mathrm{disk}}\,\mathrm{sech}^2(Z/Z_\mathrm{disk})$	0.2	2.5	5.3	...
Postnov & Prokhorov (1998)	$e^{-R/R_\mathrm{disk}} e^{-(Z/Z_\mathrm{disk})^2}$	4.2	5	7.2	...
Lipunov et al. (1995)	$e^{-R/R_\mathrm{disk}} e^{-Z/Z_\mathrm{disk}}$	0.46	4.6	5.7	See paper
Lipunov & Postnov (1987)	$e^{-R/R_\mathrm{disk}} e^{-(Z/Z_\mathrm{disk})^2}$	4.2	10	9.6	...
Lipunov et al. (1987)
Giampieri & Polnarev (1997)	$e^{-R/R_\mathrm{disk}} e^{-(Z/Z_\mathrm{disk})^2}$	0.1	3.5	4.5	ns
Evans et al. (1987)	Constant density	0.2	10		...

Notes. ns: not specified. Distance to the Galactic center is assumed to be 8.5 kpc for all models. Exponential squared bulge unless otherwise specified: $e^{-(r/r_{\rm bulge})^2}$ with $\mathcal {D}_{\rm char}({\rm bulge}) = 8.4 \, \mathrm{kpc}$. ^aNumerical fit, DDWD systems only. ^b$\mathcal {D}_{\rm char}$ computed for disk + bulge for this case.

Download table as:  ASCIITypeset image

We now present a semi-analytical derivation of the amplitude and shape of the GW continuum. Similar models have been used extensively in GW data analysis, see, for example, works by Hughes (2002) and Barack & Cutler (2004a, 2004b). Our results show that the magnitude and spectral shape of the continuum are a robust function of the merger rate and characteristic chirp mass of Galactic-binary systems, and depend only weakly on Galactic structure.

We follow the arguments presented in Evans et al. (1987) and Giampieri & Polnarev (1997). The GW luminosity of a circular binary and the change in its orbital energy as a function of GW frequency f = 2f_orb are given by (Peters & Mathews 1963)

$$\begin{eqnarray} && {dE \over dt} = {32 \over 5} \pi ^{10/3} {G^{7/3} \over c^5} {{\cal M}_z}^{10/3} f^{10/3},\nonumber\\ && {dE \over df} = {1 \over 3} \pi ^{2/3} G^{2/3} {{\cal M}_z}^{5/3} f^{-1/3}, \end{eqnarray} \tag{ A1 }$$

where E is the energy of the system and f is the GW frequency. Combining these expressions yields the frequency evolution of the system,

$$\begin{equation} {df \over dt} = {96 \over 5} \pi ^{8/3} {G^{5/3} \over c^5} {{\cal M}_z}^{5/3} f^{11/3}. \end{equation} \tag{ A2 }$$

If $\mathcal {R}_{\rm birth}(f,{\bf r})$ is the birth rate per unit frequency and volume of systems with f_orb < f/2, then in a stationary regime the density of sources must be

$$\begin{eqnarray} {d^2 N \over df d {\bf r}} & = & {\mathcal {R}_{\rm birth}(f,{\bf r})\over {df/dt}}\nonumber\\ & = & {5 \over 96} \pi ^{-8/3} {c^5 \over G^{5/3}} {{\cal M}_z}^{-5/3}\ \mathcal {R}_{\rm birth}(f,{\bf r})f^{-11/3}. \end{eqnarray} \tag{ A3 }$$

By continuity, $\mathcal {R}_{\rm birth}(f,{\bf r})$ must also be equal to the merger rate, $\mathcal {R}_{\rm mrg}(f,{\bf r})$, of systems which merge at f_orb > f/2. DDWD systems, which generate most of the GW foreground, merge at f > 10⁻² Hz, and most of them are born below a few 10⁻⁴ Hz (although estimates vary). Thus, if we restrict Equation (A3) between those two frequencies, then we may just set $\mathcal {R}_{\rm birth}(f,{\bf r})= \mathcal {R}_{\rm mrg}({\bf r})$, defined as the total merger rate above 10⁻² Hz.

Now, the direction-averaged GW strain from a circular binary is given by

$$\begin{eqnarray} h^2 & =& \big\langle h_+^2 + h_\times ^2 \big\rangle = {1 \over \pi ^2} {1 \over d^2} {G \over c^3} f^{-2} {dE \over dt}\nonumber\\ & = & {32 \over 5} \pi ^{4/3} {G^{10/3} \over c^8} {{\cal M}_z}^{10/3} f^{4/3}, \end{eqnarray} \tag{ A4 }$$

see for instance (Giampieri & Polnarev 1997), where d is the luminosity distance to the binary; the corresponding detector response is

$$\begin{equation} h_I^2= \big\langle h_+^2 F_+^2 + h_\times ^2 F_\times ^2 \big\rangle \equiv h^2 F^2 ({\bf \hat{r}},t), \end{equation} \tag{ A5 }$$

where F₊ and F_× are the antenna-pattern functions (see, e.g., Maggiore 2007) and $F^2 ({\bf \hat{r}},t)$ is the polarization- and inclination-averaged antenna pattern (〈F²₊〉 + 〈F_ײ〉)/2 in the direction ${\bf \hat{r}}$.

We combine Equations (A3)–(A5) to obtain the one-sided detector response density $S_{h_I}$. The yearly orbital motion of LISA-like detectors modulates the amplitude of signals to individual sources; this effect averages out for isotropic source distributions, but not for the anisotropic Galactic-binary foreground (Giampieri & Polnarev 1997; Edlund et al. 2005). We therefore average over a year (while assuming for now that all binaries have the same chirp mass ${\cal M}_z$), to obtain

$$\begin{eqnarray} S_{h_I}(f) &=& \biggl \langle \int h^2(d({\bf r}), f)\ F^2 ({\bf \hat{r}},t)\ {d^2 N \over df d {\bf r}}\ d{\bf r} \biggr \rangle _t \nonumber \\ &= & {1 \over 3} \pi ^{-4/3} {G^{5/3} \over c^3} {{\cal M}_z}^{5/3} f^{-7/3} \!\! \int {{\langle F^2 ({\bf \hat{r}},t)\rangle _t} \mathcal {R}_{\rm mrg}({\bf r})\over d({\bf r})^2} d{\bf r} \nonumber \\ & \equiv & {1 \over 3} \pi ^{-4/3} {G^{5/3} \over c^3} {{\cal M}_z}^{5/3} f^{-7/3}\ \langle F^2\rangle _{{\bf \hat{r}},t}\langle \mathcal {R}_\mathrm{mrg} / d^2 \rangle _{\bf r}. \end{eqnarray} \tag{ A6 }$$

Of the two quantities defined in the last line of Equation (A6), $\langle F^2\rangle _{{\bf \hat{r}},t}$ accounts for the time-, merger-rate- and direction-averaged instrument response, and it depends on the orbit and geometry of the detector. For an isotropic distribution of sources and a triangular interferometer with 60° angles, $F^2 ({\bf \hat{r}},t)= 1/5 \times \sin ^2(60^\circ) = 3/20$, which is often taken as the average isotropic response for LISA-like detectors (Giampieri & Polnarev 1997). The factor 1/5 is the average of (〈F²₊〉 + 〈F_ײ〉)/2 over source locations and polarizations (again, e.g., in Maggiore 2007). Again for a LISA-like detector, but for typical disk populations, $F^2 ({\bf \hat{r}},t)$ is found empirically to vary by ∼10% from 3/20.

The quantity $\langle \mathcal {R}_\mathrm{mrg}/{d^2} \rangle _{\bf r}$ is the merger-rate-averaged inverse squared distance. For a Galactic disk, it depends on the scale height and radial scale length, and of course on the position of the detector within the disk. A typical disk density distribution considered in many papers is the squared-hyperbolic-secant disk model,

$$\begin{equation} \mathcal {R}_{\rm mrg}({\bf r})= \frac{\mathcal {R}_{\rm gal}}{V_{\rm gal}} \, e^{-R/R_\mathrm{disk}} {\rm sech}^2(Z/Z_\mathrm{disk}), \end{equation} \tag{ A7 }$$

where $\mathcal {R}_{\rm gal}$ is the total merger rate in the Galaxy and V_gal is the volume of the Galaxy. We find empirically that $\langle \mathcal {R}_\mathrm{mrg}/{d^2} \rangle _{\bf r}$ varies slowly for 2.5 kpc ⩽ R_disk ⩽ 10 kpc and 0.2 kpc ⩽ Z_disk ⩽ 1 kpc. Table 5 summarizes the wide range of Galactic models used in the literature. For most disk models, the characteristic distance to Galactic sources $\mathcal {D}_{\rm char}\equiv (V_{\rm gal}/\,\mathcal {R}_{\rm gal})^{-1/2} \langle \mathcal {R}_\mathrm{mrg} / d^2 \rangle ^{-1/2}$ is approximately 5–6 kpc, and $\mathcal {D}_{\rm char}+ R_\mathrm{disk}$ is close to 8.5 kpc, the assumed distance to the Galactic center.

Inserting numerical values appropriate for our simulations (cases 1–3) in Equation (A6) yields

$$\begin{eqnarray} S_{h_I}(f) & \simeq & 1.9 \times 10^{-44} (f/\mathrm{Hz})^{-7/3} \, \mathrm{Hz}^{-1} \times \langle F^2\rangle _{{\bf \hat{r}},t}\nonumber\\ &&\times \biggl (\frac{\mathcal {D}_{\rm char}}{6.4\ {\rm kpc}}\biggr)^{\!-2} \! \biggl (\frac{\mathcal {R}_{\rm gal}}{0.015\ {\rm yr}^{-1}}\biggr) \biggl (\frac{\mathcal {M}_{z, \mathrm{char}}}{0.35\, M_\odot }\biggr)^{5/3},\qquad \end{eqnarray} \tag{ A8 }$$

where ${\cal M}_{z, {\rm char}}$ is the characteristic chirp mass $\langle {\cal M}_z^{5/3}\rangle ^{-5/3}$. The power spectral density S_h(f) of the GW strain is obtained by dropping the factor $\langle F^2\rangle _{{\bf \hat{r}},t}$ (which is close to 3/20 for LISA-like detector). Equation (A8) is in good agreement with the results of our simulations. Note however that the spectra of Figure 2 are plotted for the TDI observable X (rather than for strain h or strain response h_I), so they must be compared with

$$\begin{equation} S_X(f) = \frac{16 \, x^2 \sin ^2 x}{1 + 0.6 \, x^2} S_{h_I}(f), \end{equation} \tag{ A9 }$$

where x = 2πfL/c, and L is the detector armlength.

Our analysis indicates that $S_{h_I}(f)$ and S_h(f) are relatively insensitive to the structure of the Galactic disk, consistent with the results of our simulations, and instead depend strongly on the merger rate $\mathcal {R}_\mathrm{gal}$ and the characteristic chirp mass $\mathcal {M}_{z, \rm {char}}$. By contrast, we expect that the local population density of sources, as well as the number of resolvable sources as a function of their GW strain, will depend roughly linearly on the disk's scale height.

The transition between the confused and resolved detection regimes is characterized by dN/df, which is determined by the evolution of the chirp frequency and by the absolute merger-rate density of sources. Integrating Equation (A3) over $\mathbf {r}$, we find

$$\begin{eqnarray} \frac{dN}{df} & = & 5 \times 10^{-3} \, (f/\mathrm{Hz})^{-11/3} \,{\rm Hz}^{-1} \left({\mathcal {R}_{\rm gal}\over {0.015\ {\rm yr}^{-1}}}\right)\nonumber\\ &&\times \left({M_{z, {\rm char}} }\over {0.35\, M_\odot }\right)^{-5/3}, \end{eqnarray} \tag{ A10 }$$

consistent with our simulations.

This simple analytical derivation of the GW continuum shows that the major differences between the various estimates listed in Table 4 arise primarily from two sources: the total merger rate of systems $\mathcal {R}_{\rm gal}$ and the composition of the source population, leading to a different characteristic chirp mass M_z, char. $\mathcal {D}_{\rm char}$ is also relevant, but less important. (Remember however that our derivation is only predictive for frequencies high enough that most binaries form at longer periods—for LISA-like detectors, this is true above 10⁻⁴ Hz. For lower frequencies, the shape of the continuum depends on the detailed accounting of birth rate versus period.)

LISA (Prince et al. 2009; Jennrich 2009), the space-based GW observatory studied jointly by NASA and ESA until early 2011, comprises three spacecraft orbiting the Sun in a quasi-equilateral configuration with armlengths ∼5 × 10⁶ km. LISA measures GW strain by monitoring the oscillating distance between freely falling test masses, using heterodyne laser interferometry. The spacecraft hover around the test masses to protect them from external perturbations, performing slight orbital corrections with extremely precise micro-Newton thrusters.

Although LISA is sometimes described as a large Michelson interferometer in space, the lasers exchanged by spacecraft are neither split nor reflected. Rather, two lasers on each spacecraft are used to establish six one-way links; each link provides a test-mass-referenced interferometric measurement y_ij of the phase (or equivalently, the frequency) of the incoming laser, compared to the local, outgoing beam. In terms of fractional laser frequency, the effect of a plane GW on the links is given by the two-pulse response (Estabrook & Wahlquist 1975)

$$\begin{equation} y_{ij}(t) = \frac{1}{2} \frac{\hat{n}_{ij} \cdot \bigl [{\hbox{\sf h}}(p_i,t-L_{ij}) - {\hbox{\sf h}}(p_j,t)\bigr ] \cdot \hat{n}_{ij}}{1 - \hat{n}_{ij} \cdot \hat{k}}, \end{equation} \tag{ B1 }$$

where $\hat{n}^i_{ij}$ is the unit vector from spacecraft i to j, $\hat{k}$ is the GW propagation vector, ${\hbox{\sf h}}(p_k,t)$ is the TT-gauge GW strain (see, e.g., Maggiore 2007) at the spacecraft position p_k and at time t, and L_ij is the photon time of flight along the arm. Thus, a GW pulse is registered twice in y_ij: once when it impinges on spacecraft j, and a time L_ij after it has impinged on spacecraft i. In the limit of long GW wavelengths λ ≫ L, the y_ij response to plane GWs is therefore proportional to L.

Each y_ij measurement is also affected by displacement noise from residual test-mass accelerations, position noise from optical path and shot noise in the heterodyne measurements, and, most importantly, by oscillations in the central laser frequency, which are many orders of magnitude stronger than expected GW imprints, and which appear as y_ij = C_ji(t) − C_ij(t − L_ij), with C_ij as the frequency noises for the six lasers. Laser frequency noise is removed by combining multiple one-way measurements with appropriate time delays so that each laser appears in canceling pairs, as in the unequal-arm Michelson observable

$$\begin{eqnarray} X(t) & = & [ y_{31}(t) + y_{13}(t - L_{31}) + y_{21}(t - L_{31} - L_{13})\nonumber\\ && + y_{12}(t - L_{31} - L_{13} - L_{21}) ]\nonumber \\ && - [ y_{21}(t) + y_{12}(t - L_{21}) + y_{31}(t - L_{21} - L_{12})\nonumber\\ && + y_{13}(t - L_{21} - L_{12} - L_{31}) ]. \end{eqnarray} \tag{ B2 }$$

This noise-suppression technique, specific to LISA-like interferometers, is known as TDI (see, e.g., Vallisneri 2005 for a review). TDI observables can be seen as synthesized interferometers, since their y_ij components trace the paths (and reproduce the phase delays) of light across recognizable interferometric configurations (hence X's "Michelson" designation). Many different observables are possible with LISA's six laser links, but all can be reconstructed from a set of three observables, such as the unequal-arm Michelson X, Y, and Z centered on the three spacecraft (Vallisneri et al. 2008).

In this paper, we assume an rms displacement noise

$$\begin{eqnarray} S^{1/2}_\mathrm{pm}(f) & = & 3 \times 10^{-15} \, \mathrm{m} \, \mathrm{s}^{-2}\, \mathrm{Hz}^{-1/2} \left(1 + \left (\frac{f}{10^{-4} \,\mathrm{Hz}}\right)^{-1}\right)^{1/2}\nonumber\\ && \times \left(1 + \left (\frac{f}{0.008 \,\mathrm{Hz}}\right)^{4}\right)^{1/2}, \end{eqnarray} \tag{ B3 }$$

for each of the proof masses, and a total rms position noise

$$\begin{equation} S^{1/2}_\mathrm{op}(f) = 18 \times 10^{-12} \, \mathrm{m}\, \mathrm{Hz}^{-1/2} \left(1 + \left (\frac{f}{0.002 \, \mathrm{Hz}}\right)^{-4}\right)^{1/2} \end{equation} \tag{ B4 }$$

for each y_ij measurement.

The GW mission currently under study as an ESA-led project (Amaro-Seoane et al. 2012) would decrease cost primarily by reducing spacecraft and propellant mass, allowing for launch on smaller rockets than envisaged for LISA. Propellant is saved by placing the spacecraft closer together and closer to Earth. At the low-frequency end, the shorter armlengths reduce the response to GWs proportionally; for the same displacement noise, sensitivity also decreases by the same ratio. At the high-frequency end, the laser power available for position measurement increases as L⁻², since beams are broadly defocused at millions of kilometers, improving shot noise (but not other optical noises) by an rms factor L⁻¹. As a consequence, the sweet spot of the LISA sensitivity shifts to higher frequencies.

Spacecraft mass may also be saved by adopting a two-arm "mother–daughters" configuration with only four laser links, two each between the middle "mother" spacecraft and one of the "daughters." With four links, only one independent TDI observable can be formed, with a loss in sensitivity averaging $\sqrt{2}$ at compact-binary frequencies. Because no trade studies have yet been performed to optimize the payload for the new orbits, in this paper we consider a "short LISA" configuration with aggressively reduced 1 × 10⁶ km armlengths, and with the LISA displacement and position noises (except for the power rescaling of the position-noise component, which we assume as five times smaller in rms). The actual sensitivity of a cost-constrained space-based GW observatory is likely to fall somewhere between this short LISA and the classic configuration.

Searches for GWs of known shape proceed by correlating detector data with theoretical waveform templates computed for the expected range of source parameters. The statistical confidence of detection for a GW signal h immersed in Gaussian additive noise n is described by the S/N of the data s = h + n after filtering by the best-fitting normalized template $\hat{t}$. For a single TDI observable X, this is given by (Cutler & Flanagan 1994)

$$\begin{equation} \rho (s;\hat{t}) = 4 \, \mathrm{Re} \int _0^{\infty } \frac{\tilde{X}_s(f)^* \tilde{X}_{\hat{t}}(f)}{S_X(f)} \,df, \end{equation} \tag{ B5 }$$

where $\tilde{X}_s(f)$ and $\tilde{X}_{\hat{t}}(f)$ are the Fourier transforms of the TDI data and of the TDI response to waveform $\hat{t}$, respectively; "*" denotes the complex conjugate; and S_X(f) is the noise power spectral density of TDI observable X. Normalization amounts to requiring $\rho (\hat{t};\hat{t}) = 1$.

With this noise weighting, the false-alarm probability of S/N exceeding a chosen detection threshold ρ_thr is proportional to exp (− ρ²_thr/2). Because the many templates used in typical searches amount to repeated independent trials, ρ_thr is routinely set to relatively large values (e.g., 5–10) to yield one false alarm per many years of data. The detectability of a signal h is characterized by its optimal S/N ρ(h; h), which is the S/N that would be obtained on average over noise realizations for a perfectly matching template.

The noise S_X(f) includes the confusion-noise foreground S^conf_X(f) from unresolved sources (see Section 5), as well as instrument noise S^inst_X. For the unequal-arm Michelson observable X in the (rather reasonable) limit of equal interferometer armlengths, S^inst_X is given by (Vallisneri et al. 2008)

$$\begin{equation} S^\mathrm{inst}_X = 16 \sin ^2(2 \pi L f) [2\bigl (1 + \cos ^2(2 \pi L f)\bigr) S^{\Delta f/f}_\mathrm{pm}(f) + S^{\Delta f/f}_\mathrm{op}(f) ]. \end{equation} \tag{ B6 }$$

In these units S^Δf/f_pm = S_pm(f)/(2π c f), S^Δf/f_op(f) = S_op(f) × (2πf/c).

To obtain the TDI response to a GW signal (5), we insert Equation (5), multiplied by the appropriate phasing function and sky-position-dependent polarization tensors, into Equation (B2). The response is modulated further by the heliocentric motion of the detector (which yields a time-dependent Doppler shift) and its yearly rotation in a plane inclined 60° from the ecliptic (which generates sidebands at frequency multiples of 1 yr⁻¹).⁷ Both effects imprint sky-position and polarization angles into the TDI response; this structure makes it possible to resolve sufficiently strong GW sources that overlap in frequency space, and to determine their approximate sky position (see, e.g., Królak et al. 2004). The TDI response to a population of sources is given simply by the linear superposition of the responses to individual sources.

Barack and Cutler (2004a, 2004b) codified the lore about the total noise of a LISA-like detector (i.e., the instrument noise plus Galactic-binary confusion noise) by providing a heuristic derivation for the well-known expression

$$\begin{eqnarray} && S^\mathrm{total,BC}_h(f) = \mathrm{min} \big [ S_h^\mathrm{inst}(f) e^{\kappa T^{-1} dN/df}{,}\;\; S_h^\mathrm{inst}(f) + S_h^\mathrm{conf,BC}(f) \big ], \nonumber \\ \end{eqnarray} \tag{ C1 }$$

where they estimate

$$\begin{equation} S_h^\mathrm{conf,BC}(f) = 1.4 \times 10^{-44} (f / \mathrm{Hz} )^{-7/3} \, \mathrm{Hz}^{-1}, \end{equation} \tag{ C2 }$$

$$\begin{equation} \frac{dN}{df} = 2 \times 10^{-3} (f / \mathrm{Hz} )^{-11/3} \, \mathrm{Hz}^{-1}, \end{equation} \tag{ C3 }$$

and κT⁻¹ ≃ 1.5 yr⁻¹. As we discuss in Appendix A, under generic assumptions, the unsubtracted Galactic-binary foreground will follow a f^−7/3 spectrum; this is confirmed in our simulations, with a coefficient of 2.6 × 10⁻⁴⁴ for cases 1–3 and 5, and five times smaller for case 4 (see Figure 3). After the subtraction of detectable sources, the total noise found in our simulations is fit better by a smaller amplitude, a slightly different exponent, and a smoother transition function:

$$\begin{equation} S^\mathrm{total,sim}_h(f) = S_h^\mathrm{inst} + S^\mathrm{conf,sim}_h \times \mathrm{tanh}^\alpha \bigl ({\textstyle \frac{1}{2}} \beta \, dN/df \bigr), \end{equation} \tag{ C4 }$$

with

$$\begin{equation} S^\mathrm{conf,sim}_h(f) = 1.4 \times 10^{-45} (f / \mathrm{Hz} )^{-8/3} \, \mathrm{Hz}^{-1}, \end{equation} \tag{ C5 }$$

$$\begin{equation} \frac{dN}{df} = 5 \times 10^{-3} (f / \mathrm{Hz} )^{-11/3} \, \mathrm{Hz}^{-1}. \end{equation} \tag{ C6 }$$

For cases 1–3 and 5, and for our reference subtraction run (one year, one interferometric observable, ρ_thr = 5), taking α = 1, β = 1 yr⁻¹ yields a good fit for a 5 Mkm LISA-like detector, while α = 0.7, β = 1.2 yr⁻¹ is appropriate for 1 Mkm arms. For the entire range of ρ_eff that we probed, good fits are obtained by setting

$$\begin{eqnarray} && \alpha = 2 - 1.2 \, (\rho _\mathrm{eff} / 5) + 0.2 \, (\rho _\mathrm{eff} / 5)^2,\nonumber\\ && \beta = [-0.2 + 1.5 \, (\rho _\mathrm{eff} / 5) - 0.3 \, (\rho _\mathrm{eff} / 5)^2] \ \mathrm{yr}^{-1} \quad (5 \, \mathrm{Mkm}),\nonumber \\ && \alpha = 1.6 - 1.2 \, (\rho _\mathrm{eff} / 5) + 0.3 \, (\rho _\mathrm{eff} / 5)^2, \nonumber\\ && \beta = [0.3 + 1.3 \, (\rho _\mathrm{eff} / 5) - 0.4 \, (\rho _\mathrm{eff} / 5)^2] \ \mathrm{yr}^{-1} \quad (1 \, \mathrm{Mkm}). \qquad \end{eqnarray} \tag{ C7 }$$

For case 4, an acceptable fit follows from reducing S^{conf, sim}_h by a factor of five.

The equivalent confusion noise in the TDI X combination (expressed as fractional-frequency fluctuations) is given by

$$\begin{equation} S^\mathrm{total,sim}_\mathrm{X}(f) = \frac{3}{20} \frac{16 \, x^2 \sin ^2 x}{1 + 0.6 \, x^2} S^\mathrm{total,sim}_h(f), \;\; \mathrm{with} \;\; x = 2 \pi f L / c. \end{equation} \tag{ C8 }$$

Here, 3/20 is the all-sky averaging $\langle F^2\rangle _{{\bf \hat{r}},t}$ discussed above, 1 + 0.6 x² is an approximation to the frequency-dependent interferometer response function, and 16 x²sin ²x is the TDI X transfer function. (In Barack and Cutler's convention (2004a), which we follow here, the average response function is included in the definition of the equivalent strain noise S^total_h, so that the sky-averaged S/N can be computed directly as $4 \, \mathrm{Re} \int |\tilde{h}|^2/S^\mathrm{total}_h(f)\, df$; therefore S^conf_h, rather than $S^\mathrm{conf}_{h_I}$, appears in S^total_h. By contrast, TDI observables already include the position-dependent instrument response, so it is appropriate to define the TDI noise S^total_X without the sky averaging.)

Due to the evolving orientation of the LISA-like formation with respect to the Galaxy, the Galactic foreground (before and after subtraction) has a strong seasonal variation in amplitude, almost constant across frequencies (Giampieri & Polnarev 1997; Seto 2004; Edlund et al. 2005). For a single TDI X combination, for both our detector configurations, for all scenarios, and for our assumptions about the relative orientation of detector and Galaxy at t = 0 (consistent with the standard MLDC LISA orbits, Arnaud et al. 2006), this variation is fit well by

$$\begin{eqnarray} S_h^\mathrm{conf,sim}(f,t) & = & S_h^\mathrm{conf,sim}(f) \times [ 1 - 0.21 \cos (2 \pi \, t/\mathrm{yr})\nonumber \\ & & - 0.45 \cos (4 \pi \, t/\mathrm{yr}) - 0.05 \sin (2 \pi \, t/\mathrm{yr})\nonumber \\ & & + 0.01 \sin (4 \pi \, t/\mathrm{yr}) ], \end{eqnarray} \tag{ C9 }$$

where S^{conf, sim}_h(f, t) denotes the power spectral density computed for data segments shorter than one year, but longer than the typical measurement timescale of a few hours. An equation analogous to Equation (C9), with the same coefficients, holds also for the unsubtracted background S_h(f, t).