University of Amsterdam
Faculteit der Natuurwetenschappen, Wiskunde en Informatica
Report Bachelor Project Physics and Astronomy, size 15 EC, conducted between 29-03-2016 and 07-07-2016.
Determining the Mass Distribution of Neutron
Stars in Binaries Using Gravitational Waves
Author:
Anniek Gloudemans 10631356
Supervisor: Dr. Chris van den Broeck Second assessor: Prof. dr. ir. Paul de Jong
Abstract
Recently, research has been done to determine if the neutron star equation of state can be constrained using gravitational wave signals from coalescing binary neutron stars.3 To determine the necessary parameters from the gravitational wave signal without any biases, the neutron star mass distribution has to be well-known. Unfortunately, there is currently no accurate knowledge of the neutron star mass distribution.14 In this work it is investigated if the correct mass distribution can be determined using gravitational wave signals. A scaling factor will be introduced to correct for systematics due to the mass dependent brightness of the sources. Using the definition of the signal-to-noise ratio, the scaling factor is derived to be s =
Mref
M 5/2
. The effectiveness of the scaling factor is tested with simulated gravitational waves. The outcome suggests that the proposed scaling factor is indeed useful in retrieving the correct neutron star mass distribution from gravitational wave signals.
Contents
1 Introduction 3
2 Theory 6
2.1 Gravitational wave detector . . . 6
2.2 Signal-to-noise ratio . . . 6
2.3 TaylorF2 waveform model . . . 9
2.4 Parameter estimation . . . 10 2.4.1 Bayesian inference . . . 10 2.4.2 Nested Sampling . . . 12 3 Method 14 3.1 Scaling factor . . . 14 3.2 Data analysis . . . 15 4 Results 17 5 Discussion 21 6 Conclusion 22 7 Acknowledgements 22 References 23
1
Introduction
In 1916, Albert Einstein predicted the existence of gravitational waves in his theory of general relativity.10 Gravitational waves (GW) are ripples in spacetime produced by accelerating or decelerating objects. The strongest gravitational waves originate from the most violent events in the universe, such as coalescing binary neutron stars and/or black holes. September 14th 2015, the first gravitational wave was detected with the LIGO detectors, by the LIGO and Virgo collaboration. This might prove to be one of the biggest discoveries of the century.2 Gravitational waves originating from binary
neutron star systems (BNS) might also soon be discovered. The detection rate of coa-lescing binary neutron stars by Advanced LIGO at final design sensitivity is estimated to range from 0.2 yr−1 - 1000 yr−1.19
Neutron stars are the most dense objects known in the universe containing matter at nuclear densities. The final goal is to determine the neutron star mass distribution. The mass distribution is not only interesting on it’s own. One of the main motivations to determine the mass distribution is to constrain the properties of matter at extremely high densities. Another motivation is the role of neutron stars in many astrophysical phenomena, such as gamma-ray bursts.14
A correctly known neutron star mass distribution is also necessary to constrain the neutron star equation of state (EOS) with gravitational waves without any biases. The equation of state gives a relation between the pressure P and the density ρ. Much, yet unsuccessful, research has been done to determine the neutron star equation of state.11, 14 Any gravitational wave signal originating from a coalescing binary neutron
star system will contain information on the equation of state and can therefore be used in future research. A few tens of gravitational wave detections could give a constraint on the equation of state.3, 9 Recently, Agathos and others demonstrated
how the neutron star equation of state can be determined using gravitational waves from coalescing binary neutron stars.3 In this research it is stated that the gravitational wave signal is mainly altered by the equation of state through tidal deformation. When the neutron stars in a binary system are close to each other, their tidal fields will induce a quadrupole moment in one another. This causes the neutron stars to deform. The orbital motion and hence the gravitational waveform is effected by this deformation. The contribution of this effect to the gravitational waveform is dependent on the tidal deformability λ(m). In the quadratic approximation λ(m) can be written as:
λ(m) ≈ c0+ c1 m − m0 M +1 2c2 m − m0 M 2 . (1)
In this equation m is the neutron star mass, m0 the reference mass, and cj coefficients
with j = 0, 1, 2. In practice, only the value of coefficient c0 will be measurable. The
equation of state can already be constrained by measuring the value of c0. However, if
the neutron star mass distribution is not correctly known, large biases appear in the measurement of c0. This effect is shown in Figure 1. Three different equations of state
(MS1). A reference mass of m0 = 1.35 M is used. Large biases appear as a result of a
wrongly assumed mass distribution. Consequently, the neutron star mass distribution will have to be correctly determined to constrain the equation of state.
Figure 1: Measurements of c0, c0= λ(m0), in the case of a soft EOS (SQM3), a moderate one (H4), and a stiff one (MS1). The x-axis depicts the number of sources observed. Left: The masses of the neutron stars are drawn from a Gaussian distribution, whereas a uniform mass distribution of [1, 2]M is supposed. Large systematic error appear in the estimation of c0. Right: The masses of the neutron stars are drawn from a Gaussian distribution that exactly matches the injected mass distribution. The significant biases have largely disappeared.3
Future measurements of gravitational waves originating from binary neutron stars can potentially be used to determine the neutron star mass distribution. A different kind of bias appears when determining the neutron star mass distribution with grav-itational waves, because more massive systems are relatively easier to detect by the gravitational wave detectors. This causes the neutron star mass to be overestimated. A scaling factor is necessary to lift this effect and correctly determine the neutron star mass distribution when using gravitational waves. In this work a scaling factor will be derived to determine the underlying mass distribution from the measured one. The effect of the scaling factor is tested by a simulation. Gravitational waves from coalescing binary neutron stars are simulated and Advanced LIGO detector noise is added to imitate real data. Initial values of distance range, mass distribution, and other necessary parameters are injected in the simulation. When analysing the simu-lated gravitational waves, parameters from the waveform are recovered. The recovered data can then be scaled with the scaling factor. Finally, the injected and the scaled recovered data are compared to determine if the mass distribution is correctly retrieved.
The report is structured as follows. In section 2 the simulation and waveform model are introduced and an expression for the signal-to-noise ratio is given. In section 3 the scaling factor is derived and the method used in data analysis is explained. The main results are shown in section 4. Finally, the result is discussed and concluded in section 5 and 6 respectively. Natural units (G=c=1) are used throughout the report.
2
Theory
In this section first the working principle of a gravitational wave detector is discussed. Then an expression for the signal-to-noise ratio is given. Finally, the method of param-eter estimation to dparam-etermine paramparam-eters from gravitational waves by the simulation is elaborated. This includes Bayesian inference and the nested sampling algorithm.
2.1 Gravitational wave detector
All modern gravitational wave detectors, such as LIGO and Virgo, are interferometers. The current ground-based detectors have two arms of a few kilometers long at an angle of 90 degrees with laser beams running through. Where the arms come together, the laser beams interfere destructively when no distortions are present. When a gravi-tational wave passes, the arms relatively change in length. The laser beams interfere constructively as a consequence, causing light with a certain intensity to appear (Figure 2).17 The sensitivity of the detector is dependent on the frequency of the gravitational wave. In Figure 3 the noise spectrum of Advanced LIGO at design sensitivity, also called noise power spectral density Sn, used in the simulation is shown. The noise
spectrum shows that the detector is able to detect gravitational waves up to a strain sensitivity of approximately 4 × 10−24√1
Hz with an optimal frequency range between
100 and 500 Hz.
Figure 2: Schematic overview of an interferometric gravitational wave detector. (a) The laser beams interfere destructively when no distortions are present. (b) A passing gravitational wave causes the arms to change relatively in length. This causes the laser beams to interfere constructively as a consequence.17
2.2 Signal-to-noise ratio
The relative change in distance caused by a gravitational wave is called the strain h(t) and is used to express the gravitational wave amplitude. Gravitational waves have two
Figure 3: Noise spectrum of Advanced LIGO used in the simulation. The detector is able to detect gravitational waves up to a strain sensitivity of approximately 4×10−24√1
Hz with an optimal frequency range between 100 and 500 Hz. At higher frequencies the sensitivity is limited by shot noise (quantum noise), which arises from statistical fluctuations in photon arrival time. At lower frequencies the sensitivity is limited by thermal and seismic noise.1
independent polarizations, called ’plus’ polarization and ’cross’ polarization. A circular ring of particles will in succession take the form of two ellipses when a gravitational wave is passing, forming either a plus sign or a cross sign (see Figure 4).15 The strain
h(t) can be written as a combination of these two polarizations:
h(t) = F+h+(t) + F×h×(t). (2)
In this equation F+ and F× are called ’beam pattern’ functions, which depend on the
direction of propagation (θ, ϕ) and the polarization angle ψ.17 The two polarizations
h+ and h× from a binary with a circular orbit is to lowest order given by:15
h+= 2ηM D v 2 (1 + cos2ι) cos(2Φ(t)), h×= 4ηM D v 2cos ι sin(2Φ(t)). (3)
In this equation M is the total mass M = m1+ m2, η the symmetric mass ratio
de-fined by equation 6, D the distance, v the velocity, Φ(t) the orbital phase function, and ι the inclination of the orbit with respect to the line of sight. The velocity is defined as
v(t) = Rω, where R is the radius of the orbit and ω the angular velocity of the orbit. If the binary system is observed edge-on (ι = π/2), then h× will be zero. If the system
is observed face-on (ι = 0), then h+and h×are equal in amplitude and π/2 out of phase.
Figure 4: Deformation of a ring of particles caused by a passing gravitational wave. The left figure shows deformation caused by a plus-polarized gravitational wave and the right figure deformation caused by a cross-polarized gravitational wave. The circular ring will sequentially deform in both ellipses when the wave passes.15
The orbit of a coalescing binary neutron star will shrink in size due to the energy loss caused by the emission of gravitational waves. Both the frequency and the amplitude of the emitted gravitational wave will increase in time as the size of the orbit shrinks. This is called ’chirping’.18 Finally, when the size of the orbit reaches a critical point,
the innermost stable orbit, the neutron stars will merge and most likely form a black hole. The mass of a binary neutron star is expressed in the chirp mass M, as defined by equation 4.12 Other useful variables are the mass ratio q and the symmetric mass
ratio η. These are defined as:
M = (m1m2) 3/5 (m1+ m2)1/5 , (4) q = m2 m1 , (5) η = m1m2 (m1+ m2)2 . (6)
The signal-to-noise ratio (SNR) is used as a measure of the strength of the gravita-tional wave signal compared to the background. The SNR is defined as:
(SN R)2 = 4 Z |˜
h(f )|2 Sn(f )
In this equation ˜h(f ) is the waveform and Sn(f ) the noise power spectral density
(PSD) of the detector. In the case of a coalescing binary the waveform is given by equation 11 and the SNR can be written as13 :
(SN R)2 = 5 6 1 π4/3 c2 r2 GM c3 5/3 | Q(θ, φ; ι) |2 fmax Z 0 dff −7/3 Sn(f ) , (8) where Q(θ, φ; ι) = F+(θ, φ) 1 + cos2ι 2 + iF×(θ, φ) cos ι. (9)
In this equation r is the distance from the detector to the source, M the chirp mass, and f the gravitational wave frequency. The integral boundary fmax is approximately
2fs,ISCO, which is the source frequency of the inner most stable orbit where the in-spiral
phase ends. It is defined as:
fmax = 2fs,ISCO =
2 6√6(2π)
1
M. (10)
In this equation M is the total mass. The orientation function Q is given by equa-tion 9. In this equaequa-tion F+ and F× denote the beam pattern functions as described
earlier and ι the inclination of the binary neutron star orbit. If the direction and or-bit inclination is optimal then Q will be equal to 1, otherwise it will be less than 1. The strength of the signal will thus depend on the mass of the neutron stars, the dis-tance to the object, the inclination of the orbit, the polarization, the frequency, and the noise curve of the detector. The noise curve used in the simulation is given by Figure 3.
Heavy neutron stars will produce stronger gravitational waves than lighter neutron stars and can therefore be detected in a greater volume. Consequently, the mass distri-bution of neutron stars directly derived from gravitational waves will be biased towards higher mass and will therefore be incorrect. To neutralize this effect a scaling factor is necessary. As mentioned before, a correct mass distribution is necessary to be able to constrain the neutron star equation of state.3 The scaling factor will be derived from
the equation of the SNR in section 3.1.
2.3 TaylorF2 waveform model
The gravitational waveform produced by a BNS can be accurately described by the post-Newtonian approximation.5 In this simulation the simple frequency domain Tay-lorF2 waveform model for systems with no spin is used. In our case of coalescing binary neutron stars neglecting spin is a good approximation.12 The TaylorF2 waveform is
waveform can be derived from the Einstein field equations. The derivation of the Tay-lorF2 waveform will not be given here. The expression for the TayTay-lorF2 waveform used in the simulation can be written as6, 7, 8 :
˜ h(f ) = 1 r A(θ, φ, ι, ψ, M, η) q ˙ F (M, η; f ) f2/3eiΨ(M,η;f ), (11) where Ψ(M, η; f ) = 2πf tc− φc− π 4 + 7 X i=0 [ϕi+ ϕ (l) i ln(f )]f (i−5)/3. (12)
In equation 11 is r the distance and f the frequency. The amplitude A depends on the sky position (θ, φ), the orientation of the orbital plane to the line of sight (ι, ψ), the chirp mass M, and the symmetric mass ratio η. The frequency sweep ˙F is an expan-sion in powers of frequency arising from the approximation of the Fourier transform. An explicit expression for the frequency sweep ˙F is given by Arun and others.4 The
phase Ψ is given by equation 12 where tc is the time at coalescence, φc the phase at
coalescence, and ϕi and ϕ (l)
i the phase coefficients. The phase coefficients depend on
the component masses of the neutron stars as set by general relativity.
2.4 Parameter estimation
The simulation uses the characteristics of the Advanced LIGO detector to detect the injected gravitational waves and determine the parameters. Bayesian inference is used to discriminate between noise and gravitational wave signals and estimate the param-eters such as the mass and distance of the neutron stars. An algorithm called ’nested sampling’, developed by John Skilling, is used to determine these parameters of the recovered signal.16 The simulation used in this work is in good agreement with the one
described by Veitch and Vecchio in 2010. Their article is used in this section to give a summary of Bayesian inference and nested sampling.19
2.4.1 Bayesian inference
Bayes’ theorem gives the probability of a model Hi being correct, given background
information I and data ~d. This can be expressed as:
P (Hi| ~d, I) =
P (Hi|I)P ( ~d|Hi, I)
P ( ~d|I) . (13)
In this equation P (Hi|I) is the prior probability of the model, Hi, P ( ~d|Hi, I) the
likelihood function, and P ( ~d|I) the marginal probability of the data ~d. Parameter es-timation is used to obtain the unknown parameters in the data and model selection to
detect gravitational waves.
In model selection a comparison between two models is made by calculating the ratio of their probabilities. This ratio of probabilities is called the odds ratio Oi,j and is given
by: Oi,j = P (Hi| ~d, I) P (Hj| ~d, I) = P (Hi|I) P (Hj|I) P ( ~d|Hi, I) P ( ~d|Hj, I) = P (Hi|I) P (Hj|I) Bi,j. (14)
Here Bi,j is called the Bayes factor and is equal to the ratio of the likelihoods.
The gravitational wave form depends on a set of parameters ~θ. The likelihood p( ~d|H, ~θ, I) of both models must therefore be marginalised. This means obtaining the probability distributions of variables contained in the set of parameters ~θ. The marginalisation is done over all parameters and weighted by their prior probability dis-tribution p(~θ|H, I). This procedure results in an integral over parameter space called ’the marginal likelihood’ or ’evidence’ Z, given by:
Z = Z
Θ
p(~θ|H, I)p( ~d|H, ~θ, I)d~θ. (15)
The evidence integral is determined by nested sampling as will be discussed in the next section.
In parameter estimation the posterior density function (PDF) of the parameters is calculated by:
P (~θ| ~d, H, I) = P (~θ|H, I)P ( ~d|~θ, H, I)
P ( ~d|H, I) . (16)
The PDF can be marginalised on a subset of parameters ~θA by integrating over the
other parameters in the set ~θB, given by:
p(~θA| ~d, H, I) =
Z
ΘB
p(~θ| ~d, H, I) d~θB. (17)
From this the posterior mean can be computed:
h~θAi =
Z
ΘA
~
θAp(~θA| ~d, H, I) d~θA (18)
The parameters are estimated by calculating the marginalised posterior density func-tions (equation 17). Besides the evidence integral, nested sampling also calculates the marginalised posterior density functions as a side product. The nested sampling algo-rithm will now be discussed.
2.4.2 Nested Sampling
The goal of the nested sampling algorithm is to calculate the high-dimensional evidence integral with the marginalised posterior density functions as a by-product. In the nested sampling algorithm a number of N samples, called ’live points’ ~θi with i =
1, . . . N , are taken from the prior distribution. The integral in equation 15 can then be approximated as the sum of the weight wi and the likelihood Li of each sample,
resulting in Z ≈ m X i=1 Liwi. (19)
In this equation the weight wi = p(~θi|H, I)d~θ, the likelihood Li = p( ~d|H, ~θi, I), and
m is the number of steps in the algorithm.
After drawing N samples from the prior distribution all the corresponding likelihoods are calculated. The point with the lowest likelihood will then be replaced by a new point with higher likelihood sampled from the prior distribution. This process goes on until the integral does not increase any more than a certain fraction, determined by the termination condition. Each point can be thought of as lying an a contour line of equal likelihood as shown in Figure 5. The enclosed prior mass by the i-th contour surface is denoted as Xi and defined as the fraction of the total prior volume. The largest volume
is enclosed by the lowest likelihood contour line and the smallest volume is enclosed by the maximum likelihood contour line. The weight of each sample is calculated by taking the enclosed prior mass by the contour surface, Xi, and subtracting that of the
previous sample Xi−1.16
Figure 5: Each point can be thought of as lying on a contour line of equal likelihood in parameter space. The enclosed prior mass X decreases when the likelihood increases.16
Finally, all the live points with corresponding likelihood and weight are added in the sum of equation 19 to estimate the evidence integral.
Approximately a number of eH live points are needed to certainly find the maximum
of the likelihood function, where H is the so-called ’information’ measured in nats (1 nat ≈ 1.44 bits). From this is can be derived that approximately N H ±√N H iter-ations will be necessary to reach the zone of high likelihood. As a consequence, the computation time is dependent on the amount of information in the data, determined by the SNR of the signal and the number of live points.
The posterior PDF can be obtained by post-processing the nested sampling output. Samples can be taken from the area Z under the curve L(X) as shown in Figure 6. This area was already decomposed into P
iLiwi by nested sampling. Therefore, we
can just use the live points ~θi with i = 1, . . . N , since the live points already give a
set of posterior representatives. The appropriate weights Liwi are already assigned to
each live point ~θi. The posterior PDF can therefore be written as:
p(~θi| ~d, H, I) ∝ Liwi (20)
Integrating the posterior PDF over the parameters space yields:
1 = Z
p(~θ| ~d, H, I)d~θ discrete≈ AX
i
Liwi = AZ (21)
From this it becomes clear that the factor A = 1/Z. So eventually we obtain:
p(~θi| ~d, H, I) =
Liwi
Z (22)
Figure 6: Samples from posterior PDF are obtained by sampling randomly from the area under the curve of the likelihood L(X).16
3
Method
The injected neutron star mass distribution is uniform between [1, 3] M . In reality we
do not expect the mass distribution to be uniform, but it is useful to clearly demon-strate the effect of the scaling factor. The value of m1 will always be bigger than m2
by convention. In the simulation 5000 sources with a uniform mass distribution are injected within a distance range of [0, 500] Mpc. The noise curve used in the simulation is shown in Figure 3.
In this section an expression for the scaling factor will be derived and the important aspects of the data analysis are highlighted.
3.1 Scaling factor
The scaling factor can be derived from the expression of the SNR (equation 8). The scaling factor has to take into account the neutron star mass and the distance. The integral in equation 8 depends on fmaxas defined by equation 10, which in turn depends
on the neutron star mass. However, the change in value of the integral when varying the neutron star mass is negligible and therefore the integral term in equation 8 does not contribute to the scaling factor. The SNR also depends on the orientation function Q(θ, φ; ι), but since the neutron stars are injected with random angles, the contribution of Q is also negligible. The dominant terms in equation 8 are therefore the chirp mass and the distance, leading to:
SN R ∝ M
5/6
d . (23)
Below a certain value of the SNR no more sources will be measurable. This is called the SNR threshold and is denoted as ρth. The maximum number of sources N (M) the
detector is able to detect is given by:
N (M) ∝ d3max ∝ M
5/2
ρ3th , (24)
since the detection volume V is proportional to d3. Taking a reference mass M ref: N (M) N (Mref) = M Mref 5/2 . (25)
To eliminate the bias, the number of sources the detector is able to detect must not depend on the mass of the source. The scaling factor therefore becomes:
s = Mref M
5/2
. (26)
In this work a reference mass of Mref = 1.38 M is used. Note that the choice of
reference mass does not have any influence on the result, since it will disappear when normalizing. The normalized intrinsic mass distribution is given by:
Pintrinsic(M) = M ref M 5/2 Pmeasured(M) R Mref M 5/2 Pmeasured(M)dM . (27)
Since Mref is just a constant, it can be taken out of the integral. The reference mass
in the numerator and denominator then cancel each other out.
The scaling factor is used to rescale the recovered neutron star mass distribution. By comparing the scaled and unscaled recovered mass distribution to the injected mass distribution, the influence of the scaling factor can be determined. Details of the anal-ysis of the recovered data are given in the next section.
3.2 Data analysis
From the simulation, posterior density functions (PDF’s) of the chirp mass M and mass ratio q are retrieved for every detected source. The posterior density functions give the probability for a variable to be a certain value. Samples are drawn from the PDF’s of the chirp mass and mass ratio. We start the data analysis with these samples.
The effect of the scaling factor can be shown by comparing the injected and recov-ered chirp mass. First, the median value of the posterior density function of the chirp mass is calculated for each source. Then, a histogram is made of both the injected and recovered chirp mass. The recovered histogram is scaled using the scaling factor and finally all histograms are normalized to unity. The result is shown in Figure 9.
Another way to show the effect of the scaling factor is to directly use the PDF’s to determine the distribution of the component masses m1 and m2. The first step is
performing a coordinate transformation. Each point in the posterior density functions of M and q, pdf (M, q), is converted to the component masses m1and m2. The resulting
posterior density function pdf (m1, m2) will thus depend only on the component masses.
Using the definition of q and M,
m1 = M
(1 + q)1/5
q3/5 , (28)
m2 = M(1 + q)1/5q2/5. (29)
By convention, m1 is always bigger than m2, since q is always smaller than 1.
After coordinate transformation, a kernel density estimation (or KDE) is used to estimate the probability density functions of the component masses. The next step is marginalising and rescaling the KDE’s of m1 and m2. In probability theory the
posterior density function pdf (m1, m2) can be marginalised by integrating over m1 or
pdf (m1) = m2,max Z m2,min s(m1, m2)pdf (m1, m2)dm2, pdf (m2) = m1,max Z m1,min s(m1, m2)pdf (m2, m1)dm1. (30)
The scaling factor is calculated using the definition of the chirp mass and equation 26. Equation 30 will thus be used to rescale the KDE’s of m1 and m2. In the code
we assume that pdf (m1, m2) = pdf (m1)pdf (m2). This is in fact only true if masses m1
and m2 are uncorrelated. This is not the case however, due to the definition of the
chirp mass and mass ratio. For instance, if component mass m1 is small then m2 will
have to be large to obtain the same chirp mass value. Consequently, the component masses m1 and m2 are correlated. The assumption pdf (m1, m2) = pdf (m1)pdf (m2) is
therefore not mathematically correct, but still useful to illustrate the efficiency of the scaling factor.
Finally, the scaled KDE’s of m1 and m2 of each source are separately added to obtain
the mass distribution of m1 and m2 (Figure 10). The total mass distribution mtot is
calculated by adding both the KDE’s of m1 and m2 (Figure 11). After adding the
KDE’s for each source, the total KDE is normalized to unity. This method of adding KDE’s is useful for now to illustrate the effect of scaling, but not completely valid. This will be discussed in section 5. One last remark: in this method the scaling factor is used to rescale the individual KDE’s instead of the total distribution. Although this is not the correct way of applying the scaling factor, it can be used in this case since the KDE’s are added in the end.
4
Results
In the simulation 5000 sources within a distance range of [0, 500] Mpc and uniform mass distribution were injected. A number of 613 sources were recovered from the simulation and used in the analysis. In this section the results are given.
The injected and recovered number of sources at each distance is displayed in Figure 7. Sources are recovered within a distance range of [47,411] Mpc. At higher distances the SNR is to low to detect any gravitational wave from neutron stars using the Ad-vanced LIGO detector.
Figure 7: Distance plot of injected and recovered sources. A number of 5000 sources are injected in a distance range of [0,500] Mpc and 613 sources are recovered in a distance range of [47,411] Mpc. Approximately 12 percent of the sources is recovered by the simulation.
The median chirp mass M of each source is plotted against its recovered distance to determine if low mass neutron stars are indeed harder to be found at large distances. The result is shown in Figure 8. The plot is divided into four boxes: [0,110] Mpc, [110,220] Mpc, [220,330] Mpc, and [330,440] Mpc. For each box the median chirp mass and standard deviation of the sources is calculated. The results are shown in Table 1. Figure 8 and Table 1 show that the higher the distance, the larger the typical chirp mass of the recovered sources.
Figure 8: Scatter plot of the chirp mass and recovered distance. The plot is divided in four boxes: [0,110] Mpc, [110,220] Mpc, [220,330] Mpc, and [330,440] Mpc. For each box the median and standard deviation of the chirp mass of the sources found in that particular box is calculated. Fewer low mass neutron stars are found at a distance above 250 Mpc.
Distance (Mpc) [0,110] [110,220] [220,330] [330,440]
Median M 1.65 1.70 1.80 2.04
Standard deviation 0.32 0.38 0.33 0.30
Number sources in box 14 118 256 225
Table 1: Median chirp mass and standard deviation of sources in a distance range of [0,110] Mpc, [110,220] Mpc, [220,330] Mpc, and [330,440] Mpc. The typical median chirp mass of the recovered sources increases towards higher distances.
Figure 9 shows a histogram of the injected and recovered median chirp mass. The recovered chirp mass is, as expected, biased towards higher values. However, scaling the recovered mass distribution with the derived scaling factor result in a distribution very similar to the injected distribution. This result suggests that the scaling factor indeed eliminates the bias.
As explained in section 3.2, the posterior density functions of the chirp mass M and mass ratio q can be used to determine the distribution of m1 and m2. The estimates
of the posterior density functions (KDE’s) of m1 and m2 of each source are added
and eventually normalized. When adding both the KDE’s of m1 and m2 the total
Figure 9: Histogram of the injected, recovered, and recovered scaled median chirp mass. When recov-ering the data the chirp mass is slightly overestimated. Using the scaling factor the mass distribution is retrieved fairly well.
left figure shows the recovered mass distribution of the component masses when no scaling factor is applied and the right figure when it is applied. In Figure 11 the total mass distributions are plotted together to make a comparison. The injected uniform mass distribution of [1,3] M is not correctly retrieved by the simulation without the
scaling factor and we do indeed see a bias appearing towards larger mass. On the other hand, the uniform mass distribution seems to be correctly retrieved when applying the derived scaling factor.
Figure 10: Left: Unscaled recovered mass distribution of p(m1), p(m2), and p(m1) + p(m2). The uniform mass distribution is not correctly retrieved, since there is a bias towards larger mass values. Right: Scaled recovered mass distribution of p(m1), p(m2), and p(m1) + p(m2). The uniform mass distribution seems to be correctly retrieved. In both plots m1 is bigger than m2 by convention. The distributions are normalized to unity.
Figure 11: The neutron star mass distribution with and without scaling. The uniform mass distribu-tion [1,3] M is to good estimation retrieved when applying the scaling factor.
5
Discussion
For now, it seems justifiable to assume that the derived scaling factor is correct. Since this research has never been done before, the results can not be compared. In this section the reliability of the results and the possible next steps are discussed.
The method used to test the scaling factor is a simplification of reality. Firstly, an uniform mass distribution is injected in the simulation. In reality the neutron star mass distribution is not expected to be uniform, but more likely some Gaussian distri-bution. Secondly, the number of recovered sources is unrealistic high when comparing to a realistic situation. For the upcoming years, we do not expect to observe this many gravitational waves signals originating from coalescing binary neutron stars. Finally, only one detector is used in the simulation, which differs from reality. The uniform mass distribution, the high number of sources, and the single detector are merely used to demonstrate the efficiency of the scaling factor and do not make the results less reliable.
In Figure 10 and 11 the estimates of the posterior density functions (KDE’s) of each source are added to obtain the mass distribution from the recovered data. This method is not completely valid, since adding KDE’s does not have any direct statistical mean-ing. A single KDE gives the probability of a variable having a certain value for a given detection. When adding these probabilities it does not give you the overall probability distribution for masses of neutron stars in binaries. Even though the method does not have any statistical meaning, it is still useful in showing the effect of the scaling factor. A next step in this research could be combining the KDE’s in a more appropriate way. Other possible next steps could be injecting a Gaussian mass distribution or including several detectors in the simulation.
When the first gravitational waves from coalescing binary neutron stars are indeed detected, we have to keep in mind that we can only determine both the mass distribu-tion and EOS parameters from the data if they are uncorrelated. To test this, half of the gravitational wave detections can be used to calculate the mass distribution and the other half to calculate the EOS parameters. When also done the other way around the result has to be consistent.
6
Conclusion
To conclude, it seems that the neutron star mass distribution can be correctly deter-mined using gravitational waves. A scaling factor is needed to correct for the relatively higher amount of heavy neutron stars measured by the detector. The scaling factor can be derived from the signal-to-noise ratio and is equal to s =
Mref
M
5/2
. In this work the scaling factor is tested through simulation and the results look promising. A well determined mass distribution is necessary to eliminate the bias in the tidal de-formability λ(EOS; m), which can be used to constrain the neutron star equation of state.
The first gravitational wave measurement of coalescing binary neutron stars will in all probability happen in the near future. A few tens of gravitational wave detections could give a constraint on the equation of state.
In follow-up research one could also derive the scaling factor needed to determine the mass distribution of binary black hole systems by measuring gravitational waves generated in the merge and ring-down. The SNR in this stage will not be proportional to M5/6, so the scaling factor needed to determine the mass distribution will differ
from the one derived in this work.
7
Acknowledgements
I would like to express my great appreciation to the people of the Nikhef Virgo group. Especially, I would like to thank Laura van der Schaaf for her excellent daily supervision and Chris van den Broeck for his invaluable and constructive suggestions.
References
[1] B. Abbott et al. GW150914: The Advanced LIGO Detectors in the Era of First Discoveries. 2016.
[2] B. Abbott et al. Observation of Gravitational Waves from a Binary Black Hole merger. Phys. Rev. Lett., 116:061102, 2016.
[3] M. Agathos, J. Meidam, W. D. Pozzo, et al. Constraining the neutron star equation of state with gravitational wave signals from coalescing binary neutron stars. 2015. [4] K. G. Arun, B. R. Iyer, B. S. Sathyaprakash, and P. A. Sundararajan. Parameter estimation of inspiralling compact binaries using 3.5 post-newtonian gravitational wave phasing: The nonspinning case. Phys. Rev. D, 71:084008, 2005.
[5] L. Blanchet. Gravitational radiation from post-newtonian sources and inspiralling compact binaries. Living Reviews in Relativity, 17(2), 2014.
[6] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet. Gravitational-wave inspiral of compact binary systems to 7/2 post-newtonian order. Phys. Rev. D, 65:061501, 2002.
[7] T. Damour, B. R. Iyer, and B. S. Sathyaprakash. Comparison of search templates for gravitational waves from binary inspiral. Phys. Rev. D, 63:044023, 2001. [8] T. Damour, B. R. Iyer, and B. S. Sathyaprakash. Comparison of search
tem-plates for gravitational waves from binary inspiral: 3.5pn update. Phys. Rev. D, 66:027502, 2002.
[9] W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale. Demon-strating the feasibility of probing the neutron-star equation of state with second-generation gravitational-wave detectors. Phys. Rev. Lett., 111:071101, 2013. [10] A. Einstein. Die Grundlage der allgemeinen Relativit¨atstheorie. Annalen der
Physik, 354:769–822, 1916.
[11] J. M. Lattimer. The nuclear equation of state and neutron star masses. 2013. [12] T. G. F. Li, W. D. Pozzo, S. Vitale, et al. Towards a generic test of the strong
field dynamics of general relativity using compact binary coalescence: Further investigations. 2011.
[13] M. Maggiore. Gravitational Waves. Oxford University Press, 2008.
[14] F. Ozel and P. Freire. Masses, Radii, and Equation of State of Neutron Stars. 2016.
[15] B. S. Sathyaprakash and B. F. Schutz. Physics, Astrophysics and Cosmology with Gravitational Waves. 2009.
[16] J. Skilling. Nested sampling for general Bayesian computation. Bayesian Anal., 1(4):833–859, 2006.
[17] G. F. L. Tjonnie. Extracting Physics from Gravitational Waves. PhD thesis, Nikhef, 2013.
[18] C. van den Broeck. Gravitational wave searches with Advanced LIGO and Ad-vanced Virgo. 2015.
[19] J. Veitch and A. Vecchio. Bayesian coherent analysis of in-spiral gravitational wave signals with a detector network. Phys. Rev., D81:062003, 2010.
Populair wetenschappelijke samenvatting
In 1916 formuleerde Albert Einstein de algemene relativiteitstheorie, waarin onder andere werd gesteld dat versnelde objecten rimpelingen in de ruimte veroorzaken, genaamd gravitatiegolven. Als een gravitatiegolf langs de aarde komt dan zal de ruimte een klein beetje vervormd worden, waardoor objecten verder of minder ver van elkaar af komen te staan. Een detector met armen van een paar kilometer lang kan de vervorming van de ruimte meten als lengte verandering van de armen. Deze lengte verandering is ontzettend klein en staat gelijk aan het veranderen van de afstand tot de dichtstbijzijn-ste dichtstbijzijn-ster buiten ons zonnedichtstbijzijn-stelsel met de dikte van een haar. De afgelopen 40 jaar hebben wetenschappers gezocht naar deze gravitatiegolven. Op 14 September 2015, 100 jaar na de publicatie van Einstein, is voor het eerst een gravitatiegolf waargenomen met de LIGO detectoren door de LIGO en Virgo collaboratie. Dit signaal was afkomstig van twee om elkaar heen draaiende zwarte gaten die op het punt stonden samen te smelten. Deze ontdekking is niet alleen een goede test voor de algemene relativiteitstheorie van Einstein, maar geeft ons ook veel nieuwe mogelijkheden om het heelal te bestuderen. Door het meten van gravitatiegolven kunnen we namelijk een stuk verder in het heelal kijken dan met telescopen. ”Up until now, we’ve been deaf to gravitational waves, but today we are able to hear them.” zei de uitvoerend directeur van LIGO tijdens de ceremonie, ”What’s going to come now is we’re going to hear more things, and no doubt we’ll hear things that we expected to hear. . . but we will also hear things that we never expected.”
Artist impression van gravitatiegolven afkomstig van dubbele neutronensterren.
In de nabije toekomst zullen we waarschi-jnlijk ook gravitatiegolven meten afkomstig van twee neutronensterren die dicht om elkaar heen draaien. Deze meting zou mogelijk kun-nen bijdragen aan het bepalen van de struc-tuur van neutronensterren. Neutronensterren zijn de kleinste en meest compacte sterren in het heelal die onstaan door het instorten van een ster met een massa 8 keer zo groot als onze zon. Om de structuur te kunnen bepalen is het belangrijk te weten hoe zwaar de ver-schillende neutronensterren zijn. Een prob-leem is echter dat telescopen te weinig neutro-nensterren waarnemen om hiervan een goede schatting te kunnen maken. We moeten dus
proberen om gravitatiegolven te gebruiken om de ’massaverdeling’ van neutronenster-ren nauwkeurig te bepalen. In dit onderzoek heb ik mij bezig gehouden met het bepalen van de massaverdeling van neutronensterren. Hierbij heb ik data gebruikt van ges-imuleerde gravitatiegolven afkomstig van neutronsterren, aangezien echte metingen van deze systemen nog niet zijn gedaan. De gravitatiegolf detectoren zullen zware neutronsterren makkelijker meten dan lichte neutronensterren, omdat deze sterkere
golven uitzenden. Het gevolg hiervan is dat neutronensterren waarschijnlijk gemiddeld zwaarder worden geschat dan ze daadwerkelijk zijn. De hoofdvraag in mijn onderzoek is welke correctie er gedaan moet worden om de juiste massaverdeling te bepalen met behulp van gravitatiegolven. Wellicht kan met deze informatie de structuur van deze fascinerende objecten bepaald worden in de nabije toekomst.
