The gravitational-wave signal generated by a galactic population of double neutron-star binaries

1.2 MB·English

by  Shenghua Yu

Mon.Not.R.Astron.Soc.000,1–24(2011) Printed13January2015 (MNLATEXstylefilev2.2)
The gravitational-wave signal generated by a galactic
population of double neutron-star binaries
5
1 Shenghua Yu1,2⋆, C. Simon Jeffery3,4†
0
1National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China
2
2The KeyLaboratory of Radio Astronomy, Chinese Academy of Sciences

n 3Armagh Observatory, College Hill, Armagh BT61 9DG, N. Ireland
a 4School of Physics, TrinityCollege Dublin, Dublin 2, Ireland
J

0
1 Accepted .Received;inoriginalform

]
R
ABSTRACT
S
Weinvestigatethegravitationalwave(GW)signalgeneratedbyapopulationofdouble
.
h neutron-star binaries (DNS) with eccentric orbits caused by kicks during supernova
p collapse and binary evolution. The DNS population of a standard Milky-Way type
- galaxy has been studied as a function of star formation history, initial mass function
o
(IMF) and metallicity and of the binary-star common-envelope ejection process. The
r
t model provides birth rates, merger rates and total numbers of DNS as a function of
s
time. The GW signal produced by this population has been computed and expressed
a
[ in terms of a hypothetical space GW detector (eLISA) by calculating the number of
discrete GW signals at different confidence levels, where ‘signal’ refers to detectable
1
GW strain in a given frequency-resolution element. In terms of the parameter space
v
explored, the number of DNS-originating GW signals is greatest in regions of recent
4
star formation, and is significantly increased if metallicity is reduced from 0.02 to
1
0.001, consistent with Belczynski et al. (2010). Increasing the IMF power-law index
3
2 (from–2.5to–1.5)increasesthenumberofGWsignalsbyalargefactor.Thisnumber
0 is also much higher for models where the common-envelope ejection is treated using
. theα−mechanism(energyconservation)thanwhenusingtheγ−mechanism(angular-
1
momentumconservation).WehaveestimatedthetotalnumberofdetectableDNSGW
0
signals from the Galaxy by combining contributions from thin disc, thick disc, bulge
5
1 and halo. The most probable numbers for an eLISA-type experiment are 0−1600
: signals per year at S/N>1, 0−900 signals per year at S/N>3, and 0−570 at S/N>5,
v
coming from about 0−65,0−60 and 0−50 resolved DNS respectively.
i
X
Key words: Gravitational waves - neutron stars - Stars: binaries: close - Galaxy:
r structure - Galaxy: stellar content
a
1 INTRODUCTION binary. Evidence that such systems do form is provided by
thedoublepulsarPSRJ0737-3039A/B (Burgay et al.2003;
Most stars are members of binary or multiple star systems.
Lyneet al. 2004; Kramer & Stairs 2008). Such compact bi-
Over 70% of massive stars (O-type stars) have a nearby
nariesareexpectedtobeasignificantsourceofgravitational
companion which will affect their evolution, with over one
wave(GW)radiation (Barish & Weiss1999;Ricci & Brillet
half doing so before they leave the main sequence (MS)
1997; Postnov & Yungelson 2006). They are amongst the
(Sana et al.2012;Langer2012).Ifbothmembersofabinary
sources most likely to be detected by a gravitational wave
are massive enough to end their evolution as core-collapse
detector in the frequency range 10−5 100 Hz. Aasi et
supernovae1, after one or two explosions, the final prod- −
al. (2014) measured upper limits of the GW strain ampli-
uct could become a double neutron-star binary (DNS), a
tudesfrom hundredsof pulsarsusing datafrom recent runs
neutron-star plus black-hole binary, or a double black-hole
of the ground-based GW observatories - LIGO, Virgo and
GEO600, and showed that there are good prospects for de-
tectionsinthe10 -1000 Hzrangewith theadvancedLIGO
⋆ [email protected] and Virgo detectors.
† [email protected]
1 NeutronstarscanalsobeformedbyaccretionofOxygen-Neon Double neutron-star and black-hole binaries play an
whitedwarfs important role in testing the theory of General Relativ-
2
double compact objects using the binary-star population-
Table1.Orbitalperiods(Porb)andfrequencies(forb),eccentric-
ities (e), GW timescales (τGW) and whether masses have been synthesis (BSPS) method. They concluded that only a few
establishedforknownandsuspectedDNSbinaries;afterLorimer (2-4) NS-NS binaries in the Galaxy would have been de-
(2008). tectablebythecancelled spaceobservatoryLISA.However,
approximationsfor(i)thecalculationoftheGWsignalfrom
individualDNSbinaries,and(ii)theemploymentofcrucial
DNS Porb e M log10 forb initial conditions and the treatment of important physical
PSR (d) τGWyr−1 Hz processes in the BSPS method, may result in quite large
J0737–3039 0.102 0.09 Yes 7.9 1.1×10−4 uncertainties.
J1906+0746 0.17 0.09 Yes 8.5 6.8×10−5 In order to understand the signal detected by suffi-
B1913+16 0.3 0.62 Yes 8.5 3.9×10−5 cientlysensitiveGWdetectors,itisnecessarytocharacterize
B2127+11C 0.3 0.68 Yes 8.3 3.9×10−5 the radiation from all GW-emitting populations, including
J1756–2251 0.32 0.18 Yes 10.2 3.6×10−5
DNS. This paper presents a study of the GW signal from
B1534+12 0.4 0.27 Yes 9.4 2.9×10−5 theGalacticpopulationofsteady-state DNSs2 includingthe
J1829+2456 1.18 0.14 No 10.8 9.8×10−6
Galactic star-formation history, the initial-mass function,
J1518+4904 8.6 0.25 No 12.4 1.3×10−6
J1811–1736 18.8 0.83 Yes 13.0 6.2×10−7 metallicity, and the physics of common-envelope evolution.
B1820–11 357.8 0.79 No 15.8 3.2×10−8 ItexamineshowthisDNSpopulationdiffersfromtheGalac-
tic double-white-dwarf (DWD) population, and establishes
Inorderofdiscovery:PSRB1913+16: Hulse&Taylor(1975), a basis for computing the GW signal due to extragalactic
PSRB1820–11: Lyne&McKenna(1989),PSRB1534+12:
DNS populations. We describe the methods used to model
Wolszczan(1991),PSRB2127+11C: Princeetal.(1991),PSR
the DNS population in the Galactic disc, the emission and
J1518+4904: Niceetal.(1996),PSRJ1811–1736: Lyneetal.
superposition of GW signals, and thereduction of the data
(2000),PSRJ0737–3039A/B: Burgayetal.(2003);Lyneetal.
in terms of a conceptual GW experiment (eLISA) in 2.
(2004),PSRJ1829+2456: Championetal.(2004),PSR
§
J1756–2251: Faulkneretal.(2005),andPSRJ1906+0746: Majorresultsarepresentedin 3.Implications,observations
§
Lorimeretal.(2006). andpreviousworkarediscussedin 4.Themainconclusions
§
are reviewed in 5.
§
ity,whilst double pulsars provideprobes of magnetospheric
physics.Theaccuratetimingofpulsarsorbitingablackhole 2 METHODS
can be used to constrain the strain amplitude of gravita-
2.1 Binary-star population synthesis
tional waves and the physical properties of the black hole.
(Kramer et al. 2004). In this section, we review the important physical processes
At present (2014) 7 DNSs, have been confirmed and 3 inbinarystarevolutionandtheinitialandboundarycondi-
more are suspected (Table 1). Of these, half have merger tions for population synthesis to obtain a sample of double
timescales shorter than a Hubble time. Half also show ec- compact objects. Population synthesis, including individ-
centricorbits,evenatrelativelyshortperiods.Bayesiansta- ual stellar-evolution tracks and initial conditions, was car-
tisticalanalysesbasedontheseobservationsindicatethatan riedoutusingthemethoddescribedbyYu& Jeffery(2010,
optimisticGalacticDNSmergerratemaybeupto1.810−4 2011)andHurley et al.(2000,2002)inaninitialstudyofthe
yr−1, implying that their number should be approximately present Galactic double degenerates population. Note that
a few million (Kalogera et al. 2001, 2004) if we assume the wetake some standard parameters for theGalactic disc (or
age of the Galaxy is 14 Gyr. This number is roughly 1-2 theGalaxy) to represent a Milky-way typegalaxy.
∼
ordersofmagnitudehigherthanestimatedfrombinarystar
population synthesis (Nelemans et al. 2001; Os lowski et al.
2011;Dominik et al. 2012),although theuncertainty in the 2.1.1 Common-envelope evolution
Bayesian analysis can exceed one order of magnitude.
When one star in a binary fills its Roche lobe either by
TheoreticalinvestigationsoftheGWsignalfromGalac-
evolutionary expansion or by orbital shrinkage, Roche lobe
ticDNSshavebeencarriedoutby(e.g.)Allen et al.(1999);
overflow (RLOF) occurs. The Roche radius of the primary
Belczynski et al. (2010a); Rosado (2011); Allen et al.
is given by
(2012). The stochastic GW background produced by DNS
mergers in low-redshift galaxies (up to z 5) has been in- R 0.49q2/3
∼ L1 = 1 (1)
vestigated by Regimbau & de Freitas Pacheco (2006) (also a 0.6q2/3+ln(1+q1/3)
1 1
see Rosado (2011)), who found that the signal should be
whereaismoregenerallythesemimajoraxisoftheorbitand
detectable by the new generation of ground-based inter-
q (= m /m ) is the mass ratio of primary and secondary
ferometers. Zhu et al. (2013) studied the GW background 1 1 2
(Eggleton1983).Hurley et al.(2000,2002)haveshownthat
from compact binary mergers, and showed that, below 100
thereisacriticalmassratioq ofabinarystarwhichcanbe
Hz, the background depends primarily on the local merger c
used todistinguish thestable RLOFandcommon envelope
rate and the average chirp mass and is independent of the
chirp mass distribution. In addition, the effects of cosmic
star formation rates and delay times between the forma- 2 A DNS merger would be clearly indicated by a strong GW
tion and merger of binaries are linear below 100 Hz in signalatrapidlyrisingfrequencies>0.1Hz.Themostoptimistic
theirmodel.Belczynski et al.(2010a)studiedtheGWback- GalacticDNSmergerratesare<10−3yr−1–see§3.asafunction
∼
groundandforegroundsignalfromaGalacticpopulationof ofseveralkeyparameters,
Gravitational waves from double neutron stars 3
(CE)phases,whereq isafunctionofprimarymassm ,its spiral-instageinordertoejecttheenvelope(ifweonlycon-
c 1
coremassm ,andthemass-transferefficiencyofthedonor. sider the orbital energy as the main engine) whilst the γ
1c
−
Weadopt mechanism does not, implying that the orbital separation
can be larger after CE ejection in thelatter case. In partic-
q = 1.67 x+2 m1c 5 /2.13, (2) ular,undertheassumption of noexternalmoment imposed
c − (cid:18)m1 (cid:19) ! onaconservativemass-transferbinary,theangularmomen-
tumofthebinarymustbeconservative,andthefinalorbital
where x = 0.3 is the exponent of the mass-radius relation
separation can be written as
at constant luminosity for giant stars (Hurley et al. 2000,
2002). a (m ∆m)(m +∆m) 2
f = 1− 2 , (7)
The calculation of the orbital parameters (e.g. orbital a m m
i (cid:20) 1 2 (cid:21)
separation) of a binary after CE ejection in our model is
where ∆m is the fraction of mass transferred from the pri-
basedononeoftwoassumptions:either1)angularmomen-
marytothesecondary.Thismeansthatinconservativeevo-
tumconservation(γ mechanism)or2)energyconservation
(α mechanism). − lution, if m1 >m2 prior to mass transfer, the orbital sepa-
− ration incfreases after mass transfer.
Inthefirstcase,weconsidertheangularmomentumlost
AlthoughbothCEejection formulationscanreproduce
byabinarysystemundergoingnon-conservativemasstrans-
observations (Nelemans & Tout 2005; Webbink 2008) via a
fer to be described by the decrease of primary mass times
variation of free parameters, considering both conservation
a factor γ (Paczyn´ski & Zi´o lkowski 1967; Nelemans et al.
lawsmaybeabetterapproachtoconstraintheCEevolution
2000):
andfinalorbitofabinary.Notethatinthispaperweneglect
Ji−Jf =γm1−m1c, (3) viscosity,friction betweentheCEandthestellar cores, and
J m +m thepotential nuclear and chemical energy in thesystem.
i 1 2
where J is the orbital angular momentum of the pre-mass
i
transfer binary; J is the final orbital angular momentum
f 2.1.2 Gravitational radiation, magnetic braking, and tidal
after CE ejection; m and m is the primary mass and its
1 1c interaction
core mass respectively; m is thesecondary mass.
2
Combining the above equation with thefraction of an- Othermechanisms which reducethe orbital separation of a
gularmomentumlost duringthemasstransfer,J J ,and binary system include gravitational radiation and magnetic
i f
−
Kepler’slaw,wehavetheratiooffinaltoinitialorbitalsep- braking. A close compact binary system driven by gravi-
aration tational radiation may eventually undergo a mass transfer
phase, ultimately leading to coalescence. Using the average
a m 2 m +m m m 2
f = 1 1c 2 1 γ 1− 1c , (4) energy (E) and angular momentum (Jorb) loss during one
ai (cid:18)m1c(cid:19) (cid:16) M (cid:17)(cid:16) − M (cid:17) orbitalperiod,wededucethedecayoforbitalseparationand
wherea anda are theorbital separations before and after eccentricity with respect to time to be
i f
the CE phase; M = m1 +m2 is the sum of the primary da 64G3m m (m +m )1+ 73e2+ 37e4
andsecondarymassbeforetheCEphase.Nelemans & Tout = 1 2 1 2 24 96 , (8)
dt − 5 c5a3 (1 e2)7/2
(2005) investigated the mass-transfer phase of the progeni- −
torsofwhitedwarfsinbinariesemployingtheγ-mechanism
based on 10 observed systems and deduced a value of γ in de = G3m1m2(m1+m2) 31054e+ 11251e3. (9)
the range of 1.4 – 1.7. In order to investigate the influence dt − c5a4 (1 e2)5/2
−
of angular momentum loss on the rates of DCOs, we here
ThiscalculationisconsistentwithourcalculationoftheGW
adopt γ =1.3 and 1.5.
signal from DNS described in 2.2.
Inthesecondcase,CEejectionofabinarystarrequires §
Gravitational radiation could explain the formation of
that the envelope binding energy, including gravitational-
cataclysmic variables (CVs) with orbital periods less than
binding and recombination energies, must represent a sig-
3h, while magnetic braking of the tidally coupled primary
nificant fraction of the orbital energy (Webbink1984).
byitsownmagneticwindwouldaccountfororbitalangular-
momentumlossfromCVswithperiodsupto10h(Faulkner
G(m m )m Gm m Gm m
1− 1c 1 =α 1c 2 1 2 , (5) 1971; Zangrilli et al. 1997).Weusetheformula for therate
λrL1 (cid:18) 2af − 2ai (cid:19) ofangular-momentumloss duetomagneticbrakingderived
where λ is a structure parameter depending on the evolu- byRappaport et al. (1983) and Skumanich (1972):
tionarystateofthedonor,α istheCEejection efficiency
representinghow muchorbitCalEenergy was usedto eject the J˙ = 5.83 10−16menv rωspin 3 M R2yr−2, (10)
CE,rL1 istheRocheloberadiusofthedonorattheonsetof mb − × m (cid:18)R⊙yr−1(cid:19) ⊙ ⊙
masstransfer,andGisthegravitationalconstant(Webbink
wherer,m andmaretheradius,envelopemassandmass
env
1984).Rearranging in theform of Eq.4, we have
of a star with a convective envelope, and ω is the spin
spin
a m 2(m m )a −1 angular velocity of the star.
f = 1c 1+ 1− 1c i . (6) Tidalinteractioncausedbythegravitydifferentialplays
a m αλm r
i 1 (cid:18) 2 L1 (cid:19) an important role in the synchronization of stellar rota-
In this paper, we adopt λ = 1.0, and α = 0.5 and 1.0. tion and orbital motion, and the circularization of the or-
A major difference between the two CE ejection for- bit. Relatively complete descriptions of the tidal evolution
mulations is that energy conservation implies a significant havebeengivenbyHurley et al.(2002)andBelczynski et al.
4
(2008). In this paper, we adopt the same formulae and whileitisintheHertzsprunggaporontheredgiantbranch,
procedures to deal with the tidal evolution as by Hut then RLOFwill occur.
(1981); Zahn (1977); Campbell (1984); Rasio et al. (1996) If the adiabatic response of the radius of the mass
and Hurley et al. (2002). donor is less than the change of its Roche lobe radius
with respect to a change of mass, i.e. ∂lnRdonor <
∂lnMdonor ad
2.1.3 Formation of double neutron stars ∂∂llnnRMRdLonOoFr RLOF,masstransferwillbeun(cid:16)stableand(cid:17)acom-
(cid:16)monenvelop(cid:17)e(CE)willform.Interaction(friction)between
Our simulations assume three routes for the formation of
the compact cores and the CE will convert orbital energy
neutronstars: i) if astar has coremass ofm .2.25M at
c ⊙ into kinetic energy, heating and expanding the CE. If the
shellheliumignition,itevolvethroughdouble-shellthermal-
energyconversionmechanismissufficientlyefficient,theCE
pulsesuptheasymptoticgiantbranch.Thestarmaybecome
will be expelled and a compact binary with a short orbital
aneutronstarifitscoremassgrowsandeventuallyexceeds
period will result (see 2.1.1).
2.25M⊙; ii) if the core mass of the star on the thermal- The above format§ion channels for DNSs represent the
pulsing asymptotic giant branch does not exceed 2.25M
⊙ combinationofRLOFandCEprocesses. Notethatchannel
butitisheavyenough(m &1.6M )tobecomeanelectron-
c ⊙ Voccurswhentheenvelopeofamassiveprimaryisremoved
degenerate oxygen-neon white dwarf which may become a
byastellarwindratherthanafirstCEejection.CEejection
neutronstarviaaccretion-inducedcollapse; iii) ifastarhas
following evolution of the secondary may also give rise to a
acoremassofm &2.25M atthestartoftheearlyasymp-
c ⊙ compact binary.
totic giant (or red giant) branch, it will become a neutron
starwithoutascendingthethermal-pulsingasymptoticgiant
branch. If the core mass of a star at the time of supernova
explosion is sufficiently high (& 7M ), it will most likely 2.1.4 Neutron star kicks
⊙
becomeablackholeunlesssignificantmasslosstakesplace.
Observationsofpropermotionsindicatethatpulsarshavean
These criteria are consistent with Hurley et al. (2000).
extraordinary natal velocity higher than their nominal pro-
Thegravitationalmassofneutronstarsiscalculatedby
genitors(Minkowski1970;Lyneet al.1982;Lyne& Lorimer
m =1.17+0.09m , (11) 1994; Hansen & Phinney 1997; Fryer& Kalogera 1997).
ns α
This may result from binary evolution (Iben& Tutukov
where mα represents either the mass of the carbon-oxygen 1996) andasymmetric collapse and explosion of supernovae
core at the time of supernova explosion or the mass of the (Lai et al. 1995, 2001; Nordhauset al. 2012). In this paper,
oxygen-neoncore,estimatedbymα=max mCh,0.773mc weassumethatbothcasescancontributetotheacquisition
{ −
0.35 withmCh beingtheChandrasekharmass.Sincemα is of kick velocities by neutron stars.
}
in the range of 1.4 7M⊙, the masses of neutron stars Although Fryer& Kalogera (1997) and
≈ −
are in the range of 1.3-1.8 M⊙. This is consistent with ob- Arzoumanian et al. (2002) suggested a possible bimodal
servationalandtheoreticalconstraints.Lattimer & Prakash distribution for the kick velocities, we here simply assume
(2007)showthatabout83%ofobservedneutronstarshave that the kick velocities have a Maxwellian distribution
mass in the range 1 2M⊙, while 100% of observed neu- following the best estimate of Hansen & Phinney (1997),
−
tron stars have mass in the range 0.8 2.5M⊙. For this takingselection effects into account, with
−
paper, we assume the radius of a neutron star to be 10
km(Lattimer & Prakash(2007)giveempiricalvaluesinthe dN =(2/π)1/2vk2e−vk2/2σk2, (13)
range 9 – 15 km). Ndvk σk3
For the formation of DNS, we assume seven evolution
wherev isthekickvelocity andσ isitsdispersion, dN/N
channels: k k
is the normalized number in a kick velocity bin dv . We
I: CE ejection + CE ejection; k
take σ = 190 km s−1, so that the probable kick velocity
II: Stable RLOF+ CE ejection; k
v 268kms−1.Aswithotherparametersinourpopula-
III:CE ejection + stable RLOF; kp ≈
tion synthesis,aMonteCarlo procedureisused togenerate
IV: Stable RLOF+ stable RLOF;
the individual kick velocities for the neutron stars. Other
V: Exposed core + CE ejection;
parameters associated with the kick velocity (e.g. the di-
VI: Solo CE ejection;
rection) are assumed to follow a uniform distribution. For
VII: Solo RLOF.
comparison, Hansen & Phinney (1997) gives a kickvelocity
A binary with mass ratio q = m /m less than some
1 2 distributionofpulsarswithameanvalueofabout250 300
critical value qc will experience dynamically stable mass km s−1 and σ =190 km s−1. Without considering an−y se-
transfer if the primary fills its Roche lobe while the star k
lection effects, Lyne& Lorimer (1994) found amean pulsar
is in the Hertzsprung gap or on the red giant branch. The
kick velocity of 450 90 km s−1. Hobbset al. (2005) sug-
primary will become a compact object and the orbital sep- ±
gested that there is lack of evidence for the bimodal distri-
aration will change as
bution of kick velocities and a σ =265 km s−1.
k
dlna=2dlnm +2α dlnm +dln(m +m ) (12) Akickaddsacomponenttotheorbitalvelocityofneu-
2 RLOF 1 1 2
−
tron stars imposes, leading to a change of the binary orbit.
whereα isthemass-transferefficiencyforstableRoche-
RLOF Wecorrectthevaluesoforbitalparametersofneutronstars
lobe overflow (RLOF) (Han et al. 1995). Here, we take using Kepler′s laws and the binary dynamics. Wehave
α = 0.5 (Paczyn´ski & Zi´o lkowski 1967; Refsdal et al.
RLOF
1974). Subsequently, if the secondary fills its Roche lobe a(1 e2)= d v 2 /(GM), (14)
o
− | × |
Gravitational waves from double neutron stars 5
where G is gravitational constant, M the total mass of bi- (1993) and Kroupa (2001) constrained by the observations
nary, a the semi-major axis of the orbit, e the eccentricity, of Wielen et al. (1983), Popper (1980) and the Hipparcos
dthedistancevectorbetween thetwostars, andv theor- mission (Creze et al. 1998; Jahreiß & Wielen 1997).
o
bitalvelocityofneutronstar.Theorbitalvelocityv canbe In this paper, we adopt the frequently used power-law
o
expressed as IMF
1v2 =GM 1 1 . (15) ξ(m)= Amσ, 0.16m6100M⊙ (18)
2 o d − 2a
(cid:18)| | (cid:19) where m is the primary mass; ξ(m)dm is the number of
The velocity of the mass centre of the binary vc relative to starsinthemassintervalmtom+dm;Aisthenormaliza-
theold mass centrecan be simplified to tioncoefficientdeterminedbyA 100ξ(m)dm=1.Sincethe
0.1
(M′−∆m′1)vc =m1cvk+∆m′1Mm′2′vo, (16) rtehsautltσs oisfiKnrothueparaentgael.o(f1−9913.3) atRnod−K2.r7o,uwpae(t2a0k0e1σ) i=nd−ica1t.5e
and 2.5 for comparison.
where M′ is the total mass of the binary before supernova −
(iii) Wehaveadopted a metallicity Z =0.02 (Population I)
collapse, ∆m′ themasslossofthe(exploding)primarydue
1 and 0.001 (Population II).
to the supernova collapse and explosion, and m′ the mass
2 (iv) We assume a constant mass-ratio distribution
of secondary. From Eqs.14-15 and Kepler’s laws, the new
(Mazeh et al. 1992; Goldberg & Mazeh 1994),
orbital parameters are determined. We use the new orbital
parameters to calculate the GW emission from DNS. Note n(1/q)=1,0.001<1/q<1. (19)
thatourmodelforpost-kickneutronstarsisconsistentwith
The inverse mass ratio has a minimum value of 0.001 since
the dynamic model in Brandt & Podsiadlowski (1995) and
the maximum and minimum mass of MS stars is 100 and
Hurley et al. (2002).
0.1 M in our simulations.
⊙
(v)Weemploythedistributionofinitialorbitalseparations
used by Han (1998) and Han et al. (2003), where they as-
2.1.5 Initial conditions for population synthesis
sumethatallstarsaremembersofbinarysystemsandthat
InordertoobtainasampleofcompactbinariesintheGalac- the distribution of separations is constant in loga (a is the
ticthindiscwhichiscomparablewithobservations,wehave separation) for wide binaries and falls off smoothly at close
performed a Monte-Carlo simulation in which we need the separations:
six physical inputs described below. In this study, we only
da α ( a )k,a6a ,
investigate the effect of the first three. We use the differ- = sep a0 0 (20)
dn α ,a <a<a .
ent cases in our simulations to obtain information on the (cid:26) sep 0 1
population of DNS and their GW radiation in otherGalac- where α 0.070, a = 10R , a = 5.75 106R =
sep 0 ⊙ 1 ⊙
tic components via mass-scaling. The Galactic structure is 0.13 pc, k ≈1.2. This distribution implies that×the number
described in 2.1.6, and the results are presented in 3.4. of binary s≈ystems per logarithmic interval is constant. In
§ §
(i) We adopt three star formation (SF) models in our bi- addition,approximately50percentofallsystemsarebinary
nary population synthesistosee theinfluenceof SFhistory starswithorbitalperiodsoflessthan100yr.Thesebinaries
on thepopulation of compact binaries. These are: areexcludedwhen,duringevolution,theconditionthatthe
Instantaneous SF:asinglestarburstattheformationofthe sum of their radii exceeds their initial orbital separation is
thin disc with a constant SF rate of 132.9M⊙ yr−1 from 0 satisfied.
to 391 Myr.followed by noSF from 391 Myr to10 Gyr; (vi) The distribution of initial eccentricities of binaries fol-
Constant SF: SF occurs at a constant rate of 5.2M⊙ yr−1 lows Pe =2e (Heggie 1975; Nelemans et al. 2001).
from 0 to 10 Gyr;
Quasi-exponential SF:thestarformation rateS isthecom-
bination of a major star-forming process (the first term of 2.1.6 Galactic mass distribution - potential, and
thefollowingfunction)andaminorstarformation (thesec- rotational velocity
ond term of the function),
We consider the Galaxy to comprise three components,
S(t )=7.92e−(tsf)/τ +0.09(t ) M yr−1 (17) namelythebulge,disc,andadark matterhalo. Weassume
sf sf ⊙
thatthepositionoftheSunisgivenbyitsdistanceftomthe
where tsf is the time of star formation, and τ = 9 Gyr Galactic centreRsun =8.5kpcand height abovetheGalac-
(Yu& Jeffery 2010), which produces ≈3.5 M⊙yr−1 at the tic plane zsun = 16.5pc (Freudenreich 1998). Our approxi-
current epoch. mation for the Galactic density distribution is summarized
All three models produce a thin-disc star mass of inTable2.Thedetailedexpressionsaredescribedasfollows.
5.2 1010 M⊙ at the thin-disc age t = 10 Gyr. Obse∼r- (1) We adopt a normal density distribution for
×
vations place the current thin-disk SF rate in the range the spherical bulge with a cut-off radius of 3.5kpc
3 5M⊙yr−1 (Smith et al. 1978; Timmes et al. 1997; (Nelemans et al. 2004),
≈ −
Diehl et al.2006)andimplythatInstantaneous SFishighly
implausible for the SFH in thethin disc. It is retained here ρ (r)= Mb e−(r/r0)2 M pc−3, (21)
b 4πr3 ⊙
for comparison. 0
(ii) The initial mass function (IMF) can be constrained where r is the radius from the center of the Galaxy, r =
0
by the local luminosity function, stellar density and po- 0.5kpcisthebulgescalelength,andM =2.0 1010 isthe
b
×
tential. The IMF for the Galactic components may be dif- massofbulge.Robin et al.(2003)suggestthatthestructure
ferent as indicated by Robin et al. (2003), Kroupa et al. of the inner bulge (< 1◦ from the Galactic center) is not
6
Table2.Densitylawsandassociatedparameters.risthespher-
ical radius from the center of the Galaxy and r0 is bulge scale
length;Randz arethenaturalcylindricalcoordinatesoftheax-
isymmetricdisc,hR isthescalelengthofthedisc,hz isthescale 350
heightofthethindisc,h′ isthescaleheightofthethickdisc;a′ observational estimate
z
is the radius of the halo and a′0 is a constant; ρc is the central 300 Brand & Blitz (1993)
massdensity.
m/s) 250 total
k
densitylaw co(nksptacn)ts (M⊙ρpcc−3) city ( 200
Bulge e−(r/r0)2 r0=0.5 4Mπrb03 =12.73 ar velo 150 halo with dark matter
ul
Thindisc e−R/hRsech2(−z/hz) hR=2.5 4πMh2Rtnhz =1.881 Circ 100 thin+thick disc
hz=0.352
Thickdisc e−R/hRe−z/h′z hh′zR==1.21.558 4πMh2Rtkh′z =0.0286 50 bulge
0
Halo [(1+(aa′0′)2)]−1 a′0=2.7 0.108 0 10 20 R (kpc) 30 40 50
Figure 1. Rotational velocity as a function of galactocentric
yet well constrained observationally. Consequently we here distanceRfromtheGalaticmodel,showingthecontributiondue
focus on the outer bulge and make no allowance for any to different components, i.e. the bulge, thin disc + thick disc,
additional contribution to the compact-binary population and halo including dark matter. The dashed line indicates the
from the central region. observational estimate by Brand&Blitz (1993). The spheroidal
We use the potential proposed by Miyamoto & Nagai component due to the interstellar medium was not considered
(1975) in cylindrical coordinates to calculate the rotational separately.
velocity of stars in thebulge.
(2)Wemodelthethinandthickdisccomponentsofthe
Galaxy using a squared hyperbolic secant plus exponential Caldwell & Ostriker (1981). The observational estimate by
distribution expressed as: Brand & Blitz (1993) is included for comparison.
Fig.1demonstratestheinfluenceoftheGalacticmodel
M
ρ (R,z)= d e−R/hRρ(z) M pc−3, (22) on the rotation curveof theMilky Way.
d 4πh2Rh ⊙ From the Galactic model, the total mass of the halo
including dark matter is 4.5 1011M inside a sphere of
whereRandz arethenaturalcylindricalcoordinatesofthe ⊙
×
axisymmetric disc, and h = 2.5kpc is the scale length of radius50kpc.Weonlyfocusonthebaryonicmassinthehalo
R
thedisc, h=hz for thethin disc, h=h′z for the thick disc, which is considered to be 5×1010M⊙, constrained by the
andM =M =5.2 1010 M isthemassofthethindisc; density of double white dwarfs (Yu & Jeffery 2010). With
d tn ⊙
M = M = 2.6 1×09 M is the mass of the thick disc. the SFR adopted here, the baryonic mass in the bulge and
d tk ⊙
ρ(z) is thedistribu×tion in z,with: disc is at least 2 1010M⊙ and 5.5 1010M⊙ respectively,
× ×
implyingthatourmodelassumesnodarkmattercomponent
ρ(z)=sech2( z/hz) (thin disc) (23) localized to thebulge or the thin disc.
−
Combining the Galactic model and the mass of the
and
Galactic components, the stellar density in the solar neigh-
ρ(z)=e−z/h′z (thick disc), (24) bourhood is 0.064M⊙pc−3, of which 6.27 10−2M⊙pc−3 is
×
inthethindisc,9.4 10−4M pc−3 isin thethickdisc, and
⊙
where hz = 0.352kpc is the scale height of the thin disc 2.18 10−5M⊙pc−3×is in the halo. This is consistent with
and h′z =1.158kpc is the scale height of the thick disc. We the H×ipparcos result, 0.076 0.015M⊙pc−3 (Creze et al.
neglect the age and mass dependenceof the scale height. 1998). The local dark matter±density in our model is about
The Miyamoto &Nagai potential in cylindrical coordi- 0.01M pc−3.
⊙
natesisalsousedtocalculatetherotationalvelocityofstars
in thedisk.
(3)For the halo, we employ a relatively simple density
2.2 Gravitational wavesfrom double neutron stars
distribution which is consistent with Caldwell & Ostriker
(1981), Paczynski (1990) and Robin et al. (2003): We calculate the GW strain amplitude from a single NS-
NS binary by linearizing the equations of general rela-
a′ 2 −1 tivity (Peters & Mathews 1963; Landau & Lifshitz 1975).
ρ (a′)=ρ 1+ , (25)
h ch (cid:18)a′0(cid:19) ! We assume the following properties of a binary. 1) The
masses of the two components are m and m respectively.
1 2
where a′ is the radius of the halo, ρ = 0.108 M pc−3 Therefore the total mass M = m +m , and the reduced
ch ⊙ 1 2
and a′ = 2.7 kpc. mass µ = m m /M. 2) The semi-major axis of the or-
0 1 2
For the dark matter halo, we adopted the potential of bit is a. 3) The eccentricity is e. 4) From Kepler, the or-
Gravitational waves from double neutron stars 7
bital separation d = a(1−e2) and the angular velocity g(n,e)= n4 [J (ne) 2eJ (ne)+ 2J (ne)
1+ecos(ϕ) 32{ n−2 − n−1 n n
ϕ˙ = (GMa3)1//22(1[1−+ee2c)o3/s(2ϕ)]2, where ϕ is the angle between the +2eJn+1(ne)−Jn+2(ne)]2 (32)
orbitalseparationandthex-directionofanarbitraryCarte-
+(1 e2)[J (ne) 2J (ne)+J (ne)]2
sian coordinate in the plane of the binary orbit, and G is − n−2 − n n+2
4
thegravitationalconstant.Wetakethez directionperpen- + [J (ne)]2 ,
dicular to the orbit plane and the origin−is at the center of 3n2 n }
mass.Clearly, 02πdϕ= 0Porbϕ˙dt,wherePorb istheorbital whereJn(ne)areBessell functionsofthefirstkindandn=
period. 1,2,3,.... Since the sum of the power in each harmonic is
R R
The components of the strain amplitude can be ex- equalto the total power emitted from the binary,we have
pressed as (Landau & Lifshitz 1975)
∞
g(n,e)=z(e). (33)
 hhhxyzxxx hhhxyzyyy hhhxyzzzz =−3c24GRb  AAA¨¨¨yxzxxx AAA¨¨¨yxzyyy AAA¨¨¨yxzzzz , equaAtifotnesraabomvea(tNheemlneXam=t1aicnaslettraaln.s2fo0r0m1;atYioun&bJaesffeedryo2n01t0h)e,
   (26)
weobtain thestrain amplitudeh(n,e)inthevicinityof the
where A¨ represents the second order differential of the
αβ Earth at GW frequency f in thenth harmonic as
mass-quadrupole tensor with respect to time, suffix αβ de- n
notesthedirection,Rbisthedistanceformtheobserver,and h(n,e)≡hn
c is thespeed of light. Sincethe mass quadrupoletensor is G5/3 g(n,e) 1/2
=4√2(2π)2/3 M2/3µP−2/3R−1
A A A c4 orb b n2
xx xy xz (cid:18) (cid:19)
Ayx Ayy Ayz =1.14 10−21
  ×
A A A
zx zy zz g(n,e) 1/2 5/3 P −2/3 R −1
 µd2(3cos2ϕ 1) 3µd2cosϕsinϕ 0 M orb b ,
= 3µd2sinϕco−sϕ µd2(3sin2ϕ 1) 0 , ×(cid:18) n2 (cid:19) (cid:18)M⊙(cid:19) (cid:18) h (cid:19) (cid:18)kpc(cid:19)
 −  (34)
0 0 µd2
−
 (27)
f =n/P , (35)
n orb
we have where µ3/5M2/5 is the so-called chirp mass. Eqs.31,
M ≡
A¨ = C 6(cos2ϕ+ecos3ϕ)+A¨ , 32, 34, and 35 are themain equations used to calculate the
xx − 2 zz power and strain amplitude of the GW signal from one in-
A¨yy =C26(e2+ecosϕ+cos2ϕ+ecos3ϕ)+A¨zz, dividual NS-NS pair in frequency space. These equations
A¨ = C 2(e2+ecosϕ), also tell usthat thepower and strain amplitudeof theGW
zz 2
− signal from a binary consisting of two point-masses can be
A¨xy = C26(sin2ϕ+2esinϕcos2ϕ+esin3ϕ), (28) determinedbyfourparameters -thechirpmass, orbitalpe-
−
A¨ =A¨ , riod,eccentricityanddistance.Inthispaper,werefertothe
yx xy
first three as orbital parameters, and use the mass of each
GMµ
C2 = a(1 e2). component instead of thechirp mass.
− The energy fluxof GW waves can be expressed as
From Eqs.26, 27and 28,it can beshown that theGW c3Ω2 c3πf2
from a binary with circular orbit is monochromatic with F = h2 = h2, (36)
16πG 4G
frequency2/P . Whentheorbit iseccentric, thewavebe-
orb
whereΩ=2πf istheangularfrequency.Wedefinethespec-
comespolychromaticandthereisafrequencybroadeningin
tral function of theenergy fluxas
theGW signal.
The average power of theGW radiated from two point 1 dF 2π2 d(f2h2)
S = = , (37)
masses over one orbital period can be obtained by solv- ρ c3 df 3H2 df
c 0
ing the third order differential of the mass-quadrupole
tensor with respect to time (Peters & Mathews 1963; where ρc = 3H02/8πG is the critical mass density of the
Landau & Lifshitz 1975). Wequotetheresult present universe with H0 73 km s−1 Mpc−1 being the
≈
Hubbleconstant (Freedman & Madore 2010). We show the
32 G4µ2M3 GW energy flux spectrum for selected cases to assist com-
L = z(e), (29)
GW 5 c5 a5 parison with other work which uses this metric.
(cid:18) (cid:19)
1+(73/24)e2+(37/96)e4 2.3 The sensitivity of eLISA
z(e)= . (30)
(1 e2)7/2
− Evolved-LISA(eLISA)isdesignedas anupdatedversionof
AfterFourieranalysisofKeplermotion,weobtainthepower LISA, a space-based GW detector, consisting of 1 mother
in thenth harmonic (Peters & Mathews 1963) and 2 daughter satellites flying in formation to form a
Michelson-type Laser interferometer with an arm length
32 G4µ2M3 of 1 106 km (see http://www.elisa-ngo.org/). Noise arises
Ln = g(n,e), (31) ×
GW 5 c5 a5 mainly from the displacement noise (including the noise
(cid:18) (cid:19)
8
caused bylasertrackingsystemand otherfactors) andpar-
Table3.Mainparametersandtheirvaluesinoursimulation.See
asitic forceson theproofmass ofanaccelerometer (acceler- §2.1foranexplanation oftheparameters.
ationnoise) (Larson et al.2000)3.Wecanconvertthenoise
signal to an equivalent GW signal in frequency space by
h =2 Sn, (38) CEE(αλ=) CEA(γ=)
f R 1.0 0.5 1.3 1.5
r
where Sn is the total strain noise spectral density,hf is the SFH IMF(σ=)
root spectral density and R is the GW transfer function Con −1.5 C1 C3 C25 C27
givenbyLarson et al.(2000).Inthesimulations,wetakethe −2.5 C2 C4 C26 C28
displacementnoisetobe1.110−11 mHz−1/2at10mHz,and Z=0.02 Exp −1.5 C5 C7 C29 C31
theaccelerationnoisetobe310−15 ms−2Hz−1/2 at10mHz. −2.5 C6 C8 C30 C32
By comparison, the arm length of LISA would have been Inst −1.5 C9 C11 C33 C35
5 106 km. The displacement noise and acceleration noise −2.5 C10 C12 C34 C36
wouldhavebeen410−11mHz−1/2and310−15ms−2Hz−1/2. SFH IMF(σ=)
For a continuous monochromatic source, such as a NS- Con −1.5 C13 C15 C37 C39
NS binary with a circular orbit, which is observed over a −2.5 C14 C16 C38 C40
time T, the root spectral density will appear in a Fourier Z=0.001 Exp −1.5 C17 C19 C41 C43
spectrum asa single spectral linein theform (Larson et al. −2.5 C18 C20 C42 C44
2000) Inst −1.5 C21 C23 C45 C47
−2.5 C22 C24 C46 C48
h =h√T. (39)
f
So,foran observation timeT =1yr,theroot strain ampli-
tudespectral density h =5.62 103h Hz−1/2. (i) Foreachstarformationepochinthedisc,wecalculate
f
To demonstrate the detecta×bility of the predicted GW asampledistributionofk coevalMSbinarieshavingatotal
signal duetothegalactic DNSpopulation,we show theex- massmp andgeneratedbythefourMonte-Carlosimulation
pected eLISA sensitivity for S/N=1 in thefigures in 3. parameters m, q, a and e.
§ (ii) Wefollowtheevolutionofeachprimordialbinaryina
time grid consisting of many time intervals to establish the
2.4 Data reduction propertiesofDNSsformedfromtheaboveMSbinariesupto
t .FromthetimescalesfromMSbinaryformationtoDNS
We reduce the simulated GW signal by using its mean in- disc
formation andfrom DNSformation toDNSmerger,wecan
tegrated value h . To do this, we first choose a frequency
interval∆f′ whhichi isgreaterthantheintervalusedtocom- obtain the contribution function from each star formation
pute the simulations (i.e. ∆f = 1yr−1). We then calculate epoch to the number of new-born and merged DNSs in a
giventimeinterval,aswellasthetotalnumberofaliveDNSs
themeanvalue h ofthestrainamplitudeanditsstandard
h i duringthat interval.
deviation σ in this large frequency intervalusing
hhi
(iii) By summing the contributions from all star forma-
h = ji=1hi,σ2 = ji=1(hi−hhi)2, (40) tion epochs in the thin disc, we can obtain all the physi-
h i j hhi j cal information of theDNSpopulation, including birth and
P P
where j represents the number of small-frequency inter- merger rates, present number, and distributions of orbital
vals in the large frequency interval. In this paper, we take parameters.
∆logf′ =0.03, soj isalsoafunctionoffrequency.Weplot (iv) Since most of the DNSs have eccentric orbits, we
h asafunctionofGWfrequencyinallfiguresunlessspec- compute the Fourier transform of each orbit to obtain the
hifieid otherwise. In each panel, we also show the maximum valueof GW strain emitted ateach harmonic frequency(as
standard deviation which represents the maximum uncer- indicated by Eqs.32 and 34), and then sort the DNSs by
tainty in each large frequency interval in different cases. the harmonics of orbital frequency (equivalent to the GW
frequency).
(v) Wethencalculate thetotal strain amplitudeh2 from
2.5 Computation procedure
thenumberanddistanceofDNSsineachfrequencybin(for
InordertocomputethesuperpositionoftheGWsignalfrom oneyearobservationofe-LISA,see 2.3).Thisyieldsaraw
§
theentireDNSpopulationintheGalacticdiscandcompare data of strain amplitude against theGW frequency.
thesignalwiththesensitivityoftheproposeddetectors,we (vi) Finally, we use the method described in 2.4 to re-
§
needtoknowtheirbirthrates,mergerrates,presentnumber, duce the raw data, producing a reduced relation between
space, mass, eccentricity and orbital distributions. we have thestrain amplitude and GW frequency.
adoptedthefollowingproceduretoobtaintheabovephysical
Inourpopulationsynthesissimulation,westartedwith
properties. For theGalactic thin disc with age tdisc, having 16 107 primordial MS binaries, distributed equally between
a total mass M (t )= tdiscS(t )dt 4:
tn disc 0 sf sf all48cases,andassumethattdisc=10Gyr.Theparameters
R in the present study are SFH (Instantaneous, Continuous
3 The values of parameters to calculate the total noise can be and Quasi-exponential), IMF (σ = −1.5 and −2.5), metal-
licity (Z =0.02 and 0.001), and two different CE evolution
foundonhttp://www.elisa-ngo.org/.
4 Wehereneglecttheinterstellarmediumwhichmakesupabout processes (α and γ). For each CE formalism we adopt two
20%ofthethindiscmass(Robinetal.2003). parameters (αλ = 1.0 and 0.5, γ = 1.5 and 1.3). The 48
Gravitational waves from double neutron stars 9
Figure 2.BirthratesofDNSsintheGalacticdiscindifferentcases.
10
Figure 3.MergerratesofDNSsintheGalacticdiscindifferentcases.