LIGO GRAVITATIONAL WAVE DETECTION, PRIMORDIAL BLACK HOLES, AND THE NEAR-IR COSMIC INFRARED BACKGROUND ANISOTROPIES

LIGO's recent discovery of the gravitational wave (GW) from an inspiraling binary black hole (BH) system of essentially equal mass BHs (∼30M_⊙) at z ∼ 0.1 (Abbott et al. 2016b) has led to the suggestion that all or at least a significant part of the dark matter (DM) is made up of primordial BHs (PBHs; Bird et al. 2016; Clesse & García-Bellido 2016). In particular, Bird et al. (2016) argue that this PBH mass range is not ruled out by astronomical observations and the observed rate at ∼(a few) Gpc⁻³ yr⁻¹ (Abbott et al. 2016a) can be accounted for if DM PBHs are distributed in dense, low velocity-dispersion concentrations that escaped merging. There is abundant motivation for PBHs forming in the very early universe (Hawking 1971; Carr & Hawking 1974), e.g., during phase transition at the QCD epoch when horizon mass is of the right magnitude (Jedamzik 1997); see a nice overview of possible mechanisms in Mack et al. (2007). If PBHs indeed constitute the bulk or all of DM, they would contribute an additional Poissonian component to the power spectrum of the mass distribution from adiabatic fluctuations from the earlier inflationary era (Afshordi et al. 2003). If so, this component would dominate small scales leading to significant modification of the history of collapse (and possibly formation of the first sources), resulting in greater rates of cosmic infrared background (CIB) production at z > 10.

CIB contains emissions produced over the entire history of the universe including from sources inaccessible to direct telescopic studies (see the review by Kashlinsky 2005) with early stars and BHs contributing to its near-IR (∼1–5 μm) fluctuation component (Cooray et al. 2004; Kashlinsky et al. 2004). In this context, Kashlinsky (2005) and Kashlinsky et al. (2007b) have identified from deep Spitzer data significant CIB fluctuations remaining on sub-degree scales after subtracting individual galaxies to faint levels. The measurement was later extended to degree scales (Kashlinsky et al. 2012) and confirmed in subsequent analyses of Akari and Spitzer data (Matsumoto et al. 2011; Cooray et al. 2012b). It is now established that these fluctuations cannot arise from remaining known galaxy populations (Kashlinsky et al. 2005; Helgason et al. 2012) and it has been suggested that they arise from new populations at early epochs (Kashlinsky et al. 2005, 2007c; Yue et al. 2013a). This proposition is currently a subject of debate (Cooray et al. 2012b; Gong et al. 2015), although the CIB fluctuations on relevant scales appear to be uncorrelated with the diffuse light in the visible produced by sources down to AB mag >28 (Kashlinsky et al. 2007a).³ It was established that the CIB fluctuations are coherent with unresolved soft X-ray cosmic X-ray background (CXB) (Cappelluti et al. 2013) at levels greater than expected from remaining known populations (Helgason et al. 2014) and that the measured coherence levels require much higher proportions of BHs among the CIB sources than in the known populations. At the same time, Helgason et al. (2016) have argued that if early populations were to produce the measured CIB signal that would require higher than expected efficiencies of early star formation (cf. Kashlinsky et al. 2015b).

In this Letter, we point out that if indeed the LIGO discovery is indicative of PBHs making up the DM, the extra Poissonian isocurvature component of the fluctuations would lead to much greater rates of collapse at early times, which would naturally produce the observed levels of the CIB fluctuations. We briefly revisit the required near-IR CIB energetics in Section 2 and the effects of the extra power component from PBHs on the collapse of the first halos in Section 3. We discuss the effects PBHs and the extra power may have on the formation of first populations in Section 4. The discussion below adopts cosmology with (h, Ω_tot, Ω_CDM, Ω_bar, σ₈) = (0.7, 1, 0.23, 0.05, 0.9).

The observed CIB fluctuations reflect several aspects of the sources producing them: (1) the measured shot-noise power characterizes the typical flux-magnitude of the sources producing the large-scale (clustering) power; the fact that the arcminute fluctuations arise at very low shot-noise levels means that the individual sources must be very faint consistent with their high-z origin (Kashlinsky et al. 2007c). (2) The angular shape of the CIB power spectrum on arcminute scales reflects the epochs spanned by the sources; the shape from the Spitzer data is consistent with high-z origin within the current errors, and the upcoming Euclid all-sky survey would further probe the epochs and history of emissions much more accurately (Kashlinsky et al. 2015a). (3) Given the angular power template, the amplitude of the fluctuations spectrum at some fiducial scale (we adopt 5' below) reflects the overall abundance of the sources with fluxes constrained by (1) and (2) via the corresponding mean CIB flux.

The integrated/bolometric CIB flux levels depend on three efficiency parameters: (1) the efficiency of collapse of halos suitable for forming luminous sources, or the mass-fraction of the universe in these halos, denoted f_Halo; (2) the formation efficiency of conversion of baryons inside each halo into luminous sources, f_*; and (3) the radiation efficiency of converting the rest mass into radiation for the luminous sources inside the collapsed halos, .

We now briefly revisit the arguments in Kashlinsky et al. (2015b, their Section 2) for a general set of efficiency requirements for sources at high z to reproduce the observed CIB anisotropies between 2 and 5 μm.⁴ The integrated CIB fluctuation at 5' between 2 and 5 μm from the AKARI to Spitzer bands is $\delta {F}_{2\mbox{--}5\mu {\rm{m}}}({5}^{\prime })\simeq 0.09\ \mathrm{nW}\ {{\rm{m}}}^{-2}\ {{\rm{sr}}}^{-1};$ this arises as excess over that from known galaxies remaining in the data (Helgason et al. 2016). Populations at high z are strongly biased, span a short period of cosmic time, and would be expected to produce CIB with a relative fluctuation amplitude of Δ_5' ∼ 10% on arcminute scales, which would then require producing ${F}_{{\rm{CIB}}}(2\mbox{--}5\;\mu {\rm{m}})\;=\delta {F}_{2\mbox{--}5\mu {\rm{m}}}({5}^{\prime })$ $/{{\rm{\Delta }}}_{{5}^{\prime }}\sim 1\ \mathrm{nW}\ {{\rm{m}}}^{-2}\ {\mathrm{sr}}^{-1}$ (Kashlinsky et al. 2007c) in the integrated flux at near-IR wavelengths ∼(2–5) μm. The Lyman cutoff would cut the emissions below the observer wavelengths ∼(1 + z)/10 μm.

Let us assume that a fraction f_Halo of all matter in the universe collapses, in halos capable of producing luminous sources, at a given redshift converting on average a fraction f_* of the halo baryons into luminous sources. The bolometric diffuse flux produced by these populations, after they have converted their mass–energy into radiation with radiation efficiency , is ${F}_{{\rm{tot}}}\simeq {f}_{{\rm{Halo}}}{f}_{\ast }\left(\tfrac{c}{4\pi }\epsilon {\rho }_{{\rm{bar}}}{c}^{2}\right){z}_{{\rm{eff}}}^{-1}$ $\simeq 9.1\times {10}^{5}\epsilon \;{f}_{{\rm{Halo}}}$ ${f}_{\ast }{z}_{{\rm{eff}}}^{-1}$ $\tfrac{{{\rm{\Omega }}}_{{\rm{bar}}}{h}^{2}}{0.0227}\;\mathrm{nW}\ {{\rm{m}}}^{-2}\ {{\rm{sr}}}^{-1}$ where ${z}_{{\rm{eff}}}\equiv 1/\langle {(1+z)}^{-1}\rangle$ is a suitably averaged effective redshift factor that accounts for the radiation energy density decreasing with expansion as ∝(1 + z)⁻⁴ versus the matter density ∝(1 + z)⁻³. The overall fraction of universe's baryons needed to explain the CIB is f_Halof_* (see Section 2.3.2 in Kashlinsky et al. 2015b). For massive stars, which are fully convective and radiate close to the Eddington limit, ≃ 0.007 for the H-burning phase of a few megayears per star. Accreting BHs can reach electromagnetic radiation efficiencies = 0.4 for a maximally rotating Kerr BH. If the integrated CIB fluctuation approximates the bolometric flux produced by these sources, the fraction of baryons that on average go into the sources inside each halo is

$$\begin{eqnarray}{f}_{\ast } & = & 0.1\times {\left(\displaystyle \frac{{f}_{{\rm{Halo}}}}{0.01}\right)}^{-1}\ {\left(\displaystyle \frac{\epsilon }{0.01}\right)}^{-1}\ \left(\displaystyle \frac{{z}_{{\rm{eff}}}}{10}\right){\left(\displaystyle \frac{{{\rm{\Delta }}}_{{5}^{\prime }}}{0.1}\right)}^{-1}\\ & & \times {\left[\displaystyle \frac{{F}_{{\rm{CIB}}}(2\mbox{--}5\mu {\rm{m}})}{{F}_{{\rm{tot}}}}\right]}^{-1}.\end{eqnarray} \tag{ 1 }$$

Thus, in order to produce the measure CIB at z > 10 with "reasonable" formation efficiencies (f_* < 10%), one requires a large fraction of matter in collapsed halos capable of producing luminous sources (see the next section).

Helgason et al. (2016) discuss the requirements of high-z sources to produce the observed CIB fluctuations within the conventional, if necessarily simplified, framework of gravitational clustering and spherical collapse of adiabatic ΛCDM fluctuations. They conclude that (1) first galaxies, if extrapolated to z > 8 from known UV luminosity functions would produce much less CIB fluctuation power than observed (cf. Cooray et al. 2012a; Yue et al. 2013b), and (2) at still higher z (first) stars would have to form inside the collapsed halos at substantial formation efficiencies (converting f_* ≳ 5% of the available baryons in collapsing halos) and be very massive (∼500M_⊙) if they are to explain the observed CIB anisotropies. Kashlinsky et al. (2015a) reproduce the observed Spitzer signal with massive early stars forming at the mean formation efficiency f_* ≃ 4% out to z = 10.

The "high-mean-formation-efficiency" difficulty can ultimately be traced to a relative paucity of high-z collapsed halos—with the parameters appropriate for star formation—due to the limited amount of power on the relevant scales set by the adiabatic ΛCDM matter fluctuations, which arose from the period of inflation. The next section discusses how the abundances of the high-z collapsed halos are dramatically increased if PBHs constitute the DM, and reduce—by large factors—the efficiencies required to produce the observed CIB anisotropies.

LIGO's GW 150914 originated at z = 0.09 from the merger of two BHs of essentially identical masses at ${36}_{-4}^{+5}$ and 29 ± 4 M_⊙ (Abbott et al. 2016b). This mass range lies within the horizon mass-scale at ∼0.01–0.1 Gev where various mechanisms for generating PBHs in the very early universe operate, such as discussed by, e.g., Jedamzik (1997). Bird et al. (2016) discuss how the observed detection rate, inferred from the so far single published event, can be made consistent with that expected from the PBHs making up the DM such that their comoving mean mass density, assumed constant since their formation until at least their possible later evolution (discussed in Section 4 and references therein), is given by

$$\begin{eqnarray}{n}_{{\rm{PBH}}} & = & \displaystyle \frac{1}{{M}_{{\rm{PBH}}}}{{\rm{\Omega }}}_{{\rm{CDM}}}\displaystyle \frac{3{H}_{0}^{2}}{8\pi G}\\ & \simeq & {10}^{9}{\left(\displaystyle \frac{{M}_{{\rm{PBH}}}}{30{M}_{\odot }}\right)}^{-1}\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{CDM}}}{h}^{2}}{0.1}\right)\;{{\rm{Mpc}}}^{-3}.\end{eqnarray} \tag{ 2 }$$

Below we will assume, for simplicity, that all PBHs have identical mass. The arguments that follow can be generalized to a PBH mass distribution, such as in, e.g., Carr (1975) and Choptuik (1993), with M_PBH being the effective mass leading to the overall n_PBH comoving number density. We note that this mass range is allowed by, although close to, the limits from the MACHO microlensing surveys (Alcock et al. 2001). Ricotti (2007) and Ricotti et al. (2008) have argued that if PBHs are of this mass range, accretion onto them may violate COBE/FIRAS constraints on the CMB blackbody energy spectrum, but as Bird et al. (2016) discuss, such arguments are model-dependent and subject to complex physics assumptions. Afshordi et al. (2003) limit M_PBH < 4 × 10⁴M_⊙ from Lyα forest data.

As pointed out by Afshordi et al. (2003), the DM from PBHs will contain an extra (isocurvature) component due to Poissonian fluctuations with the power component at the time of the PBH formation being ${P}_{{\rm{PBH,}}{\rm{initial}}}={n}_{\mathrm{PBH}}^{-1}$ in comoving units. From their formation to today (z = 0), these isocurvature fluctuations would grow, at wavelengths below the horizon at matter-radiation equality z_eq, by a scale-independent factor of $\tfrac{3}{2}(1+{z}_{{\rm{eq}}})$, so the extra power component at redshift z is given by (Afshordi et al. 2003)

$$\begin{eqnarray}{P}_{{\rm{PBH}}}(z) & = & \displaystyle \frac{9}{4}{(1+{z}_{{\rm{eq}}})}^{2}{n}_{{\rm{PBH}}}^{-1}{[g(z)]}^{-2}\\ & \simeq & 2\times {10}^{-2}\left(\displaystyle \frac{{M}_{{\rm{PBH}}}}{30{M}_{\odot }}\right)\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{CDM}}}{h}^{2}}{0.13}\right)\left(\displaystyle \frac{1}{{g}^{2}(z)}\right)\ {{\rm{Mpc}}}^{3}\end{eqnarray} \tag{ 3 }$$

where g(z) is the linear growth factor of fluctuations from z to today, with g(0) = 1. Figure 1 shows the extra power component for M_PBH = 30M_⊙ compared to the ΛCDM power spectrum from the purely adiabatic fluctuation component. The power is plotted versus the mass contained in wavelength 2π/k, which is $M{(r)=1.15\times {10}^{12}(r/1{\rm{Mpc}})}^{3}{M}_{\odot }$ for the adopted cosmological parameters. This extra power is ∝ M_PBH and for M_PBH > 1M_⊙ dominates the small scales relevant for collapse of the first halos at z > 10. This isocurvature power component dominates very small scales and has no impact on the observed CMB anisotropies or baryonic-acoustic oscillations (Eisenstein & Hu 1999) that appear in CIB fluctuations on arcminute scales and can be probed with Lyman-tomography of CIB from the upcoming Euclid survey (Kashlinsky et al. 2015a). Furthermore, unlike the part of the power from clustering, white-noise power contributions to the angular CIB power spectrum are not affected by the biasing amplification (Kashlinsky et al. 2004).

Zoom In Zoom Out Reset image size

The net power spectrum would be given by ${P}_{{\rm{tot}}}(k,z)={P}_{{\rm{\Lambda }}{\rm{CDM}}}(k,z)+{P}_{{\rm{PBH}}}(z)$, which we use to evaluate the rms density contrast at z over a sphere of comoving radius r_M containing mass M(r_M) as ${\sigma }_{M}(z)={\left[\tfrac{1}{2{\pi }^{2}}\int {P}_{{\rm{tot}}}(k,z){W}_{{\rm{TH}}}({{kr}}_{M}){k}^{2}{dk}\right]}^{1/2}$ after normalizing to σ₈ over 8h⁻¹ Mpc at z = 0. (W_TH is the top-hat function.) Assuming spherical collapse, masses with density contrast >δ_col = 1.68 at that epoch will have collapsed by z. In general, collapse and subsequent formation of compact objects is driven by balance between pressure and gravity, which in turn is determined by cooling in the collapsing gas. Two modes of halo collapse are relevant here: if enough H₂ forms, the gas will have T ≃ 10³ K, and in the absence of metals, the H cooling will in any event keep the gas at T ≃ 10⁴ K (see the review by Bromm & Larson 2004). Figure 2 shows the resultant rms density fluctuation versus mass at various z relevant here for M_PBH = (0, 15, 30) M_⊙ with the vertical lines demarcating where the halo virial temperatures exceed these limits. The strong increase in the rms density contrast, over that in the absence of the PBHs, at masses of the first halos capable of producing luminous objects, is obvious. This increase will lead to substantially more collapsed halos capable of forming luminous sources at z > 10.

Zoom In Zoom Out Reset image size

We use the Press–Schechter formalism (Press & Schechter 1974) to compute the fraction of collapsed halos as the probability of a density field region with virial temperature T_vir having overdensity >δ_col. For Gaussian-distributed density fluctuations, the fraction of the halos that collapsed by redshift z is

$$\begin{eqnarray}&&{f}_{{\rm{Halo}}}(M({T}_{{\rm{vir}}}),z)=\displaystyle \frac{1}{2}{\rm{erfc}}\left(\displaystyle \frac{{\delta }_{{\rm{col}}}}{\sqrt{2}{\sigma }_{M({T}_{{\rm{vir}}})}(z)}\right).\end{eqnarray} \tag{ 4 }$$

Figure 3 shows the fraction in halos that collapsed by z having T_vir > 10⁴ (right) and 10³ K for M_PBH = (0, 15, 30) M_⊙. The increase in f_Halo is large enough to produce the required CIB flux with very modest baryon conversion efficiencies (Equation (1)) of well below ∼1% for the H-burning = 0.007 even by z ∼ 20 in halos with H₂-cooling. Even in halos with T_vir > 10⁴ K, the required f_* remains at a modest few percent level at z ≃ 12–15 for M_PBH = 30M_⊙. If the bulk of the CIB comes from BH accretion, the values of the required f_* drop by over an order of magnitude. Thus, to account for the observed near-IR CIB fluctuation signal with high-z emissions, very few baryons would need to be converted into luminous sources inside first collapsed halos at z > 10–15 if the DM is made of PBHs such as discovered by LIGO.

Zoom In Zoom Out Reset image size

If PBHs make up the DM, luminous sources within the much more abundant early collapsed halos would reproduce the observed Spitzer/Akari source-subtracted CIB fluctuations with modest formation efficiency requirements. This can be demonstrated by taking population models from Helgason et al. (2016) and rescaling them by the collapse-efficiency ratio from Figure 3. Thus, Figure 2 (upper right) from Kashlinsky et al. (2015a) would now reproduce the observed CIB signal with only f_* < 0.5% forming out to z ≳ 15 (instead of 4% with formation continuing to z ≃ 10), and the lines in Figure 5 of Helgason et al. (2016) need to be rescaled down by the corresponding factors. Additionally the measured CIB-CXB coherence (Cappelluti et al. 2013) would require that at least ≳(10%–15%) of the luminous CIB-producing sources are accreting BHs, broadly consistent with this discussion.

We now outline briefly the possible modifications in the early collapse and source formation that the PBHs may require. Two temperature regimes are relevant for the description of emitting sources: (1) minihalos where H₂ formation is efficient evolve at T ≲ 10³ K and the gas converges toward density of n_gas ∼ 10⁴ cm⁻³ (Bromm & Larson 2004 and references therein); (2) in the absence of H₂, the metal-free gas will be able to cool to 10⁴ K and collapse in halos with larger virial temperature will proceed isothermally. Feedback effects from first sources would affect H₂ formation via a resulting Lyman–Werner (LW) radiation at [11.2–13.6] eV (see the review by Bromm 2013). Gas collapse/evolution in the PBH minihalos may affect the subsequent emitting source formation inside them.

PBHs will accrete the minihalo gas resulting in both the additional source of emission from PBH accretion as well as the LW radiation feedback. The gas at sound speed c_s within the halo of velocity dispersion v_d will be accreted within the typical radius r_acc = GM_PBH/u² with ${u}^{2}={v}_{d}^{2}+{c}_{s}^{2}$. The total accretion mass will be ${M}_{{\rm{acc}}}=2{({n}_{{\rm{gas}}}/{10}^{4}{{\rm{cm}}}^{-3})({M}_{{\rm{PBH}}}/30{M}_{\odot })}^{3}$ ${(u/1\mathrm{km}{{\rm{s}}}^{-1})}^{-6}{M}_{\odot }$. For typical parameters, this may be a non-negligible fraction of the minihalo baryons at $\sim {M}_{{\rm{acc}}}/{M}_{{\rm{PBH}}}\times {{\rm{\Omega }}}_{{\rm{CDM}}}/{{\rm{\Omega }}}_{{\rm{bar}}}\propto {M}_{{\rm{PBH}}}^{2}{u}^{-6}$ up to a few percent, but will not increase the PBH mass dramatically. (Note the sensitive dependence on u, so M_acc is rapidly decreased when u ≫ 1 km s⁻¹.) The spectrum of the resultant emission may be modeled after Yue et al. (2013a): (1) the multicolor blackbody from different parts of the accretion disk with temperatures up to ${T}_{{\rm{max}}}\simeq 0.4{({M}_{{\rm{PBH}}}/30{M}_{\odot })}^{-1/4}\;\mathrm{keV}$ shifted mainly into the observer's near-IR, after reprocessing by the surrounding medium, and (2) hot corona and reflection emissions, which leave their mark in the observer soft X-rays. The coherence between the near-IR and X-ray emissions would be strong for the PBHs because of the larger value of T_max than for DCBHs. This mode of evolution, inevitable if PBHs make up the DM, may influence adjacent star formation and DCBH collapse and evolution as discussed, e.g., in Bromm & Loeb (2003), Agarwal et al. (2012), and Yue et al. (2014).

The PBHs in minihalos will evolve via stellar dynamical effects similar to that discussed in Kashlinsky & Rees (1983) and by loss of energy to GW emissions (Bird et al. 2016). Stellar evaporation will lead to a core–halo structure with the isothermal core of radius r_c and N_PBH PBHs evolving on gigayear-timescales ${t}_{{\rm{evap}}}\sim {N}_{{\rm{PBH}}}/\mathrm{ln}{N}_{{\rm{PBH}}}\times {r}_{{\rm{c}}}/{v}_{d}$, at constant binding energy, or ${v}_{d}\propto {N}_{{\rm{PBH}}}^{-1/2}$, because evaporating PBHs carry zero energy. At the same time, a fraction of PBHs will become binary when GW emission exceeds their kinetic energy ($\sim {v}_{d}^{2}$); the cross-section for this process being ${\sigma }_{{\rm{GW}}}\simeq {10}^{-8}{({M}_{{\rm{PBH}}}/30{M}_{\odot })}^{2}{({v}_{d}/1\mathrm{km}{{\rm{s}}}^{-1})}^{-18/7}\;{{\rm{pc}}}^{2}$ (Bird et al. 2016). The fraction of PBHs that will form binaries before evaporation is then

$$\begin{eqnarray}&&{f}_{{\rm{PBH,}}{\rm{binary}}}\sim \displaystyle \frac{{N}_{{\rm{PBH}}}^{2}}{\mathrm{ln}{N}_{{\rm{PBH}}}}\displaystyle \frac{{10}^{-8}\;{{\rm{pc}}}^{2}}{{r}_{c}^{2}}{\left(\displaystyle \frac{{M}_{{\rm{PBH}}}}{30{M}_{\odot }}\right)}^{-2}{\left(\displaystyle \frac{{v}_{d}}{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-18/7}.\end{eqnarray} \tag{ 5 }$$

Instead of evaporating, the resultant binaries will spiral in to the center due to dynamical friction possibly forming a central large BH contributing to the massive BH formation in early universe.

Finally, we note further constraints from reionization by both first stars and BHs as discussed in this context in Atrio-Barandela & Kashlinsky (2014) and Helgason et al. (2016). While reionization would be complicated in the presence of X-rays (Ricotti & Ostriker 2004; Ricotti et al. 2005), the Thomson optical depth of τ ≲ 0.1 may imply that the ionizing photons at rest <0.0912 μm are mostly absorbed in their paternal minihalos, although Atrio-Barandela & Kashlinsky (2014) recover τ ≲ 0.05 from hot gaseous bubbles reionized by first stars producing the observed CIB levels and forming to z = 9. The situation, while important, is clearly model-dependent: emissions probed at observer λ ≳ 2 μm, where current CIB data appear established, would translate directly into ionizing, Lyman-continuum photons only at $1+{z}_{{\rm{ion}}}\geqslant 22\left(\tfrac{\lambda }{2\;\mu {\rm{m}}}\right);$ at these epochs, the baryon density is high and harder to reionize and only a small fraction of the CIB is expected to be produced. A possibility, discussed in Yue et al. (2013a) for the DCBH model, whereby the gaseous collapsed halos are Compton thick so the ionizing photons are absorbed and reprocessed into a two-photon continuum, may also appear relevant here.

This work was supported by NASA/12-EUCLID11-0003 "LIBRAE: Looking at Infrared Background Radiation Anisotropies with Euclid" project (http://librae.ssaihq.com). I thank my LIBRAE colleagues for comments on the draft of this paper.