January 29, 2010 The search for spinning black hole binaries in mock LISA data using a genetic algorithm Antoine Petiteau, Yu Shang and Stanislav Babak Max-Planck-Institut fuer Gravitationsphysik, Albert-Einstein-Institut, Am Muchlenberg 1, D-14476 Golm bei Potsdam, Germany∗ Farhan Feroz Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK Coalescing massive Black Hole binaries are the strongest and probably the most important gra- vitational wave sources in the LISA band. The spin and orbital precessions bring complexity in 0 the waveform and make the likelihood surface richer in structure as compared to the non-spinning 1 case. We introduce an extended multimodal genetic algorithm which utilizes the properties of 0 the signal and the detector response function to analyze the data from the third round of mock 2 LISA data challenge (MLDC 3.2). The performance of this method is comparable, if not better, n to already existing algorithms. We have found all five sources present in MLDC 3.2 and recovered a the coalescence time, chirp mass, mass ratio and sky location with reasonable accuracy. As for J the orbital angular momentum and two spins of the Black Holes, we have found a large number of 9 widely separated modes in theparameter space with similar maximum likelihood values. 2 PACSnumbers: 04.30.Tv,04.80.Nn ] c q - r g [ 1 v 0 8 3 5 . 1 0 0 1 : v i X r a ∗[email protected], [email protected],[email protected],[email protected] 2 I. INTRODUCTION The existence of massive black holes (MBH) with masses ranging from 105M to 109M in majority of galactic ⊙ ⊙ centers is confirmed by several observations (see [1] and reference therein) with the MBH in our Galaxy, SgrA* [2] being the best example. Mergers of galaxies are common events in the Universe, it is believed that each galaxy has had at least one merger event during its life. During these mergers, the MBH of each galaxy is driven to the center of the remnant full of stars and gas via dynamical friction [3]. The pairs of MBH separated by about kiloparsec are observed in some active galactic nuclei such as NGC6240 [4], Arp299 [5], ESO509-IG066 [6], Mrk463 [7] and J100043.15+020637.2 [8]. The interaction with the gas disc can bring the binary on a tighter orbit down to a few parsecs in a reasonable amount of time (few Myr) [9]. There are few candidates of MBH binaries on the sub-parsec scale: the quasars OJ287 [10] ( 0.05 pc) and SDSSJ092712.65+294344.0 [11, 12]( 0.1-0.3 pc). To overcome the ∼ ∼ last parsec separating the MBHs and bring them to the efficient gravitational wave (GW) driven inspiral several scenario have been proposed. Here are few possibilities: rotation of the merging galaxies and triaxial potential [13], processes involving gas [14], resonant relaxation [15], massive perturber [16], young compact stars cluster [17], effect from IMBH [18], etc. When the separation is less than 10−3 pc, the binary evolution is efficiently driven by the gravitationalradiation and can reach the coalescence in less than 109 years. The GWs emitted by the binary at the end of the inspiral phase followed by the merger and the ringdown will be detected by the future space born mission LISA with a very high signal-to-noise ratio (SNR). The MBH binaries are the strongestsourcefor LISA andseveralsuchevents per yearare expected[19,20]. It is believedthat almostallthe MBHs are spinning. However the predictions for the magnitudes and directions of the spins of MBHs in the binary systems differ largely depending on the models, the environment of the binary and the physical processes involved (coherentaccretion,alignmentofspins with the gasdisc [21–23], sequenceof randomlyorientedaccretionevents [24], etc). Severaltechniques to measurethe spins using electromagnetic radiation[25, 26] havebeen proposedto measure thespinsofMBHbinaries,butuncertaintiesarestillverylarge. ThegravitationalwaveobservationsofMBHbinaries with LISA should enable us to measure masses and spins of MBHs in the binary with unprecedented accuracy [27]. The knowledge of spins could give us a lot of information about the kick velocity of remnant MBH, the engines of active galactic nuclei, the mechanisms involved in galactic centers, etc. Finally, the spin measurements combined with a precise estimation of masses and sky positionmade with LISA will increase our understanding of the originof MBHs, the galactic evolution, the galactic center, cosmology, etc. Severalalgorithmsfordetectingnon-spinningMBHbinariesinsimulatedLISAdatahavealreadybeendemonstrated [28–31]. In this paper we consider inspiralling spinning MBH binaries and we present a particular adaptation of the Genetic Algorithm (GA) to search for GW signals from those systems. Genetic algorithms belong to the family of optimization methods, i.e. they look for extrema. The first application of GA in LISA data analysis was proposed in [32]forGalacticbinaries. Inthis methodthewaveformtemplateisassociatedwithanorganism,andparametersplay the role ofthe setofgenesdefining this organism. The logarithmoflikelihoodobtainedwith a giventemplate defines the quality of the organism. A set (colony) of organisms is then evolved through breeding, mutation and custom designedacceleratorswiththe aimoffinding the genotypewiththe highestquality. This correspondstothe standard Darwin’s principle: weak perishes, strong survives, or, translated into the conventional data analysis language: by evolving a set of templates, we are searching for the parameter set that maximizes the likelihood. We have applied the GA to the analysis of the third round of mock LISA data challenge. The mock data set 3.2 consistedoftheGaussianinstrumentalnoise,Galacticbackgroundandbetweenfourtosixsignalsfromtheinspiralling spinning MBH binaries in a quasi-circular orbit [33]. We have found several almost equal in value maxima of the likelihood which are widely separated in the parameter space. We search for each such strong maximum, which we call mode, and then explore it by a designated set of organisms. We refer to this extension of the standard GA as a multimodal GA. The mutlimodal GA applied to the blind search has shown an excellent performance: we have detected all present signals with a very accurate estimation of the parameters (were possible). The structure of this paper is as follows. In the next Section II we describe a model of the signal and the search template we utilized. Thenwe discuss the detection statistics we associatewith the quality ofthe organismandtheir efficientmaximizationoverasubsetofparameters. InSectionIIIweintroduceaGeneticAlgorithmanditsparticular implementation in GW data analysis. Then, in the Section IV we introduce accelerators for the rapid evolution – mechanisms which allow efficiently explore the parameter space. Besides few standard accelerators often used in GA we have introduced few new ones, specific to the MBH binary searchproblem. We introduce multimodal GA (MGA) in Section V and describe how to identify the modes and to follow their evolution. Our complete algorithm as part of the search pipeline is presented in Section VI which describes the search pipeline. The results are presented and discussed in the Section VII and finally we give a short summary in the concluding Section VIII. 3 II. FORMULATION OF THE PROBLEM The mock LISA data challenges are organized in order to encourage the development of efficient algorithms for gravitational wave data analysis and to evaluate their performance. The third round of MLDC consisted of five challenges but in this work we focus our attention here on the data set 3.2 which contained GW signals from 4-6 binaries of spinning MBHs (exact number was not revealed to the participants), on top of the confusion Galactic binariesbackgroundandthe instrumentalnoise. Thesedatawasanimprovementuponthe MLDC challenges1.2and 2.2 by adding spins to MBHs. The spin-spin and spin-orbital coupling causes the orbital and spins precession which results in the modulation of the amplitude and phase of the GW signal. The prior range on the parameters and the detailed set up of the challenge can be found in [33]. A. Model of the template The signal used in MLDC is modeled as the amplitude-restricted waveform (i.e. only dominant harmonic at the leading order is used) with the phase taken up to the second Post-Newtonian (PN) order with the leading order contributions from the spin-orbital and spin-spin coupling. The binary evolution is described as a quasi-circular adiabatic inspiral. The waveformis describedby fifteen parameterswhich are: the two massesm and m , the time atcoalescencet , 1 2 c the sky location of the source in ecliptic coordinates, latitude β and longitude λ, the dimensionless spin parameters, χ and χ , the initial direction of the spins, polar angles θ and θ and azimuthal angles φ and φ , the initial 1 2 S1 S2 S1 S2 direction of the orbital angularmomentum, polar angle θ and azimuthal angle φ , the phase at coalescence Φ , and L L c the luminosity distance D . L We denote the unit vector in the direction of the orbital angular momentum as Lˆ and two spins are S~ =χ m2Sˆ , 1 1 1 1 S~ =χ m2Sˆ , where Sˆ are unit vectors and 0<χ <1. The precession equations are given in [34] 2 2 2 2 1,2 1,2 S~˙ = (Mω)2 η(Mω)−1/3 4+ 3m2 Lˆ+ 1 S~ 3(S~ Lˆ)Lˆ S~ , (1) 1 2M m M2 2− 2· × 1 (cid:26) (cid:18) 1 (cid:19) h i(cid:27) S~˙ = (Mω)2 η(Mω)−1/3 4+ 3m1 Lˆ+ 1 S~ 3(S~ Lˆ)Lˆ S~ , (2) 2 2M m M2 1− 1· × 2 (cid:26) (cid:18) 2 (cid:19) h i(cid:27) Lˆ˙ = (Mω)1/3 S~˙ +S~˙ = ω2 4+ 3m2 S~ + 4+ 3m1 S~ Lˆ − ηM2 (cid:16) 1 2(cid:17) 2M((cid:20)(cid:18) m1 (cid:19) 1 (cid:18) m2 (cid:19) 2(cid:21)× 3ω1/3 S~ Lˆ S~ + S~ Lˆ S~ Lˆ . (3) −ηM5/3 2· 1 1· 2 × ) h(cid:16) (cid:17) (cid:16) (cid:17) i The modulation of the waveformdue to the presence of spins is taken at the leading order. The orbital angular frequency with spin effect up to 2 PN order is given by 1 743 11 3 β Mω = τ−3/8 1+ + η τ−1/4 π τ−3/8 8 " (cid:18)2688 32 (cid:19) − 10(cid:18) − 4(cid:19) 1855099 56975 371 3 + + η+ η2 σ τ−1/2 , (4) 14450688 258048 2048 − 64 # (cid:18) (cid:19) where M =m +m is total mass, η = m1m2 is the symmetric mass ratio and 1 2 M2 η τ = (t t) (5) c 5M − 1 m2 β = χ Lˆ Sˆ 113 i +75η (6) 12 i · i M2 iX=1,2(cid:20) (cid:16) (cid:17)(cid:18) (cid:19)(cid:21) 1 σ = ηχ χ 247 Sˆ Sˆ 721 Lˆ Sˆ Lˆ Sˆ . (7) 1 2 1 2 1 2 −48 · − · · h (cid:16) (cid:17) (cid:16) (cid:17)(cid:16) (cid:17)i 4 In our following analysis, we use η and the chirp mass M = (m1m2)3/5 as independent parameters instead of m c (m1+m2)1/5 1 and m . 2 The intrinsic phase is τ5/8 3715 55 3 Φ =Φ 1+ + η τ−1/4 (4π β)τ−3/8 orb C − η " (cid:18)8064 96 (cid:19) − 16 − 9275495 284875 1855 15 + + η+ η2 σ τ−1/2 . (8) (cid:18)14450688 258048 2048 − 64 (cid:19) # The phase is defined with respect to the ascending node, however the spin-orbital coupling causes precession of the orbit, therefore we need to introduce precessional correction to the phase, δΦ(t). It depends on the unit vector nˆ pointing from the solar system barycenter to the source: Φ(t)=Φ (t)+δΦ(t), (9) orb (Lˆ nˆ)(Lˆ nˆ) Lˆ˙ Φ˙(t)=ω+δΦ˙ =ω+ · × · , (10) 1 (Lˆ nˆ)2 − · δΦ(t)= tc Lˆ·nˆ (Lˆ nˆ) Lˆ˙dt. (11) −Zt 1−(Lˆ·nˆ)! × · The gravitationalwave polarization components in the source frame are given by 2Mη h =h cos2Φ= (1+cos2ι)(Mω)2/3cos2Φ + +0 − D L 4Mη h =h sin2Φ= cosι(Mω)2/3sin2Φ, (12) × ×0 D L where cosι=Lˆ nˆ. · The strain h(t) that the GW exerts on the LISA detector is the following linear combination of h and h + × h(t)=F (β,λ)hS(t)+F (β,λ)hS(t), (13) + + × × where F and F are “detector beam-pattern” coefficients. The polarization components in the radiation frame, hS + × + and hS, are expressed as × hS = h cos2ψ h sin2ψ, + − + − × hS =h sin2ψ h cos2ψ, (14) × + − × where the polarization angle ψ is defined by sinβcos(λ φ )sinθ cosθ cosβ L L L tanψ = − − . (15) sinθ sin(λ φ ) L L − Note, due to the precession of the orbital plane, this polarization angle varies during the evolution of the binary. The data sets in MLDC are the TDI (time delay interferometry)variables. These are the combinations ofthe time delayed measurements which drastically reduce the laser frequency noise [35, 36]. In our search, we adopted the two orthogonal(i.e. with uncorrelated noise) TDI channels, A and E, in the phase domain (i.e. strain). In our template, weconsideralongwavelengthapproximationtothesesignals[37,38]. Thisapproximation(Lω 1,whereLisLISA’ ≪ armlength and ω is an instantaneous frequency of GW) works pretty well below approximately 5 mHz. Assuming rigid LISA with equal arms, the waveformsampled at discrete times takes the following form [37, 38] h (t ) 2L sin∆φ (t ) I k 2L k ≃ × h (t )[cos(2ψ(t ))F (t) sin(2ψ(t ))F (t)]sinφ′(t ) +0 k k +I k ×I k {− − +h (t )[sin(2ψ(t ))F (t)+cos(2ψ(t ))F (t)]cosφ′(t ) , (16) ×0 k k +I k ×I k } 5 where I = A,E , ∆φ (t) =(φ(t ) φ(t 2L))/2, φ′(t) = (φ(t )+φ(t 2L))/2 with φ(t) being the phase of 2L k k k k { } − − − GW and tk =t+n.−→r0 is the time in LISA frame with −→r0 the vector from the Sun to LISA barycenter. The antenna pattern functions F and F corresponding to the TDI channel, have the following expressions +I ×I b 1 F (θ ,λ ;t,Ω)= [6sin(2θ )(3sin(Φ (t)+λ +Ω) sin(3Φ (t) λ +Ω)) + d d d T d T d 32 − − 18√3sin2θ sin(2Φ (t)+Ω) √3 1+cos2θ d T d − − × (sin(4Φ (t) 2λ +Ω)+9sin(2λ +Ω))], (17) T − d d (cid:0) (cid:1) 1 F (θ ,λ ;t,Ω)= √3cosθ (cos(4Φ (t) 2λ +Ω) × d d d T d 16 − − 9cohs(2λ +Ω))+6sinθ (cos(3Φ (t) λ +Ω)+ d d T d − 3cos(Φ (t)+λ +Ω))] (18) T d with θ = β +π/2, λ = λ+π, Φ (t) = 2πt/Year and Ω = 0, π/2 for channels A and E respectively. Note that d d T − this is the long wavelength approximationto the signal injected in the simulate data, we found it to be a reasonably accurate representation until the last 1 1/2 cycles before the merger. The end of the signal is discussed in more − detail later. B. quality estimation The signal from one detector is s (t)=h (t,θˆ)+n (t), (19) i i i where h (t,θˆ) is a signal described by a set of parameters θˆand n (t) is the stationary Gaussian noise characterized i i by the power spectral density (PSD) S (f), n S (f)δ(f f′)=2n˜(f)n˜(f′). n − Heretheover-barmeanstheaverageoverensembleofthenoiserealizationsandthetildedenotestheFouriertransform defined as ∞ n(f)= dt n(t)e2πif. (20) Z−∞ We use the maximum likelihood method (MeL) [39–41] for estimating the parameters of the waveformθ. Assuming thatthe noise n(t)isa normalprocesswith zeromean, thenthe likelihood(probability)ofthe presenceofsignalh(θˆ) in the detector output is given by b p(sθˆ)=Ce−hs−h(θˆ)|s−h(θˆ)i/2, (21) | where h s is the inner product of the signal given by h | i ∞ h(f)s∗(f)+h∗(f)s(f) h s =2 df . (22) h | i S (f) Z0 n ThenoiseisthesumofinstrumentalnoiseSinst.(f)anedtheeGWconefusionenoisefromGalacticbinariesSGal. Bin.(f). n n Instraindata(i.e. phasemeasurements),theinstrumentalnoiseforTDIvariablesAandEisdescribedbythefollowing PSD Sinst.(f)=16sin2(φ ) SOPN(f)+(cos(2φ ))Sacc.(f) n L n L n −4sin(2φL)(cid:0)sin(φL) 4SnOPN(f)+Snacc.(f) (cid:1), (23) nwohiesreeatrheeSaOcPceNl(efra)t=ion3.n6o7i5seis10S−na4c2c.(Hf)z−=1.5.75×10−53(f−4+1(cid:0)0−8f−6)Hz−1 andth(cid:1)e opticalpathnoise andthe shot n × The Galactic GW confusion noise is a combination of the unresolved signals from 30 millions of white dwarf ∼ binaries. This noise is modeled by the following function, in units of Hz−1, [42, 43] 10−44.62f−2.3 10−4 f 10−3 ≤ ≤ 10−50.92f−4.4 10−3 f 10−2.7 SnGal. Bin.(f)= 10−62.8f−8.8 10−2.7≤ f≤ 10−2.4 (24)  10−89.68f−20 10−2.4 ≤f ≤10−2 ≤ ≤  6 1. Maximization of the likelihood The goal of the method presented in this paper is to find the maximum of the likelihood in the 15-dimensional parameter space, and, thus, obtain the maximum likelihood estimation of the parameters of the signal. The value of the likelihood tells us also about the statistical significance of the detected event. In case of LISA data, the signals usually have high signal-to-noise ratio (SNR), so the probability of the false detection is rather low. However, the data is signal dominated and several GW signals of one type (say, Galactic binaries) could conspire and produce significantlyhighSNR atthe outputofthe matchedfiltering during the searchfor anothertype ofsignal(say,SMBH binary) [44]. It is possible to maximize the likelihood analytically over two parameters and we will call the resulting function MaximizedLikelihood(orquality). Theprocedureofmaximizationissimilartotheoneusedtoproducethe -statistic F [29, 45]. Due to amplitude modulation we can only maximize over two parameters: the luminosity distance D and L the phase at coalescence Φ . c We will be working with the logarithm of the likelihood (getting rid of the normalization factors which does not depend on parameters by adopting the relative likelihood): 1 ln s h h h , (25) I I I I L≃ h | i− 2 h | i I I X X here the sum is overthe independent detectors (TDI streams, I = A,E ). The GW template (16) canbe express as { } h = a h (26) I k kI k X and the extrema of ln over a are found by solving the coupled set of equations k L ∂ln L =0 (27) ∂a k and it turns out that a =X (M )−1 withX = sh , M = h h . (28) k k kl k kI kl kI lI h | i h | i I I X X The log likelihood maximized over a is called the -statistic in the case where a are four functions of (constant) k k F amplitude, polarization angle, inclination and initial phase: 1 X (M )−1X . (29) k kl l F ≃ 2 InthecaseofspinningMBHbinary,theanalyticmaximizationispossibleonlyovertheluminositydistanceD and L the phase at coalescence φ . Consequently, the dimension of the search parameter space is reduced to 13. Following c (26) the template (16) can be expressed as h (t)=a h (t)+a h (t) (30) I 1 1I 2 2I with a =cos(2φ )/D h =h (D =1Gpc; φ =0) 1 c L , 1I I L c (31) a =sin(2φ )/D h =h (D =1Gpc; φ =π/4) 2 c L 2I I L c (cid:26) (cid:26) Using the above expressions and the orthogonality h˜ i h˜ we obtain the maximized log likelihood. 2I 1I ≃ ( <s h >)2+( <s h >)2 I I| 1I I I| 2I (32) F ≃ ( <h h >)2 P I 1I|P1I P 7 2. Maximization over the time of coalescence In order to efficiently find the time of coalescence, we use correlation in place of the inner products. Given a template h which is constructed with the initial value (usually taken at the edge of the prior) t and using the c,0 inverse Fourier transform, we find the value of τ which maximizes (32) or which is equivalent to maximizing max ∞ h(f)s∗(f)+h∗(f)s(f) c(τ)=2 df e−2iπfτ. (33) Sn(f) Z0 e e e e Note that the amplitude of the signal depends on the choice of t via annual modulation caused by LISA’s orbital c motion, therefore the new value t = t +τ is not necessarily the final answer. The time of coalescence which c,1 c,0 max maximizesthequality (32)forgivenotherparametersshouldcorrespondtomaximumof (33)atzero(oralmostzero) lag. Using the new value of t we repeat the maximization, and we stop iterations when the difference t t c c,i c,i−1 | − | is sufficiently small. Usually few iterations are sufficient to find t which maximizes the quality. c 3. The waveform termination The signal from MBH binaries is band limited, the lower frequency limit is defined approximately by twice the orbital phase at t=0. The upper frequency is introduced somewhat arbitrarily. To terminate both the signal and the template smoothly anexponentialtaper isapplied. The taperaffectsthe data whentwoblackholesareseparatedbyadistance R=7M and kills the signal completely around R = 6M (which is the last stable orbit for the test mass in Schwarzschild space-time). Therefore, in computing the overlaps, we use the maximum frequency in the integration corresponding to the orbital separation 6M: 1 η3/5 f = = . (34) max πM(R/M)3/2 π(R/M)3/2M c The exponential taper causes problems for the long-wavelength approximation, and our template deviates from the signal during the last cycle. Unfortunately these small deviations fall in the most sensitive part of the LISA band and are further enhanced by high SNR. This causes a significant problem: the bias caused by this deviation is unacceptably large because there is a large region of the parameter space that produces templates which fit the end part of the signal perfectly (using incorrect parameters) but fail to reproduce the low frequency part of the signal. Inordertosolvethisproblemweterminatethetemplatewaveformfewcyclesearlierbyfixingcutofffrequencywhich corresponds to the orbital separation R > 7M. Our approximation becomes better as we go to lower frequencies, howeverwestartlosingpowerofthesignal(SNR)whichishighlyundesirable. Weautomaticallyreadjustthefrequency cut-off if the SNR drops below a certain threshold. We want to emphasize a very important feature which accompany the earlier termination of the waveform. The mapofthe quality changes: in the Figure 1we showthe mapofthe quality inthe “chirpmass”- “eta”plane keeping other parameters fixed to their true values. On the left panel we show (we use no frequency cut off other that full F introduced by the taper), and, on the right panel, we plot with template cut at f = f = 0.26mHz. One cut max cut F canseemultiplemaximainbothplots,but(!) thepositionofthesecondarymaximaaredifferentwhereasthelocation of the true (global) maximum (indicated by an arrow)is the same. It can also be seen that the size of the secondary maxima on the right panel is smaller. We will use these features later in our search. 4. A-statistic Chopping the template at lower frequency solves the problems mentioned above but is not completely satisfactory. We lose some SNR and consequently some accuracy in the parameter estimation, we also lose information stored at the end of the signal which is especially important to recover spin-related parameters. In order to reduce the impact of the coalescence part, without killing it completely. For that, we introduce a new function, called A-statistic which is simply a geometricalmean of the Maximized Likelihood of the cut waveformand the Maximized Likelihood of the full waveform: = . (35) cut full A F ×F A-statistic is not log likelihood anymore, but one ofpits advantages is that it keeps the information from the full waveformincluding the coalescence but at the same time it enhances the information coming from the low-frequency 8 2.5e+06 100000 1.62e+06 1.62e+06 90000 1.6e+06 2e+06 1.6e+06 80000 70000 1.58e+06 1.58e+06 Chirp Mass 1.56e+06 11e.5+e0+606 Chirp Mass 1.56e+06 456000000000000 1.54e+06 1.54e+06 30000 500000 20000 1.52e+06 1.52e+06 10000 1.5e+06 0 1.5e+06 0 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 Eta Eta FIG.1. DistributionoverMc andη,oftheMaximized Likelihood (quality)computedwiththefullwaveform onleft paneland withthewaveformcutatfmax =0.26 ontherightpanel. Thisexamplecorrespondstoasignalwiththefollowing parameters: β = −0.38896 rad, λ = 3.28992 rad, tc = 19706568.3273 sec, Mc = 1589213.34 M⊙, η = 0.23647, θL = 2.78243 rad, φL =1.53286 rad, χ1 =0.24115, χ2 =0.16145, θS1 =1.20839 rad, φS1 =5.61808 rad,θS2 =0.39487 rad, φS2 =5.82937 rad, DL =6856164697.8 parsec,φc =4.96746 rad . The arrow points to thetrue parameters. part. A-statistic also reduces the number of local maxima as can be seen in the Figure 2. In this example we have reduced the size and number of maxima from five to three. 500000 1.62e+06 450000 1.6e+06 400000 350000 1.58e+06 Chirp Mass 1.56e+06 223050000000000000 1.54e+06 150000 100000 1.52e+06 50000 1.5e+06 0 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 0.25 Eta FIG. 2. Distribution of A-statistic over Mc and η. This example corresponds to the same signal as in Figure 1. The arrow points to thetruelocation of parameters of thesignal. III. GENETIC ALGORITHM A. The basic principle In order to find all the parameters of the signal, we need an effective algorithm to search over the 13 dimensional parameterspace. Buildingthegridinthemulti-dimensionalparameterspaceisahighlynon-trivialproblem. Theuse of the stochastic/randombank [46–49] is a feasible method for the template placement, however a full grid scan over the whole parameter space would be prohibitively computationally expensive. Alternative wouldbe to use variations ofthe MarkovchainMonte-Carlo[29]ornestedsampling[30]methods. Here wehavechosento usegeneticalgorithm (GA) (adjusted to our needs) to search for the global maximum of the likelihood in multi-dimensional parameter space. The GA is derived from the computer simulations of the biological system, which were originally introduced by ProfessorHollandandhisstudentsinMichiganUniversity. Itisamethodfortheglobalsearch(optimizationmethod) based on the natural selection principle – the basis for the evolution theory established by C. Darwin. In the nature, organisms adapt themselves to their environment: the smartest/strongest/healthiest organisms are more likely to survive and participate in the breeding to produce the offsprings. These two processes, selection and breeding, are used in genetic algorithms to produce a new generation of organisms. Since the best organisms are more likely to participate in breeding, the new generation should be better than the previous one (at least no worse). So this procedure induces the evolution of the organism, just like in the nature, the good qualities of the parents can be transferred to their offsprings. In the biological world, besides these two basic operations, among every generation, 9 there are always few individuals which have better characteristics to adopt to the environment, produced as a result of a positive mutation. By introducing the new genotype into the population, mutation can potentially improve the forthcoming generations and consequently accelerate the evolution towards the global maximum. Some measure of “goodness” needs to be associated with each organism. In the case of gravitational wave search, it is natural to associate the logarithm of the likelihood (or any other equivalent detection statistic e.g. Maximized LikelihoodorA-statisticdiscussedintheprevioussection)withthe“goodness”whichneedstobe“improved”through the evolution of the organisms. We will call the chosen measure of “goodnees”, the quality of an organism Q. Following is a brief description of how a typical GA works. We start with a randomly chosen group of organisms (templates), we evaluate the quality of each organism (log likelihood). We select set of pairs (parents) based on their quality, the organisms with better quality (templates with higher likelihood) are chosen more often than weak organisms. We combine genotype of two parents to produce a child (we combine parametersof two chosentemplates to produce a new one). Number of produced children is equal to the number of parents (we keep number of evolving organisms (generation) fixed). Next we allow with a certain probability a random mutation in the children’s genes (with some probability we randomly change the parameters of the new templates, exploring a larger area of the parameter space). The parents are discarded and the resulting children form a new generation. We repeat the procedure until we reachsteady state (maximum in the quality). In this simple example we keep only one generation active (one group of templates). A list of (biological) GA terms with the equivalent terms in GW data analysis is given in the Table I. Genetic algorithm GW search organism ⇐⇒ template gene (of an organism) ⇐⇒ parameter (of a template) allele (of a gene) ⇐⇒ bits(of thevalueof theparameter) quality Q ⇐⇒ Maximized Likelihood or A-statistic colony of organisms ⇐⇒ evolving group of templates n-th generation ⇐⇒ thestate of colony at n-th step of evolution (selection + breeding) + mutation ⇐⇒ way of exploring theparameter space TABLE I.Relation between GA and GW notions. In the following subsections we give a detailed description of each element of the basic GA and then we introduce the specific modifications to speed up the search. B. Code of the gene As we have discussedabove,everyorganismis associatedwith a template and the parametersof the template play the role of genes. So each organism is described by 15 genes, two of them are chosen optimally (maximization of the log likelihood, see Sec. IIB1) and the gene corresponding to the time of coalescence is efficiently found using correlation (see Sec. IIB2). We imitate the DNA structure by describing the gene (parameter value) by a set of alleles. Inourimplementationweadoptabinaryrepresentationofthe gene (parameter)whichmeansthateachallele (bit) has two possible values: 0 or 1. In practice we first fix the precision of each parameter (by fixing the number of significant digits in the decimal format) and then we translate it to standard binary and/or in Gray form. In our method we use both representations,the reason will be explained later when we discuss quantization issue. Let us show how this is done in practice. Consider a parameter θ with the uniform prior range [x ,x ]. k k,min k,max First we convert a value xk of θk into an integer ck =(xk xk,min)/∆xk where ∆xk =(xk,max xk,min)/2Nk is the − − resolution of θ and N is the number of bits. Then, we convert c into the set of bits b [i] using the coding rule k k k k of the chosen representation. As we see, the resolution for each parameter depends on the number of bits N used k for describing it and is the same for both representations. The importance of the bit is determined by its position. A change of a bit in a higher position (significant bit) corresponds to a big change in the parameter value. In our convention, the first bit, b [0], is the lowest significant bit and the last bit, b [N 1], is the highest significant bit. k k − Thereisacloserelationshipbetweentwogenerepresentations. WecantransformthebinaryrepresentationtoGray representation by the following procedure: given a string of binary code with N bits B[0],B[1], ,B[N 1] , we { ··· − } set B[N]=0 then the Gray code with the same N bits is G[i]=B[i+1] B[i], (36) ∧ where the operator “ ” corresponds to the XOR operator in programming languages. Other way round, by setting ∧ 10 again B[N]=0, we can get the binary representation with N bits from the Gray representation as B[i]=B[i+1] G[i]. (37) ∧ C. Selection The selection process chooses the parents for breeding. The probability of selecting an organism is defined by its quality. Organisms with higher quality have a better chance of being chosen to participate in the breeding. First the quality Q , i.e. the Maximized Likelihood or A-statistic, for all organisms is computed (index i refers to the i-th i organism). Then each organism is assigned the probability of being chosen for breeding as p = Q / NQ . The i i j j selection is made by the roulette selection method: we choose a random number uniformly from [0,1]; if it is bigger P thanp andsmallerthanp ,thentheith organismisselected. Thisselectionensuresthatthe“good”organismsare i i+1 chosen more often than the “bad” ones and guarantees that the genotype responsible for a high quality propagates in generations approaching the optimal value. In our approach we do not take into account the geographicalproximity between parents (in other words possible correlation between templates in the same generation). By forbidding the breeding between the correlated parents, it might be more efficient to explore a larger region of parameter space, but the overallresolution of the method will be reduced. We therefore do the selection based only on the quality. D. Breeding After selecting the parents, we need to produce the new generation, this can be achieved through “breeding”. Breeding is the rule according to which a child is produced from the selected parents. The genes of the child are constructed by mixing the corresponding genes of each parent. We take one part from the first parent and the other part from the second one. Depending on which parts are chosen, there are several types of breeding. We usually use three different types: cross-over one random point, cross-over two fixed points, and random. For the cross-over one point, we choose one bit (denoted by i) randomly as the cross-over point and the child’s genes are created by combiningthefirstibitsofthegenesfirstparentwiththelastN ibitsofthegenesofthesecondparent(seetheleft − panel of the Figure 3). For the cross-overtwo fixed points, the genes of the child are built from three equal parental parts (see the middle panel of the Figure 3). In the random breeding, each child’s bit is chosen randomly from the corresponding bits of the parents (see the right panel of the Figure 3). parent 1 0 1 1 1 0 0 1 0 1 0 0 1 parent 1 0 1 1 1 0 0 1 0 1 0 0 1 parent 1 0 1 1 1 0 0 1 0 1 0 0 1 parent 2 1 0 0 0 0 1 1 1 0 0 1 0 parent 2 1 0 0 0 0 1 1 1 0 0 1 0 parent 2 1 0 0 0 0 1 1 1 0 0 1 0 child 0 1 1 1 0 0 1 1 0 0 1 0 child 0 1 1 1 0 1 1 1 1 0 0 1 child 0 1 0 0 0 0 1 1 1 0 0 0 FIG. 3. Examples of used breeding: cross-over one random point on the left, cross-over two fixed points in the middle and random on the right panels. E. Mutation The first generation is chosen randomly by drawing parameters uniformly within the priors specified in [33]. The chosen selection implies that the quality of our organisms is likely to be increased with each generation. But, if we use only these two processes,the range of resulting genes is quite restricted: it totally depends on the initial random state and is just a combination of the parts from the first generation. The combination of genes and therefore the exploration of the parameter space is very limited and completely dependent on the initial choice. This undesirable feature can be cured by introducing mutation. MutationinGAworksinawaysimilartohowitoperatesinthenature. Mutationisarandomchangeoffewalleles inageneofanorganism;inouralgorithmitcorrespondstochangingfewbitsinarepresentationofaparametervalue of a template. The probability of mutation is called the probability mutation rate (PMR). We mutate each gene of each child independently and there are several types of mutation. First we need to decide whether we mutate a gene or not, and, if yes, we need to decide on the mutation rule (how we do it). The first possibility is that we always