UNIVERSITÀ DI URBINODIPARTIMENTO DI SCIENZE PURE E APPLICATE

Corso di Dottorato in Scienza della Complessità

Detectors characterization and lowlatency search of gravitational waves

from binary neutron stars

Francesco Piergiovanni

Tutore: Dott. Gianluca Maria Guidi

Settore scientifico disciplinare FIS/01 FISICA SPERIMENTALE

ANNO ACCADEMICO 2014 / 2015

A Benedetta, Chiara e Matteo

Acknowledgements

Ringrazio per primo Gianluca Guidi per il paziente e puntuale sostegno fornito come tutore,per la indispensabile collaborazione prestata come collega e ancora di più per la franchezzadimostrata come amico.

Un grazie al Prof. Flavio Vetrano per avermi consentito di svolgere da ormai più di diecianni un lavoro davvero stimolante in una grande collaborazione internazionale.

Grazie al gruppo Virgo di Annecy LAPP per l’enorme lavoro di sviluppo della pipelineMBTA e in particolare a Frédérique Marion e Benoît Mours per avermi aiutato ad orientarmitra i meandri del codice.

Un rigraziamento a tutti i colleghi della Sezione di Fisica del DiSPeA dell’Università diUrbino, in particolare a Matteo Montani per le illuminanti discussioni sui punti più osticidi questa ricerca ed anche per la indispensabile quotidiana collaborazione nel lavoro noncompreso in questa tesi.

Ai miei genitori, per avermi sostenuto anche in questa mia seconda vita da studente, ungrande abbraccio e grazie di cuore.

Infine a Benedetta, per essere stata come sempre la compagnia più preziosa, va il ringrazi-amento più grande.

Abstract

One hundred years ago, in 1916, Albert Einstein published the first comprehensive overviewof the General Relativity theory. The theory predicted the existence of gravitational waves:ripples of the space-time traveling at the speed of light, produced by accelerating masses innon-axisymmetric motion. Due to the intrinsic weakness of the gravitational field, experi-mentally detectable gravitational waves can only be generated by astrophysical objects or canhave some cosmological origin and even in these cases the expected amplitude is extremelyfaint. A gravitational wave can be thought as an oscillation of the space-time producing adeformation of the distance among objects. The strain (∆d/d)produced by a gravitationalwave which arrives on the earth is predicted to be around ∼ 10−20. This means that, forexample, the distance between the Earth and the Sun would be subjected to a variation ofthe order of a picometer. This is the reason why gravitational waves have not been observedup to now. Nevertheless, their existence has been proved by Hulse and Taylor observing theorbital decay of the PSR 1913+16 pulsar.

Presently we have the instruments to observe gravitational waves with a significantexpected rate: the advanced interferometric detectors. The first observational run of advancedLIGO has recently starts, with some hope of a first detection, while advanced Virgo isexpected to be on-line for the end of this year. To put this ultra-sophisticated instrumentsin operation, improving their performances and increasing the confidence in a possibledetection, an accurate characterization of the detector have to be performed and a great effortin unfolding noise sources is needed.

In this thesis it is described a tool for understanding noise sources: SILeNTe (SystemIdentification Linear et Non-linear Technique). SILeNTe is specifically designed to uncoverlinear and non-linear relationships between the gravitational wave channel and the auxiliarychannels which record the detector behavior and monitor the surrounding environment.The tool is essentially a fast system identification algorithm based on forward-regressionorthogonal estimator. It can analyze and rank linear and non-linear couplings between agiven noise structure and hundreds of auxiliary channels. The tool has been implementedand tested on some noise structure that afflicted the initial LIGO and Virgo detectors. Theresults are encouraging, giving hints on non-linear coupling mechanisms.

viii

One of the most promising sources of gravitational wave is the coalescing of binarycompact objects as neutron stars and black holes. The physics of these sources is wellunderstood and the expected waveforms are accurately modeled. Different data analysispipelines, using matched filtering technique, have been implemented to drive out signalscoming from these sources, buried in the detector noise. MBTA (Multi Band TemplateAnalysis) is a pipeline devoted to find this kind of signals with a latency of just some tens ofseconds. MBTA has already run in previous initial detector runs performing matched filteringwith bank of theoretical waveforms which considered only the two masses of the binarysystem objects. Since the spins of the individual objects give a significant contribution to thewaveform models, considering the spin on the template generation improves the matchedfiltering searches.

In this work, some functions of MBTA has been upgraded and new features has beenadded. The most relevant development concerns making the pipeline able to handle spinningtemplates and spinning template banks. Two analysis, on binary neutron star coalescencesimulated signals, have been carried on: one with no-spinning templates and the otherwith templates having spin aligned with the orbital momentum. The performance of thetwo analysis have been extensively compared. This comparison has also represented acomprehensive test in order to check the sanity of the new feature implementations.

This is the outline of the thesis: the basic concepts about general relativity and gravita-tional waves is briefly introduced in chapter 1. In chapter 2, a description of the gravitationalwave detectors with the most important noise sources that afflict them, is given. Chapter 3is devoted to a detailed description of SILeNTe and the results of some preliminary test arepresented. Coalescence of binary objects as source of gravitational waves and the techniqueto analyze these kind of signals are described in chapter 4 and 5. An overview of MBTApipeline is given in chapter 6 with a particular attention to the new introduced features.The results of the comparative test on simulated binary neutron star coalescence signals arereported in chapter 7.

Table of contents

1 The Gravitational waves 11.1 From the General Relativity to the Einstein equations . . . . . . . . . . . . 11.2 Gravitational waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2.1 Gravitational wave interactions . . . . . . . . . . . . . . . . . . . . 41.2.2 Gravitational wave emission . . . . . . . . . . . . . . . . . . . . . 51.2.3 Radiated energy . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

1.3 Sources of gravitational waves . . . . . . . . . . . . . . . . . . . . . . . . 81.3.1 Transient signals . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3.2 Long-lived signals . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.4 The binary pulsar PSR1913+16 . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Gravitational waves interferometric detectors 132.1 Basic concepts of Virgo and LIGO advanced detectors . . . . . . . . . . . 14

2.1.1 Optical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.1.2 Mirrors: the test masses . . . . . . . . . . . . . . . . . . . . . . . 162.1.3 Mirror suspensions . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1.4 Thermal compensation system . . . . . . . . . . . . . . . . . . . . 18

2.2 Sensitivity and noise sources . . . . . . . . . . . . . . . . . . . . . . . . . 192.2.1 Seismic noise and gravity gradients . . . . . . . . . . . . . . . . . 202.2.2 Quantum noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.3 Thermal noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.3 Detector characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3 Non-linear system identification: an application to detector characterization 253.1 The Non Linear Identification problem . . . . . . . . . . . . . . . . . . . . 263.2 Model structure determination with OLS algorithms . . . . . . . . . . . . . 273.3 The fast orthogonal search for system identification . . . . . . . . . . . . . 293.4 Model identification with error filtering . . . . . . . . . . . . . . . . . . . 30

x Table of contents

3.5 Applications of system identification OLS technique to Virgo data . . . . . 313.5.1 Calibration lines modulation noise . . . . . . . . . . . . . . . . . . 323.5.2 The Vela case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.6 Applications of OLS on LIGO data . . . . . . . . . . . . . . . . . . . . . . 363.7 Applications of SFOS technique on Virgo data . . . . . . . . . . . . . . . . 373.8 Relative phase and time lags . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.8.1 Cross correlation normalization . . . . . . . . . . . . . . . . . . . 413.9 Virgo noise modeling: a blind search approach . . . . . . . . . . . . . . . 42

3.9.1 Computation time and parallelization . . . . . . . . . . . . . . . . 443.10 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4 Gravitational waves from compact binary coalescence 474.1 CBC formation scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.2 Coalescence rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.3 The signal from a compact binaries coalescence during the inspiral phase . 52

4.3.1 The radiated power . . . . . . . . . . . . . . . . . . . . . . . . . . 544.3.2 The phase evolution . . . . . . . . . . . . . . . . . . . . . . . . . 544.3.3 The post-newtonian approximation . . . . . . . . . . . . . . . . . 574.3.4 Waveform model . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

4.4 Detector response to the inspiral waveform . . . . . . . . . . . . . . . . . 594.5 Merging phase emission . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

5 Detecting CBC signals 635.1 The matched filtering technique . . . . . . . . . . . . . . . . . . . . . . . 635.2 Generating the template bank . . . . . . . . . . . . . . . . . . . . . . . . . 655.3 The aligned-spin template banks . . . . . . . . . . . . . . . . . . . . . . . 67

6 MBTA: a CBC low latency pipeline 736.1 Filtering data in multiple bands . . . . . . . . . . . . . . . . . . . . . . . . 73

6.1.1 Real and Virtual templates . . . . . . . . . . . . . . . . . . . . . . 746.1.2 Hierarchical search . . . . . . . . . . . . . . . . . . . . . . . . . . 78

6.2 The spin template banks case . . . . . . . . . . . . . . . . . . . . . . . . . 796.2.1 Splitting the large banks . . . . . . . . . . . . . . . . . . . . . . . 81

6.3 The χ2 cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 826.4 The matched filter output shape cut . . . . . . . . . . . . . . . . . . . . . . 836.5 Searching the coincidences . . . . . . . . . . . . . . . . . . . . . . . . . . 856.6 Clustering the triggers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

Table of contents xi

6.7 Data quality: science segments and vetoed time periods . . . . . . . . . . . 906.8 The false alarm rate estimation . . . . . . . . . . . . . . . . . . . . . . . . 92

7 A search for BNS with spin 957.1 The analyzed data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 957.2 The template banks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 967.3 Bank splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1007.4 Running MBTA on the GRID . . . . . . . . . . . . . . . . . . . . . . . . . 1017.5 Filtering step and the single detector triggers . . . . . . . . . . . . . . . . . 1027.6 The coincident search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.7 Analysis results and comparison . . . . . . . . . . . . . . . . . . . . . . . 103

7.7.1 Efficiency comparison . . . . . . . . . . . . . . . . . . . . . . . . 1047.7.2 Accuracy in injection reconstruction . . . . . . . . . . . . . . . . . 1067.7.3 The False Alarm Rate and the significance assessment . . . . . . . 110

7.8 Final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

8 Conclusions 115

References 117

Chapter 1

The Gravitational waves

The General Relativity theory, formulated by Einstein in 1915, is based on the EquivalencePrinciple of Gravitation and Inertia, which introduces a geometric approach to the gravita-tional interactions, and the Principle of General Covariance, which states that all the physicallaws preserve their form under an arbitrary change of coordinate system, in other words theequations must be covariant and expressed in tensor form [1].

1.1 From the General Relativity to the Einstein equations

According to the Equivalence Principle in a free falling frame, the distance between twopoints on the spacetime can be computed as:

ds2 = ηµνdξµdξ

ν (1.1)

where ξ µ are the locally euclidean coordinates and ηµν = diag(−1,1,1,1) is the metrictensor of the locally flat spacetime. On that frame the motion of a free falling particle (i. e.with no other forces than the gravitational field acting on it) is given by:

d2ξ α

dτ2 = 0 (1.2)

with τ the particle proper time, chosen to be the time coordinate of the local frame. Changingthe frame with a coordinate transformation ξ α → xα = xα(ξ α) the distance between twoneighboring points is:

ds2 = gµνdξµdξ

ν (1.3)

2 The Gravitational waves

where gµν is the metric tensor of the new frame:

gµν =∂ξ α∂ξ β

∂xµ∂xνηαβ (1.4)

Defining the Christoffel symbol, also named affine connection, as:

Γαµν =

12

gαδ

(∂gδγ

∂xβ+

∂gδβ

∂xγ−

∂gβγ

∂xδ

)=

∂xα

∂ξ λ

∂ 2ξ λ

∂xµ∂xν(1.5)

in the new frame the equation 1.2 will change in this way:

∂ 2xα

∂τ2 =−Γαµν(

∂xµ

∂τ

∂xν

∂τ) (1.6)

This equation, known as the geodesic equation, describes the motion of a free falling particlein an arbitrary coordinate system. Comparing it with the equation 1.2 it is evident the roleof the affine connection that acts like a generalization of the Newtonian gravitational field,carrying also the effects of the apparent forces in case of a not inertial frame.

The stress-energy tensor (Tµν ) contains the energy density and the flux of energy andmomentum, representing the General Relativity generalization of the gravitational fieldsources. The Einstein field equations establish the relation between the metric tensor andthe affine connection, on one hand, and the metric tensor on the other [2]. Writing the Riccitensor Rαβ and the scalar curvature R as a contraction of the Riemann tensor Rα

βγδ:

Rα

βγδ=

∂Γα

βδ

∂xγ−

∂Γα

βγ

∂xδ−Γ

α

δεΓ

ε

βγ+Γ

αγεΓ

ε

βδ(1.7)

Rαβ = Rγ

αγβ(1.8)

R = Rαα (1.9)

and defining the Einstein tensor Gµν :

Gµν = Rµν −12

gµνR (1.10)

the Einstein field equation can be written in a very compact form:

Gµν =8πGc4 Tµν (1.11)

1.2 Gravitational waves 3

1.2 Gravitational waves

The Einstein field equations are intrinsically non linear. Any radiative solution carriesenergy and momentum which are thereby sources of the gravitational field and appear on theright side of the equations. Since the gravitational radiation have a very low intensity, theWeak Field Approximation can be assumed far enough from the source and the perturbativeapproach can be followed [3]:

gµν = ηµν +hµν |hµν |<< 1 (1.12)

where hµν is a small perturbation to the flat Minkowski metric tensor ηµν . Hence, consideringa point far from sources and in vacuum, where Tµν = 0, and using the approximation up tothe first order on h, the Einstein field equations can be read as:

hµν −∂ 2hλ

ν

∂xλ ∂xµ−

∂ 2hλµ

∂xλ ∂xν+

∂ 2hλ

λ

∂xµ∂xν= 0 (1.13)

The equation 1.13 contains a set of 10 equations for the 10 independent components of thesymmetric tensor hµν . However the symmetric properties of the Riemann tensor, given bythe Bianchi identities, limit the number of the independent equations to six, leaving fourdegrees of freedom. This is what expected according to the Principle of General Covariance,which implies that under a coordinate transform, the covariant transform of the solution mustbe a new solution, representing the same geometry on the new reference frame. Therefore theEinstein equations do not determine a unique solution but leave a gauge freedom. Obviouslythe equation 1.13 remains valid only if the condition |hµν | << 1 is satisfied in the newcoordinate system. This is what happens in the harmonic gauge that is the reference frame inwhich a free falling body is at rest and the following equation is satisfied:

gµνΓ

λµν = 0 (1.14)

that for h becomes:∂hµ

ν

∂xµ=

12

∂hµ

µ

∂xν(1.15)

The harmonic gauge condition does not determine the gauge uniquely and it is alwayspossible to find a reference frame where equation 1.15 is satisfied, the tensor hµν is tracelesshµ

µ = 0 and the terms h0ν = 0. This gauge is called transverse traceless (T T ). In that framethe equation 1.13 assume a very simple form:

hµν = 0 (1.16)

4 The Gravitational waves

Equation 1.16 admits the general waveform solution:

hµν = Aµνeıkλ xλ (1.17)

where kλ is the wave vector satisfying the condition:

kλ kλ = 0 (1.18)

This means that the wave vector is a null vector and consequently the wave propagates at thespeed of light. From equation 1.16 another condition on the wave vector follows:

kµAµ

ν =12

kνAµ

µ (1.19)

from which the tensor hµν for a monochromatic wave, propagating along the z axis, can bewritten as follow:

hµν =

0 0 0 00 h+ h× 00 h× −h+ 00 0 0 0

(1.20)

where h+ and h× are defined as follow:

h+ = A+cos[ω(t − zc)+φ0]

h× = A×cos[ω(t − zc)+φ0]

(1.21)

being A+ and A×, the complex amplitudes of the two components h+ and h×, referred as theplus (+) and the cross (×) components and φ0 a phase constant.

1.2.1 Gravitational wave interactions

To see how a gravitational wave interacts with matter, let us consider two neighboring freefalling particles initially at rest on the plane xy and a gravitational wave moving along the zaxis an impinging orthogonally on the plane. In a T T gauge reference frame the two particlesare assumed to remain at rest also when the wave arrives. So the position of the two particleswill remain constant but, since the metric change, the proper distance between them willchange as well. For example for two particles placed on x axis the proper distance (∆l) will

1.2 Gravitational waves 5

Fig. 1.1 Effect of the two polarizations of a monochromatic gravitational wave on a circularmass distribution.

be:

∆l =∫

ds =∫ x2

x1

|gxx|12 dx =

∫ x2

x1

|1+hxx(t)|12 dx ≈ (1+

12

hxx(t))(x2 − x1)≈

≈ [1+12

A+cos(ω(t +φ0)](x2 − x1)

(1.22)

where x2 and x1 are close enough to assume hxx independent from x. Hence, the properdistance between the two particles oscillates at the wave frequency with an amplitude whichdepends on the plus polarization of it (A+). It is easy to show that for two particles thatlies on the bisector of the plane xy, the amplitude of the proper distance oscillation will bedependent from the cross polarization (A×) and the phase will be rotate by 45. The effect ofthe gravitational wave on a ring of particles is graphically shown on figure 1.1

1.2.2 Gravitational wave emission

The intensity of a gravitational wave far from the source can be computed starting again fromthe Einstein field equation 1.11. In T T gauge and weak field conditions, the equation can berewritten as:

hνµ =−16πGc4 Tµν (1.23)

Considering source mass-energy distribution ρ(x) placed around the origin of the referenceframe and a point in n at a distance D = |n|>> |x| where |x| is the typical dimension of the

6 The Gravitational waves

source, a general solution is:

hµν =4πGc4

∫ Tµν(t − |x−n|c ,x)

|x−n|d3x (1.24)

Far from the sources the term |x−n| can be approximate as |x−n| ≈ D(1+O( |x|D )) and thesolution becomes:

hµν =4πGc4D

∫Tµν(t −

Dc,x)d3x (1.25)

Since the conservation law of the stress energy tensor implies that T µν

;ν = 0, it can be shownthat hµν depends on the second time derivative of the mass quadrupole momentum associatedto the mass-energy distribution of the source:

Mi j =∫

ρ(x)xix jd3x (1.26)

where ρ = T 00

c2 . The equation for h in a T T gauge far from the sources will be:

hµ0 = 0 µ = 0,3

hi j =8πG3c4D

Mi j(t −Dc) i = 1,3, j = 1,3

(1.27)

This means that gravitational wave comes from acceleration of systems with time dependentmass quadrupole moment, while system with spherical or cylindrical symmetry cannot besources of gravitational waves. The tensor Mi j in equation 1.27 is computed in the T T framewhere the wave is propagating along the z axis. The expression for the two wave polarizationson this frame will be:

h+ =4πG3c4D

(Mxx − Myy)

h× =8πG3c4D

(Mxy)

(1.28)

This frame is commonly named "radiation frame". In case of sources with particular geomet-rical characteristics, as the case of two binary objects, is simpler to compute the quadrupolemoment in a different frame where the orbit lies in the x′y′ plane. This frame is generallycalled the "source frame". Since the rotation of the reference frame around the directionof propagation (z axis) is free, two angles are enough to define a rotation from a frame tothe other, as in figure 1.2 [4]. Naming θ and φ the two angles and placing the x axis of the

1.2 Gravitational waves 7

Fig. 1.2 Relationship between source frame and radiation frame [5].

radiation frame on the plane x′− y′, the rotation can be written as:x′

y′

z′

= R

xyz

(1.29)

where R is the product of the two rotations R = Rz(φ)Rx(θ):

Rx(θ) =

1 0 00 cosθ −sinθ

0 sinθ cosθ

(1.30)

Rz(φ) =

cosφ −sinφ 0sinφ cosφ 0

0 0 1

(1.31)

The quadrupole moment will thus transform accordingly:

M′i j = RikMklRl j (1.32)

1.2.3 Radiated energy

The power emitted by the sources as gravitational waves can be computed in a T T frame as[3]:

P =dEdt

=r2c3

32πG

∫ ⟨hi jhi j

⟩dΩ (1.33)

8 The Gravitational waves

where Ω is the solid angle and the bracket indicate the temporal average. Considering thetwo wave polarizations and integrating over Ω, the equation can be rewritten to get:

P =c3

16πG

∫ ⟨h+(t)2 + h×(t)2⟩dΩ (1.34)

1.3 Sources of gravitational waves

The gravitaional wave sources which are interesting for the current or the foreseen gravita-tional waves detectors, are brefly described. The sources are classified on the base of theduration of the emitted signal in the detector frequency bandwidth and on the knowledge ofthe expected signal.

1.3.1 Transient signals

Burst or transient signals are related to sources that release gravitational wave energy in ashort time period of the order of some seconds or minutes.

• Burst sources. Gravitational wave burst refers to waves emitted in a very short time,of the order of seconds or less, with an un-modeled waveform. The most relevant burstsources are the typeII-supernovae violent core collapse. To emit gravitational wavesthe supernovae have to exhibit some asymmetry like non-spherical motion given bydynamic instabilities of the progenitor stars. Another emitting mechanisms can be thehydrodynamical oscillation of the nascent neutron star core, as hypothesized by recentstudies[6]. The gravitational wave strain amplitude expected for the gravitational waveemission can be written as [7]:

h ≈ 6 ·10−21(

E10−7M⊙c2

)1/2(1msT

)(1kHz

f

)(10kpc

r

)(1.35)

where E is the energy emitted as gravitational wave, T and f the characteristic timeand frequency of the emission and r the source distance. If some realistic (but largelyuncertain) values are inserted in equation 1.35, it is found that advanced interferometricdetectors will be able to detect supernovae emission within the Milky Way and thenearby dwarf companions galaxies, with a expected rate of the order of 1/30 yr. Evenif simulations give important results to understand the mechanisms involved in super-novae explosions, the poor knowledge of many of the simulation parameters makesthe prediction on the waveform subject to large uncertainty. The search technique ofburst signals is based on the detection of the coherent excess power on the outputs

1.3 Sources of gravitational waves 9

of different detectors [8]. No model of the emitted signal is generally assumed. Thismeans that signals coming from unknown sources or emission mechanisms can be alsodetected.

• Compact binary coalescence. A binary system composed by two compact objects asneutron stars or black holes, rotating around their center of mass, loses its energy asgravitational waves, shrinking the orbital radius. If the two objects are tight enough tocoalesce in a Hubble time, GW emission can fall inside the interferometer bandwidth.This kind of sources can give a detectable transient signal during the last part (minutesor seconds) of the inspiral orbital motion, the merger of the two objects and the possibleoscillations of the remnant body. The emission mechanisms are very well definedand accurate model of the expected waveform are available. Thus, matched filteringtechnique can be used in powerful searches for this kind of sources, with suitabletemplate banks, covering the parameter space for the different kind of orbiting objects(black holes or neutron stars). The rates of compact binary coalescence which isexpected to be detected by advanced interferometers, are summarized on table 1.1.

Since the searches of the gravitational waves from compact binary objects is one of themain arguments of this thesis, it will be discussed in more details in later chapters.

Source NLow [yr−1] NReal [yr−1] NHigh [yr−1]

NS-NS 0.4 40 400NS-BH 0.2 10 300BH-BH 0.4 20 1000

Table 1.1 Expected detection rates for coalescence of binary neutron stars (NS), black holes(BH) and NS-BH systems. Higher, lower and realistic detection rate estimations are reportedconsidering the advanced interferometer sensitivity curve[9].

1.3.2 Long-lived signals

Continuous gravitational waves refer to signals that last for long period within the detectorbandwidth.

• Periodic sources. Gravitational wave signals at well defined frequency, can be emittedby rapidly rotating no-axisymmetric systems [7] [10]. Typically a rotating neutron starwhich presents some spherical asymmetric deformation, can be a source of periodic

10 The Gravitational waves

signals. A precessing motions of its rotation axis, can be also a source of monochro-matic waves as well as the normal mode oscillations of the stars. Rapidly rotatingneutron stars have been observed by their radio emission (radio-pulsar). For theseobjects spindown rate produced by energy emission can be observed. Most of theradio-pulsar energy is emitted as a wind of energetic particles but a significant amountof energy can be possibly emitted as gravitational waves. Equation 1.35 and spindownmeasurements, can give an upper limit to the gravitational wave amplitude supposingthat all energy is emitted as gravitational waves. The spindown limit has been beatenby initial detector searches, for two known radio-pulsars: Vela and Crab. The lackof gravitational waves in their data have set the upper limit of the fraction of energyemitted as GW to ≈ 1% for Crab, and ≈ 10% for Vela pulsar[11].

• Stochastic background. The superposition of incoherent emissions, far unresolvedtransient and periodic sources can give a floor of confused signals which contribute to agravitational wave background [12]. Another stochastic background source consists onwaves generated at the very beginning of the life of the Universe, or in other transitionphases during the Universe evolution [13] [14]. This radiation can be really primordialbecause, while the cosmic microwave radiation come to us from ∼ 105 yr after theBig Bang, the present-day GW background can be emitted in the first 10−35 s afterthe Big Bang. It is very arduous to distinguish this signal from the detector noise [15][7]. Using a single interferometer, a very accurate detector characterization, whichunfolds very precisely every noise source and coupling, could permit to discriminatethe GW background. A correlation between two detectors could be useful to canceluncorrelated noise, but with the condition that the two detectors must have the identicalresponse to the signal. This is not the case for terrestrial advanced GW detectors thatare separated by significant distances and with different orientations, whereby thesignals will be uncorrelated. The 3rd generation of GW detectors will have co-locatedmultiple interferometers, increasing the possibility to detect background signal.

1.4 The binary pulsar PSR1913+16

The PSR 1913+16, or Hulse Taylor binary system [16], is a binary neutron star discoveredby Joseph Taylor and Russell Hulse in 1974, using the Arecibo antenna. The system iscomposed by a pulsar and a companion neutron star rotating around their center of mass.During this motion the system loose rotational energy emitting gravitational waves [17].Hence, the orbital radius will decrease over the time and, according to the third Kepler law,also the rotating period will decrease.

1.4 The binary pulsar PSR1913+16 11

From 1975 to the present, the rotating period has been measured and compared with thetheoretical expectation. The results shows a extremely good agreement between measure-ments and expectation, with a ratio between the predicted and observed energy loss quotedas 0.997±0.02 [18]. The figure 1.3 shows graphically the measured cumulative decreasingof the binary period and the General Relativity prediction.

12 The Gravitational waves

Fig. 1.3 Orbital period decay of the binary system PSR1913+16: dots are measurements andthe line is the General Relativity prediction[18].

Chapter 2

Gravitational waves interferometricdetectors

The first experiments to detect gravitational waves go back to the 1960s, with the buildingof a network of resonant mass detectors working at cryogenic temperatures [19]. Ten yearslater, the next step, was the proposal of detectors based on L-shape Michelson interferometrictechniques using arms of the order of kilometers long with an strongly improved sensitivity.Three large interferometric detectors have been build up and put in working condition in therecent years: the two detectors of the Laser Interferometer Gravitational wave Observatory(LIGO) and the detector Virgo [20] [21]. These instruments have constituted a networkthat have taken and shared scientific data from 2005 to 2011. Other less sensitive or withnarrower band has been realized: GEO600 [22] and TAMA300 [23] that have been usedto search for specific sources and to be a test bench to develop new detection techniques.During the last five years, major upgrading of the initial LIGO and Virgo interferometersmarks the beginning of the Advanced Detector Era (ADE) with a second generation ofinstruments: the aLIGO and the AdVirgo interferometers[24] [25]. The expected sensitivityof the advanced detectors is roughly ten times better than the initial interferometers. Whilethe end of AdVirgo construction, in an early configuration, is foreseen for the very beginningof 2016 to be ready to take science data for the fall 2016, the first observational run (O1)of the early configuration aLIGO has covered the period from September 2015 to January2016. With those instruments in their final configuration, the expected rate based on themore realistic models will be of the order of tens per year. Another detector, KAGRA[26],will join to the network in the next few years with extremely innovative techniques. It isan underground and cryogenic detector and it will represent a transitional model for a thirdgeneration of instruments. In this chapter the basic concepts of the GW interferometricdetectors will be described, with particular attention to the improvements between the initial

14 Gravitational waves interferometric detectors

Fig. 2.1 Advanced LIGO optical configuration. ITM and ETM are the test masses at theinput/end of the FP cavities; the power recycling mirror (PRM), the beam splitter (BS), thesignal recycling mirror (SRM) are also shown. PD is the output photodetector. The laserpower numbers referring to the ITF in final configuration [24].

and the advanced version of LIGO and Virgo. The fundamental noise contributions thatdetermine the designed sensitivity curve will be described, as well as the technical noisewhich affects the real sensitivity curve.

2.1 Basic concepts of Virgo and LIGO advanced detectors

The conceptual design of Virgo and LIGO interferometers (ITFs) is very similar and it issketched in figure 2.1. The two detectors are basically two Michelson interferometers whereeach arm is a Fabry-Perot resonant cavity, which convert an arm length change on a laserlight phase shift, producing a signal on the output port. The dimensionless amplitude of thespace-time perturbation h produced by a gravitational wave is extremely faint, of the order of10−21 −10−22 for solar mass sources at few tens of Mpc. Since the length change producedby an impinging gravitational wave is proportional to the arm length itself, ∆L = h(t)L/2,

2.1 Basic concepts of Virgo and LIGO advanced detectors 15

the gravitational interferometers need to be very big instruments [27]. The Virgo arms size3 km, while the LIGO ones are even longer: 4 km. The advanced interferometers have adual-recycled configuration: beside the power recycling, already implemented in the initialITF, a signal recycling is inserted at the detector output. To beat the noises sources manychallenging issues have to be faced. The long arms must be kept under ultra-high vacuum toprevent effects due to fluctuations of the refractive index along the path. To remove the effectsof the seismic motion, the mirrors must be isolated from the ground by very sophisticatedsuspension systems. The thermal motion produces a noise proportional to the temperatureand the energy dissipation. Working at room temperature, the minimization of the mechanicallosses has been a major goal of the the detector systems design. The uncertainty related tothe photon counting at the output and the fluctuation in radiation pressure, determines theso-called quantum noise, that is faced increasing the effective power inside the interferometer.

2.1.1 Optical scheme

The basic optical scheme is a Michelson interferometer, enlightened by a 1064 nm laserlight, with orthogonal arms and mirrors which are suspended and isolated from the groundwith a suspension resonance frequency below the GW detection band. In this way, mirrorscan be considered like free falling test masses. When a gravitational wave approaches theinstrument, the travel time of the light on the two arms changes. For frequency much smallerthan 1/τ0, where τ0 = 2L/c, the difference is approximately equal to δτ ≈ 2L

c h(t). Therelative dephasing is δφ ≈ 4πL

λh(t) with λ the laser wavelength. To maximize the output

contrast the arms length are locked to keep the dark fringe in absence of GW signal [28].Fabry-Perot cavities are realized placing two semitransparent mirrors at the beginning of thetwo arms to enhance the detection performances. The photons that enters in the cavity arereflected back and forth storing light power on the arms that results in a detection performanceimprovement. The cavity finesse F depending on the mirror optical characteristics givesthe average number of the back and forth travel done by the photons: N = 2F/π . For theadvanced detector a finesse of about 450 has been chosen increasing of about a factor 9 theinitial detector value. The increment on the dephasing produced by a GW is proportional toF : δφ ≈ 4πL

λ

2Fπ

h(t).To increase the power circulating inside the interferometer a power recycle system is

implemented. When the output is locked on the dark fringe all of the power injected insidethe interferometer by the laser is reflected back, this means that the optical impedance ofthe ITF does not match the laser beam. Inserting a power recycling mirror it is possible tochange the optical impedance chosen correctly the transmission properties. The reflectedback light can be minimized building a new cavity which involved the ITF and the power

16 Gravitational waves interferometric detectors

recycling mirror, allowing the instrument to store a huge amount of power. For advancedLIGO and Virgo detector in their final configuration the power will be of the order of 100 Wat the input, 5.2 kW on the beam splitter and up to 750 kW inside the arm cavities. Signalrecycling cavity [29] has been added in the advanced detectors, this adds degrees of freedomto the optical configuration making the instruments more flexible in tuning the shape of thesensitivity curve. The requirements on the arm cavity finesse are also relaxed and the samelevel of sensitivity with a larger bandwidth can be obtained.

The readout of the interferometer is based on hom*odyne technique (DC readout). Thistechnique consists in introducing a carrier signal offsetting the cavities by some picometers.Then all the light that cannot contain gravitational wave signal is removed by an output modecleaner cavity and the demodulation is provided by the interference between the carrier-modeand the GW-induced sidebands[30].

2.1.2 Mirrors: the test masses

The design of the mirror substrates and coatings is a crucial point for the instrument per-formance, hence all the optical and mechanical characteristic must meet very stringentspecification. The mirror substrates material (high purity synthetic fused silica Suprasil3001/3002) is chosen to have ultra low absorption to minimize the thermal effects, very highquality factor to reduce the thermal noise and a low level of inhom*ogeneity. To deal withradiation pressure and thermal fluctuations, the mirror must be large (∼ 30 cm in diameter)and heavy (∼ 40 kg, two time heavier than the initial detector test masses). Super-polishingtechniques and ion-beam milling provide a flatness better than 0.5 nm rms in the central area.The most challenging issue was represented by the mirror coating that determines both thetotal mirror thermal noise and the optical quality. The optimized coating used for advancedLIGO and Virgo [31][32] represents the state of the art in this field, with extremely smallabsorption losses (of the order of 0.4 ppm at 1064 nm) and low dissipation (and consequentlylow thermal noise) achieved by the use of Ti doped Ta2O5 and SiO2 layers. To optimizethe cavity symmetry the input mirrors have been coated together, as well as the end mirrorsand to not spoil the flatness of the substrate a new planetary motion ion-beam sputteringtechnique has been used.

2.1.3 Mirror suspensions

To isolate the detector test masses from the ground motion a large attenuation, of the order of1010 above 10 Hz, must be provided by the mirror suspension system. In the last stage ofLIGO and Virgo advanced detectors suspension, the main mirrors (those at the input and at

2.1 Basic concepts of Virgo and LIGO advanced detectors 17

Fig. 2.2 LIGO suspension system. The first scheme on the left shows the three steel filterstages and the last stage composed by the penultimate silica mass and the test mass. Thescheme on the middle the last stage is visible in more detail. The test mass mirror is hung byfour silica fibers to an identical silica mass. On the right some details of the mirror-fiberscoupling is visible. The fibers are welded to silica pieces called ears which are hydroxide-catalysis (silicate) bonded to the silica masses [24].

the end of the FP cavities), are hold by means of four fused silica fibers. This solution wasalready tested in GEO600 [33] and in the first major upgrade of Virgo named Virgo+[34].This monolithic configuration, with fused silica test masses suspended by fused silica fibers,permits to achieve a very low level of dissipation. The fibers are 400 µm thick to providethe proper violin and bouncing mode frequencies, but they have two thicker end parts, in theregion where the bonding energy is mostly stored, to lowers the thermoelastic dissipations[35].The contribution of the suspension thermal noise to the sensitivity is lowered by more thanan order of magnitude with respect the steel wire suspensions of the initial detectors [36] .The upper part of the isolation system is very different for LIGO and Virgo, even though thetwo isolation strategies can be both considered an hybrid mix of passive and active damping.LIGO uses mainly an active damping strategy with quadrupole pendula, supported by anactive isolation platform [37] [38] (figure 2.2). The Virgo Super Attenuator approach is more

18 Gravitational waves interferometric detectors

Fig. 2.3 Virgo suspension system. On the left the Super Attenuator with the inverted pendulum(IP) legs and the four suspended steel filters. On the middle a magnification of the last stagewith the test mass suspended to a steel plate called "marionetta". All around the mirrorother many parts of the thermal compensation system and baffle to prevent scattered lightare visible. On the left a naked view of the last stage with the fibers welded to silica piecessilicate bonded to the "marionetta" and the mirror [25].

oriented to a passive damping (figure 2.3): a chain of four mechanical seismic filters arehung to an inverted pendulum consisting in three 6 m long aluminum legs connected to theground through flexible joints. In both LIGO and Virgo, vertical attenuation is provided bysteel blade spring mounted on the suspension filters to achieve a vertical resonance under 1Hz at the test mass.

2.1.4 Thermal compensation system

Although the mirrors used by advanced detectors have extremely low absorption, the powercirculating in the cavities is so high that can generate thermal deformation of the mirrorsurfaces. This can produce aberrations spoiling the detector performances. A very effectivethermal compensation system (TCS) permits to remove the thermal effects and also to cure

2.2 Sensitivity and noise sources 19

some residual inhom*ogeneities on the mirror coating and substrate. The system is also usedto tune the radius of curvature of the cavity mirrors [39].

The TCS is composed by a sensing system and actuators. The sensing is provided bya phase camera that measures the intensity and phase distribution of the fields inside therecycling cavity and an Hartmann Wavefront Sensor that looks to the thermal aberration ofa probe beam. The actuators are: a resistive ring-heater that acts directly on the mirror tomodify the radius of curvature and a CO2 laser that shines a compensation plate (CP), that isa trasmissive optic facing the mirror. The presence of a CP is needed to isolate the mirrorfrom the CO2 laser noise. The heat is applied to the CP with optimized pattern in order tominimize aberrations due to thermal lensing and other possible defects.

2.2 Sensitivity and noise sources

The advanced detectors have been designed to improve the sensitivity with respect to theinitial ones, in each part of the bandwidth. Over a broad band (the most important for manyof the relevant GW sources) the improvement is of the order of a factor ten. To reach this goalevery detector subsystem has been upgraded to break down all the possible noise sources.The design sensitivity curve is limited by the noise floor determined by the fundamentalnoises: seismic noise, Newtonian noise, quantum noise and thermal noise. Other noisesare classified like "technical", such as: laser fluctuations in frequency, amplitude and phase,interferometer control noise, photodetectors noise, scattered light etc.

To evaluate the performance of the detector, beside the sensitivity curve, the standardfigure of merit is considered to be the BNS range. Considering a binary neutron star systemcomposed by two objects of masses 1.4 M⊙ each, the BNS range is defined as the distanceto which the interferometer detects the merger of this system with a signal to noise ratioequal to 8, averaging over all sky location and orbital plane inclination [40]. The range ofadvanced detectors are one order of magnitude larger than the initial ones; this means avolume of Universe 103 times bigger is reachable. The figure 2.4a shows the aLIGO (finalconfiguration) sensitivity curve and the contribution of the main noise sources. The BNSrange for that curve is 190 Mpc. The sensitivity curve and the main noises of AdVirgo areshown in figure 2.4b. The designed BNS range for AdVirgo in final configuration is 140Mpc.

20 Gravitational waves interferometric detectors

(a) (b)

Fig. 2.4 Designed sensitivity curve and principal noise contribution of aLIGO (a) and AdVirgo(b). The plots refer to the detectors in their final configurations [24] [25].

2.2.1 Seismic noise and gravity gradients

Ground motion can be produced by micro-seismicity and by anthropic activities. The spectralamplitude density have the following empirical law:

x( f )≈ Af 2 (2.1)

with A of the order of 10−7 mHz3/2.Mirror suspension systems provides to a huge attenuation making the seismic noise

negligible above ∼ 10 Hz, while below it represents a severe limit for the detection.The ground motion transmitted to the mirrors by the suspension is not the only effect

related to the seismic activity. The seismic waves perturb the distribution of matter, andconsequently the gravity field that acts directly on the mirrors[41]. This produces theso-called gravity-gradient (or Newtonian) noise which sets the limit of the sensitivity atlow frequency. It is not possible to shield or attenuate this kind of noise for an advancedinterferometric detector, at least without building it underground as for KAGRA detector.However, recent accurate models of terrestrial gravity perturbations show the possibility ofan effective active mitigation[42].

2.2.2 Quantum noise

The quantum noise is the sum of two effects related to the corpuscular nature of the light: theshot noise and the radiation pressure fluctuation.

2.2 Sensitivity and noise sources 21

Shot noise is generated by Poisson statistical fluctuation of the number of photons in theoutput signal [29]. If the light modulation is taken into account [43], the contribution to thenoise strain is:

hshot( f ) =

√32

18LF

√2hP

λcηCPin

√1+(

ffFP

)2

(2.2)

where L is the arm length, C the recycling factor, F the FP cavities finesse, η the photodetectorefficiency , Pin the injected light power and fFP the cavity characteristic frequency fFP = c

4LF .

The radiation pressure is the fluctuation of the force that photons apply to the mirrors. IfM is the mass of the mirror, the noise contribution is given as follows [27]:

hRP =F

ML f 2

√hPCPin

2π4cλ(2.3)

Hence the radiation pressure noise increase as the square root of injected power Pin while theopposite happens for the shot noise.

2.2.3 Thermal noise

Thermal noise is the main noise source for a large part of the detector bandwidth, fromsome tens to few hundreds of Hz. Thermodynamical fluctuation around the equilibriumis a general topic related to every degrees of freedom of a system at a given temperature.For GW detectors the relevant thermal effects are associated to the mirror and the laststage of suspension degrees of freedom. From the Fluctuation-Dissipation (F-D) theorem[44] the thermal fluctuations at the equilibrium are completely determined by the energydissipation of a mechanical system. The F-D theorem states that the power spectral densityof the fluctuation is proportional to the temperature T and the real part of the inverse of themechanical impedance Z(ω):

⟨x(ω)⟩2 =4kbTω2 ℜ

[1

Z(ω

](2.4)

where ℜ[Z(ω)−1] is proportional to the mechanical damping. Thus all the materials used in

mirror substrates and coatings, as well as the mirror suspensions are chosen to minimize theenergy losses. Moreover all the coupling between the different parts of the suspension aredesigned to lower the dumping mechanisms as much as possible.

The total thermal noise is given by the sum of the suspension last stage and the mirrorcoating losses. The suspension contribution, in the region of tens of Hz, is the sum of intrinsicdissipation of the fused silica fibers (mostly due to the surface losses [45]), thermoelastic

22 Gravitational waves interferometric detectors

losses (due to thermal cycles induced by the bending [46] [35]) and clamping losses, due tothe coupling between the fibers and the other components of the suspension [47].

The mirror coating determines the thermal noise contribution above few tens of Hz.Different mechanisms play a role in coating thermal noise [48] [49]:

• intrinsic dissipation of the SiO2 −Ta2O5 layers, mostly due to the Ta2O5 which havethe poorest mechanical characteristics;

• thermoelastic dissipation due to heat flux produced by elastic vibrations and by laserphotons absorption;

• thermorefractive fluctuation produced by the coupling between the temperature fluctu-ations and the variation of the refraction index.

2.3 Detector characterization

In figure 2.4a and 2.4b the design sensitivity of the detectors are shown. However, thesensitivity does not tell the full story about the detector performances. The output signal canbe disturbed by a large variety of noise sources originating inside the detector, or coming fromthe environment around it. These disturbances produce effects that restrict the gravitationalwave sensitivity. This effects can be persistent, like narrow bandwidth spectral bumps orlines, or they can have short duration (s or ms), like glitches.

Searches for transient gravitational wave sources can be affected and spoiled by glitches.As will be discussed in section 6.7, when this kind of glitches appears, the data are flaggedand vetoed, i.e. short period of time surrounding these events are excluded from the analysis.The result is a reduction of the effective duty cycle.

The persistent disturbances give elevated noise at given frequencies. The detector ispartially or even completely blinded to that frequencies, affecting the searches of continuousgravitational waves for emissions in frequency regions close to that disturbances and thestochastic background search can be also spoiled.

During the last run of the initial interferometers, many source of disturbances have beenwidely studied both in LIGO [50] and Virgo [51] detectors. In the following a list of sourcesthat are coupled with disturbances, are reported.

• Glitch sources.

– Seismic glitches. Low frequency ground motion is effectively attenuated by themirror suspension systems, nevertheless the seismic noise can be re-injected and

2.3 Detector characterization 23

up-converted on the interferometer by means of several mechanisms [52] [53].These include light that is scattered by tower wall and suspended baffles movingwith the micro-seismic motion of the ground. Bad weather conditions, strongwind against the detector building, storms and also heavy human activities canincrease the rate and the amplitude of the seismic glitches.

– Acoustic glitches. The detector systems can be mechanically affected by acousticwaves produced around the detector. For example, airplanes and helicopter flyingnear the detector site give transient signals on the detector output. Most of thetime the source are mechanical devices operating too close to the interferometer.Acoustic enclosures have been realized to isolate the detector optical benches toattenuate these effects.

– Electrical cables and connections. Especially air condition, heating systemsand electronic cooling fans can produce big electric transients when they areswitched on and off. Thus all the detector subsystems have to be isolated andelectrically decoupled from the environment. However the detectors are affectedby the electric main frequency (50Hz for Virgo in Europe and 60Hz for LIGOin the USA) that can couple with magnetic glitches. Another source of glitchescan be the connection between electronics devices like electronic board solderingand faulty connectors.

– Radio frequency glitches. High-frequency electromagnetic glitches close to thelaser modulation frequency can be piked up and then demodulated. The sourcesof RF noise are typically electronic devices, or high order harmonics of the powersupply switching rate or digital devices.

• Persistence sources:

– Vibration noise. Vibration is re-injected in the interferometer mainly by means ofscattered light which transmits to the detector the motion of pipes, tower walls,optical windows etc. If the motions are periodic, lines appear on the output powerspectrum. The scattered light, together with other mechanisms, can transmit themotion of the electric and electronic devices, e.g. cooling fans and pumps. Inthe initial detectors the main laser beam jitter was coupled with vibration of theoptical benches driven by the cooling fans. In advanced interferometers, to isolatethe detector from possible vibration pollution, all the main optical benches willbe suspended in vacuum and the optical windows between tower and pipes havebeen replaced with cryogenic vacuum traps [25].

24 Gravitational waves interferometric detectors

– Magnetic noise. While low frequency magnetic noise can affect directly thesensitivity curve, high frequency electromagnetic field can couple with the lasermodulation frequency thus the bandwidth around that frequency must be cleanedup. The main magnetic noise sources consists in the mirror local control system,where coils and magnets can pick up environmental magnetic field (and itsgradient). All the electronic, that can be source of magnetic noise, must be placedon a safe distance from the mirrors, while magnet and coils must be as weak aspossible.

– Sideband line. The strongest noise lines, for instance electric main frequencyand suspension violin mode, are surrounded by several sideband lines, probablyrelated to non-linear coupling with suspension mechanical modes.

The detector characterization consists in following up each transient or persistent noisefeature to understand the source and the coupling mechanism. Then the next step is topossibly remove or mitigate the noise source and second to reduce the coupling with thedetector output.

In order to detect noise sources, environmental monitor devices are installed: seismome-ters, accelerometers, microphones, magnetometers, radio frequency receivers. Channelscoming to these devices are recorded together with hundreds of channels that record thedetector signals from the interferometer subsystems and the global control that can helpto identify noise generating within the detector. All of these channels (named auxiliarychannels) are analyzed to search correlations with noise features in the detector output.

Several software have been designed to understand persistent noise and glitches. One ofthem is the Omega tool, a detector characterization software which is able to scan a hugenumber of auxiliary channels beside the output channel, to search for excess of energy in atime-frequency space[54]. This search permits to identify coincidences between glitches inthe output and the auxiliary channels and similarities in the time-frequency behavior, helpingto locate the noise source. NoEMi (Noise Event Miner)[55] is another tool which analyzesthe output channel and a large number of auxiliary channels, but it is devoted to spectral lineidentification. This tool is developed to run on time following the non-stationary of the noiselines looking for coincidence in frequency and in the time they appear or change. MoreoverNoEMi gives alarms in case a noise line appears close to some known pulsar candidate asa continuous wave source. Nevertheless many non-linear mechanisms can produce up ordown conversion of the noise making the search of correlation between output and auxiliarychannels more challenging. Various tools mostly based on computing learning techniquesare developed to search for non-linear coupling mechanism. In chapter 3 a tool based onnon-linear system identification is described in details.

Chapter 3

Non-linear system identification: anapplication to detector characterization

The GW interferometer is in general locked in a dark fringe to maximize the sensitivity[27]. Hence the detector output signal is mainly the power of the dark fringe that is sensitiveto a GW event. However a great number of other signals are recorded to control theinterferometer in all of its parts, to check the status of the subsystems, to register theenvironmental conditions and so on. All those signals are called auxiliary channels. A GWdetector presents many kinds of non linear physical processes. Noise structures present inauxiliary channels, including narrow spectral features, may be converted linearly and nonlinearly into noise polluting the gravitational channel in a wide frequency region. Uncoveringsuch relationships between auxiliary channels and the gravitational wave channel can bevery useful to characterize the detector and possibly also to improve the confidence of theGW searches. Here we describe a non-linear system identification techniques developed toidentify linear and non-linear noise in gravitational waves detector output. The algorithmsperform the system identification purposes modeling the system with a non linear functionexpanded in a Volterra series. Series coefficients are found by means of Gram-Schmidtorthogonalization techniques. The main feature of the method is that it can compute thespecific contribution to the model of a channel, or a combination of channels, through theprocess of orthogonalization with respect to any other channels contribution. The magnitudeof the specific contribution is then given by a parameter called error reduction ratio (Q). Inthis way the redundancy produced by many inputs correlated to the same particular featureof the model, is prevented.

In section 3.1 the system identification problem is sketched. Two different algorithmsare described in section 3.2 and 3.3: the Orthogonal Least Square technique (OLS) and theSorted Fast Orthogonal Search (SFOS). A method to limit the system identification problem

26 Non-linear system identification: an application to detector characterization

in a desired region of the noise spectrum is described in section 3.4. Sections 3.5, 3.6 and 3.7are devoted to the application of OLS and SFOS to the Virgo and LIGO GW detectors noise,highlighting the advantages of the SFOS method in terms of computation speed and memoryusage. Problems related to the number of time lags of the Volterra series are discussed insection 3.8 while in section 3.9 a blind search approach to some Virgo noise structure ispresented. Finally in section 3.10 some conclusions are drawn.

3.1 The Non Linear Identification problem

The non linear system identification problem is shown in figure 3.1 in MISO (Multiple Inputsand Single Output) case.

Fig. 3.1 System identification. Basic scheme for MISO systems. u1(t) . . .ur(t) are the inputs,y(t) the output and e(t) the error function.

The detector is modeled as a box, with GW channel as the output and auxiliary channels,and non linear combination of them, as the inputs (figure 3.1). The system can be modeledas non linear function of r inputs u1(t) . . .ur(t) and one output y(t), where t is the index ofthe discrete time serie [56]:

y(t) = F(u1(t −1),u1(t −nu1), ...,ur(t −1), ...,ur(t −nur)) (3.1)

F is the function to be determined; inputs with different time lags are considered andnu1 . . .nur are the maximum numbers of lags for each input. The search of the function F isbased on the minimization of the error function (also called cost function):

∥Ξ∥= 1N

N

∑t=1

ξ (t)2 (3.2)

where ξ (t) = y(t)− y(t) is the difference between the observed and the modeled output.

3.2 Model structure determination with OLS algorithms 27

The function F can be expanded in a series of function. Using a polynomial expansion aVolterra series can be written as:

y(t) = θ0 +K

∑i=1

θixi(t)+K

∑i1=1

K

∑i2≥i1

θi1i2xi1(t)xi2(t)+K

∑i1=1

...K

∑in≥in−1

θi1...inxi1(t)... xin(t) (3.3)

where: x1(t) = u1(t −1),x2(t) = u1(t −2)...,xnu1(t) = u1(t −nu1),xnu1+1(t) = u2(t −1) are

the inputs at every considered time lag, K = nu1 +nu2 · · ·+nur is the sum of the time lagsand θi are the series coefficients. In a more compact form the model can be written as:

y(t) =M

∑i=1

pi(t)θi (3.4)

where pi are monomial of xi up to degree n henceforth called regressors and M is the numberof regressors. The equation 3.4 remains valid also if y and pi are vectors of length N; whereN times the sample frequency is the length in seconds of the analyzed time slot. In this way amatrix form for equation 3.4 can be written:

Y = PΘ (3.5)

with P the regressor matrix with regressors as columns and Θ is the parameter vector.The issue is to estimate the parameters which minimize the norm of the difference

between the output and the model ∥Y −PΘ∥. This is a linear least square problem. Theexpansion as in equation 3.3 produce a huge number of regressors proportional to the numberof the inputs, the order of the non linearity and the number of the time lags. Many regressorswill result negligible or redundant, hence a subset selection is necessary. This is a combinedproblem of structure selection and parameter estimation: find the subset Ps of P and findthe θs which fit the observed data. A useful technique is the orthogonal least square (OLS)algorithm. The orthogonal decomposition of the regressors matrix P may be carried out in anumerically stable manner with the modified Gram-Schmidt (MGS) procedure.

3.2 Model structure determination with OLS algorithms

The orthogonal least square technique allows to estimate the parameter of the model alsoevaluating the contribution of each regressor to the model. Moreover thanks to the orthogo-nalization procedure the evaluated contributions of each regressor is orthogonal with respectto the others. Therefore each regressor of the final list of most significant ones, modelsa different feature of the noise structure. The regressor matrix P in equation 3.5 can be

28 Non-linear system identification: an application to detector characterization

factorized by mean of Gram-Schmidt procedures in the following way:

P =WA (3.6)

where:

A =

1 a1,2 a1,3 · · · a1,M

1 a2,3 · · · a2,M. . . . . . ...

1 aM−1,M

1

(3.7)

is an M×M upper triangular matrix with M the number of the regressors, W is an N ×Mmatrix with orthogonal columns and N is the time series length. W satisfies the propertiesW TW = D with D a positive diagonal matrix. Defining g as follows:

g = D−1W TY (3.8)

the parameters θ1,··· ,M can easily be computed from:

AΘ = g (3.9)

using backward substitution.

To factorize the matrix P a classical Gram-Schmidt (CGS) and modified Gram-Schmidt(MGS) procedure can be considered [57]. In CGS the matrix A is computed one column at atime and W is evaluated at each step taking a new column of P and making it orthogonal withrespect to the previous ones. On the contrary in MGS the matrix A is evaluated one row at atime and W is computed at each step making all the remaining columns of P orthogonal tothe column considered in that step. From numerical point of view MGS is much less affectedby rounding errors that limit CGS in practical cases. On the other hand in MGS procedurethe storage of the entire regressor matrix P is required and this can be an issue if the numberof the regressors is quite high. The second parallel issue of the combined problem of systemselection and parameters estimation, is then the selection of the most significant regressorsand ultimately the determination of the model structure. Assuming the choice of a subset ofMs regressors in a matrix Ps from the M regressors of the matrix P, the error function can bewritten:

Ξ = PsΘs −Y = PsA−1s AsΘs −Y =Wsgs −Y (3.10)

3.3 The fast orthogonal search for system identification 29

The square norm of the vector Ξ is:

∥Ξ∥2 =< Ξ,Ξ >=Ms

∑i=1

g2i < wi,wi >−< Y,Y > (3.11)

The reduction of the mean square error Q, also called error reduction ratio, produced by thei-th regressor is:

Qi =g2

i < wi,wi >

< Y,Y >(3.12)

The values of the Q coefficients give the magnitude of the contribution of each regressor tothe model. Hence Q values can drives the regressors selection and the structure determination.Moreover the sum of the Q values represents a good estimation of the model accuracy.

3.3 The fast orthogonal search for system identification

Chosing the candidates with Q greatest of a specified threshold and rejecting the others canallow to build up an accurate and parsimonious model of real systems. But the evaluationof Q values in OLS algorithm, requires the computation of the orthogonal vectors wi thatis expensive in computation time. Moreover in MSG the entire regressor matrix P must bestored to perform the orthogonalization and to evaluate the Q values. The number of theregressors may be huge also for low order nonlinearities, because of the high number ofpossible candidates in the determination of a total unknown system, so the memory usage isan important issue that can be vexing in OLS-MGS system identification method. The socalled fast orthogonal search (FOS) can be a solution with some trade off to asses [58]. FOSis a system identification technique based on Gram-Schmidt factorization method. Basicallythe strategy of FOS is to create a set of recursive difference equation, that produces theelements of the matrix A and the vector Θ of equation 3.9, avoiding the creation of theorthogonal base vectors wi. The robust fast orthogonal search (RFOS) uses the same strategyof FOS but the factorization in equation 3.6 is carried out using MGS technique that is morenumerically reliable. The candidate regressors (Pi) are added to the model one at a timeand the respective Qi value can be evaluated without the computation of wi [59]. Then ifQi is greatest than a certain threshold the candidate Pi is added to the model, otherwise itis discarded. As mentioned before, in MGS procedure the matrix A is computed one rowat a time hence it is evident that the matrix A has to be totally recomputed whenever a newregressor is added. This may look like a trouble, but often in physical system only fewregressors (some tens or less) bring a significant contribution and will be chosen, althoughthe initial number of candidates can be of the order of thousands.

30 Non-linear system identification: an application to detector characterization

The RFOS can be further improved arranging the list of candidates by some characteristic,increasing the probability to treat first the most relevant terms with a sorted fast orthogonalsearch (SFOS) [60]. For example first we can select the linear regressors then the bilinearones and so on, taking into account the nonlinearity degree. The process can be stopped whena model with the required accuracy is produced. Another possibility is obtained by sortingthe list on the base of the regressors correlation with the output. In fact it can be provedthat the candidates with the highest contribution in the model, are the one with the highestcorrelation with the input and it was also shown that this kind of arrangement produces aneffective speed up in the achievement of a desired level of model accuracy.

3.4 Model identification with error filtering

The output of an GW interferometer presents a complex behavior. It is normally impossibleto consider GW output only as a function of the copious auxiliary signals that are recorded.A huge number of random and deterministic processes determine both stationary noisestructures and narrow time glitches. Therefore, the claim to consider the GW detector asa black box and find a good model of that is hopeless. However it is possible to confinethe problem around a region of frequency spectrum where some kind of noise structure ispresent. Possibly, only few noise sources will be in charge for that structure, hence if theinputs list contains auxiliary channels which record the effects coming from that sources,a model can be obtained. OLS or FOS techniques allow to evaluate the magnitude of thecontribution of the single auxiliary channel, or of a combination of them, to the model, that isto understand how much an auxiliary channel is involved in a certain noise feature, hopefullygiving good hints to drive out the related noise sources.

The confinement of the system identification problem cannot be achieved simply filteringall the signals in the interested band. Non linear phenomena can shifts the frequency of anoise source in a completely different region of the spectrum. For instance the very lowfrequency (< 1Hz) noise related to the pendula motion in a GW detector can be shifted by amodulation with some carrier at high frequency. At the same time the beating between twohigh frequency lines can produce a low frequency signal. A solution is the error filteringapproach [61]:

∥ΞF∥= ∥(Y −PΘ)F∥= ∥YF −PFΘ∥ (3.13)

where F means a linear filter application. The important point is that filtering is applied tothe regressors and not to channels, as the frequency shifts produced by channels combinationare preserved. Some attention has to be dedicated to search a suitable filter for the purpose.In order to preserve the relative phase among the signals, zero-phase filters can be used. Zero-

3.5 Applications of system identification OLS technique to Virgo data 31

phase filters work by processing the input data both in the forward and reverse directions:after filtering the data in the forward direction, the filtered sequence is reversed and runs backthrough the filter. This provides zero-phase distortion and a filter order that is double theorder of the original one. A 4th order IIR Butterworth has been used to meet the requirementsin terms of cut-off slope and computation speed.

The choice of the sampling frequency is also a crucial point. The channels are sampledat different rate but a common sampling frequency must be set to perform the systemidentification procedure. According to the Nyquist theorem, the natural choice might seemsto be to double the characteristic frequency of the considered noise structure. But in thatcase, non linear effects, such as frequency down-conversion, would be discarded. On theother hands long time slots are necessary to catch up-conversions of low frequency signals.High sampling rate and long time slots means big amounts of data to manage, therefore atrade-off solution must be carefully assessed.

3.5 Applications of system identification OLS technique toVirgo data

System identification techniques using OLS and SFOS methods have been applied to Virgodata. The gravitational wave interferometric detector Virgo [21] has been consider has ablack box, with a number of auxiliary channels as the inputs and the dark fringe signal as theoutput in a MISO (Multiple Inputs and Single Output) configuration. The dark fringe signalpresents a lot of peaks and noise structure as shown in figure 3.2.

Fig. 3.2 A view of the Virgo interferometer dark fringe signal (Pr_B1_ACp), during VSR4.

The system identification procedure can be divided in steps. The first step is to select thenoise structure and to choose a number of auxiliary channels that are suspected to be related

32 Non-linear system identification: an application to detector characterization

with the noise. Then the signals have to be read choosing a suitable time slot and a commonsampling rate. The problem must be confined in a band around the noise structure (± somehertz) performing error filtering.

In the analysis of Virgo data, OLS technique with both CGS and MGS has been applied.To avoid the analysis of a too high number of regressors, only linear and bi-linear combinationhave been taken into account. A number from 3 to 5 time lags has been considered, takingonly combination of channels with the same time lag.

To test the goodness of the procedures, an artificial noise structure has been injected inthe output in a relatively quite region of the spectrum. MATLAB environment has been usedfor the algorithm coding and the native MATLAB function linsolve using LU factorizationhas been used to test the result of OLS procedures with a few inputs. The test has shown thatCGS fails because of round-off errors even with only few inputs (> 10) are involved. Thus,only MSG technique will be used from now on.

The number of possible regressor and the computing time, for a given non-linearitydegree, grows quadratically as a function of the number of channels. A non negligible part ofcomputing time is devoted to the signals reading and decimation. The same channel signal isused several time in different combinations. Three possible options can be chosen in the code.The channels signals can be read, decimated and stored in memory (1) or written in the localdisk (2) or else they can be read and decimated every time they are used (3). These solutionshave a respectively increasing computation time, while the memory usage decreases. A tradeoff has to be assessed pondering the memory capability of the used machine.

3.5.1 Calibration lines modulation noise

Calibration lines, pure sine signals at different frequencies, are injected in various pointsof the interferometric system to meet the control procedures and to calibrate the sensitivitycurve. Because of detector’s nonlinearities some noise structure combines with calibrationlines in an amplitude modulation process. This produces characteristic sidebands around thecalibration peaks. OLS with MSG have been applied on those noise structures around twocalibration lines at 113 Hz and 444 Hz and a list of 32 auxiliary channels have been selectedas inputs. That list has been filled following the indication of the commissioning crew aboutthe channels which could be possibly involved in that kind of noise. Since the calibrationsignal was not recorded, a pure sine has been added to the list of the selected channels. Atime slot of 600 s of the fourth Virgo Science Run (VSR4) data was analyzed with a samplingfrequency of 2 kHz and a bandwidth of ± 2 Hz around the calibration line frequencies.As mentioned before OLS with MSG procedure needs to store the whole regressor matrix.Taking into account linear and bilinear combinations, N inputs yield M = N(N+3)

2 regressors

3.5 Applications of system identification OLS technique to Virgo data 33

(a) (b) (c)

Fig. 3.3 Noise modulation of calibration line at 113Hz. In (a) the Pr_B1_ACp signal (thedark fringe) is shown. The noise lines around 113Hz are visible. The model output is thegreen line in (b) while the red line in (c) is the Pr_B1_ACp signal cleaned by the modelednoise.

per time lag, thus N=32 correspond to a 1200000x560-per-time-lag regressor matrix size,which poses a serious matter in terms of memory usage. To pick up the most significantregressors, a preselection has been made on the base of the correlations of each regressorwith the output. A preselection has been made on the base of the correlation between eachregressor and the GW output. A subset of the Mp << M highest correlation regressors havebeen selected, with Mp depending on time slot length, sampling rate and memory capabilityof the computing machine. However, this operation is quite tricky. In fact, it may be thatthe regressors having the higher correlation with the output can be very similar and theycan contain contributions for the same part of the noise structure, such as a carrier for amodulation-like noise, and the regressors related to the rest of the noise structure, for instancethe sidebands, can be discarded. This makes the preselection the most knotty problem inOLS MSG system identification approach.

In figure 3.3 the result of the OLS MSG procedure are shown and in figure 3.4 there arezoomed views of the same result. To check the goodness of the procedure, a subtracted signalis also shown. The output is well reconstructed and a list of the main regressor used in themodel is in table 3.1 where sin113 means the sine at 113 Hz for the calibration line.

Regressor name Time lag Qsin113 2 0.476

sin113*V1:Gx_OB_ty 2 0.200sin113 1 0.161

sin113*V1:Gx_PR_PSDm_PosX 3 0.067sin113*V1:Gc_Al_DiffEnd_ty 1 0.051

Table 3.1 A list of the most significant regressors with time lag and error reduction ratio (Q)in OLS MSG model of 113 Hz calibration line modulation noise.

The most important regressor is the pure sine related to the central peak while a seconddegree regressor settles the sidebands. In figure 3.5 the results of OLS MSG identification

34 Non-linear system identification: an application to detector characterization

(a) (b)

Fig. 3.4 A zoomed view of the noise modulation of calibration line at 113Hz. In (a) thePr_B1_ACp signal (the dark fringe) in blue and the reconstructed signal in green are shown.In (b) the red line is the cleaned output signal. Analysis parameters: GPS start 991157550,duration 600s, filter range 111-115 Hz, N. of time lags 5, N. of input channels 32.

procedure applied to 444 Hz calibration line modulation noise is shown. Some peaks are notreconstructed, this means that the noise accountable to these peaks are not contained in theconsidered input channels or in their combination up to the second degree. The list of themost important regressors are shown in table 3.2

(a) (b)

Fig. 3.5 A zoomed view of the noise modulation of calibration line at 444Hz. In (a) thePr_B1_ACp signal (the dark fringe) in blue and the reconstructed signal in green are shown.In (b) the red line is the cleaned output signal. Analysis parameters: GPS start 991157550,duration 600s, filter range 442-446 Hz, N. of time lags 5, N. of input channels 32.

3.5.2 The Vela case

During VSR4 a series of noise peaks affected the Virgo sensitivity in the frequency regionof the expected gravitational wave signal coming from Vela pulsar. Previous studies proved

3.5 Applications of system identification OLS technique to Virgo data 35

Regressor name Time lag Qsin444 3 0.693

sin444*V1:Gc_Al_DiffEnd_ty 4 0.150sin444 1 0.161

sin444*V1:Gc_Al_PR_ty 4 0.010sin444*V1:Gc_Al_DiffEnd_ty 1 0.006

sin444*V1:Gc_Al_CommEnd_tx 4 0.004Table 3.2 A list of the most significant regressors with time lag and error reduction ratio (Q)in OLS MSG model of 444 Hz calibration line modulation noise.

that the origin of that noise was a third degree modulation of two calibration lines and theresidual motion of the interferometer pendula due to a non linearities in an ADC behavior[51]. Because we already know the non linear coupling, we checked that our procedurecould effectly find the combination of channels involved. OLS with MSG procedure hasbeen applied. The channels related to the known noise source have been inserted in a list of30 channels randomly chosen. Non linear combinations up to degree 3 have been taken intoaccount. The constrain that in the third degree combination at least one channel related to theresidual pendula motion should be contained, has been applied to avoid a very large numberof regressors. Indeed this constrain is not very restrictive because the relationship betweenthe noise and the pendula motion is evident from the spacing of 0.2 Hz among the carrier andthe sidebands, the well known fundamental frequency of the Virgo suspensions. In figure 3.6the result is shown and in table 3.3 the most significant regressors are listed. The proceduresucceeds to recognize the channels related to the noise. The first regressor of table 3.3 is thethird degree combination of a signal with the 0.2 Hz pendula motion with the two calibrationline identified in the previous analysis. The small values of error reduction ratio are due tothe small amplitude of the noise structure with respect to the background noise that remainunmodeled.

Regressor name Time lag QV1:Sc_NE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Sc_NE_zCorr356Hz 2 0.048V1:Sc_NE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Sc_NE_zCorr356Hz 3 0.043V1:Sc_NE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Sc_NE_zCorr356Hz 1 0.005V1:Sc_WE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Gc_DARM_AddExc 2 0.003V1:Sc_WE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Sc_WE_zCorr356Hz 2 0.003V1:Sc_WE_zCorr0.2Hz*V1:Gc_DARM_AddExc*V1:Sc_WE_zCorr356Hz 3 0.003

Table 3.3 A list of the most significant regressors with time lag and error reduction ratio (Q)in OLS MSG model of noise in Vela pulsar frequency region.

36 Non-linear system identification: an application to detector characterization

(a) (b)

Fig. 3.6 The noise in the region of Vela pulsar GW signal. In (a) the Pr_B1_ACp signal (thedark fringe) in blue and the reconstructed signal in green are shown. In (b) the red line is thecleaned output signal. Analysis parameters: GPS start 992721500, duration 600s, filter range21-24 Hz, N. of time lags 5, N. of input channels 8.

3.6 Applications of OLS on LIGO data

The OLS with MSG technique has been also applied to model the modulation noise of thepower line in LIGO H1 detector data. LIGO (Large Interferometer Gravitational Observatory)is composed of two interferometric GW detectors very similar to Virgo in terms of workingprinciple [20]. H1 is a LIGO detector set in Hanford WA USA (the other is placed inLivingstone LA USA). The power lines are noise lines due to the power distribution systemtransmitting AC power at typical frequencies of 60Hz in USA (50 Hz in Europe). So anoise line at 60 Hz and the related harmonics (120 Hz, 180 Hz... ) affect all the electricand electronic systems. In LIGO H1 detector the power line presents a structure composedof side modulation peaks produced by some non linear phenomena. A completely blindsearch has been tried to model this noise structure. Namely no assumptions have been donein auxiliary channels choice. A list of 630 channels has been considered and the procedureinvolved all of them. Since the possible linear and bilinear regressors are 199395 per time lagthat is a number that cannot be handled, the preselection phase was crucial. As mentionedbefore the preselection is based on the correlation between each regressor and the output.Nevertheless there is the risk to choose only regressors which contribute for the same partof the noise structure. The problem is amplified for power lines, because all signals containthese noise lines. Moreover the power line is much larger than the side peaks. Hence a largenumber of regressors show a high correlation with the output. For this reason the procedurehas been divided in two steps. The first step consists in the construction of a model of thepower line central peak only with linear regressors . A second step is the subtraction ofthe reconstructed signal from the output and then a new model identification procedure is

3.7 Applications of SFOS technique on Virgo data 37

made with bilinear regressors to model side peaks. Before beginning the bilinear modelconstruction the preselection is applied, to reduce the amount of the bilinear regressor to anhandleable number. To make hundreds of thousands of correlations in an allowable timelength, the correlation process has been parallelized in 315 jobs. Each of them computes allthe correlation of 2 of the 630 channels. The code was wrote in MATLAB environment andsubmitted with Condor framework for distributed computing in LIGO computing cluster. In

(a) (b)

Fig. 3.7 Subtraction of the reconstructed noise in power lines bands from the GW outputsignal of LIGO H1 detector. In red and green, the GW channel after the removals ofrespectively linear and bilinear noise model. For 120 Hz (a) and 180 Hz (b) power lineharmonics. Analysis parameters: GPS start 863635679, duration 600 s, filter range 118-122 /178-182 Hz, N. of time lags 5, N. of input channels 630.

figure 3.7 the result of the subtraction of linear and bilinear noise model from the output isshown for the first power line harmonics. It can be seen that the subtraction of the linearreconstruction cleans the output from the central peak while the bilinear works in side peaksreconstruction. A list of the most significant contribution to 180 Hz power line harmonicmodel is in table 3.4. In the upper part there are the regressors involved in linear model,below there are the regressors that model the output cleaned from the linear contributions.

3.7 Applications of SFOS technique on Virgo data

The OLS with MSG system identification technique is limited by the need to store the entireregressor matrix to determine the model structure. A preselection based on correlationbetween regressors and output can be useful to reduce the regressor matrix size. However, thepossibility to have a large number of regressors with high correlation but containing the samecontribution, is a problem to be faced. As explained in section 3.3 the fast orthogonal search

38 Non-linear system identification: an application to detector characterization

Linear ModelRegressor name Time lag Q [%]

H2:SUS-SM1_SENSOR_SIDE 3 81H1:SUS-MMT3_SENSOR_UL 5 0.01H2:SUS-FMY_SENSOR_LL 0 0.01

Bilinear ModelRegressor name Time lag Q [%]

H1:ASC-WFS2_IP*H1:SUS-MC1_SENSOR_UL 3 33H1:ASC-WFS2_IP*H1:SUS-MC1_SENSOR_SIDE 3 0.14H1:ASC-WFS2_IP*H1:SUS-MMT2_SENSOR_LL 3 0.02

Table 3.4 The most significant regressors for 180 Hz power line harmonic. Results of thelinear modeling of the GW output channel are reported on the top. On the bottom part thebilinear model of the output cleaned from the linear part

(FOS) is a system identification technique which allows to evaluate each regressor contributebefore the addition of the regressor to the model. Nevertheless correlation values can beuseful to sort the regressor list to speed up the construction of an effective model in a so-calledsorted FOS (SFOS). SFOS has been first applied in a case already studied with OLS to verifythe consistency of the results. In figure 3.8 a comparison among OLS and SFOS techniquesis done for the power line modulation at 50 Hz in Virgo data. The same list of 15 channel waschosen with a second order non-linearity, 5 time lags and no preselection done before OLS.Results are in good agreement in terms of signal reconstruction and regressor contributionestimation, with a time saving of about a factor 5 for SFOS technique. SFOS technique hasbeen applied also for modulation noise in 113 Hz and 444 Hz calibration lines. Figure 3.9 and3.10 show the results. The same list of 32 auxiliary channels has been used and 5 time lags.No preselection has been done for the 2800 regressors. The contribution of each regressorhas been evaluated setting Q = 10−4 as the acceptance threshold. The computing time hasbeen of the order of 50 minutes running on a quad-core 8Gb RAM machine. Comparingthis results with the ones find in section 3.5 some difference in the selected channels list areevident. This can happen because the same noise source could affect in a similar way severaldifferent channels. Hence, in some condition, a different reconstruction method or even samevariation in recontruction parameters can leads to a different channels selection. This meansthat the results must be interpreted by commissioning experts that can finally possibly linkthe selected channels to the noise source.

3.7 Applications of SFOS technique on Virgo data 39

Fig. 3.8 A comparison between OLS and FOS technique is shown in modeling of power linenoise at 50 Hz in Virgo VSR4. In blue the original output signal, the OLS model is in greenand the FOS model in dotted brown.

(a) (b)

Fig. 3.9 SFOS model of the 113 HZ noise structure in Virgo data. In (a), the original outputsignal (blue), the model (black) and the output signal after the model subtraction (red). Aplot bar of the most significant regressors is shown in (b). In ordinate axis the contribution tothe model construction is reported in percent. Analysis parameters: GPS start 991157550,duration 600s, filter range 111-115 Hz, N. of time lags 5, N. of input channels 32.

40 Non-linear system identification: an application to detector characterization

(a) (b)

Fig. 3.10 SFOS model of the 444 HZ noise structure in Virgo data. In (a) the original outputsignal (blue), the model (black) and the output signal after the model subtraction (red). Aplot bar with the most significant regressors is shown in (b). In ordinate axis the contributionto the model construction is reported in percent. Analysis parameters: GPS start 991157550,duration 600s, filter range 442-446 Hz, N. of lags 5, N. of input channels 32.

3.8 Relative phase and time lags

The choice of the values and the number of time lags is an important point. GW detectorsare a very complex machine in which many physical process have complicated interactionwith the system. Therefore a markovian hypothesis results too reductive. Moreover, even ifa noise source affect at the same time different detector subsystems, sensor time responseand signal transfer delay can change the relative phase of the signals. Thus the same noisefeature can be recorded in different inputs not at the same time and again it can affects theoutput at another different time. This makes the markovian hypothesis worth even less. It iseven possible that a noise feature appears in the output channels before than in the auxiliarychannels. That is the inputs can be in principle delayed with respect to the output. Hence,even the causality could not be respected. This last possibility has not been taken into accountin this work, so no negative time lags have been considered, but it could be interesting forfuture improvements. An empiric strategy has been followed for the time lags choice with theaim to have enough time lags to fulfill the model without producing a number of regressorstoo high to handle.

Cross-correlation function can help in the selection of the best time lags. The cross-correlation between regressors and output can be written as follows:

(y⋆ pi)(m) =N−m

∑j=0

y⋆( j+m)pi( j) (3.14)

where y is the output, pi the regressors and N the signals length. In general definition ofcross correlation m belongs to the range ]−N, N[. However,in this work it was decided to

3.8 Relative phase and time lags 41

preserve the causality and therefore m is taken greater than 0. Moreover a maximum timelag can be imposed, m

samplerate ∈ [0,T ]. A maximum of T = 1s has been chosen. The ideais to select the time lags with the highest values of cross-correlation. Considering a singletime lag, with the maximum cross-correlation value, the model results as good as the onebuilt with 5 time lags randomly chosen. Cross-correlation is more time expensive than thesimple correlation, but the number of time lags can be reduced. Accordingly computing timeof model building can be sped up. For the next analysis, to increase the performance in orderof noise phase reconstruction, three time lags are considered: the time lags of the maximumand the minimum cross correlation and a third time lag which is the average of the first two.

3.8.1 Cross correlation normalization

A cross correlation normalization is necessary to compare contribution of different regressors.Normalize the cross correlation means to force (x⋆ x)(0) = 1. To get this, signals have beennormalized subtracting the means value (in general zero because of filtering) and dividingby their standard deviation: y = y−y

σyand pi =

pi−piσpi

. Furthermore the equation 3.14 has tobe modified dividing the sum by the length of the vectors N. This kind of normalizationis called biased normalization because for m > 0 arguments vectors are cutted at N −m’ssample. Thus correlation with m > 0 are systematically depressed. A different normalizationcan be considered:

(y⋆ pi)(m) =1

N −m

N−m

∑j=0

y⋆( j+m)pi( j) (3.15)

This normalization is named unbiased normalization. Since T is always much smaller thanthe analyzed time slot, the difference between biased and unbiased normalization are lessthan 1%. Nevertheless this small effect can determine a different choice in the time lags of theregressors. When a lot of input channels are taken into account in the system identificationprocedure, many combinations present very similar values of correlation with the outputand even a small difference in cross correlation definition can produce a different choice ofregressors time lags, different regressor sorting in SFOS procedure and finally some changesin model structure selection. As already said at the end of section 3.7, this can happen ifsame noise source affect different channels. In this case several regressors can be highlycorrelated to each other in the interested frequency band and all of them will have a similarcorrelation with the output.

The unbiased definition of cross correlation has been chosen for our purpose.

42 Non-linear system identification: an application to detector characterization

3.9 Virgo noise modeling: a blind search approach

For OLS algorithm the preselection phase based for example on the simple correlation, asdescribed in section 3.5 and 3.6, has a crucial role to limit the number of regressors to analyzewhen a high number of input channels is considered. SFOS technique is a fast orthogonalalgorithm for system identification. Moreover it allows to estimate the contribution magnitudeof all the regressors just storing those actually involved in the model. Thus hundreds inputschannel case can be handled without encountering any memory problems. A list of 560channels has been drawn removing from the whole set of Virgo auxiliary channels thoseconsidered not safe, where not safe means channels that can contain the gravitational wavessignal. The SFOS algorithm has been runned with that list in input for some modulationnoise of Virgo calibration lines. Cross correlation has been used to choose the regressorstime lags as explained in section 3.8. Noise at calibration line 113 Hz, 444 Hz and 379 Hzhas been chosen to be reconstructed. A slot of 150 s of VSR4 data has been considered,sampled at 1 kHz. A zero-phase 4th order Butterworth filter has been performed with a bandof ±2 Hz around the calibration line.

(a) (b)

Fig. 3.11 Reconstruction of the 113 HZ noise structure using a blind search with a 560channels list. In (a), the original output signal (blue), the model (black) and the output signalafter the model subtraction (red). A plot bar of the most significant regressors is shown in (b).In ordinate axis the contribution to the model construction is reported in percent. Analysisparameters: GPS start 991157550, duration 150s, filter range 111-115 Hz, N. of lags 3, N. ofinput channels 560.

In figure 3.11 and 3.12 results for 113 Hz and 444 Hz are shown. It can be notedthat reconstruction and subtraction are better than in the case of section 3.7 which areobtained with a smaller subset of input channels which were already known to be involvedin calibration lines modulation. Regarding the most significant regressors in part b) of thefigures, it can be seen that no artificial pure sine are added to channels list to simulate thecalibration line. Anyway the calibration lines signals are picked up from channels directly

3.9 Virgo noise modeling: a blind search approach 43

(a) (b)

Fig. 3.12 Reconstruction of the 444 HZ noise structure using a blind search with a 560channels list. In (a), the fft’s of the original output signal (blue), the model (black) and theoutput signal after the model subtraction (red). A plot bar of the most significant regressorsis shown in (b). In ordinate axis the contribution to the model construction is reported inpercent. Analysis parameters: GPS start 991157550, duration 150s, filter range 442-446 Hz,N. of lags 3, N. of input channels 560.

affected by them. The channel recognized as the most involved in the main modulationprocess is the Gc_Al_DiffEnd_ty which is the same already found in section 3.7.

The 379 Hz case presents a huge amplitude of the central peak with respect to themodulations ones. To face this, a two step analysis has been performed. In the first step onlylinear regressors has been chosen. A reconstruction of the central peak is obtained and it hasbeen subtracted from the output. Then a reconstruction of the subtracted output with bilinearregressor is done to reconstruct the modulation peaks. Results for the linear reconstruction ofthe central peak are shown in figure 3.13. In this case the calibration signal is recorded in

(a) (b)

Fig. 3.13 Linear reconstruction of the 379 Hz noise structure with a 560 channels list usinglinear regressors only. In (a) the original output signal (blue), the model (black) and theoutput signal after the model subtraction (red). A plot bar of the most significant regressorsis shown in (b). In ordinate axis the contribution to the model construction is reported inpercent. Analysis parameters: GPS start 991157550, duration 600s, filter range 377-381 Hz,N. of lags 3, N. of input channels 560.

Gc_DARM_AddExc channel and this allows a good reconstruction of the central peak. Infigure 3.14 the bilinear reconstruction is shown.

44 Non-linear system identification: an application to detector characterization

(a) (b)

Fig. 3.14 Bilinear reconstruction of the 379 HZ noise structure after the subtraction of thelinear model. In (a), the fft’s of the original output signal (blue), the model (black) and theoutput signal after the model subtraction (red). A plot bar of the most significant regressorsis shown in (b). In ordinate axis the contribution to the model construction is reported inpercent. Analysis parameters: GPS start 991157550, duration 600s, filter range 377-381 Hz,N. of lags 3, N. of input channels 560.

The time slot analyzed in this case is 600 s long to have a better resolution to reconstructsome peaks which are very close to the central line.

3.9.1 Computation time and parallelization

The computation time for blind searches is of the order of 8 hours for a data slot of 150 sin a completely dedicated quad-core 8Gb machine. The time increase quite linearly withthe analyzed time hence for 600s of data more than 1 day is needed. Anyway it is possibleto parallelized a part of the SFOS algorithm to gain computation time. About 40% of thecomputation time is dedicated to cross correlation estimation for the whole set of regressorsand this phase can be completely parallelized. In figure 3.15 the gain in computing time withrespect to the number of cores is shown. It can be seen that the saved time with 8 cores isvery close to the maximum and is quite useless to go further in core number.

3.10 Final remarks

The search of noise source that spoils the sensitivity curve of a GW detector is the firstcommissioning task. In a few years new advanced GW detectors will be ready to run. Newinstruments will be affected by lots of new and unknown noise sources. A tool whichgives informations and hints about noise source localization, can be very useful. Thiswork is devoted to develop that kind of tool basing the procedure on nonlinear systemidentification techniques. The procedure as been coded in MATLAB language and showsgood performance. Two different system identification techniques, using Volterra regression,

3.10 Final remarks 45

Fig. 3.15 Time saving vs core numbers for a partially parallelized code.

have been applied to some noise structure in LIGO and Virgo data. Sorted fast orthogonalsearch technique has shown to be effective to perform blind searches among several hundredsof auxiliary channels in linear and bilinear combination. Cross correlation works well in theinputs - output dephasing identification, allowing to minimize time lags necessary to noisereconstruction.

However a further important step will be to test the code with the commissioning crew. Afield-test will tell us if the tool is effective, or there is something to change or to improve.If the usefulness will be proved, the possibility to integrate the code with the other existingdetector characterization tools will be taken into account.

On the other hand we plan to investigate the potentiality of system identification techniquein glitches source analysis. Glitches are transient excesses of power in the output signals dueto something happened inside or around the detector, that can be mistaken for a GW signal.So another big issues for commissioning crew is to eliminate as much as possible the glitchsources. The idea is to use the system identification techniques also for glitches analysis.The aim is to identified channels that are involved in false alarm events obtaining informationabout glitch sources, in a similar way to what developed for persistent noise.

Chapter 4

Gravitational waves from compactbinary coalescence

Binary systems consisting of two compact stellar remnants like white dwarf, neutron stars andblack holes, lose their orbital energy by emission of gravitational waves; this leads eventually,to the coalescence of the system and the merging of the components. The last part of the lifeof binaries composed by neutron stars (NS) and black holes (BH), represents the brightestsource of gravitational waves in the frequency band of the ground based interferometricdetectors. Moreover the GW emission can be accurately modeled. Indeed considering thesystem close to the coalescence, the process can be divided in three phases:

• the inspiral phase: where the two objects orbit around the center of mass loosing energyand approaching each other;

• the merging phase: when the two objects fall one on the other collapsing in a singleremaining object;

• the ring down phase: if the final object is a black hole it will be characterized byquasi-normal mode oscillation towards the equilibrium.

The gravitational wave emission on the inspiral phase can be analytically modeled as adominant Newtonian term plus a series of post-Newtonian terms to improve the prediction ofthe expected signal. Numerical relativity can also model the merge, while reliable analyticmodel can predict the ring down phase, thus a complete knowledge of the signal produced bycompact binary coalescence (CBC) is available. Matched filtering technique with accuratetemplates can be used to extend the observable horizon of the advanced detector to hundredsof Mpc.

48 Gravitational waves from compact binary coalescence

A big amount of electromagnetic energy can be emitted during the merger of the twocompact objects, making this kind of process probably accountable for the so called shortgamma-ray bursts (GRBs). The join analysis between gravitational and electromagneticobservations consists both in electromagnetic trigger for a focused gravitational wave search(external triggered search) and in a electromagnetic follow-up of the gravitational wave event.

4.1 CBC formation scenarios

To be detectable by ground based interferometric detector, the gravitational wave sourceshave to emit on a frequency band that goes from ∼ 10 Hz to some kHz. This means that theyhave to be in rapid rotation in very thigh orbit, with a time to coalescence on sub-cosmologicaltime scale. The evolution scenario of this kind of systems, depicted in figure 4.1, has beenconfirmed by several years of astronomical observations [62], and is now considered asa standard. The steps to the formation of these systems, also called "relativistic compactbinaries", are sketched in the following:

• The process begins with a binary system composed by two OB main-sequence starswell inside their Roche lobes. The Roche lobe is the gravitational equipotential surfacebetween the two stars, it have a teardrop shape with a typical size Rlobe ≈ a

( m1m1+m2

)1/3,

where the masses of the heavier star and the companion are respectively m1 and m2,and a is the orbital separation. The tidal force prevents the filling of the lobe by themore massive star, removing a possible initial eccentricity. The system remains inthis state until the core hydrogen fuel of the heavier star is exhausted. This time isrelatively long compared to the time scale of the CBC formation and can take up to107 years. Thus a lot of such binaries are expected in our galaxy: ≈ 104.

• After the core hydrogen of the heavier star is run out the star begins to expand,overflowing its Roche lobe. Then matter begins flowing into the less massive starincreasing its mass and consequently decreasing the orbital separation. Indeed fromthe orbital momentum conservation law J = µ

√G(m1 +m2)a and considering the

total mass (M) conservation, when the mass m2 is increased by the in-falling matter,the reduced mass µ = m1m2

m1+m2goes up and the orbital separation a goes down. However,

mass lost by stellar wind can unbound the system in particular for higher mass and highmetallicity stars. Low metallicity environments, are more transparent to the photonsand this means that they are less effective in pushing away the external matter and thestellar wind is depressed. For this reason, low metallicity environments are consideredthe best place for the birth of "relativitic compact binaries" in particular for higher

4.1 CBC formation scenarios 49

Fig. 4.1 Formation of neutron stars or black holes in close binaries: evolution scenario. Foreach evolutionary stage the typical time scale (T) and the estimated number of objects in ourGalaxy (N) are given [62].

50 Gravitational waves from compact binary coalescence

mass systems . After a typical time of the order of 104 years, the matter transfer stopswhen all the hydrogen envelope of the donor star is lost, leaving a Wolf-Rayet (WR)star with an helium naked core. The number of these systems in the Galaxy is of theorder of tens.

• The system can remain in this stage, a WR star with a heavy OB companion, can lastfor about 105 years. Hundreds of such binaries are expected to be in the Galaxy.

• At the end of the thermonuclear evolution of the WR star, if the remaining mass isenough (at least 8M⊙), a Ib (or Ic) supernova explosion occurs and a neutron star (NS)or a black hole (BH) if formed. The explosion can break the binary system and the twoobjects can be kicked away in this stage.

• If the binary survives to the supernova explosion a result should be a compact objectwith a rapidly spinning Be star (a B star surrounded by a ionized disk) as a companion.The orbital eccentricity is enhanced by the supernova explosion and this increases thestellar wind that brings matter to the compact object with powerful X-ray emission. Theduration of this stage depends on many parameters but ends when the core hydrogenrun out (≈ 104 years).

• The helium star expands and overflows the Roche lobe, until it creates a commonenvelope if the mass ratio is high enough. At this point the binary system can be foiledby the merging of the helium core and the compact object. In this case a recentlyhypothesized Thorne-Zytkow (TZ) object can be formed [63]. Possibly, single heavyneutron stars and black holes could descend by TZ objects.

• If the common envelope still contains the compact object and the helium star, a secondsupernova explosion can occur. Again the binary can be unbound and the two objectscan be thrown away by the supernova kick. Otherwise a close binary system of twocompact objects survives with an evolution completely determined by the gravitationalwave emission. This emission decreases the energy of the system, thus the orbitseparation will decrease until the merging of the two objects in a single compactremnant. If at least one of the two objects is a neutron star, the merging can be followedby a powerful electromagnetic emission as a short γ-ray burst. All the compact binaryobserved up to now are expected to merge not before 106 years.

Other formation processes are possible, as the dynamical capture of a second compact objectby the first one, moreover even the standard evolution scenario contains a lot of poorlyconstrained parameters, that makes the estimation of the formation rate largely uncertain.

4.2 Coalescence rate 51

4.2 Coalescence rate

To estimate the detection rate of the gravitational waves, it is necessary to predict thecoalescence rate density, that is the number of expected coalescences per unit of time andobservable volume. Unfortunately, this prediction is still affected by large uncertainties.Indeed, as already said in the previous section, the evolution model contains many parametersthat are not well constrained and the efficiency of the different formation model in creatingcompact binaries is highly uncertain. Estimate on the rate has been obtained by assumingphysical parameter distributions and simulating a large number of binary system formation[64]. Important constrains on binary neutron stars (BNS) rate, can be derived by astrophysicalobservations of binary system where one of the two objects is a pulsar. An estimate of thegalactic merger rate can be done evaluating the ratio between the number of observed binarypulsar (N), taking into account the selection effects [65], and their characteristic lifetime,being the lifetime the sum of the pulsar characteristic age at the moment of the observation(τc) and the time of the orbit decay, due the gravitational wave emission, up to the coalescence(τGW ): RG ≈ N

τc+τGW.

This statistic is largely dominated by short orbital period pulsars: in particular the nearbydouble pulsar system J0737-3039 with a short orbital period determines the value presentlyestimate of ≈ 10−5 mergers per years in our galaxy [66] in agreement with the theoreticalestimate from the formation models [64]. To extrapolate the rate beyond the Galaxy a scalefactor can be used. This scale factor can be evaluated as the ratio between the blue luminositydensity of the Galaxy and of the local Universe [67], being the blue luminosity correlatedwith star formation rate, or from the direct evaluation of the star formation rate between theGalaxy and the local Universe that gives a similar result [68]: Rv = 0.01RG[Mpc−3] whereRv is the merger rate in the local Universe.

For the black hole binaries things are more complicated since there are not direct ob-servation. Upper limits to the binary black holes are given by the observation of the highmass x-ray binaries that are considered as direct progenitors [69], while for NSBH systemsan upper limit can be extracted by the studies on the nucleosyntesis of heavy neutron-richelements by means of rapid neutron-capture process (r-process), in the matter resulting fromdisruption of the neutron star [70]. On the other hand, a lower limit for NSBH rate is fixedby the observation of the neutron star-black hole binary Cygnus X-1 that sets the Galacticlower limit to 10−9 per year [69].

52 Gravitational waves from compact binary coalescence

4.3 The signal from a compact binaries coalescence duringthe inspiral phase

The waveform of a gravitational wave emitted by a compact binary coalescence during theinspiral phase will be discussed in this section. Some assumption is made: the componentshave no spin and the eccentricity is negligible and thus the motion will be quasi-circular [71].The effect of the spin will be treated in later sections. Before entering in some details ofthe computation, it is useful to define some combination of the masses (m1,m2) of the twobodies, which are commonly used when dealing with CBC systems:

• the total mass: M = m1 +m2;

• the symmetric mass ratio: η = m1m2(m1+m2)2 ;

• the reduced mass: µ = m1m2m1+m2

;

• the chirp mass: Mc =(m1m2)

3/5

(m1+m2)1/5 .

The gravitational wave equations in a region far from the sources are written in section 1.2.2.In particular, choosing a transverse-traceless (T T ) gauge, the equation 1.28 can be written,while equation 1.32 shows how to rotate the reference system in a new one, where the axisx and y lie on the orbital plane of the two rotating object. In absence of spins the angularmomentum LN will not change, hence the orbital plane will not change as well. In this framecalled "source frame", the kinematic equations will be:

x(t) = r(t)cos(∫ t

0ω(τ)dτ

)y(t) = r(t)sin

(∫ t

0ω(τ)dτ

)z(t) = 0

(4.1)

choosing the direction of x-axis so that x(t = 0) = 0. In the source frame the quadrupolemomentum is easily computed as [3]:

Mi j(t) = µxi(t)x j(t) (4.2)

The rotations described in equations 1.30, 1.31 have to be applied to the equation 1.32 towrite it in the source frame. To simplify the equation, the Geometrized units where c = G = 1

4.3 The signal from a compact binaries coalescence during the inspiral phase 53

are chosen:

h+ =1D

[Mxx(cos2

φ − sin2φ cos2

θ)+ Myy(sin2φ − cos2

φ cos2θ)+

− Mzz sin2θ − Mxy sin2φ(1+ cos2

θ)+ Mxz sinφ cos2θ+

+Myz cosφ sin2θ]

h× =1D

[(Mxx − Myy)sin2φ cosθ + Mxy cos2φ cosθ −2Mxz cosφ sinθ

+2Myz sinφ sinθ]

(4.3)

The second derivative of the quadrupole of equation 4.2 must be computed, under thequasi-circular motion assumption rω >> r, to get:

Mxx =−Myy = 2µr2ω

2 cos(

2∫ t

0ω(τ)dτ

)Mxy =−2µr2

ω2 sin

(2∫ t

0ω(τ)dτ

)Mxz = Myz = Mzz = 0

(4.4)

The equation for the two polarizations of h is obtained plugging the upper equations inequation 4.3:

h+(t) =2µω(t)2r(t)2

D

(1+ cos2

θ)

cos(

2∫ t

0ω(τ)dτ +2φ

)h×(t) =−2µω(t)2r(t)2

D2cosθ sin

(2∫ t

0ω(τ)dτ +2φ

) (4.5)

Comparing these equations with the equations of motion 4.1, it is evident that the gravitationalwave phase is twice the orbital phase:

Φ(t) = 2∫ t

0ω(τ)dτ. (4.6)

The third Kepler law sets a first order relation between ω and r: ω2 = Mr3 . Hence, the

equations of h become:

h+(t) =2D

M5/3c ω

2/3 (1+ cos2θ)

cos(Φ(t)+2φ)

h×(t) =− 2D

M5/3c ω

2/32cosθ sin(Φ(t)+2φ) .

(4.7)

54 Gravitational waves from compact binary coalescence

4.3.1 The radiated power

The radiated power can be obtained plugging equation 4.7 for the two h polarizations in the1.34 which gives the radiated power of a gravitational wave. The temporal average gives:⟨cos(Φ(t)+2φ)⟩ = ⟨sin(Φ(t)+2φ)⟩ = 1

2 hence the radiated flux per unit of solid angledoes not depend on φ and can be written as:

dPdΩ

=2π

(Mcω

2

)10/3[(

1+ cos2 θ

2

)2

+ cos2θ

]. (4.8)

The integration over the solid angle gives the total radiated flux:

P =325(Mcω)10/3 . (4.9)

4.3.2 The phase evolution

The gravitational phase evolution can be found starting from the radiated energy. The energyfor a quasi-circular motion can be written as the sum of the kinetic and the potential energy:

Eorbit = Ekin +Epot =m1m2

2r− m1m2

r=−m1m2

2r=−1

2M5/3

c ω2/3 (4.10)

The gravitational wave radiation decreases the energy which becomes more and more negative.This means that the orbital radius will decrease and accordingly with the third Kepler law thefrequency increases. Looking to the equation 4.9, an increment of the frequency producesan enhancement of the radiated power in a positive feedback process that leads the binarysystem to the coalescence. The orbital energy can be rewritten in terms of the chirp mass,using the Kepler law to substitute r by ω:

Eorbit =−12

M5/3c ω

2/3 (4.11)

and the time derivative have to match the radiated flux of equation 4.9:

Eorbit =−13

M5/3c ω

1/3ω =−P. (4.12)

to give:

ω =965

M5/3c ω

11/3 (4.13)

According to the quasi-circular motion condition, the angular acceleration is much lowerthan the centripetal one: ω << rω2, hence the right-hand member of the equation above

4.3 The signal from a compact binaries coalescence during the inspiral phase 55

can be treated as constant in time. The time integration will give the orbital frequency timebehavior:

ω(t) =(−256

5M5/3

c (t − t0)+ω8/30

)−3/8

(4.14)

where ω0 is the fiducial frequency at the fiducial time t0. If ω0 is taken to be infinite, that isthe orbital radius collapses to zero at the coalescence, and the time is substituted by a "timeto the coalescence" τ = tc − t, with tc the coalescence time, the previous equation can bewritten in an easier way:

ω(τ) =18

(τ

5

)−3/8M−5/8

c (4.15)

As already stated, the characteristic frequency of the gravitational wave is twice the orbitalfrequency:

ωGW =dΦ

dt= 2ω (4.16)

Putting some reference values helps to have an idea of the typical frequency and time of thecoalescence process. For example if the reference chirp mass is taken to be 1.21 M⊙, that isthe chirp mass of two 1.4 M⊙ neutron stars, we have:

fGW =ωGW

2π≈ 134Hz

(1.21M⊙

Mc

)5/8(1sτ

)3/8

τ ≈ 2.18s(

1.21M⊙Mc

)5/3(100HzfGW

)8/3

.

(4.17)

The second of the above equations states that the signal of a binary system composed by two1.4 solar masses, enters in the detector band, with a minimum frequency of 10 Hz, about 17minutes before the coalescence.

The phase Ω(t) can be found integrating the equation 4.15:

Φ(t) = 2∫ t

0ω(t ′)dt ′ = 2

∫τ

tcω(τ ′)dτ

′ = 2

[(τ

5Mc

)5/8

−(

tc5Mc

)5/8]

(4.18)

Finally if the equation 4.18 is put into the equations 4.7, the waveform is obtained.

To find the phase equation in frequency domain, the time to coalescence τ as a functionof ω can be extracted by rearranging the equation 4.15. Plugging the result in equation 4.18gives:

Φ( f ) =1

16f−5/3M−5/3

c −2(

tc5Mc

)5/8

(4.19)

being f the gravitational wave frequency: f = fGW = ω/π .

56 Gravitational waves from compact binary coalescence

To obtain the Fourier transform of a gravitational wave in useful to consider a compactform of the equation 4.7:

h+(t) = A(t)cos(Φ(t)). (4.20)

Since the wave amplitude has a much slower time dependence than the phase:

1A

dAdt

<<dΦ

dt, (4.21)

the stationary phase approximation can be applied. Applying the Eulero formula the Fouriertransform can be written as:

h+( f ) =∫

A(t)cos(Φ(t))ei2π f tdt =12

∫A(t)

(eiΦ(t)+ e−iΦ(t)

)ei2π f tdt

=12

∫A(t)ei(2π f t−Φ(t))dt +

12

∫A(t)ei(2π f t+Φ(t))dt

(4.22)

The second integral of the sum contains a rapidly oscillating function, hence it can beneglected in stationary phase approximation, while the exponent in the first integral isexpanded around the stationary point t∗ in which:

Φ(t∗) = 2π f . (4.23)

This is the point where the variable f of the Fourier transform is actually the wave frequency.The Taylor expansion read:

Ψ(t) = 2π f t −Φ(t) = 2π f t∗−Φ(t∗)− Φ(t∗)(t − t∗)2

2+O

((t − t∗)2) (4.24)

The Fresnel integral can help to integrate equation 4.22:

∫e−iαx2

=

√π

αe−iπ/4 (4.25)

applying this result we find:

h+( f ) = A(t∗)(

π

2Φ(t∗)

)1/2

ei(2π f t∗−Φ(t∗)−π/4) (4.26)

Equations 4.7 and 4.19 can be used to obtain a more explicit formula for h+( f ) to thedominant order and, using the same procedure, an equivalent equation for h×( f ) can be

4.3 The signal from a compact binaries coalescence during the inspiral phase 57

found:

h+( f ) =1D

(5

96

)1/2

M5/6c π

−2/3 f−7/6 (1+ cos2θ)

exp [i(Ψ+( f )−2φ)]

h×( f ) =− 1D

(5

96

)1/2

M5/6c π

−2/3 f−7/6 (2cosθ)exp [i(Ψ×( f )−2φ)]

Ψ+( f ) = 2π f t∗−Φ( f )− π

4Ψ×( f ) = 2π f t∗−Φ( f )+

π

4

(4.27)

4.3.3 The post-newtonian approximation

The gravitational waveform for a coalescing binary in inspiral phase in the Newtonianapproximation is in equation 4.7. This approximation can be used in weak fields, with thebackground space-time assumed to be flat, and for non-relativistic velocities. However, apost- Newtonian approach is used to refine the accuracy of the estimated waveform increasingthe knowledge of the signal that have to be extracted from the detector noise.

The post-newtonian (PN) approximation is a formalism to calculate the waveform phaseΩ beyond the dominant newtonian order [72], using a perturbative expansion in term of thecharacteristic binary velocity which at the leading term, following the Kepler law, can bewritten in G = c = 1 units, as:

v = (Mcω)1/3 (4.28)

From the relation above and the two equations 4.11 and 4.12, is easy to obtain a couple ofdifferential equations that determine the phase behavior Φ(t):

dΦ(t)dt

− v3

Mc= 0

dv(t)dt

− E(v)McE ′(v)

= 0(4.29)

where the prime symbol means the v derivative and the integration gives:

t(v) = tre f +Mc

∫ vre f

v

E ′(v)E(v)

dv

Φ(v) = Φre f +∫ vre f

vv3 E ′(v)

E(v)

(4.30)

58 Gravitational waves from compact binary coalescence

A series that gives E(v) and E ′(v) is obtained by the PN formalism up to a given power of v:

E(v) = ENv2

(1+

nE

∑i=2

Eivi

)

E(v) = FNv10

(1+

nE

∑i=2

1

∑j=0

Fi jvi log j v

) (4.31)

Conventionally the PN order is equal to the half of the power of v in the expansion, so theexpansion above has an order nE/2 on the energy and nE/2 on the flux. Currently the termswhich involve the orbital contribution are known up to 3PN on the energy and 3.5PN on theflux, while the terms involving the spin contributions are fully known up to 2.5PN order.

For a given order, different PN models expand and resum the terms of the fraction E(v)E ′(v) ,

in different ways giving different approximation. The PN approach is valid under the quasi-circular condition, ω << ω , so this approximation fails in the very neighborhood of themerger phase.

4.3.4 Waveform model

Different PN models [73] arise with different approach to the PN expansion of the quantity:

B(v) =−E ′(v)E(v)

(4.32)

It follows a summary of the main model families of PN approximation.

In TaylorT1 model, B(v) is computed by the ratio between the Taylor series of the energyand the flux as they appear in equations 4.31, and than the differential equation 4.29 aresolved keeping the polynomial ratio B(v) up to the relative PN order.

Also in TaylorT2 model the term B(v) is treated as a ratio of polynomial, but, looking tothe integrals equation 4.30, a pair of parametric equation for t(v) and Φ(v) is written:

Φ(T 2)n/2 (v) = Φ

(T 2)re f +ΦN(v)

n

∑j=0

φkvk

t(T 2)n/2 (v) = t(T 2)

re f + tN(v)n

∑j=0

tkvk(4.33)

In TaylorT3 model the parametric equations for Φ(v) and t(v) are built as in TaylorT2,then inverting the second equation in 4.33, the expression v(t) is obtained and used to find anexpression for Φ(v(t)).

4.4 Detector response to the inspiral waveform 59

Fig. 4.2 Relationship between source frame and detector frame [5].

The ratio B(v) is instead expanded to the relative PN order in TaylorT4 model and thenthe differential equations 4.29 are solved to find the relation Φ(t).

The TaylorF2 method is equivalent to TaylorT2, but in the frequency domain. The phasein equation 4.27 is expanded in series as up to the order n to obtain the n

2PN order:

Ψ( f ) = 2π f tc −Φc +n

∑j=0

1

∑k=0

λ jk f ( j−5)/3 logk( f ) (4.34)

where tc is the coalescence time and Φc a constant phase offset.

These methods, that use different perturbative approaches with different technique toexpand and rearrange the expansion terms, can give different results for the waveform phasethat can not perfectly overlap each other. This mismatch decreases at higher PN order, hencethe requirement is to use an order high enough to make the difference negligible.

4.4 Detector response to the inspiral waveform

The response of the detector to an incoming gravitational waves depends on the geometryof the detector, the source position in the sky and the gravitational wave polarization. Toevaluate the detector response it is useful to refer to a frame where the x and y directionspoint along the arms of the detector: the detector frame. The angles between the detectorframe and the source frame is depicted in figure 4.2. Two antenna pattern parameters can be

60 Gravitational waves from compact binary coalescence

defined as follows [4]:

F+(α,δ ,ψ) =−12(1+ cos2

α)

cos2δ cos2ψ − cosα sin2δ sin2ψ

F×(α,δ ,ψ) =12(1+ cos2

α)

cos2δ sin2ψ − cosα sin2δ cos2ψ

(4.35)

Hence the induced strain on the detector is given by:

h(τ) = F+(α,δ ,ψ)h+(τ)+F×(α,δ ,ψ)h×(τ) (4.36)

Plugging equations 4.7 into equation 4.36 the expression for the strain becomes:

h(τ) =(

1MpcDeff

)A(Mc)ω

2/3 cos(Φ(τ)−2φ0(α,δ ,θ ,φ)) (4.37)

or using the equation 4.27 in frequency domain :

h(τ) =(

1MpcDeff

)A′(Mc) f−7/6 exp [−i(Ψ( f )−2φ0(α,δ ,θ ,φ))] (4.38)

where A and A′ are constant factors that correspond to the amplitude for a 1Mpc distantsource, Deff is an effective distance that takes into account the actual distance of the sourceand its direction with respect to the detector:

Deff = D

√F2+

(1+ cos2 θ

2

)2

+F2× cos2 θ (4.39)

and φ0 is a constant phase offset.

As an example, the strain in a detector for an optimal-oriented source, placed directlyover the detector with the orbit lying on a plane parallel to the detector one, that is for(α,δ ,θ ,φ) = (0,0,0,0) and 1Mpc away, is plotted in figure 4.3 where the increase of thewave amplitude and frequency is evident producing a characteristic waveform commonlycalled "chirp".

4.5 Merging phase emission

When the two orbiting objects come to be very close one to the other, the quasi-circular orbitassumption fails, and the orbital radius starts to decrease more and more quickly as the tidaleffects between the two objects become not negligible. The two objects start to collapse

4.5 Merging phase emission 61

Fig. 4.3 A gravitational time domain waveform produced by binary neutron star system(1.4−1.4 M⊙ and non spinning) optimal-oriented, placed 1 Mpc over the detector in inspiralphase.

and end up to give a single oscillating black hole. During this merger phase the amplitudeof the emitted gravitational wave is larger than the inspiral emission even if it occurs for ashorter time. Up to now, there is no analytical approach to describe the merging phase, butin the recent years the big progresses in numerical relativity have been produced reliablewaveform for gravitational wave emission during the very late inspiral, the merging and thering-down oscillation of the resulting object [74]. Combining the numerical results with theNewtonian and post-Newtonian expansions, very accurate and complete waveform can beobtained for all the lifetime of the binary system from the inspiral phase, to the merging andthe ring-down (IMR waveform) [75].

Chapter 5

Detecting CBC signals

The matched filter is the optimal linear filter technique to extract a know signal from a gaus-sian stochastic noise. Since the gravitational wave emission from binary coalescence is prettywell modeled, the matched filtering is the natural choice [3]. However some complicationsarise: the detector background cannot be described as a stationary and Gaussian process andthe knowledge of the signal is not perfect because the waveforms are approximated up to agiven order and the model parameters of the incoming wave are unkwown. To reduce theeffects of the non-stationarities, vetoes can be defined to reject triggers on the base of thedetector behavior which is analyzed by thousands of monitoring channels. To overcome theissues related to the unknown parameters, template banks are built covering the parameterspace with a given accuracy and the filter output is computed summing in quadrature thefilter outputs that involved two orthogonal-phase templates. Different approximations can bechosen to obtain waveforms as close as possible to the expected signals, with different PNorders and involving more parameters, for instance, the spin of the single components andthe orbital eccentricity. Taking into account an higher number of parameters corresponds toan increment of the parameter space dimension and a growing up of the number of templatesrequired to cover effectively the space. At the end, more templates generate more falsealarms, therefore the pursuit of the best compromise between the efficiency in extractingsignals and the false alarm rate is the final goal of a search.

5.1 The matched filtering technique

Matched filtering technique is the standard filter choice when the expected signal is known.The technique is based on evaluating the correlation between the data coming from thedetector and a bank of template that must reproduce the expected signal [76]. Each templateof the bank is produced starting from a set of waveform parameters, i.e. the masses and the

64 Detecting CBC signals

spins of the binary system objects. Therefore, once the waveform equation is chosen, with agiven approximation and up to a given order, the bank can be considered as a collection ofpoints that span the parameter space. The filter output x j(t) for the template Tj is computedin the following way [77]:

x j(t) =(s;Tj

)(t) = 4ℜ

∫∞

s( f )T ∗j ( f )

Sn( f )e2πi f td f (5.1)

where s(t) is the detector output, that can be pure noise or noise plus signal, and s( f ) isits Fourier transform, Tj( f ) the Fourier transform of the template j and Sn( f ) is the noiseone-side power spectral density that weights the convolution giving more importance tothe region of the bandwidth where the detector is more sensible. The filter is applied bothin phase and in quadrature, in this way the maximization φ0 is inherently achieved. Thequadrature output is obtained multiplying the template in frequency domain by eiπ/2:

Tj( f ,2φ0 −π

2) = Tj( f ,2φ0)eiπ/2 = iTj( f ,2φ0) (5.2)

where, according to equation 4.7, 2φ0 is the waveform initial phase. Hence the phase andquadrature filtering can be performed considering a complex output of the matched filter:

z j(t) = x j(t)+ iy j(t) (5.3)

where y j(t) is the matched filter output for the dephased template iTj. To define the signal tonoise ratio, the filter output has to be normalized by the amplitude of the template. Usuallythe template waveforms are built to imitate a signal from a source with an effective distanceDe f f = 1Mpc, being the effective distance defined as in equation 4.39. The normalizationfactor will be:

σ2j =

(Tj;Tj

)(0) =

(Tj;Tj

)= 4

∫∞

|Tj( f )|2

Sn( f )d f (5.4)

In the case of purely stationary and Gaussian noise, the quantity σ j is a measure of thematched filtering output noise for the template j: σ2

j =⟨x j(t)2⟩= ⟨y j(t)2⟩. The signal to

noise ratio time series for the template j (ρ j(t)) is defined as:

ρ j(t) =

√√√√z2j(t)

σ2j

(5.5)

From the definition, for a pure Gaussian signal s(t), the signal to noise ratio (SNR) has a nullexpect value and a variance equal to one. Hence, when ρ j(t) presents peaks over a given

5.2 Generating the template bank 65

threshold (larger than one), the search will be triggered and the event will be collected as agravitational wave candidate, ready to be processed with other quality checks looking, forinstance, to the distribution of the SNR over the bandwidth (the χ2test) and to the coherencewith other detector events. The SNR ρ associated to an event will be typically the maximumpeak value of ρ j(t) over the template bank. From equations 4.38 and 5.5 an estimate of theeffective distance of the source for a detected signal can be achieved by the ratio:

Deff =σ j

ρ jMpc. (5.6)

5.2 Generating the template bank

To make matched filter technique effective, accurate templates have to mimic as best aspossible the waveform emitted by a compact binary coalescence. In particular the matchedfiltering is very sensitive to the phasing of the template which depends on the physicalparameters of the system. Thus, the astrophysics gives the boundaries of the parameterspace that must be covered generating a set of templates. Since the number of templateshave to be finite, the parameter space will be covered with a discrete lattice. Even if thegravitational waveform is accurately simulated, it is characterized by a set of parameterswhich are different from those of all the templates and thus the matched filter will loose someSNR. To evaluate the maximum SNR that it is lost due to the discretization of the parameterspace, the minimum match (or maximum mismatch) among templates has to be computed,defining the match between two templates as [78]:

M (Tj,Tk) = maxtc,φ0

(Tj;Tk

)√(Tj;Tj

)(Tk;Tk)

(5.7)

where Tj = T (ξ j) and Tk = T (ξk) are the templates generate by starting from the two sets ofparameter ξ j and ξk, tc and φ0 are the coalescence time.

The fitting factor (FF) for a template Tj is defined as the maximum match between thattemplate and all the template of the bank:

FF j = maxk

M(Tj;Tk

)(5.8)

The correspondent template mismatch will be 1−FF. The minimum fitting factor in atemplate banks is defined as [79]:

FFmin = minj

FF j (5.9)

66 Detecting CBC signals

that is basically a measure of the sparsity of the template coverage. Generating a templatebank, the minimum fitting factor typically accepted is 0.97; this means that for each templateTj there is at least another template Tk which has a match higher than 0.97. Looking toequation 5.5 and 5.7 it is easy to notice that the maximum mismatch (relative to the minimumfitting factor), represents the maximum relative loss of SNR due to the parameter spacediscretization that for a FFmin = 0.97 will be equal to 3 %.

There are mainly two methods to build a lattice of discrete points placing the templatesto cover the parameter space with a given minimum match: the stochastic placement and thegeometric placement. Using the stochastic template placement, the templates are generatedrandomly exploring the parameter space. If a new generated template has a minimum matchlower than the accepted threshold (usually 0.97) it is kept, otherwise it is discarded [80].This method is relatively easy to implement even for complex geometries of the parameterspace and complicated waveforms, the drawback is that the result is typically not optimal.This means that a number of templates larger than the minimum is required to cover theparameter space. On the other hand, the geometric template placement method aims to createa locally flat metric on the parameter space, placing the templates on a regular grid over thecoordinates of the defined metric space [81] [82]. The metric is defined setting a distancebetween templates that has to be related to the mismatch:

d(T (ξ j),T (ξk)

)=√

1−M(T (ξ j),T (ξk)

). (5.10)

This distance is expressed by means of a metric tensor g jk(ξ ) in a quadratic form:

d2(ξ ,ξ +δ ) = 1−M (T (ξ ),T (ξ +δ )) = g jk(ξ )δjδ

k (5.11)

where the metric tensor will be derived as:

g jk(ξ ) =−12

∂ 2M (T (ξ ),T (ξ +δ ))

∂δ j∂δk(5.12)

The coordinate system has to be well chosen to minimize the curvature of the metric space,so that the metric tensor g jk(ξ ) is constant in a large scale. This permits to use the defineddistance not only for infinitesimally separated points but also for finite neighborhood. Inthis coordinate system the space is covered by an hexagonal placement of the templatesconsidering for each template a circle with a radius d =

√1−FFmin. This method is efficient

in achieving the minimal number of templates for a given minimum fitting factor, which can

5.3 The aligned-spin template banks 67

be roughly estimated as follows:

N ≈ 1∆V

∫ √∥g jk∥dξ

p (5.13)

where p is the space dimension and ∆V is the volume of the element of the space covered bya single point which depends on the minimum fitting factor and can be expressed by:

∆V =

(2

√1−FFmin

p

)p

(5.14)

For non-spinning templates the intrinsic parameter space is one dimensional for thedominant Newtonian order, but it becomes two dimensional (on m1 and m2 or two differentcombinations of the masses) when PN terms are added. Up to 2PN, a good set of coordinatesto obtain a nearly flat metric space is derived transforming the physical masses in the socalled chirp times [82]:

τ0 =5

256π flowη(πMc flow)

−5/3

τ1.5 =1

8 flowη(πMc flow)

−2/3(5.15)

where flow is the low frequency cut off of the waveform, η and Mc the masses combinationdefined in section 4.3.

5.3 The aligned-spin template banks

Adding parameters to the waveform means to produce more accurate templates that canmake, on one hand, the matched filter more effective, but on the other hand it increases thedimensionality of the parameter space that have to be covered with the template bank [83].According to the observation the binary neutron stars, whose masses range from 1 to 3 solarmasses, have very small spin of the component stars, with a maximum value of χ equal to0.05 [84], being the dimensionless spin χ:

χ1,2 =S1,2

m21,2

(5.16)

with S1,2 the spins (the rotation angular momenta) of the neutron stars and m1, m2 the masses.However, isolated neutron stars have been observed to have a value of χ up to 0.4 and inthis case the contributes of the spin interactions can make the difference. Moreover, a wide

68 Detecting CBC signals

Fig. 5.1 Fitting factor distribution for signals from BNS systems with isotropic oriented spinsand a non-spinning template bank. The spin of the signals go up to 0.4 (blue solid line), 0.2(green dashed line) and 0.05 (red dotted line). The sensitivity of advanced LIGO in its finalconfiguration is considered. A flow of 15Hz is taken [86].

class of equations of state for neutron stars sets a limit of 0.7 [85] beyond which the staris disrupted. Figure 5.1 [86] shows a fitting factor distribution for a non-spinning templatebank, generated with geometric placement using TaylorF2 approximation up to 3.5PN order.The FF are computed for about 105 signals with masses in BNS region (1−3M⊙) and spinsuniformly distributed from 0 to 0.4 with isotropic orientation. An effective spin quantityχe f f can be introduced to encompass the effects of the two individual spins on the waveform:

χe f f =m1χ1 +m2χ2

m1 +m2(5.17)

For χe f f up to 0.4 more than 59% of the signals have a FF, smaller than 0.97, i. e. smallerthan the minimum fitting factor used to produce the bank. In other words, this means aloose in SNR larger than 3% for the 59% of the signals. The spin effects become even moreimportant when at least one of the two objects of the system is a black hole forming a binaryblack hole system (BBH) or a neutron star black hole system (NSBH). The spin of a blackhole can be of the order of one thus the effective spin can go up to one in BBH systemsbut also in NSBH systems, since the mass of the black hole is generally much larger than

5.3 The aligned-spin template banks 69

Fig. 5.2 Distribution of the fitting factor of signals from BBH systems with aligned spinsup to 1 and a non-spinning template bank. The sensitivity of advanced LIGO in initialconfiguration is considered with a low frequency cut-off of 40Hz [87].

the neutron star mass, so that the waveforms depend drastically on the spin values. Figure5.2 [87] shows the fitting factor distribution for binary black hole no-spinning templatebank, produced with stochastic placement and ispiral-merger-ringdown (IMR) waveforms.The bank is checked with BBH signals with spin aligned/antialigned with the orbit angularmomentum, with a magnitude up to 1. The figure points out that the condition FF > 0.97 isachieved only for absolute value of χe f f lower than 0.2.

If the spin values are considered as additional parameters, the number of the parameterspace dimension goes from two to eight, with a drastic growth of complexity for the templateplacement and an increment of the overall number of templates needed to cover the space.Roughly, the false alarm rate of a search (triggers generated by noise glitches) is directlyproportional to the number of templates, thus the improvement in SNR is swept out by anincreased number of false alarms that spoils the search performance. However, startingfrom the most reliable compact binary system formation scenario which foresees the binarygeneration from a common envelope, a state with object spins aligned with respect to the orbitangular momentum is more favored than a whatever spin orientation [62]. Considering onlythe case of aligned spins, with possible spin orientation aligned to the angular momentum,the number of extrinsic parameters of the waveform are limited to four: the two masses and

70 Detecting CBC signals

the two spins along the z-coordinate of the source frame. In this paper [86] an effectivesolution for generating the aligned-spin banks for BNS system is described and the improvingof the fitting factor for aligned-spin sources and also for isotropic spin sources is shown.The performance of aligned-spin banks for neutron star-black hole binary (NSBH) has beenalso proved [88]. To construct the aligned-spin bank with a geometrical approach, 3.5PNwaveform are taken into account with TaylorF2 model including the leading order of spin-orbit contribution at 1.5PN and spin-spin contribution at 2PN. According to equation 4.34the waveform phase in frequency domain can be written as:

Ψ(x) = 2π f0xtc −φc+

+ x−5/3(

λ00 +λ20x2/3 +λ30x+λ40x4/3 +λ51x5/3 logx+λ60x2 +λ61x2 logx+λ70x7/3)

(5.18)

where x = f/ f0, f is the frequency, f0 a fiducial frequency, tc the coalescence time and φc

the constant phase offset. To find suitable metric, the metric tensor is defined as in equation5.12. Hence in the phase equation, 10 parameters appear to give a 10 dimension parameterspace. The dependence of g jk on the constant phase can be removed maximizing the fittingfactor over φ0 [76]. It is convenient to write the parameters in equation 5.18 as λ (i) with ifrom 0 to 8, where λ (0) = tc. To remove the other extrinsic parameter λ (0) a preliminarymetric is defined as:

γmn =−12

∂ 2M (T (λ );T (λ +δ ))

∂δ (m)∂δ (n)(5.19)

where m and n go from 0 to 8. Then, minimizing the distance γmnδ (m)δ (n) over δ (0), thatmeans to project the tensor γmn onto a subspace orthogonal to λ (0), the tensor g jk is obtained:

g jk = γ jk −γ0iγ0 j

γ00(5.20)

where j and k go from 1 to 8. The tensor g jk gives the metric in a 8-dimensional space that ischaracterized by the eight parameters in equation 5.18, with the important property that itdoes not depend on the coordinates. Hence, this metric defines a globally flat space at theprice of increasing the number of dimensions from four to eight.

Rotating and rescaling the λ -coordinates it is possible to settle a orthonormal Cartesianµ-coordinate system :

µ j =8

∑i=1

V ij√

E jλi (5.21)

5.3 The aligned-spin template banks 71

where V ij is the ith component of the jth normalized eigenvector of g jk and E j is the corre-

spondent eigenvalue. In these new coordinates the metric tensor is simply an 8x8 identitymatrix. The range of physical sensible parameters is transformed in a variety in the newcoordinate system µ .

Looking to the correlation among the coordinates it is possible to find some principalcoordinates to reduce the number of the uncorrelated dimensions. The rotation of thecoordinate set according to the principal coordinates is obtained by the following transform:

ξ j =8

∑i=1

Cijµi (5.22)

where Cij is the ith component of the jth eigenvector of the covariant matrix from a set of

µ j which corresponds to a number of points placed in the physically sensible region of thespace. The result is that for points within the BNS range of masses and aligned spin up to 0.4the dominant directions are only two, along two coordinates which will be named ξ1 and ξ2.The extension of the space along the third and the fourth most relevant coordinates is of theorder of 10−3 and 10−4 smaller than the first, where the extension is measured in terms oftemplate size, being the template size the area covered by the template within the specifiedminimum fitting factor. Hence, in these coordinates the space is bi-dimensional, the metric isglobally flat and Euclidean. This makes very easy to cover the space with a desired minimumfitting factor.

The FF distribution of BNS bank with aligned spin up to 0.4, is shown in figure 5.3a,using BNS signal with isotropic spins up to 0.4.

Even if the spin of the testing signals are not aligned, producing a precession in the binarysystem orbit, the fitting factors for an aligned-spin bank show to be much larger than for anon-spinning bank. In this case the percentage of the signal with a fitting factor smaller than0.97 falls from 59% to 9%, figure 5.3b.

Even though the aligned-spin space results to be bi-dimensional like the non-spinningspace, the number of templates needed to cover the same range of masses without spin isabout ten times smaller than the number of templates with aligned spin up to 0.4. Thismeans that the improvement of SNR due to the better fitting factor, which reflects in abetter sensitivity of the search, has to be compared to the false alarm rate that is roughlyproportional to the number of templates. This has been done for BNS with aligned spin up to0.4, using the MBTA pipeline. The details of the search and the results are in chapter 7.

72 Detecting CBC signals

(a)

(b)

Fig. 5.3 Fitting factor distribution for BNS signals and aligned-spin template bank with spinmagnitudes up to 0.4. In figure (a) the signals come from sources with aligned spins andmaximum spin equal to 0.4. In figure (b) the signals have spins up to 0.4 (blue solid line),0.2 (green dashed line) and 0.05 (red dotted line) with isotropic orientation. Advanced LIGOin final configuration is considered with a low frequency cut-off of 15Hz [86].

Chapter 6

MBTA: a CBC low latency pipeline

The Multi Band Template Analysis, MBTA, is a pipeline devoted to the search of gravitationalwave signals from compact binary coalescences [89] [90]. In particular MBTA has beenspecifically designed to run in low latency, to flush out candidate gravitational wave eventsfew minutes after the detector data acquisitions. The low latency makes MBTA suitableto triggering electromagnetic follow-up observations[91]. Basically MBTA is a completepipeline that takes the data distributed by the detectors and gives the candidate events andtheir significance. The candidates which overcames a given significance are submitted to theGravitational Candidate Event Database: GraCEDb.

MBTA pipeline is composed by four main applications:

• MbtaRT: to search triggers in the data from the single detectors, applying the dataquality prescriptions;

• MbCoinc: to perform the time coincidence of triggers coming from different interfer-ometers and to estimate the significance of the candidates;

• MbClustering: to clusterize in time the triggers of a single detector or time coincidenttriggers of different detectors;

• MbAlert: to select interesting events and submit them to GraCEDb;

These applications will be shown in more detail in the following sections.

6.1 Filtering data in multiple bands

The MBTA application involved in triggers search from the data detector is named MbtaRT(MBTA Run Time). MbtaRT search engine is an implementation of the standard matched

74 MBTA: a CBC low latency pipeline

filter, which is commonly used in CBC searches, where the detector output is cross correlatedwith a bank of templates that estimates the predicted signal for a certain parameters space.As already discussed in chapter 5 the template banks is typically built to cover the parameterspace, with a minimal match of 97%. This means that less than 3% of SNR is lost bythe discretization of the parameter space. The template banks can be loaded in MBTA, orgenerated by it, during the initialization phase. For non-spinning template banks MBTA callsa LSC Algorithm Library (LAL) [92] function to produce the banks during the initialization.On the contrary, aligned-spin banks are not produced by MBTA, but they are generated inadvanced and then loaded when the initialization starts.

The basics feature of MBTA, which also gives the reason of its name, is that the matchedfilter is divided in two (or even more) different banks [93]. This separation has two evidentbenefits. The first one is that the chirp signal phase is tracked over fewer cycles, allowing tocover the same parameter space with less templates. The second benefit is that the samplingrate can be reduced for the lower frequency band, thereby shorter Fast Fourier Transform canbe used in filtering computation. The filtering is independently performed in the differentbanks, then the full band output is coherently computed. Basically with the same principleof the integration splits in sum of different parts, the combination can be sketched for twobands as:

S(t,M) =∫ fmax

fmin

h( f )T (M, f )d f =∫ fc

fmin

h( f )TLF(M, f )d f +∫ fmax

fch( f )THF(M, f )d f

(6.1)where fmin and fmax are the minimum and the maximum frequency of the full band, fc isthe separation frequency which divides the two bands, h( f ) is the frequency representationof the detector output and T (M, f ) is the template generated with the parameter set M, infrequency domain. The separation frequency fc is chosen to obtain a balancing of the SNRover the separated banks.

6.1.1 Real and Virtual templates

For each frequency band the filtering is done independently from the others, hence a templatebank must be provided for each band. The templates of these banks are called "Real", becausethey are really used in the filtering process. Another bank, containing the "Virtual templates",are produced for the full frequency band. Virtual templates are not used to filter the signal,but their parameters are associated to the triggers on the MBTA output.

The algorithm performs a coherent combination of the separated filtered signals, as shownin chapter 5. Both the phase and the quadrature matched filtering time series are computedfor each band:

6.1 Filtering data in multiple bands 75

• (h;RTLF)P (t): is the phase output of the low frequency real template RT;

• (h;RTLF)Q (t): is the quadrature output of the low frequency real template RT;

• (h;RTHF)P (t): is the phase output of the high frequency real template RT;

• (h;RTHF)Q (t): is the quadrature output of the high frequency real template RT;

where accordingly with the equation 5.1 :

(h;T )P (t) = 4∫ h( f )T ∗( f )

Sn( f )ei2π f td f

(h;T )Q (t) = 4∫ h( f )T ∗( f )

Sn( f )ei(2π f t−π/2)d f

(6.2)

with Sn( f ) the Power Spectral Density (PSD) of the detector noise that weights the filteringoutput. Those time series have to be translated and rotated to be aligned in time and phase tobe coherently summed as follows:

(h;V T )P (t) = (h;RTLF)P (t)+(h;RTHF)P (t +∆t)cos(∆φ)− (h;RTHF)Q (t +∆t)sin(∆φ)

(h;V T )Q (t) = (h;RTLF)Q (t)+(h;RTHF)P (t +∆t)cos(∆φ)+(h;RTHF)Q (t +∆t)sin(∆φ)

(6.3)

where V T is the virtual template in the full band, while RTLF and RTHF are the real templatesassociated to that V T , which means the RT s over the RT banks , which show the highestmatch with that V T in each frequency band. The match between VT and RT is computed asin equation 5.7:

M (V T,RT ) = maxtc,φ0

(V T ;RT )√(V T ;V T )(RT ;RT )

(6.4)

and basically tells us how much two signals are similar one to the other, being M = 1 forthe match of two identical signals. In figure 6.1 the signal recombination over the frequencybands is shown, while the output of the matched filtering is shown in figure 6.2.

The association parameters ∆t and ∆φ are respectively the offsets of the high frequencyreal templates in time and phase, that is the time the signal takes to chirp from the minimumfrequency fmin to the separation frequency fc and the phase in that point. The two param-eters are practically computed from the relationship between the virtual templates and theassociated real templates.

The matched filtering in the LF band is computed at a down-sampled rate; however, torecombine the signal an up-sampling of this part is needed. While the down-sampling is

76 MBTA: a CBC low latency pipeline

Fig. 6.1 Example of matched filtering outputs for a pure signal without noise. Top plots: theraw low frequency band output (blue), the interpolated low frequency band output (purple)and the high frequency band output (green). Bottom plot: the usual single band matchedfiltered output compared to the recombined MBTA matched filter output. The left plots arefor the template in-phase, right for in-quadrature [93].

6.1 Filtering data in multiple bands 77

Fig. 6.2 Quadratic sum of the in-phase and quadrature outputs for each bands. The highand low frequency bands are plotted in purple and green. The interpolated low frequency isin blue. In red the recombined signal. An usual single band analysis output is in black forcomparison. This is for the same example as Figure 1 [93].

78 MBTA: a CBC low latency pipeline

done in the frequency domain shortening the FFT, the up-sampling is done in time domainusing a quadratic interpolation.

The association among the virtual templates and the real templates for each band, aswell as the association parameter computation, is performed during the initialization phase.Finding the match among each V T and all the RT s has a quadratic computational complexity.To strongly reduce the association time for a V T , a preliminary selection of the RT s is doneon the base of the distance between the V T and all the RT of the bank. As in section 5.2,the distance (d), typically defined during the geometric generation of a bank, is related tothe match (M ) as M = 1−d2. The closest RT is identified showing a distance dmin, thenthe actual match, equation 6.4, is computed for all the RT s that show a distance d > 1.3dmin.Finally the RT with the highest match is chosen to be associated to the V T . The associationmust be done for all the frequency band considering the distance defined in that band.Furthermore the evaluation of the match, equation 6.4, must be done in the same conditionas in the event search, with the same FFT length and sampling, which requires that the V Tmust be re-sampled to be matched with the RTLFs and it must be truncated, to fit the FFT ofthe RTHFs.

6.1.2 Hierarchical search

The output of the matched filter is analyzed in an iterative hierarchical way. The SNR timeseries (ρi(t)) for each frequency band is computed as follows:

SNRi(t) =

√(h;RTi)

2P (t)+(h;RTi)

2Q (t)

(RTi;RTi)

i = 1,Nbands

(6.5)

where RTi is the real template associated with the V T on the band i. The global SNR of thefull band will be:

SNR(t) =

√(h;V T )2

P (t)+(h;V T )2Q (t)

(V T ;V T )(6.6)

For each V T the average fraction (αi) of SNR2 expected in each band i is estimated as:

αi =(RTi;RTi)

(V T ;V T )

i = 1,Nbands

(6.7)

6.2 The spin template banks case 79

The αi are computed during the association phase for the RT s associated with the V T in thedifferent bands. A global threshold for the full SNR is set to the value given in input. Theindividual threshold on the frequency band i is computed as the global threshold multipliedby the square root of the fraction αi of SNR2 expected on that band. The factor αi is chosenfor each band i as the minimum of the αi over all the V T associated with the same RTi. Thefilter is performed independently on each band and when the SNR of a band exceed theindividual threshold, the output recombination with the other band is done (equation 6.3).Then the combined output is analyzed to look for an excess of the global SNR comparedto the global threshold. On figure 6.3 the sharing of SNR among the bands is shown, fora two bands analysis. In this example the global SNR threshold is set to 5, thus in case ofSNR2 equally divided over the two bands (α1 = α2 = 1/2), the individual threshold will be5/√

2. As it is visible in figure 6.3, the sharing of the SNR between the bands is much moreasymmetric if the trigger is produced by background noise than for the triggers produced bysignals software injected in the detector data stream.

6.2 The spin template banks case

As already said in the previous sections, non-spinning template banks can be generated byMbtaRT application using the LAL library functions. To avoid any unnecessary complication,the generation of aligned-spin template bank has not been implemented. This means thatthe banks have to be produced by an external code. Typically pyCBC [94] applications isused to produce banks with a geometric placement, as already discussed on section 5.3. Thebank generation output is a file with the list of the four physical parameters of the templates:the masses of the two objects and the component of their spins aligned with the orbitalmomentum (the z-components on the source frame).

To run MBTA pipeline Nbands+1 template banks are needed, where Nbands is the numberof the frequency bands of the multi-band analysis: the VT bank for the full frequency range,and one RT bank for each individual frequency band. Although those banks must coverthe same parameter space, they are different because of the different frequency range. Thenumber of the VT bank templates will be larger than the number of the RT banks templatestypically by one order of magnitude.

As discussed on section 6.1.1, making the association among the virtual templates andthe real ones requires information on the metric used in the geometric bank generation. Inparticular, the definition of the distance among templates is necessary and have to be passedto MbtaRT. The pyCBC geometric bank generation performs a coordinate transform fromthe physical parameters (m1,m2,s1z,s2z) to a set of coordinate ξ(1...8) in which the space is

80 MBTA: a CBC low latency pipeline

Fig. 6.3 SNR sharing between the low and high frequency bands for single detector triggers.The line in gray is the combined SNR cut set to 5. The black lines refer to the individualband thresholds. The blue dots are for injections which usually share uniformly the SNRin each bands, while for noise triggers (red dots), the SNR sharing between bands could bemuch more asymmetric.

6.2 The spin template banks case 81

locally euclidean. To simplify the implementation, a C function named MbChisFinder isadded to the MBTA applications, transforming the physical parameters in euclidean spacecoordinates, in the same way of pyCBC. MbChisFinder adds the coordinates ξ (1...8) tothe list of the template parameters. Since the distances among V T s and RT s have to beevaluated into Nbands different metric spaces defined on the different frequency bands, theVT coordinates ξ must be computed for each band. In short, for Nbands = 2, five banks haveto be prepared:

• RTLF bank with masses, spins and related ξ computed in the low frequency band;

• RTHF bank with masses, spins and related ξ computed in the high frequency band;

• V T bank with masses, spins and related ξ computed in the full frequency band;

• V T bank with masses, spins and related ξ computed in the high frequency band;

• V T bank with masses, spins and related ξ computed in the low frequency band;

Once loaded these banks, the measure of the distance can be easily evaluate:

d =

√√√√ 8

∑i=1

(ξV T (i)−ξRT (i))2 (6.8)

Actually the V T coordinates ξ in the full frequency band metric space, are not necessary forthe association, but they can be used to quickly measure the distance between the triggeredtemplate in the output of different detectors. Indeed, the similarity, or even the equality, ofthe triggers coming from the detectors, as well as the time coincidence, can be used for thecandidates selection.

The values of ξ of the V T s are computed by pyCBC only for the full band metric spacewhere the template placement is done. For this reason MbChisFinder has been implementedin the MBTA package, even if it consists essentially in a C-language translation of the pythoncoded analogue pyCBC function.

The waveform generator function have to be suitable to handle the spins. Currently theTaylorT4 approximant and the LAL library functions are used to produce waveform on timedomain. Up to now, the 3.5PN approximation order is the only supported.

6.2.1 Splitting the large banks

Since MBTA is a low latency pipeline, the computing time is crucial. The template numberof an aligned-spin bank is tipically ten time larger than the number of template of a non-

82 MBTA: a CBC low latency pipeline

spinning bank. This makes the matched filter procedure, much more time expensive. Anyhowthe matched filter can be easily parallelized, splitting the template banks and running anindependent trigger search for each sub-bank. In MBTA the parallelization is not completelystraightforward, because the search is performed in the same time in different banks and thenthe outputs must be recombined. Hence, each parallel search is done on one virtual sub-bankand Nbands real sub-banks. The templates which have to be recombined must be containedin the same set of sub-banks running on the same search. The strategy to properly split thebanks is the following:

• the V T −RT association is done with the whole virtual and real banks;

• the virtual template bank is split in a chosen number of sub-banks;

• for each VT sub-bank the relative RTLF/HF sub-banks are created selecting the RTsassociated with the VTs of this sub-bank.

At this point, each search can run with a sub-bank of virtual templates and the relative realtemplate sub-banks.

Different strategies can be adopted to group the virtual templates to optimized the latencyand the computing requirements. Since different virtual template can be associated with thesame real template, the filter of a single real template can be duplicated in different parallelsearches. This effect can be large if the VT sub-banks are built randomly choosing the virtualtemplates over the parameter space. To minimize the duplication of the real templates, thevirtual template can be collected geometrically dividing the parameter space in regions. Inthis case for most of the times, the virtual templates associated with the same real template,will lie in the same set. Another possible choice is to gather the virtual templates accordinglywith their duration. This is useful because the length of the FFTs used in a search is taken tobe as long as the duration of the longest template. This means that if a search is performedover templates with very different duration a useless long FFT is used to filter the signalwith short templates. The two choices give similar results because templates that are close induration will be not so far on the parameter space.

6.3 The χ2 cut

The matched filtering is a very effective technique to recognize the desired signal buried innoise. Nevertheless the matched filtering output can be driven at large values also by spuriousnoise, since high amplitude glitches can contain enough power to exceed the SNR thresholdeven with low values of the match. A properly defined χ2 can be used to better discriminate

6.4 The matched filter output shape cut 83

"good" signals from noise, mitigating this problem. Within the CBC community the χ2

definition is done dividing the frequency band into several smaller bands and looking to thediscrepancy between what is found and what is expected from signal for each band[95]:

χ2 =

Nbands

∑i=1

[(h;RTi)P − (h;V T )P αi cos(φ0 +∆φi−1)]2+

Nbands

∑i=1

[(h;RTi)Q − (h;V T )Q αi cos(φ0 +∆φi−1)

]2(6.9)

where αi is the fraction of the expected SNR2 of the band i (equation 6.7) and ∆φ1..(Nbands−1)

are the phase shifts of the signal at the frequency separation between the bands i and i−1with ∆φ0 = 0.

Although the number of degrees of freedom of this χ2 definition is limited to the smallnumber of frequency band (typically two), the evaluation of the equation 6.9 comes for free,being all the equation parameters computed during the multi band output recombination.Taking the signal as a sum of a perfectly matched waveform plus a Gaussian noise it canbe shown ([95]) that the quantity defined in 6.9 follows a χ2-probability distribution withan expected value equal to 2(Nbands −1). Taking into account the approximation due to theanalytical post-Newtonian expansion, as well the discrete spanning of the parameter space,the match between the templates and the expected signal is not perfect. This mismatch makesthe probability distribution a non-central χ2-distribution with a non-centrality parameterB ·SNR2, where B depends on the mismatch and on the number of the degrees of freedom.

In the MBTA pipeline a threshold following this trend has been chosen:

χ2 < A

(2+B SNR2) (6.10)

The parameters A and B are assessed in an empirical way and typically are set to 3 and 0.025.On figure 6.4 the distributions of the χ2 values of events from H1 and L1 detectors, withrespect the SNR are shown. It is visible that the tighter cut throws away a large number oftriggers coming by the noise background, while all the triggers linked to injected signals,survive the cut.

6.4 The matched filter output shape cut

The response of a matched filter for an input signal perfectly matching the template is anarrow peak. Since the templates is only an approximation of the expected signal on adiscrete lattice of the parameter space, the filter output will be broader but still gathered

84 MBTA: a CBC low latency pipeline

Fig. 6.4 Distribution of the χ2 value. The little dots are noise triggers, the line is the typicalthreshold curve define in equation 6.10 with A = 3 and B = 0.025. The circles with crossesare triggers associated to the injections. It can be noticed that all the injections meet the cut.

6.5 Searching the coincidences 85

Fig. 6.5 SNR2 time series for a software injection of SNR = 5 in presence of noise. The totaltime is 0.1s. The brackets are illustrating the time windows on which the SNR2 is integrated.

around the peak. In case of triggers produced by spurious non Gaussian glitches, the outputis expected to present a much wider distribution with multiple maximums in addition tothe main peak [96] [97] [98]. This feature can be used to discriminate signals from noise,defining a quantity linked to the shape of the matched filtering output (MFO). This quantityhas been defined as the ratio R of the MFO integrated in a region close to the maximum (thecentral part), over the integral computed around this region (side part). The figure 6.5 showsgraphically the definition of the MFO ratio R. The length of the considered SNR2 time serieis typically set to 0.1 s, the central part is taken to be the 15 ms around the main peak. Thedifference in MFO shape is quite evident for a trigger produced by a signal injected into thedata stream on the plot, figure 6.6a, and a spurious glitch, figure 6.6b.

In figure 6.7 the value of the MFO ratio versus the SNR is shown together with the curveused to discriminate signals from noise glitches and used as a cut to veto the backgroundnoise triggers:

Rmax =A

SNR2 +B (6.11)

The quantity A and B are typically set to 65 and 0.4 respectively and they are empiricallyassessed, where the curve decreasing with SNR2 because the higher is the signal the morethe matched filtering output approaches the ideal single narrow peak behavior.

6.5 Searching the coincidences

The non-Gaussian detector noise background can produce loud triggers which survive toall the applied cuts and the vetoes. The time coincident search between different detectors,placed in different site around the world, can cut down the rate of false alarm rising up theconfidence that a trigger is a real gravitational wave event. Considering the three detectorswhich should run on the next observational runs (the two LIGOs and Virgo) only the triggers

86 MBTA: a CBC low latency pipeline

(a) Matched filtering output for a signal injected inthe noise stream of the detectors. The plot shows anarrow distribution of the signal energy around thecentral peak.

(b) Matched filtering output for a background trig-ger. Most of the signal energy is spread in a widedistribution.

6.5 Searching the coincidences 87

Fig. 6.7 Distribution of the matched filtering output ratio. The red dots are the backgroundnoise while the purple line is the threshold curve in equation 6.11 with A = 65 and B = 0.4.In black the triggers related to gravitational wave signal injections.

88 MBTA: a CBC low latency pipeline

found in at least two detectors within a defined time window are not discarded. The timewindows are chosen to account for maximum time of flight between the site and the timeresolution of the matched filtering search. Typical values of the time window are 15 ms fortriggers from LIGO Hanford and Livingstone and 40 ms for the triggers from the LIGOdetectors and Virgo. Beside reducing the background, finding triggers from multiple detectorprovide information to extract the gravitational wave source position over the sky.

Nevertheless, if the search is extended for weak events setting a low threshold on theminimum SNR, the number of triggers produced by the noise is so high that many of thatsurvived also to the time coincidence cut. A further cut can be performed on the base of aconsistency test among the templates which trigger in time coincidence in different detectors.The consistency test consists on looking to how much the triggered templates are similar,where the similarity is quantified by the match, and defining a veto as a threshold on thisvalue. To avoid the waste of computing time, the template consistency test can be speedup estimating the match with the distance on the metric space used to produce the bank.Running with non-spinning templates, a chirp mass consistency test is enough to evaluatethe distance between templates. Hence, only the triggers in time coincidence on at leasttwo detectors, with an absolute difference on the chirp masses of the triggered templateslower than a threshold, are considered. Using aligned-spin template banks the mass chirp isreplaced by the euclidean distance in the ξ coordinates space (equation 6.8).

A more drastic solution is to consider only the events, within the time coincident window,that trigger the same template on different detectors. This procedure, named exact match, hasshown to be the most effective solution in order to have, at the the same level of detectionefficiency, the best performance on rejecting the noise events. The price is an incrementof computing time and the amount of temporary data to storage, because with the exactmatch it is not possible to make a clustering of the triggers before applying the coincidencetest. Indeed the clustering procedure, explained in more details in section 6.6, groups singledetector cluster of triggers which are very close in time. For each cluster, typically producedby an event which triggers many different templates, a single trigger is considered, choosingthe highest in SNR among all the triggers of the cluster. The clustering is a sort of datacompression saving computing time and memory, but since the secondary triggers of thecluster are discarded, the efficiency can be cut down searching for exact match betweentriggered templates. For this reason the clustering can be done only after the coincidencesearch. Another drawback of the exact match option is the need to run for all detector withthe same template banks, which must be produced using the sensitivity curve of the mostsensible detector.

6.6 Clustering the triggers 89

Fig. 6.8 Event clustering procedure. Triggers are displayed like vertical bars proportional tothe event amplitude. The twelve input events are turned to three clusters by the clusteringalgorithm with a minimum time distance between two cluster lager than minGap.

MbCoinc function is the MBTA application which performs the coincidence searches.Another duty accomplished by this function is the significance estimation of a found eventin terms of the inverse of the false alarm rate (IFAR). Details on the IFAR and the eventsignificance are treated in section 6.8.

6.6 Clustering the triggers

A noise glitch or a "real" signal, coming from a GW source or which has been artificiallyinjected, usually involves many more than only one template. Hence most of the triggersproduced by the match filtering appears to be gathered in time clusters where each cluster iscaused by one single event. The clustering algorithm compresses the data output mergingtogether all the triggers that are separated in time less than a given minimum gap value,typically of the order of tens of milliseconds. Performing the clustering right after the filteringreduces the memory occupation and allows the coincidence search process to run only overthe clusters. As already stated in section 6.5, this is not the case when the exact match optionis applied. Indeed, the clustering data compression yields to loose the information on thesingle template triggers, cutting down the efficiency of the multiple detector time coincidencesearch with exact match between templates. Nevertheless, when the exact match is applied,the clustering can be done after the coincidence search. In this case, since the backgroundhas been already lowered by the coincidence search, the minimum gap can be chosen largerby one order of magnitude: the default value is 0.5 s. In figure 6.8 the clustering process issketched. For each cluster of triggers only the one with the higher SNR is chosen and allits parameters are kept to be representative of the clustered event, including vector showing

90 MBTA: a CBC low latency pipeline

the match filtering output. Only a few parameter, as the SNR, the arrival time, the templatemasses and spins of the other triggers absorbed in the cluster are stored. The two parameterstime-after and time-before give the dimension of the cluster on the time scale and the minGapis the chosen minimum time distance between two different clusters.

The MBTA application performing the clustering is named MbClustering. This functionis also used to rearrange the triggers over the frames. The detectors data streams are infact organized and handled in frames. Due to the length of the FFT s, the matched filteringtrigger arrives up to tens of seconds after the corresponding data chunk has been provided. Ifthe frame size used by the filtering process (MbtaRT) is smaller (typically one second), thetrigger is stored in a frame with a timespan which does not include the trigger time. Dealingwith events across the frame boundaries the clustering process has to create buffers of triggersand, if the buffer length (which is a parameter of the MbClustering process) is long enough,the triggers are moved in the right frame. Since the coincidence search handles data frame byframe, it is mandatory to rearrange the triggers before running MbCoinc even if the exactmatch is applied. In this case MbClustering can run with a minimum gap parameter set tozero, to move the trigger in the right frame without doing any clustering.

6.7 Data quality: science segments and vetoed time peri-ods

When a gravitational wave detector is locked and in operation, data are recorded and labeledas "science" data. The estimation of the science data quality is crucial to make scientificsentences about detection or upper limits. Triggers in detection pipelines can be producedby a series of disturbances which can be grouped in two broad categories: instrumental andenvironmental noises. Thanks to the works on the detectors characterization many of thesources of this transient noise, can be identified. Analyzing the detector behavior, a qualityflag can be assigned to the data in presence of some noise source and spurious glitches canbe vetoed to reduce the false alarm rate of the searches. The most common instrumentstransient noise sources are linked to some analog or digital part of the control system or laserlight fluctuations, due to brief mirror misalignments or light scattering. The environmentalnoise can be due, for instance, to electromagnetic disturbances or to seismic activity that canbe up-conversed by some non-linear coupling up to the detection band. A huge quantity ofauxiliary channels, which control both the detector systems and the environmental behavior,can help to locate problems and identify specific issues producing vetoed intervals. Thesearch of transient or strange features on the auxiliary channels is used both for producing

6.7 Data quality: science segments and vetoed time periods 91

vetoes and for finding clues about not understood problems. Defining vetoes has the goal toreduce the number of false alarms in order to increase the confidence of gravitational wavecandidates. Their effectiveness essentially depends on three parameters[51]:

• the efficiency (E) defined as the percentage of the triggers which are vetoed;

• the deadtime (D) that is the percentage of the science mode time which is vetoed;

• the used percentage (U) defined as the percentage of the vetoed intervals which containat least one trigger.

Two figures of merit can be defined using the upper parameters:

• the ratio between the efficiency and the deadtime: RED = ED , that is large for effective

vetoes;

• the expected used percentage obtained comparing the used percentage for a veto towhat is expected to be if the vetoed intervals are not correlated with the triggers. Theexpected used percentage is defined as U over the average trigger rate and the typicallength of the veto intervals, RU = Nwt/Nt

Tv/T where Nwt is the number of interval with atleas one trigger, Nt is the number of triggers, Tv is the vetoed time and T is the totalscience time. Effective vetoes have large expected used percentage, that indicates thatvetoed triggers are more than what is expected by random chance.

The safety of a veto is evaluated by their correlation with the hardware injections, thatare artificial signals injected into the detector control system to simulate the gravitationalwave effects commonly used to test the searching pipeline. Vetoes are classified into fourcategories accordingly to the level of the problem understanding. Well understood vetoeshave a large efficiency to reduce the false alarms and low probability to accidentally vetoinga gravitational wave signal. Decreasing the understanding, the vetoes could be useful todecrease the background but also can be less safe in dismissing gravitational waves. Fourcategories are used in decreasing order of understanding:

• Category 1 (CAT1): the detector are taking data not in the design configuration. Usuallythis interval are mistakenly marked as science interval and have not to be analyzed bythe searching pipeline.

• Category 2 (CAT2): well understood noise in confined time with a known couplingwith the gravitational wave channel. Data in CAT2 intervals are analyzed but events inthose intervals are discarded.

92 MBTA: a CBC low latency pipeline

• Category 3 (CAT3): vetoed intervals correlated with transients in auxiliary channelsbut without a full knowledge of the coupling mechanism. Usually the figures of meritare worser than CAT2 vetoes. CAT3 vetoes include also the periods of time just beforea loss of lock in the interferometer, because triggers in those intervals can be probablydue to the instabilities that contributes to the lock loss.

• Category 4 (CAT4): vetoes, with low figure of merit, related with something thathappens in the detector, or in the environment around it, without a clear possibility topollute the data. The issue must be analyzed in more details during the follow-up of apossible detection.

6.8 The false alarm rate estimation

A level of significance is associated to each candidate event by the coincidence search. Thissignificance is evaluated as the inverse of the false alarm rate (IFAR = 1/FAR), being theFAR the rate of random coincidences with SNR equal or greater than the given candidate.The FAR computation is based on the assumption that triggers produced by the noise indifferent detectors are totally uncorrelated. In this case the false alarm rate for a coincidenceis simply the product of the probability to have a coincidence among the detectors andthe rates of the detectors triggers. The rate Ri in the detector i is computed in a givenamount of the last data time Ti, just as the number of triggers Ni during this time overthe time: Ri =

NiTi

. The probability of random coincidences between two detectors can beevaluated multiplying the rate of one detector by the rate of the other one and the timewindow allowed for the coincidence dti j: Ri j = RiR jdti j. For three detectors the probabilitywill be: Ri jk = RiR jdti jRkdtik. Computing the FAR, the probability of random coincidencesmust be multiplied by the probability that a combination of triggers fulfilling the followingconditions:

• present a combined SNR larger than the given candidate; this means that the combinedSNR is the parameter used to rank the trigger significance locally inside the timeinterval;

• meet the template consistency condition, typically as the exact match is applied onlythe combination between equal templates are counted.

This probability is computed just counting the combinations that meet the conditions dividedby the total number of combinations M

Ni·N j(or M

Ni·N j·Nkfor three detectors). Multiplying this

6.8 The false alarm rate estimation 93

probability by the coincidence rates, the FAR is obtained:

FAR = Ri jM

NiN j=

Mdti j

TiTjfor double coincidence

FAR = Ri jkM

NiN jNk=

Mdti jdtikTiTjTk

for triple coincidence(6.12)

It is evident from equation 6.12 that the larger is the time intervals Ti and Tj the lower will bethe minimum estimable FAR. For example, time intervals of the order of hours are needed toestimate a FAR, in the two detectors case, of the order of one per hundreds of years. Thisdefinition of the false alarm rate is local, taking into account only the last period of data,and it can accomplish the non-stationarity of the background. Enlarging the time intervalspermits a more accurate FAR computation and a lower value of the minimum estimable FAR,but the computation starts to be less local. Hence, a trade off have to be found.

The FAR calculation is performed by MbCoinc process with a typical time interval of3 hours to have a minimum estimable inverse false alarm rate IFAR of 300 years for thedouble coincidences. Information on vetoed intervals have to be given to MbCoinc to excludetriggers in those time segments. If injections are present on the data, some seconds aroundthe injection time are excluded in collecting the triggers for the FAR evaluation. Avoiding tohave a bias in the FAR due to "satellite" triggers of a loud candidate, some seconds before anevent can be excluded as well.

When running multiple searches the FAR computed in equation 6.12 have to be multipliedby the trial factor. If the data come from three detectors as in the case of a network searchesby Hanford (H) and Livingstone (L) LIGOs and Virgo (V), four possible coincidence searchesare performed: the three doubles HL, HV, LV and the triple HLV. Hence in this case thecomputed FAR has to be multiplied by four.

Chapter 7

A search for BNS with spin

The search of gravitational waves from compact binary coalescences is performed by meansof pipelines based on matched filter techniques. To make this technique effective, accuratetemplates of the expected signals have to be computed and a bank of templates must beproduced filling the space of relevant physical parameters of the binary systems. In chapter 4the effect of introducing the spin of the objects on the waveform equation has been discussed.The possibility to build affordable template bank spanning the four-dimensional spaceidentified by the masses and the spins has been addressed in the case of spin aligned with theorbital momentum. A test of the effectiveness of aligned-spin template banks compared tothe non-spinning template bank has been performed using the MBTA pipeline to simulatea real gravitational wave search of signals from binary neutron stars systems (BNS). Thedetails and the results of the search is the content of this chapter.

7.1 The analyzed data

To obtain an effective comparison between non-spinning and aligned-spin template banksearches, 8 ·105 seconds of Mock Data Challenge (MDC) have been analyzed. The MDCdata are a common set of data used to compare output of different search pipelines under sameconditions, to investigate the pipeline configurations tuning them under different assumptionson astrophysical populations.

The used data set named early recolored, consists on the noise data from the last detectorruns (S6 for LIGO and VSR2 for Virgo), rescaled with the expected sensitivity curve of theadvanced detectors in their early configurations, in short a simulation of the expected noiseduring the first advanced detector observational run (O1). Since the O1 run foresees thepresence of only the two LIGO interferometers, our analysis considers Hanford (H1) andLivingstone (L1) detector data.

96 A search for BNS with spin

The injection set used for this search was prepared to test different pipelines on the sameplayground data. The injection set parameters are characteristic of the BNS systems:

• mass range: from 1 to 3 M⊙, uniformly distributed;

• aligned spins with magnitude up to 0.4.

The simulated systems have a uniform distribution in chirp distance. To obtain statisticallysignificant results analyzing a reasonably short period of data, six injection sets have beenused with injections at an interval of about 840 s, meaning an injection every 140 s using allthe six sets together. The waveform are generated with SpinTaylorT4 approximation up to3.5PN with a low frequency cut off of 20Hz.

7.2 The template banks

The aligned-spin template banks have been produced using a pyCBC toolkit [88] application,pycbc_geom_aligned_bank. This application uses the technique described in section 5.3 withgeometric placement of the templates covering a four-dimensional parameter space definedby the physical parameters: the two masses and the two aligned spins. Since the analysis aimsto check the performance of the aligned-spin bank in BNS parameter region, the star masseswere limited from one to three solar masses with spin up to 0.4. A 3.5PN SpinTaylorF2approximation is used to generate the template waveforms. The physical parameter space canbe transformed in a two dimensional globally flat euclidean space in the coordinate (ξ1,ξ2).In this coordinate set, a grid of points is produced to guarantee a minimum fitting factorof 0.97. A coverage test has been done by mean of the pyCBC_banksim, a pyCBC toolkitapplication that computes the fitting factor of a number of injected signals evenly spread overthe physical parameter space with a given template bank. The results, in figure 7.1a, showthat, considering injection with aligned spin, the bank is successfully covered since almost allthe fitting factors exceed the 0.97 value. For signals without a preferred direction (precessingspins) in figure 7.1b, the fitting factor is larger than 0.97 for about the 90% of the injections.

A two-bands MBTA search has been performed, with low frequency cut off equal toflow = 30 Hz, a separation frequency fc = 137 Hz and a maximum frequency fupper =

2048Hz. Choosing flow = 30 Hz instead of 10Hz permits to loose less than 1% of the SNR2

having the advantage to handle much shorter signal. The value of fc balances the signalenergy (SNR2) over the two bands. The fupper is set at the Nyquist frequency.

Three banks have been produced:

• the virtual template (V T ) bank for the full band [30 - 2048] Hz;

7.2 The template banks 97

(a) (b)

Fig. 7.1 For each upper limit on FF the fraction of injected signals with a FF lower thanthat upper limit is shown. The bank covers the BNS parameter space with aligned spin up to±0.4. The plot reports the results for 104 injected signals. On the left, plot(a), the signalshave aligned spins, on the right, plot(b), the spins are isotropically oriented.

• the low frequency real template (RTLF ) bank for the band [30 - 137] Hz;

• the low frequency real template (RTHF ) bank for the band [137 - 2048] Hz.

According to what is described in section 6.2, the ξ coordinates of the bank points havebeen computed to reduce the time needed to associate virtual and real template. Eventually,five banks are produced: three virtual banks with the same physical parameters (masses andspins), but different ξi computed for the full, the low and the high frequency band and the tworeal banks in the low and the high frequency band. A summary of the banks characteristicis reported in table 7.1, where the fitting factor FF corresponds to the match between thevirtual templates and the associated real templates.

Bank N. points flow [Hz] fup [Hz] SNR2 mean FF

VT 79324 30 2048 1 -RTLF 9436 30 137 0.55 0.986RTHF 2707 137 2048 0.45 0.984

Table 7.1 For the virtual and the real aligned-spin banks the number of templates, the lowfrequency cutoff, the maximum frequency, the fraction of SNR2 with respect to the full bandand the fitting factor of the V T s over the RT banks, are reported.

The distribution of points in a m1, m2 plot is shown in figure 7.2 for the virtual bank andthe real banks. By convention, only the area where m1 > m2 is covered, while obviouslythe systems are symmetric with respect to the bisector. The distribution in chirp mass (Mc)versus effective spin (χ) is shown in figure 7.3, being the effective spin a weighted sum of

98 A search for BNS with spin

Fig. 7.2 Masses distribution for the aligned-spin template bank, for the virtual template bankon the left, the low frequency real template bank in the middle and the high frequency realtemplate bank on the right.

Fig. 7.3 Distribution of chirp masses and effective spin for the aligned-spin template banks,for the V T bank on the left, the RTLF bank in the middle and the RTHF bank on the right.

the aligned spins of the two objects:

χ =m1χ1 +m2χ2

m1 +m2(7.1)

Finally the figure 7.4 shows the point placement over the ξ1, ξ2 space. The regular grid withhexagonal point distribution is visible in the flat Cartesian space.

The non-spinning template banks has been produced as well with pyCBC tollkit. Thegeneration technique is similar to the aligned-spin placement procedure, but in this case, anearly one-dimensional space arises from the principal component analysis, with a spread ofthe second coordinate very tiny compared to the first one: about 1/1000 in unity of templatesize. In figure 7.5 the distribution of the masses of the virtual and real non-spinning banksare shown. In table 7.2 the details of the non-spinning banks are gathered. It can be noticedthat the aligned-spin virtual bank contains ten times more templates than the non-spinningvirtual bank. This means that the growth in efficiency, which is expected using aligned-spintemplates, has to deal with a larger number of false alarms.

7.2 The template banks 99

Fig. 7.4 A portion of the parameter space in ξ coordinates is magnified, with equal data unitlengths along each axis. The hexagonal placement of the templates in ξ coordinate space isvisible as a uniform distribution.

Fig. 7.5 Masses distribution for the non-spinning template bank, for the virtual template bankon the left, the low frequency real template bank in the middle and the high frequency realtemplate bank on the right.

Bank N. points flow [Hz] fup [Hz] SNR2 mean FF

VT 8026 30 2048 1 -RTLF 1459 30 131 0.49 0.981RTHF 483 131 2048 0.51 0.989

Table 7.2 For the virtual and the real non-spinning banks the number of templates, the lowfrequency cutoff, the maximum frequency, the fraction of SNR2 with respect to the full bandand the fitting factor of the V T s over the RT banks, are reported.

100 A search for BNS with spin

After the banks production and loading the VT-RT association file is generated. Foreach virtual template, this file contains all the parameters needed to reconstruct the matchedfiltering output over the full band:

• the identification numbers of the low frequency and the high frequency real templateclosest to the virtual one;

• the fraction of SNR2 of each band with respect to the full band;

• the time delay and the phase shift to recombine the waveform according to equation6.3.

This association step is time consuming, although it is performed just at the beginningduring the so called initialization phase. For about 104 templates, as in the non-spinningbanks case, this phase takes few hours, but for 105 templates, as in the case of aligned-spinbanks, the duration can be larger than one day being the computing complexity more thanlinear. Computing times are estimated running the initialization process on Xeon E5-2650 @2.60GHz CPU.

Once the association file has been produced the analysis can start.

7.3 Bank splitting

The matched filtering procedure is intrinsically parallelizable, splitting the bank in severalsub-banks which can be used to obtain completely independent filtering processes which canrun in multi-core machines. The more is the processes we are going to use, depending onthe number of available cores, the faster will be the analysis. For the non-spinning templatebanks the number of templates is quite small and the filtering can be performed with thewhole bank in a single process. Conversely, for the aligned template banks, the V T bank hasbeen split in 50 sub-banks by the MBTA feature described in section 6.2.1. The sorting onξ1 coordinate has been chosen to split the bank, minimizing the number of real templatesassociated with virtual templates belonging to different sub-banks and therefore reducingthe total number of real templates over all the sub-banks. In figure 7.6 the distribution of thepoints in different virtual sub-banks are visible with different colors in chirp mass versuseffective spin plane.

In table 7.3 the mean numbers of templates in the sub-banks are listed. Eventually, 50independent MbtaRT jobs can run in parallel, equipped with a virtual sub-bank and two realsub-banks for the high and the low frequency bands, with the real templates associated withthe virtual templates of the virtual sub-bank.

7.4 Running MBTA on the GRID 101

Fig. 7.6 Distribution of the V T s of the virtual sub-banks, in a Mc-χe f f plane. Templatesbelonging to the same sub-bank are reported with the same color.

NV T NRT LF NRT HF

Mean 1587 211 74Total 79324 10533 3710

Table 7.3 The average and the total number of templates of the aligned-spin split banks. Itcan be noticed that the total number of RT s is slightly larger than the original RT banks,because V T s, belonging to different sub-banks, can be associated with the same RT . For thisreason the RT will be duplicated in all those different RT sub-banks.

7.4 Running MBTA on the GRID

The 106 seconds long stretch of data to be analyzed has been divided in 10 chunks 105

long. This means that, for a two detectors search, the number of job running in parallel is6000: 2 (detectors) times 10 (data chunks) times 50 (sub-banks) repeated for the 6 sets ofinjections. To fully exploit the parallelization potentiality, the CNAF Grid infrastructurehas been used. CNAF is the National Center of INFN for the information technologiesapplied to high-energy physics, involved in management and developing of different Gridinfrastructures, forming a big distributed computing system. The gLite software is used tointerface the access point where the user is connected to the distributed resources of theGrid, to submit jobs, check the status of the submission and retrieve the information on therunning jobs and the final reports. When a job is submitted it is passed to a remote computing

102 A search for BNS with spin

element, in a non-interactive node. The job must be self-contained, hence all the componentsof that particular portion of the analysis must be copied on the computing element:

• the input data chunk with the detector noise that have to be analyzed by that particularjob;

• the injections data chunk;

• the files containing information on data quality;

• the filtering application MbtaRT statically compiled, to make no longer necessary tokeep the library files that the program references;

• the configuration file;

• the virtual and real template sub-banks.

When the job is successfully finished the output files have to be copied back. To transfer thedata from the user access point to the storage elements, the LCG Data Management ClientTools (lcg Util) have been used.

After the successful end of all the filtering jobs, the outputs files which refer to the sametime chunk and the same injection set are merged in 60 filtered data files containing the datafrom the all the sub-banks and the two detectors filtering. The merging processes run inparallel on the Grid as well using a statically compiled version of the Frame Library utilityFrCopy.

The coincidence searches run in parallel for all the 60 filtered files with a static version ofMbCoinc application giving 60 output files containing the triggers found in time coincidencein H1 and L1 detector, for the ten data chunks and the six injection sets.

Finally the clustering step is applied to gather the neighboring triggers running in parallel60 MbClustering jobs.

7.5 Filtering step and the single detector triggers

During the matched filtering, when the SNR in equation 6.6 given by the real templatesrecombination exceed the threshold, a trigger is stored. The parameter associated to thetrigger are referred to the virtual template associated to those real templates. In this search theSNR threshold is set to 5.0. Other selections are done to reject non-Gaussian glitches. Thefirst one is based on the χ2 estimation as explained in section 6.3, where the threshold followsthe equation 6.10: χ2 < 4

(2+0.025 ·SNR2). The second is based on the matched filtering

7.6 The coincident search 103

output shape ratio R (section 6.4) as in equation 6.11 and was applied only for the spintemplate search. With the chosen parameter the threshold function becomes: R < 65

SNR2 +0.8.The two cuts was tightened doing the coincidence between triggers, as described later.

7.6 The coincident search

After filtering the data from the two LIGOs detectors a coincident search is done as explainedin section 6.5. The time coincidence is checked considering both the time of flight of thewave between the two detectors and the uncertainty on the arrival time measurement. A timewindow size of 15 ms has been chosen. A new triggers selection was applied to the triggersbefore searching the coincidence, tightening the thresholds of the MFO ratio (R) and the χ2

cuts. The thresholds were chosen to be equal for aligned-spin and non-spinning analysis:

χ2 < 3

(2+0.025 ·SNR2)

R <65

SNR2 +0.4(7.2)

A strong coincident trigger selection was performed applying the exact match cut.To each coincident event, the associated false alarm rate (section 6.8) has been taken as

the significance estimation. Collecting the triggers for the event false alarm rate evaluation,a time window of three hour before the event has been considered in each detector, for aminimum computable FAR of 1/300 years. A time period of 20 s before the event is removedto avoid bias effects due to "satellite" triggers around the main event. To avoid bias due tothe other injections in the three hours data segments, triggers in 20 s wide windows aroundthe injection arrival times, are removed from data during the trigger collections.

The data segments labeled "category 2" have been vetoed and removed from the coinci-dent analysis as well two seconds around the hardware injections that were done during theS6 run.

Finally the coincident events have been clustered choosing a minimum gap of 0.5 samong the clusters.

7.7 Analysis results and comparison

Comparing the results from the two analysis, three main aspects are taken into account: theefficiency of the search, that is basically the number of injections recovered, the accuracy ofthe injection parameters reconstruction, and the performance about the false alarm rate, i.e.the significance associated to the recovered injections.

104 A search for BNS with spin

Fig. 7.7 The injection found by the two analysis over the analyzed time is shown. The greendots refer to injections found by both the searches. The red and the blue dots are injectionsfound respectively only by the non-spinning search and only by the aligned-spin search.

7.7.1 Efficiency comparison

In figure 7.7 the recovered injections are shown, where an injection is considered foundif a coincident event is present within ±100 ms from the injection arrival time. The totalnumber of found injections are 1061 by the aligned-spin analysis and 968 by the non-spinninganalysis. The aligned-spin analysis gains in reconstructed SNR for the injections found byboth analysis (hereafter named common injections), as it is evident in figure 7.8. The averagegain is 6.78%; this is expected because aligned-spin templates bank shows an higher fittingfactor with respect to the spinning injections than the non-spinning template bank. Thisalso explains the gain in the total number of injection found. Some injections just underthe threshold in the non-spinning search are brought up the threshold in aligned-spin searchby the highest fitting factor with the aligned-spin template bank. Anyhow, this has to beconsidered as an average effect, while for some specific injections is randomly possible tohave an higher match with a non-spinning template.

The adaptive process to compute the PSD used performing the matched filtering evalua-tion can makes some differences between the two analysis, mainly because the FFT lengthsare involved and because the lengths are related to the template durations. This affects thecomputed SNR and can produce further fluctuation on found/missed injection close to SNRthreshold.

7.7 Analysis results and comparison 105

Fig. 7.8 Histogram of the differences between the combined SNR of the injections found byaligned-spin and non-spinning searches.

An already known issue of the pipeline is the blind period after a big glitch. When a largeamplitude trigger comes, this affect the PSD adaptive updating process and some time, of theorder of some FFT lengths, is needed to adjust the PSD evaluation and to unblind the search.This issue is important especially in mock data challenge analysis where a large injectedsignal density is chosen to earn significant statistics with reasonably small amount of datatime to analyze. Being this issue related to the FFTs lengths, the behavior is different for thetwo analysis and this explains the some injections with high SNR missed by one of the twoanalysis and found by the other.

In figure 7.9 the decisive chirp distance versus chirp mass missed-found plot is shown.The decisive chirp distance in a double coincidence search is simply the maximum chirpdistance (equation 4.39) between the injection recovered in the two detectors data, being thechirp distance the effective distance normalized for the chirp mass:

Dchirp = Deff

(Mc

Mc

)5/6

(7.3)

where Mc is a fiducial chirp mass that is set to 1.2 M⊙ (the chirp mass of two 1.4 M⊙ starbinary system).

These plots are done removing injections in a period of 20 s at the beginning and at theend of each science mode and "CAT1" veto segments to avoid a bias due to an incomplete

106 A search for BNS with spin

Fig. 7.9 The injection missed and found by the two searches are reported as a function of theeffective distance and the chirp mass. The injections missed by the aligned-spin analysis isshown on the left, the injections found by the non-spinning analysis is on the middle. On theright the injections found by the aligned-spin analysis.

FFT evaluation. The plot shows that the two analysis have similar drop of efficiency from200 Mpc to 300 Mpc and there are no injection missed under 100 Mpc. Some far injectionsfound with decisive chirp distance larger than 1 Gpc are due to random associations of noisetriggers occurred at injection time. The same view is given by figure 7.10, which shows thecumulative number of injections found versus the chirp distance. The gain in the efficiencyof the aligned-spin search appears for chirp distances larger than 100 Mpc.

7.7.2 Accuracy in injection reconstruction

As already said, the aligned-spin search presents a larger recovered SNR than the non-spinning one for the common injections. This is evident in figure 7.8. Figure 7.11 showsthe distribution of the recovered SNR versus the optimal SNR, where the last is the SNRcomputed by a matched filtering with a template identical to the injected signal. It can benoticed that the difference between the optimal and the recovered SNR is smaller for thealigned-spin search. The gain in SNR using the aligned-spin template is better explainedin figure 7.12 where the relative gain in SNR is plotted versus the effective spin χe f f . Inaverage the gain is null for low value of effective spin, i.e.

∣∣χe f f∣∣< 0.1, while it becomes

7.7 Analysis results and comparison 107

Fig. 7.10 The cumulative number of injection found by the two analysis versus the decisivechirp distance limit, is reported. Below 100 Mpc the two analysis found the same number ofinjections, while above this distance the aligned-spin search shows a better efficiency.

Fig. 7.11 The comparison between recovered SNR and optimal SNR is shown for the non-spinning search (on the left) and the aligned-spin search (on the right). The bisector line isalso visible.

108 A search for BNS with spin

Fig. 7.12 The relative difference in recovered SNR is reported as a function of the effectivespin. The solid line represents the line where the gain is null.

relevant evident that for high values of the effective spin, the aligned-spin bank produces alarge recovered SNR, as expected.

For what concern the accuracy in the injection parameter reconstruction, two quantityhave been considered: the chirp mass of the system, which determines the waveform at thefirst order, and the difference between the waveform arrival time in H1 and L1, which isrelevant because it is used to point the gravitational wave source in the sky.

The discrepancies between injection and recovery in arrival times and chirp masseshave been plotted in figure 7.13. It can be noticed that for the arrival time difference theperformance of the two searches are comparable, while there is a large improvement in masschirp reconstruction accuracy using the aligned-spin analysis. In figure 7.14 the comparisonon chirp mass recovery is split for different values of the injection effective spin. An offsetin the recovery values is evident for large values of the effective spin: 0.1 < |χ|< 0.4. Inparticular, for negative effective spin, where negative means anti-aligned with respect to theorbital momentum, the non-spinning analysis tends to overestimate the chirp mass. This isconsistent to the fact that one of the effects of anti-aligned spin is to reduce the duration ofthe waveform and this can be mimic by a non-spinning larger chirp mass waveform. Forpositive aligned effective spin the effect is inverted, the spin makes the waveform longerand it is mimic by a non-spinning smaller mass chirp waveform, hence the mass chirpof the non-spinning search is underestimated. The plot of the mass chirp reconstructionaccuracy versus the optimal combined SNR (optSNR) of the injections (figure 7.15) gives an

7.7 Analysis results and comparison 109

Fig. 7.13 The performance in injection parameter reconstruction, is shown. Two parametersare considered: the discrepancy on the reconstruction of the difference between the arrivaltime in the two detectors (∆δHLt) and the reconstruction of the chirp mass (∆Mc/Mc[%]).The ∆δ t histogram is on the left and the ∆MC histogram is on the right. The red and the bluelines refer to the non-spinning search and the aligned-spin search respectively.

Fig. 7.14 The distribution of the ∆Mc/Mc[%] is reported dividing the injections in threedifferent bins. The first bin (upper plot) contains injections with −0.4 < χe f f <−0.1, thesecond bin (middle plot) has injection with

∣∣χe f f∣∣ < 0.1 and the last bin (lower plot) has

injection with 0.1 < χe f f < 0.4. Red and blue refer to the non-spinning search and thealigned-spin search respectively.

110 A search for BNS with spin

Fig. 7.15 The optimal combined SNR versus the discrepancy in chirp mass recovering isplotted for the injection found by the two analysis. Red squares refer to non-spinning analysiswhile blue dots are for aligned-spin analysis.

estimation of how many injections found are actuay random associations with backgroundtriggers, assuming that in these cases the recovered chirp mass will differ substantially fromthe chirp mass of the injected signal. Considering as "badly" recovered all the events whichhave a mass chirp discrepancy (∆Mc/Mc) larger than 10%, only one injection with optSNRequal to 4.0 is badly recovered. In aligned-spin analysis the badly recovered injection are 13with a mean optSNR equal to 3.7.

7.7.3 The False Alarm Rate and the significance assessment

The significance of a coincident event found by MBTA pipeline is taken to be the value ofthe inverse of the false alarm rate where the FAR is computed as described in section 6.8. Toassess the validity of this significance estimation, it is useful to build a plot where, for anygiven FAR upper limits, the cumulative rate of the foreground events with a FAR smaller thanthe FAR upper limit is reported. The number of foreground events are normalized for theeffective time of the search, that is the total duration of the intersection of the science timesegments between the two detectors, after removing the vetoed chunks. Obviously, injectionsare removed in counting the foreground events to avoid a bias in the measured number offalse alarms. The golden condition is when the points lie in the bisector, meaning that theestimated false alarm rate is in agreement with the number of events actually found. In

7.7 Analysis results and comparison 111

Fig. 7.16 The comparison between the estimated FAR and the measured FAR is shown. Theactual cumulative rate of triggers with a estimated FAR lower than an upper limit is reportedas a function of that upper limit. The bisector line is also visible as a reference. The plot onthe left refers to non-spinning analysis while the left plot is for the aligned-spin analysis.

figure 7.16 the significance assessment plots are reported for the non-spinning and spinninganalysis. For both the analysis, an overestimation of the FAR is visible for high values ofthe FAR (FAR > 10−3 Hz). This because the FAR is computed before the clustering of theevents while, in the computation of the cumulative rate, the clustered events are considered.To estimate the significance of one event, the triggers before the event are collected and theones that overcome the combined SNR of the considered event (let’s call it ρ) are counted. Acluster is characterized by one large main trigger and other smaller satellite triggers, thusthe probability that when the main event overcome ρ , the satellite triggers overcome ρ aswell, decrease when ρ increases. This means that the bias in FAR estimation introducedby the clustering, becomes less important for more significant events. Some procedures toremove this bias are still under test but, since generally a maximum FAR limit of ≈ 10−7 isconsidered to follow the event as potential GW candidate, the actual importance of this biasis close to be negligible.

The saturation for the highest values of the FAR is related to the threshold in SNR, thelower is the SNR cut the higher the FAR for which the curves saturate.

Considering the significance of the recovered injections by the two analysis, figure 7.17shows the distribution of the FAR versus the combined SNR. The two quantities FAR andcSNR are clearly correlated but with some jiggling due to the "locally" FAR estimation. Since

112 A search for BNS with spin

Fig. 7.17 The distribution of the estimated FAR versus the combined SNR is shown. TheFAR is anti-correlated with the SNR, as expected. A plateau at FAR equal to 2.6 ·10−10 isdue to a minimum in estimable FAR. The plot on the left refers to non-spinning analysiswhile the left plot is for the aligned-spin analysis.

the background is estimated in the three hours before the event, the results depend on howmany glitches there are in that period. In other words the non-stationarity and non-gaussianityof the noise have an impact on the FAR computation, yet the more the pipeline is effectivein removing non-gaussian glitches the more the FAR-cSNR correlation will be strong. Theplateau for FAR equal to 2.6 ·10−10 Hz is due to the minimum estimable FAR as explainedin section 6.8, while the black cross marks the events for which the time before the eventavailable for background estimation is shorter than three hours, i.e. when the event is found inthe first three hours of an analysis chunk. An evident difference between the two distributionis the maximum FAR which is the total number of noise triggers. For the spinning analysisthis number is ten times larger than the non-spinning analysis and consequently the FAR forthe faintest recovered injections is ten time higher.

Comparing the FAR values for the common injections, the FAR estimated by the non-spinning analysis is averagely lower than the spinning analysis FAR, as it can be noticedin figure 7.18. Two concurrent effects act on the FAR of the injections: the number of thebackground triggers and the magnitude of the recovered SNR. The average increment of theFAR for spinning analysis means that the gain on the recovered SNR does not compensatethe higher number of noise triggers.

7.8 Final remarks 113

Fig. 7.18 On the right the comparison between the estimated FAR of aligned-spin and non-spinning analysis. The bisector lines is shown as a reference. On the right an histogram of thedistribution of FAR for the two analysis: non-spinning in dash-dot red line and aligned-spinin blue solid line.

This impacts also on the search sensitivity curve, that shows the efficiency versus theFAR as depicted in figure 7.19. Looking to this plot, which evaluate the goodness of asearch, the two analysis performances appear to be similar. There is a little advantage inusing spinning search only for low values of the FAR limits larger than 10−4, while thenon-spinning analysis shows a slightly better behavior for lower FAR limits.

7.8 Final remarks

The performances of MBTA pipeline using the non-spinning template bank and the aligned-spin template bank were compared, using as playground data binary neutron star signalsinjected in recolored early aLIGO noise. The maximum neutron star spin of the injectedsignals was 0.4 that is also the maximum spin ever observed in a neutron star. The searchwas performed in coincidence between the two LIGO detectors for about 106 s of data.The comparison has proved the sanity of the spin template bank implementation in MBTApipeline. As expected using aligned-spin template bank, an evident gain in SNR is obtained(about +7%), especially for injection with high absolute value of spin. Moreover, theaccuracy in the mass chirp reconstruction is remarkable and the offset that is present inthe mass chirp reconstructed by the non-spinning analysis for highly negative and highly

114 A search for BNS with spin

Fig. 7.19 The cumulative number of injections with a FAR lower than an upper limit isreported as a function of that upper limit. This plot is an effective measurement of theperformances of the two searches. For every maximum accepted value of the false alarm rate,the most efficient search can be determined.

positive spin, disappears using the aligned-spin template banks. However, the larger numberof aligned-spin templates produces an increment of the FAR and it makes the non-spinninganalysis more efficient for the commonly acceptable rate of false alarms.

Chapter 8

Conclusions

The first attempts to observe gravitational waves started in the sixties. After roughly fiftyyears, instruments which are sensible enough to detect such extremely faint signals are readyto take data with a significant hope on making the first detection in the near future. Theadvanced LIGO detectors have started its first observational run in September 2015 andadvanced Virgo is expected to attend the second observational run in fall 2016.

This work deals with two aspects of the gravitational wave detection. The first one is thedetector characterization and noise "hunting", that is important for putting the detector inoperation, improving the detector performances and increasing the confidence in a possibledetection. To reach this goal the SILeNTe (System Identification Linear et Non-linearTechnique) tool has been designed and implemented. This tool can be applied to uncovernoise structures finding linear and non linear coupling between the gravitational wave detectoroutput and hundreds of auxiliary channels. The tool is based on a fast system identificationtechnique called Sorted Fast Orthogonal Search (SFOS): a time domain model of a noisestructure, existing in the output channel, is built using linear and non-linear combinations ofthe auxiliary channels; then a rank of the contributions of each channel and combination ofthem is done, possibly providing useful hints about the possible sources involved in that noisefeatures. The tool has been checked on different noise structures which plagued the initialdetector data involving persistent and narrow-bandwidth noise structures, but the tool canbe used also to entangle contributes to transient noise and glitches. Different noise featureshave been analyzed, in particular, modulations of the calibration and the electric power lines.Some of the hints, given by the nature of the involved auxiliary channels, are in agreementwith the results of extensive noise hunting campaigns performed by the commissioning crew.The translation of the code from MATLAB language to C is ongoing (thanks to the AnnecyLAPP Virgo group) to facilitate the integration of the tool with the other tools available for

116 Conclusions

the commissioning crew. Some preliminary tests on glitches seem to be promising and thepossibility to use this tool to create vetoes is currently studied.

The second aspect is the search of signals from compact binary coalescence. The lowlatency search of this kind of signals is considered compelling, because of the possibilityto promptly send alarms to telescopes to search for electromagnetic counterpart of thegravitational wave sources. My work has been focused on an upgrade of MBTA a low latencypipeline based on matched filtering technique, which has been used to search for compactbinary coalescence signals since the initial interferometer runs. The outcome is that MBTAcan now load and handle aligned-spin template banks. Banks of templates that consider spinof the individual objects aligned with the orbital momentum are the standard choice for thecurrent observational run. The possibility to manage the exact match coincidences (triggersof the same template in time coincidence between detectors) has been added, being the exactmatch another standard choice for observational run searches. The performances of spinningtemplates banks has been tested. Signals from binary neutron stars with aligned spins(maximum spin equal to 0.4) have been injected in the initial detector noise recolored withthe advanced detector sensitivity curve. Two searches with no-spinning and with aligned-spintemplate banks have been performed and the results are in agreement with the expectation:the number of injected signals, detected by the analysis with spinning templates, is largerby ∼ 10% and the signal to noise ratio is larger by an average of ∼ 7%. Moreover the mainphysical parameters of the simulated signals, that is the mass chirp (a combination of the twomasses) is more accurately recovered by spinning template analysis. These results representa sanity check of the new MBTA implementations. However, the enlarged size of the bankwith spinning templates increases the rate of the false alarms. Comparing the efficiency ofthe search for a fixed value of the false alarm rate, the non-spinning template banks stillresults the better choice for searching signals from binary neutron stars systems.

The aligned-spin bank implementation in MBTA is nonetheless very important for thesearch of signals from binary black holes and neutron star-black hole systems, since the theblack hole can reach large values of spin. For this reason the MBTA analysis which currentlyrun on-time for the observational run uses aligned-spin template banks.

References

[1] A. Einstein. The Foundation of the General Theory of Relativity. Annalen Phys.,49:769–822, 1916.

[2] B. Schutz. A First Course in General Relativity. Cambridge University Press, 2009.

[3] M. Maggiore. Gravitational Waves. Volume 1: Theory and Experiments. CambridgeUniversity Press, 2007.

[4] I W Harry. Can we hear black holes collide? PhD thesis, School of Physics andAstronomy Cardiff University, 2011.

[5] Dave McKechan. On the use of higher order waveforms in the search for gravitationalwaves emitted by compact binary coalescences. PhD thesis, Cardiff University, 2010.

[6] C D Ott, A Burrows, L Dessart, and E Livne. A New Mechanism for GravitationalWave Emission in Core-Collapse Supernovae. Phys. Rev. Lett., 96(20), 2006.

[7] K Riles. Gravitational waves: Sources, detectors and searches. Prog.Part.Nucl.Phys.,68:1–54, 2013.

[8] S Klimenko, I Yakushin, A Mercer, and G Mitselmakher. Coherent method for detectionof gravitational wave bursts. Class.Quant.Grav., 25(11), 2006.

[9] J Abadie et al. Predictions for the Rates of Compact Binary Coalescences Observableby Ground-based Gravitational-wave Detectors. Class.Quant.Grav., 27(17), 2010.

[10] J Aasi et al. Application of a Hough search for continuous gravitational waves on datafrom the fifth LIGO science run. Class.Quant.Grav, 31(8), 2014.

[11] J Aasi et al. Gravitational waves from known pulsars: results from the initial detectorera. Astrophys. J., 785(2), 2014.

[12] T Regimbau. The astrophysical gravitational wave stochastic background. Res. Astron.Astrophys., 11(4), 2011.

[13] L P Grishchuk. Amplification of gravitational waves in an isotropic universe. Sov. J.Exp. Theor. Phys., 40:409, 1975.

[14] L P Grishchuk. Gravitational wave astronomy: in anticipation of first sources to bedetected. Phys. Usp., 44(1), 2001.

118 References

[15] B Allen. The Stochastic gravity wave background: Sources and detection. In CambridgeUniv. Pr. Cambridge, editor, Proceedings of: Les Houches 1995, Relativistic gravitationand gravitational radiation, 1997.

[16] R A Hulse and J H Taylor. Discovery of a pulsar in a binary system. Astrophys. J.,195:L51–L53, 1975.

[17] J H Taylor, L A Fowler, and P M McCulloch. Measurements of general relativisticeffects in the binary pulsar PSR1913 + 16. Nature, 277:437 – 440, 1979.

[18] J M Weisberg, D J Nice, and J H Taylor. Timing Measurements of the RelativisticBinary Pulsar PSR B1913+16. Astrophys. J., 722(2), 2010.

[19] J Weber. Detection and Generation of Gravitational Waves. Phys. Rev., 117(306), 1960.

[20] B Abbot et al. LIGO: the Laser Interferometer Gravitational-Wave Observatory. Rep.Prog. Phys., 72(7), 2009.

[21] T Accadia et al. Virgo: a laser interferometer to detect gravitational waves. JINST, 7,2012.

[22] H Grote et al. Status of GEO600. Class.Quant.Grav., 27(8), 2010.

[23] K Arai et al. Recent progress of TAMA300. JPCS, 120(Part3), 2008.

[24] J Aasi et al. Advanced LIGO. Class.Quant.Grav., 32(7), 2015.

[25] F Acernese et al. Advanced virgo: a second-generation interferometric gravitationalwave detector. Class.Quant.Grav., 32(2), 2015.

[26] Y Aso et al. Interferometer design of the KAGRA gravitational wave detector. Phys.Rev. D, 88(4), 2013.

[27] P R Saulson. Terrestrial gravitational noise on a gravitational wave antenna. Phys. Rev.D, 30(4), 1984.

[28] P R Saulson. Fundamentals of Interferometric Gravitational Wave Detectors. WorldScientific Pub Co Inc, 1994.

[29] B J Meers. Recycling in laser-interferometric gravitational-wave detectors. Phys. Rev.D, 38(8), 1988.

[30] T T Fricke et al. DC readout experiment in Enhanced LIGO. Class.Quant.Grav., 29(6),2012.

[31] C Comtet et al. Reduction of tantala mechanical losses in Ta2O5/SiO2 coatings for thenext generation of VIRGO and LIGO interferometric gravitational waves detectors. InProceedings of the 42th Rencontres de Moriond, 2007.

[32] B Cimma et al. Ion beam sputtering coatings on large substrates: toward an improvementof the mechanical and optical performances. Appl Opt., 45(7), 2006.

[33] B W Barr et al. Silica research in Glasgow. Class.Quant.Grav., 19(7), 2002.

References 119

[34] T Accadia et al. Status of the Virgo project. Class.Quant.Grav., 28(11), 2011.

[35] G Cagnoli and P A Willems. Effects of nonlinear thermoelastic damping in highlystressed fibers. Phys. Rev. B, 65(17), 2002.

[36] A V Cumming et al. Design and development of the advanced LIGO monolithic fusedsilica suspension. Class.Quant.Grav., 29(3), 2012.

[37] S M Aston et al. Update on quadruple suspension design for Advanced LIGO.Class.Quant.Grav., 29(23), 2012.

[38] F Matichard et al. Seismic isolation of Advanced LIGO: Review of strategy, instrumen-tation and performance. Class.Quant.Grav., 32(18), 2015.

[39] T Accadia et al. A thermal compensation system for the gravitational wave detectorvirgo. In Proceedings of the 12th Marcell Grossmann Meeting, 2011.

[40] L S Finn and D F Chernoff. Observing binary inspiral in gravitational radiation: Oneinterferometer. Phys. Rev. D, 47(6), 1993.

[41] M Beccaria et al. Relevance of Newtonian seismic noise for the VIRGO interferometersensitivity. Class.Quant.Grav., 15(11), 1998.

[42] J Harms. Terrestrial Gravity Fluctuations. Living Rev. Relativity, 18(3), 2015.

[43] F Bondu and M Barsuglia. Laser Frequency Stabilization Topology. Virgo internalnote VIR-NOT-OCA-1390-247, 2003.

[44] H B Callen and T A Welton. Irreversibility and Generalized Noise. Phys. Rev., 83(1),1951.

[45] A M Gretarsson and G M Harry. Dissipation of mechanical energy in fused silica fibers.Rev.Sci.Instrum., 70:4081, 1999.

[46] C Zener. Internal Friction in Solids II. General Theory of Thermoelastic InternalFriction. Phys. Rev., 53(1), 1938.

[47] M Lorenzini. Suspension thermal noise issues for advanced GW interferometricdetectors. PhD thesis, Università di Firenze, Dip. Astronomia, 2007.

[48] E Campagna. Modeling some Sources of Thermal Noise through Finite ElementAnalysis: Application to the Virgo Interferometric Antenna. PhD thesis, Galileo GalileiSchool of Graduate Studies, Università di Pisa, 2004.

[49] I Martin and S Reid. Coating thermal noise, chapter Optical Coatings and ThermalNoise in Precision Measurement. Cambridge University Press, 2012.

[50] J Aasi et al. Characterization of the LIGO detectors during their sixth science run.Class.Quant.Grav., 32(11), 2015.

[51] J Aasi et al. The characterization of Virgo data and its impact on gravitational-wavesearches. Class.Quant.Grav., 29(15), 2012.

120 References

[52] D M Macleod et al. Reducing the effect of seismic noise in LIGO searches by targetedveto generation. Class.Quant.Grav., 29(5), 2012.

[53] T Accadia et al. Noise from scattered light in Virgo’s second science run data.Class.Quant.Grav., 27(19), 2010.

[54] S K Chatterji. The search for gravitational wave bursts in data from the second LIGOscience run. PhD thesis, Massachusetts Institute of Technology. Dept. of Physics.,2005.

[55] T Accadia et al. The NoEMi (Noise Frequency Event Miner) framework. JPCS, 363(1),2012.

[56] S Chen, S A Billings, and W Luo. Orthogonal least squares methods and their applica-tion to non-linear system identification. Int. J. Contr., 50(5), 1989.

[57] S A Billings, S Chen, and M J Korenberg. Identification of MIMO non-linear systemsusing a forward-regression orthogonal estimator. Int. J. Contr., 49(6), 1989.

[58] M J Korenberg. Identifying nonlinear difference equation and functional expansionrepresentations: The fast orthogonal algorithm. Ann. Biomed. Eng., 16(1), 1988.

[59] M J Korenberg. A robust orthogonal algorithm for system identification and time-seriesanalysis. Biol. Cybern., 60(4), 1989.

[60] H M Abbas. System Identification Using Optimally Designed Functional Link Networksvia a Fast Orthogonal Search Technique. Journal of Computers, 4(2), 2009.

[61] L Piroddi and M Lovera. NARX model identification with error filtering. In Proceedingsof the 17th IFAC World Congress, 2008.

[62] K A Postnov and L R Yungelson. The Evolution of Compact Binary Star Systems.Living Rev. Relativity, 17(3), 2014.

[63] K S Thorne and A N Zytkow. Stars with degenerate neutron cores. I - Structure ofequilibrium models. Astrophys. J., 212(1), 1977.

[64] M Dominik et al. Double Compact Objects I: The Significance of the Common Envelopeon Merger Rates. Astrophys. J., 759(1), 2012.

[65] C Kim. Galactic Merger Rate of Double-Neutron-Star Systems and Prospects forDetecting Gravitational Waves. In The Seventh Pacific Rim Conference on StellarAstrophysics, ASP Conference Series, volume 362, 2006.

[66] C Kim. Implications of PSR J0737-3039B for the Galactic NS-NS binary merger rate.Mon.Not.Roy.Astron.Soc., 448(1), 2015.

[67] R K Kopparapu. Host Galaxies Catalog Used in LIGO Searches for Compact BinaryCoalescence Events. Astrophys. J., 675(2), 2008.

[68] L Chomiuk and M Povich. Toward a Unification of Star Formation Rate Determinationsin the Milky Way and Other Galaxies. Astrophys. J., 142(6), 2011.

References 121

[69] K Belczynski. The fate of Cyg X-1: an empirical lower limit on BH-NS merger rate.Astrophys. J., 742(1), 2011.

[70] A Bauswein, R Ardevol Pulpillo, H T Janka, and S Goriely. Nucleosynthesis Constraintson the Neutron Star-Black Hole Merger Rate. Astrophys. J. Lett., 795(1), 2014.

[71] D A Brown and P J Zimmerman. Effect of eccentricity on searches for gravitationalwaves from coalescing compact binaries in ground-based detectors. Phys. Rev. D, 81(2),2010.

[72] L Blanchet. Gravitational Radiation from Post-Newtonian Sources and InspirallingCompact Binaries. Living Rev. Relativity, 9(4), 2006.

[73] A Buonanno, B R Iyer, E Ochsner, Y Pan, and B S Sathyaprakash. Comparison ofpost-Newtonian templates for compact binary inspiral signals in gravitational-wavedetectors. Phys. Rev. D, 80(8), 2009.

[74] J Centrella, J G Baker, B J Kelly, and J R van Meter. Black-hole binaries, gravitationalwaves, and numerical relativity. Rev. Mod. Phys., 82(4), 2010.

[75] P Ajith et al. Inspiral-Merger-Ringdown Waveforms for Black-Hole Binaries withNonprecessing Spins. Phys. Rev. Lett., 106(24), 2011.

[76] B J Owen. Search templates for gravitational waves from inspiraling binaries: Choiceof template spacing. Phys. Rev. D, 53(12), 1996.

[77] B Allen et al. FINDCHIRP: An algorithm for detection of gravitational waves frominspiraling compact binaries. Phys. Rev. D, 85(12), 2012.

[78] C Cutler and É E Flanagan. Gravitational waves from merging compact binaries: Howaccurately can one extract the binarys parameters from the inspiral waveform? Phys.Rev. D, 49(6), 1994.

[79] T A Apostolatos. Search templates for gravitational waves from precessing, inspiralingbinaries. Phys. Rev. D, 52(2), 1995.

[80] I W Harry, B Allen, and B S Sathyaprakash. Stochastic template placement algorithmfor gravitational wave data analysis. Phys. Rev. D, 80(10), 2009.

[81] F Beauville et al. Placement of templates technique in a 2-D parameter space for binaryinspiral searches. Class.Quant.Grav., 22(20), 2005.

[82] T co*kelaer. Gravitational waves from inspiralling compact binaries: Hexagonal tem-plate placement and its efficiency in detecting physical signals. Phys. Rev. D, 76(10),2007.

[83] P Ajith. Addressing the spin question in gravitational-wave searches: Waveformtemplates for inspiralling compact binaries with nonprecessing spins. Phys. Rev. D,84(8), 2011.

[84] M Burgay et al. An increased estimate of the merger rate of double neutron stars fromobservations of a highly relativistic system. Nature, 426:531–533, 2003.

122 References

[85] K Lo and L Lin. The spin parameter of uniformly rotating compact stars. Astrophys. J.,728(1), 2011.

[86] D A Brown, I Harry, A Lundgren, and A H Nitz. Detecting binary neutron star systemswith spin in advanced gravitational-wave detectors. Phys. Rev. D, 86(8), 2012.

[87] S Privitera. Improving the sensitivity of a search for coalescing binary black holes withnonprecessing spins in gravitational wave data. Phys. Rev. D, 89(2), 2014.

[88] T Dal Canton et al. Implementing a search for aligned-spin neutron star-black holesystems with advanced ground based gravitational wave detectors. Phys. Rev. D, 90(8),2014.

[89] T Adams et al. Low latency search for compact binary coalescences using MBTA. InProceedings of the 50th Rencontres de Moriond Gravitation: 100 years after GR, pages327–330, 2015.

[90] T Adams et al. Low-latency analysis pipeline for compact binary coalescences in theadvanced gravitational wave detector era. arXiv:1512.02864 [gr-qc], 2015.

[91] J Abadie et al. First low-latency LIGO + Virgo search for binary inspirals and theirelectromagnetic counterparts. Astron. Astrophys., 541, 2012.

[92] LSC Algorithm Library Suite. https://www.lsc-group.phys.uwm.edu/daswg/projects/lalsuite.html.

[93] F Marion. Multi-band search of coalescing binaries applied to Virgo CITF data. InProceedings of the Rencontres de Moriond Gravitational Waves and ExperimentalGravity, 2003.

[94] T Dal Canton. Efficient searches for spinning compact binaries with advancedgravitational-wave observatories. PhD thesis, Von der Fakultät für Mathematik undPhysik der Gottfried Wilhelm Leibniz Universität Hannover zur Erlangung des Grades,2015.

[95] B Allen. χ2 time-frequency discriminator for gravitational wave detection. Phys. Rev.D, 71(6), 2005.

[96] G M Guidi. A time-domain veto for binary inspirals search. Class.Quant.Grav., 21(20),2004.

[97] P Shawhan and E Ochsner. A new waveform consistency test for gravitational waveinspiral searches. Class.Quant.Grav., 21(20), 2004.

[98] C Hanna. Searching for Gravitational Waves From Binary Systems in Non-stationaryData. PhD thesis, Louisiana State University, Dept. of Physics., 2008.