Stochastic background of gravitational waves generated by a cosmological population of young, rapidly rotating neutron stars

We estimate the spectral properties of the stochastic background of gravitational radiation emitted by a cosmological population of hot, young, rapidly rotating neutron stars. Their formation rate as a function of redshift is deduced from an observation-based determination of the star formation history in the Universe, and the gravitational energy is assumed to be radiated during the spin-down phase associated to the newly discovered r-mode instability. We calculate the overall signal produced by the ensemble of such neutron stars, assuming various cosmological backgrounds. We find that the spectral strain amplitude has a maximum $\approx (2-4)\times 10^{-26} {Hz}^{-1/2}$, at frequencies $\approx (30-60)$ Hz, while the corresponding closure density, $h^2 \Omega_{GW}$, has a maximum amplitude plateau of $\approx (2.2-3.3) \times 10^{-8}$ in the frequency range $(500-1700)$ Hz. We compare our results with a preliminary analysis done by Owen et al. (1998), and discuss the detectability of this background.


INTRODUCTION
In the last two years, a series of investigations on the perturbations of rotating relativistic stars have revealed the existence of a class of r-modes that are unstable due to the emission of gravitational radiation (Andersson 1998;Friedman & Morsink 1998;Lindblom, Owen & Morsink 1998). This gravitationally driven instability is of considerable importance because it has a number of interesting astrophysical implications. For instance, it determines a spin down of newly born neutron stars to periods closer to the initial periods inferred for ms pulsars. Moreover, as the star spins down an energy equivalent to ≈ 1% of a solar mass is emitted in gravitational waves making the process a very interesting source for gravitational-wave detection Andersson, Kokkotas & Schutz 1998).
In this paper we study the stochastic background of gravitational waves contributed by a cosmological population of hot, young, rapidly rotating neutron stars with a small initial excitation in the l = 2 r-mode. A rough estimate of the spectrum of this background radiation has recently been given by Owen et al. (1998), assuming a comoving number density of neutron star births which is constant in the range 0 < z < 4 and zero at earlier times.
Here we present a more complete analysis, which derives the neutron star birth rate as a function of redshift from an observation-based determination of the star formation density evolution (Madau, Pozzetti & Dickinson 1997), within three different cosmological background models. Thus, the present investigation relies on a source rate which correctly accounts for the observed strong luminosity evolution of star forming galaxies between z = 0 and z = 1 − 2 (Lilly et al. 1998;Connolly et al. 1997). With a similar approach, the gravitational background emitted by an ensemble of stellar mass black holes has been recently computed by Ferrari, Matarrese & Schneider (1998). The plan of the paper is as follows. In Section 2 we calculate the universal rate of neutron star formation; in Section 3 we briefly introduce the adopted model of gravitational-wave emission by an unstable rapidly rotating neutron star, closely following the recent analysis by Owen et al. (1998); in Section 4 we estimate the overall stochastic background produced by the ensemble of such sources throughout the Universe and calculate the total duty cycle of the process; Section 5 contains a brief analysis of the detectability of such a background, while conclusions are drawn in Section 6.

THE RATE OF NEUTRON STAR FORMATION IN THE UNIVERSE
Type II and Type Ibc supernovae are fuelled by the gravitational contraction of their core (e.g. Ruiz-Lapuente 1997).
Numerical studies have shown that single stars with masses between 8 and 20M⊙ form iron cores with masses near the Chandrasekhar limit (e.g. Timmes, Woosley & Weaver 1996). Thus, there is universal agreement that these stars leave neutron star remnants. The fate of stars with masses greater than 20M⊙ is less certain and depends on the iron core masses formed and on the effect of fallback during the explosion (e.g. Woosley & Weaver 1995;.
For the calculations to follow, we assume that stars with masses between 8 and 25M⊙ give origin to neutron stars. However, we will also investigate the effect of an upper mass cutoff of ∼ 20M⊙.
The number of neutron stars formed per unit time within the comoving volume out to redshift z is given by (Ferrari, Matarrese & Schneider 1998) Hereρ * (z) is the star formation rate density (mass of gas that goes into stars per unit time and unit comoving volume), Φ(M ) is the initial mass function (IMF), and dV /dz is the comoving volume element. Following Madau, Pozzetti & Dickinson (1997), we model the star formation evolution from the recent collection of UV-optical observations of star forming galaxies out to redshifts ∼ 4 − 5, by assuming a Salpeter IMF, normalized through the relation The star formation history is plotted in Figure 1 for a flat cosmology with vanishing cosmological constant and h = 0.5 (here h is the Hubble constant, H0, in units of 100 km s −1 Mpc −1 ). The star formation rate rises sharply from its local value to a peak at z ∼ 1.5 to fall again out to z ∼ 4. This behaviour indicates that the bulk of the stellar population was assembled in the redshift range 1 < ∼ z < ∼ 2. This is consistent with the indications of QSO absorption lines and metallic clouds observations (Lanzetta, Yahil & Fernandez-Soto 1996;Pei & Fall 1995) and with the predictions of semianalytic models for a broad class of hierarchical clustering cosmologies (Baugh et al. 1998). A thorough review of the current state of observations of the star formation history of field galaxies is given in Ellis (1997). A summary of the relevant aspects for the present analysis is given in Ferrari, Matarrese & Schneider (1998).
The predictions for the rate of core-collapse supernovae are in good agreement with the observed local values (Cappellaro et al. 1997). At z > 1, the detection of Type II events must await future experiments, such as the Next Generation Space Telescope (Madau, Della Valle & Panagia 1998).
Here we shall only stress that one major source of uncertainty in the interpretation of the data points in Figure  1 is the effect of dust extinction. Even a small amount of  (Madau et al. 1996) (filled squares). The solid line represents the SFR density evolution assumed in the present analysis. The fit has been kindly provided by P. Madau. dust in primeval galaxies can significantly attenuate the UV luminosity and re-radiate it in the far IR thus leading to a potentially serious underestimate of the star formation activity. Dust correction factors are still very uncertain mainly because of the unknown shape of the dust UV-extinction law (Pettini et al. 1997). The model for the global SFR evolution presented in Figure 1 assumes an extinction law similar to that which applies to stars in the Small Magellanic Cloud (SMC) and an amount of dust which results in an upwards correction of the comoving UV luminosity by a factor 1.4 at 2800Å and 2.1 at 1500Å (Madau, Pozzetti & Dickinson 1997). Recent near infrared observations of a galaxy at z = 4.92 has given evidence that absorption by dust is important up to redshifts of ∼ 5 (Soifer et al. 1998). Similarly, analysis of the Infrared Space Observatory (ISO) data have shown significant departure from the UV-optically derived SFR evolution at z > 2 ( Rowan-Robinson et al. 1997). This may seem a serious caveat for the present investigation. However, this is not necessarily the case because the extra-galactic astrophysical gravitational-wave backgrounds prove to be mainly contributed by low-to-intermediate redshift sources. High redshift sources play a negligible role because the energy flux emitted by each source decreases as the inverse of the squared luminosity distance. On the other hand, if an upwards correction of the UV-optically derived SFR density is required at z ≈ 1.5, as suggested by a recent detection of the infrared background by the Diffuse Infrared Background Experiment (DIRBE) and the Far Infrared Absolute Spectrometer (FIRAS) on board of the Cosmic Background Explorer (COBE) satellite (Dwek et al. 1998), this would result in an upwards correction of roughly the same amplitude on our gravitational-wave background predictions.
We have evaluated equation (1) for three sensible cosmological models (e.g. Baugh et al. 1998), namely a flat geometry with vanishing cosmological constant, an open background and a flat low-density model (see Table 1). For a gen- Table 1. Cosmological parameters adopted in the three backgrounds considered.
Model eral background, the comoving volume element is related to z through, where H0 is the Hubble constant and In eq. (4) r is the comoving distance, and the function S has the following expression (Kim et al. 1997) ΩM + ΩΛ < 1 S(x) = sinh(x) Ωκ = ΩM + ΩΛ − 1, Consistently, the SFR history shown in Figure 1 needs to be properly rescaled in order to account for the cosmological background dependence of both the luminosity distance and the comoving volume within each redshift bin. By using equation (4) we have evaluated the function RNS(z) defined in eq. (1) for the three cosmological models, and the results are shown in Figure 2.
It should be mentioned that the evaluation of the rate of neutron star formation, based on the UV-optically estimated star formation rate evolutionρ * (z), and on an assumed universal IMF, is largely independent of the specific IMF adopted (Madau 1998). In fact, the stars which mainly contribute to the UV-optical luminosity are the most mas- sive ones which, in turn, are those that at the end of their evolution give origin to core-collapse supernovae. The total number of supernova explosions per unit time leaving behind a neutron star is given in Table 1 for the considered cosmological models, and for two values of the upper mass cutoff, M = 25M⊙ and M = 20M⊙. As a consequece of the largest number of stars having masses < 25M⊙, the total rate of core-collapses to neutron stars is a factor ∼ 4 larger than that of core-collapses to black holes considered in Ferrari, Matarrese & Schneider (1998).

GRAVITATIONAL EMISSION FROM RAPIDLY ROTATING NEUTRON STARS
In this section we shall briefly describe how the timescale of the spindown phase and the energy spectrum of the gravitational radiation emitted during the r-mode instability have been estimated (Lindblom, Owen & Morsink 1998;Andersson, Kokkotas & Schutz 1998;Andersson, Kokkotas, Stergioulas 1998). These quantities will be used to compute, respectively, the duty cycle and the spectral energy density of the stochastic background. The excitation of an r-mode starts as a small perturbation of the velocity field. As the amplitude of the mode grows, non-linear hydrodynamic effects become important and eventually dominate the dynamics. Under the assumption that the star is uniformly rotating with angular velocity Ω, and by parametrizing the r-mode by its amplitude α, the star can be considered as a system having only two degrees of freedom, Ω and α. In this simple model, the total angular momentum for the ℓ = 2 r-mode can be shown to be wherẽ (9) and the rotational energy is Following Owen et al. (1998) we shall consider a Newtonian model of neutron star composed of a fluid with a polytropic equation of state, p = k ρ 2 , with k chosen so that M = 1.4M⊙ and R = 12.53 km. For this model,Ĩ = 0.261 andJ = 1.635 × 10 −2 . It should be noted thatĨ andJ depend only on the chosen polytropic index.
The evolution of the angular momentum of the star is determined by the emission of gravitational waves, which couple to the r-modes through the current multipoles, primarily that with l = m = 2. For this mode, the frequency of the emitted gravitational radiation is ν = (2/3π) Ω. The evolution of Ω and α during the phase in which α is small, can be determined from the standard multipole expression for angular momentum loss, and from the energy loss due to the gravitational emission and to the dissipative effects induced by the bulk and shear viscosity. In this phase Ω is nearly constant and α grows exponentially. After a short time, that Owen et al. estimate to be of order 500 s, the amplitude of the mode becomes close to unity and non-linear effects saturate and halt a further growth of the mode. This phase lasts for approximately 1 yr, during which the star loses angular momentum radiating approximately 2/3 of its initial rotational energy in gravitational waves, up to a point where the angular velocity and the temperature are so low that the r-mode ceases to be unstable. Viscous forces and gravitational radiation damp out the energy left in the mode, and the star slowly reaches its final equilibrium configuration.
Since gravitational radiation is emitted over such a long time interval, the background contributed by neutron stars is continuous. In fact, the total duty cycle for this process is It has been suggested that if the proton fraction in the star is sufficiently large to make the direct URCA process possible, the star would cool in a much shorter time, τ NS ≈ 20 sec. However, if this were the case, the duty cycle would still be ≫ 1 (D ≈ 10 2 ) and the resulting background would still be continuous. It is known that the angular velocity of a rotating star cannot exceed the Kepler frequency ΩK at which mass shedding at the stellar equator makes the star unstable. A reasonable approximation for this limiting frequency is where C ≈ 7.8×10 3 sec −1 (Friedmann, Ipser & Parker 1989), and M and R refer to the mass and radius of the corresponding non-rotating star. Following Owen et al., we shall assume that neutron stars are born with angular velocity close to the maximum value ΩK , and have a rotational kinetic energy In the non-linear saturation phase, where most of the gravitational radiation is emitted, α ≈ 1. Since in this simple model approximately 2/3 of the stars rotational energy is radiated in gravitational waves, the expression of the energy spectrum can be approximated as follows where νmax = (2/3π)ΩK .
The r-mode instability ceases to be effective for some value of the rotational frequency Ωc, which corresponds to the final spin period, and which is associated to a gravitational emission frequency νmin = (2/3π)Ωc.
Ωc can be determined by solving the equation 1/τ (Ωc) = 0, where τ is the total dissipation timescale which can be decomposed as a sum of the damping times associated to the gravitational emission, to the shear and to the bulk viscosity. τ (Ωc) is clearly a function of the temperature of the star, and it has been shown that the r-mode instability operates only in hot neutron stars (T > 10 9 K). Below 10 9 K, superfluidity and other non-perfect fluid effects become important and the damping due to viscosity dominates with respect to the destabilizing effect of the gravitational radiation. To be consistent with the analysis by Owen et al. (1998), we shall adopt a value of Ωc ≈ 566 Hz for the final angular velocity which corresponds to a final spin period of ≈ 11 ms and to νmin ≈ 120 Hz. The qualitative picture that arises from this simple model is believed to be sufficiently reliable, even though various uncertainties and approximations might affect the quantitative results for the initial rotation of the star after collapse, for the spindown timescales as well as for the final rotation period (Andersson, Kokkotas & Schutz 1998). Thus, it will be interesting to evaluate the effects that different values of νmax, EK and νmin have on our main conclusions.

SPECTRAL PROPERTIES OF THE STOCHASTIC BACKGROUND
In order to evaluate the spectral energy density of the stochastic background produced by the spindown radiation from newly born neutron stars, we need to know how the flux emitted by a single source at redshift z would be observed today, and convolve it with the differential rate of the sources dRNS(z) as in Ferrari, Matarrese & Schneider (1998). The average energy flux per unit frequency emitted by a single neutron star at cosmological redshift z is where ν obs is the frequency at which the radiation would presently be observed. The energy flux f (ν obs ) can be evaluated by making use of eq. (14) suitably redshifted as follows where ν em min 1 + z ν obs ν em max 1 + z , and dL(z) = (1 + z) r is the luminosity distance, which depends on the cosmological model as indicated in eq. (7). Thus, the spectral energy density of the produced background is where the differential rate of sources is given by  Table 1).
[see eq. (1)]. From the spectral energy density the corresponding values of the closure energy density of gravitational waves where ρcr = 3H 2 0 / 8π G, and of the spectral strain amplitude can be derived.
In Figures 3 and 4 we plot the spectral energy density and the strain amplitude for the flat background with zero cosmological constant corresponding to model A (see Table  1). The effect of a varying cosmological background is negligible on both the spectrum and the strain amplitude. In fact, the amplification of the rate at high redshifts shown in Figure 2 for the models B and C, is mostly suppressed by the inverse squared luminosity distance dependence of the single source spectrum [see eq. (18)] for the same models. Conversely, the closure density exhibits a more evident dependence on the Hubble constant, as shown in Figure 5, where we have considered only model A and B because the dependence on the other cosmological parameters (namely ΩM and ΩΛ) is negligible.
The estimate of the intensity of the stochastic background also depends on the chosen values of νmin and νmax. We have repeated our calculations for two values of νmin, chosen so as to cover the range of pulsar rotational periods, Pin, inferred from observations with Pin between 9 and 22.8 ms, where the upper limit is the maximum rotational period at which the r-modes are unstable (Andersson, Kokkotas & Schutz 1998). Since the strain amplitude is sensitive mainly to the low frequency region of the spectrum, due to the factor ν −2 obs in eq. (23), a decrease of νmin determines an increase of the peak amplitude and shifts its location towards lower frequencies, as shown in Figure 6.    Table 1).
The value of νmax, depends on ΩK which, in turn, depends on the mass and radius of the star, as shown in eqs. (15) and (12). The plots shown in Figures 3, 4 and 5 were obtained by assuming M = 1.4M⊙ and R = 12.53 km. In reality, the possible radii for rotating neutron stars range from 10 to 15 km, therefore, since the dependence of ΩK on the mass is less significant than that on the radius, we have computed the strain amplitude and the closure density by considering M = 1.4M⊙, and changing the radius to R = 10 km and to R = 15 km. The effect of reducing νmax, i.e. of increasing the radius, is that of narrowing the interval of frequencies over which each source contributes. Since the energy flux of a single event does not depend on ΩK , [see eq. (18)], this will result in an increase of the amplitude of the peak of the spectral energy density. This effect is amplified in the strain amplitude, that is more sensitive to low frequencies. However, it should be stressed that the increase is smaller than a factor 2. Conversely, the closure density is more sensitive to the high frequency contributions, and

DETECTABILITY
The value of the optimal signal-to-noise ratio (SNR) for Lshaped interferometers can be computed according to the formula (Allen 1996) where S (i) h are the spectral noise densities of two detectors operating in coincidence, and γ(ν) is the overlap function which takes into account the difference in their location and orientation. We shall calculate SNR for several pairs of detectors, and for ΩGW computed for a flat cosmology, νmin = 120 Hz, and νmax = 1396 Hz. We shall assume an integration time of T = 1 yr. For interferometric antennas, we have computed S (i) h from the expected sensitivity curves given for VIRGO (Beccaria et al. 1996) and for ADVANCED LIGO (Flanagan & Hughes 1997). The sensitivity curve for GEO600 has been extracted from the Web site http://www.geo600.unihannover.de. For VIRGO and GEO600 we find SN R = 3.7 × 10 −3 , assuming optimal orientation, and SN R = 5.0 × 10 −3 if GEO600 is narrowbanded around 100 Hz. For two AD-VANCED LIGO detectors SN R = 0.84. Since the correlation of two interferometers at a distance d starts to be very poor around a frequency ≈ (70 Hz)(3000 km/d) (Allen & Romano 1998), two nearby detectors respond better to higher frequencies, where our predicted signal has a larger amplitude (see Figure 8). Thus, the signal to noise ratio rises to 8.5 if we consider two nearby ADVANCED LIGO (or VIRGO) interferometers at a distance of ≈ 300 km.
We have similarly investigated how interferometric and resonant detectors operating in coincidence might respond to our predicted background. The correlation of VIRGO and NAUTILUS, situated 265 km apart, gives SN R = 3.0×10 −4 . Assuming that a similar correlation is done with a truncated icosahedral antenna (TIGA) (Merkowitz & Johnson 1995) and an ADVANCED LIGO interferometer located at the same site, the resulting SN R is ≈ 1.3×10 −3 , whereas it is SN R ≈ 1.7×10 −3 if the interferometer is ADVANCED VIRGO.
Thus, the detection of this background requires two nearby interferometers with 'advanced' sensitivities, as similarly noted by Owen et al. (1998).

CONCLUSIONS
In this paper we have evaluated the gravitational background produced by a cosmological population of hot, young and rapidly rotating neutron stars, emitting gravitational radiation during the spin-down phase associated to the rmode instability. The dependence of this background on the cosmological model and on the parameters that enter into the determination of the spectrum of each single event, have been discussed. Our predictions for the spectrum, spectral strain amplitude and closure density, appear to be rather insensitive to the assumed cosmological parameters. This is because the increase of the neutron star birth rate we find in an open or a flat low-density cosmology is compensated by the increase of the luminosity distance as a function of redshift, which suppresses the contribution of the farthest sources.
For the same reason, the uncertainties that still trouble the high redshift tail of the SFR density evolution due to dust obscuration, are not severe for our signal, as this is mainly contributed by low-to-intermediate redshift sources. Moreover, if dust obscuration produces an upward correction of the SFR density at intermediate redshifts (i.e. around the location of the maximum) then the amplitude of (S h ) 1/2 and of ΩGW will be similarly amplified at ≈ 30 Hz and ≈ 600 Hz respectively.
We find that, within a reasonable range of values of the main parameters which characterize the energy spectrum of a single source (namely νmin, νmax and E k ), the spectral strain amplitude has a maximum ranging within ≈ (2 − 4) × 10 −26 Hz −1/2 at frequencies ≈ (30 − 60) Hz, whereas the closure density, h 2 ΩGW , is shown to have a long plateau extending from ≈ 300 Hz up to ≈ 1700 Hz with an amplitude of ≈ (2.2 − 3.3) × 10 −8 .
It is interesting to compare our results with those of a similar study done by Owen et al. (1998), who assume a constant comoving number density of neutron star births in the range 0 < z < 4, and zero elsewhere.
In Figure 9 we show the closure density predicted by the present analysis and that obtained by Owen et al. The two curves exhibit quite a different high frequency behaviour. Since the cosmological expansion shifts the frequencies towards smaller values, the high frequency tail corresponds to closer sources, thus the difference between the two curves can be traced back to the SFR density evolution between z = 0 and z ≈ 1 − 2. The peak at ≈ 600 Hz exhibited by our curve, corresponds to z ≈ 1.3, where we expect to have the maximum number of sources. Conversely, the maximum amplitude predicted by Owen et al. is located around 280 Hz, and corresponds to the energy emitted at maximum spin rate by a source at z ∼ 4. Thus, their signal appears to be mostly contributed by distant sources near their maxi-mum rotational velocity, and this is due to the fact that they do not consider the luminosity distance damping, which is particularly effective for high redshift sources. It is important to stress that in the region where cross-correlation between two LIGO interferometers can be accomplished, i.e. for ν < 50 − 60 Hz, the signal is entirely emitted at z > 1, where it is important to embed the gravitational sources in the correct cosmological environment.
Both in our analysis and in the study of Owen et al., all neutron stars have been assumed to be initially rotating with an angular velocity close to the Keplerian value, while there might be a fraction of stars born with lower initial spin rate. In particular, no spindown gravitational radiation would be emitted if the initial angular velocity of the star is too small for the instability to set in Spruit & Phinney (1998). However, according to Andersson, Kokkotas & Schutz (1998) this is unlikely to be the generic situation, because the angular momentum of the degenerate core required to make the star rotating with velocity ΩK is rather small. Moreover, there is clear evidence that some pulsars were born spinning rapidly: as an example, the observed spin period for the young Xray pulsar in the supernova remnant N157B is 16 ms, and its inferred initial spin period is estimated to be ≈ 6 − 9 ms (Marshall et al. 1998).

ACKNOWLEDGMENTS
We greatly acknowledge Pia Astone, Piero Madau, Cedric Lacey and Silvia Zane for useful discussions.