Likelihood Investigation Between the Time of Detection Between Two IFO's and Their SNR Ratio



One can establish preliminary detection likelihood cuts based on the detection delay between two detectors and, assuming the delay is physically possible for a gravitational wave traveling at the speed of light, the SNR ratio between the two detectors. 

For this preliminary investigation, only the polarization averaged antenna patterns for each detector are used and it is assumed that all of the detectors have an identical noise profile.  All of the values shown here are the direct measurements from functions used in the GravEn GW software simulator available in the MATAPPS repository (namely ifodelay.m [time projector for GW incident on the center of the Earth to the detector of interest] and detproj.m [antenna pattern generator]).  The codes used in the production of these results are also available in the MATAPPS repository (matapps/src/searches/inspiral/delayRatio).

A first likelihood cut is obvious - is the delay in detection between the two detectors physically possible given the maximum time of flight (TOF) of light between the two detectors?  The numbers listed below assume a zero timing error:

TOF in ms
 LHO
 LLO
 VIRGO
 GEO
 TAMA
LHO
 0
 10.013 27.288 25.106
 24.929
LLO
 10.013
 0
 26.448
 25.107
 32.299

Of course, there is never zero timing error.  Therefore, the maximum allowable delay between two detectors will be the above times +/- \sqrt{\epsilon_1^2 + \epsilon_2^2} where \epsilon_n is the timing error associated with detector n.

Assuming the candidate event passes this first cut, we can then compare the SNR ratio to the range of possible ratios for the given delay.  First, the circle on the sky for the given delay is determined and then applied to the polarization averaged antenna pattern for each detector.  The antenna pattern coefficients for each coordinate on this circle are then used to generate a ratio.  This ratio is analogous to the SNR ratio of the candidate events in the 2 detectors.

The time of flight from the center of the Earth to a detector is the negative dot product of the detector's WGS-84 coordinates and the sky location of the gravitational wave relative to the detector (i.e. Earth fixed coordinates) divided by the speed of light.  The time of flight between two detectors is simply the differences in these times.  Since this resulting expression is not separable into the sky coordinates producing a desired time of flight, a table of TOFs was produced and desired values were used to select the sky locations that can produce this delay.  This method introduces an error in the sky locations since it is unlikely to produce a continuous circle of sky locations from a finite table.  However, one can specify the tolerance in the allowable TOF from the exact desired value and this can also represent the actual timing errors from the detectors.  In the following results, the TOF table used was a 1000x1000 sample table evenly distributed over the Mercator projection of Earth's coordinates and tolerance is set to be +/- 0.01 times the maximum TOF between the 2 detectors.

Below is the TOF table between LHO and LLO:
TOF between LHO and LLO
TOF tables for other LIGO combinations:

LHO and VIRGO
 LLO and VIRGO
LHO and GEO
 LLO and GEO
LHO and TAMA
 LLO and TAMA

An illustration of the selection of source locations for a zero difference in delay time (+/- 0.01*10.013 ms) between LHO and LLO is shown below:
TOF = 0 +/- 0.10013 ms
where the black line represents the selected sky locations for the TOF = 0 +/- 0.10013 ms.

These discrete sky locations are then used to derive the polarization averaged antenna pattern coefficients from theory.   The projection of these sky locations on the antenna pattern are shown below for LHO and LLO:
TOF = 0 +/- 0.10013 ms projected on LHO antenna patternTOF = 0 +/- 0.10013 ms projected on LLO antenna pattern
Since the Mercator map projection is used, each sample does not represent an equal area of the globe; the larger the absolute value of the latitude, the smaller the area each sample represents.  Therefore, considering the distribution of calculated ratios isn't meaningful since the poles are more heavily sampled than the equator.  However, the minimum and maximum ratios set bounds for the range of likely ratios.  The ratio bounds for TOF = 0 +/- 0.10013 ms are 0.6593 to 1.5262.

Measuring the range of likely ratios for discrete TOF values between + and - the maximum time of flight will generate a table that can be referenced when weighing the likelihood of a candidate event.  The number of samples is determined so that the dt between samples is less than the 1/16384 sec (the sample period of LIGO) - this should insure an acceptable interpolation error when using these min/max ratio tables in practice.   The measure minimum and maximum ratios WRT the candidate event's TOF between LHO and LLO is shown below:

Min/max ratios vs TOF

The poles and zeros are the result of the sky location circle passing near the zenith of a detector antenna pattern zero.  Since LIGO's detector's antenna patterns are nearly collocated, the poles and zeros come in pairs for specific times.  This is not the case for the combination of LIGO and other international detectors.  An SNR ratio that falls between the minimum and maximum ratio curves is a plausible candidate gravitational wave event, while above the maximum curve or below the minimum curve are unlikely.

Minimum and maximum ratio plots are shown below for other LIGO combinations:

LHO and VIRGO (900 samples)
 LLO and VIRGO (870 samples)
LHO and GEO (826 samples)
 LLO and GEO (826 samples)
LHO and TAMA (818 samples)
 LLO and TAMA (1060 samples)

The following are text files containing the minimum and maximum ratios (1st column: time delay, 2nd column: minimum ratio, 3rd column: maximum ratio):

LHO and LLO
LHO and VIRGO
 LLO and VIRGO
LHO and GEO
 LLO and GEO
LHO and TAMA
 LLO and TAMA


It is also interesting to note that there is a band of allowable ratios for all physically possible time delays (TOFs).  These values are a useful intermediate check once it is known that the time of flight of a candidate event is possible.  These ratios are shown in the table below:

Ratio Pair
 Minimum ALWAYS Allowed
 Maximum ALWAYS Allowed
LHO/LLO 0.659
 1.526
LHO/VIRGO 0.766
 2.062
LHO/GEO 0.523
 1.460
LHO/TAMA
 0.708
 2.392
LLO/VIRGO
 0.527
 1.167
LLO/GEO
 0.506
 2.198
LLO/TAMA
 0.670
 1.657

Once it has been determined that the SNR ratio is plausible for the candidate event given its detection delay between IFOs, the analysis can be used to seek the possible sky locations of the source assuming unpolarized waves (see findskyloc.m).  Given the actual detection delay, determine the circle on the sky that produces that delay within a given tolerance (I will still use 0.01*maxDelay) and calculate the antenna pattern ratios for each point on the circle.  Find the calculated ratios that match with a minimum error (the code begins seeking ratios within +/-  0.5% and increments the error by 0.5% iteratively until suitable ratios have been found).  Separate sky locations (which are arcs of the circle produced by detection delay) are isolated and the center (half arc length point) of each location arc are output.  The black arcs are all points on the circle that matched the desired ratio within 0.5% and the green circles with an x inside is the center of the arc which is reported (shown below the plot):
Found sky locations; t = 0, SNR ratio = 1
>> findskyloc('LHO', 'LLO', 0, 1)

The ratio error used is +/- 0.5%.

There appears to be 4 possible sky location(s) for this event:
  Location 1:     -38.16 latitude,    -153 longitude
  Location 2:     39.24 latitude,     -103.68 longitude
  Location 3:     38.52 latitude,     27.72 longitude
  Location 4:     -38.88 latitude,    75.96 longitude

N.B.:
  1. Negative latitude values indicate S; Negative longitude values indicate W.
  2. These coordinates are in Earth fixed coordinates, i.e. the sky revolves around the Earth.
     Further calculation is needed to determine what part of the universe was zenith to these locations.

Note that these locations do not give error bars.  What is shown here is the just an illustration that possible sky locations can be isolated.


Application to the follow-up pipeline:

Since the calculations used to measure the plausible SNR ratios shown above are computationally intensive, it isn't feasible to make these computations on the fly.  Therefore, I suggest using the min/max ratio data files given above and interpolating ratios between TOFs for a specific candidate event.  For now, this will only be applied to two detectors (H1 and L1) due to the assumption that the detectors have identical noise profiles.

Suggested algorithm (see tof2snr.m):


Further Work:


The Effect of Polarization on the Observed SNR Ratio

When the polarization of gravitational waves are taken into account, it is expected that observed SNR ratios can fall outside the bounds determined using the polarization averaged antenna patterns.  To investigate this effect, GravEn was used to simulate 10,000 coincident GW simulations between LHO and LLO that had random polarization angle and sky locations.  The observed SNR ratio for each simulation is plotted (with an x) against the simulation's time of flight between detectors in the following illustration:


It is clear here that the effect of polarization can indeed pull a physically possible signal outside of the SNR ratio bounds determined with polarization averaged antenna patterns.  However, there does not appear to be any change in distribution around the unity ratio with respect to the time delay of the signals between detectors.  Also, the distribution is an even function around the unity SNR ratio.  Therefore, we can use the absolute difference between the unity ratio and the measured ratio if we also use the simulations in a similar manner and fit the distribution.  This is shown in the plots directly below:


The distribution is exponential with a probability density function (pdf) of:

pdf(ln(SNR)) = 0.726479*exp(-|ln(SNR)|*0.726479)

Therefore, the probability (P) of observing a gravitational wave between the LIGO detectors an absolute distance away from the unity SNR ratio is:

P(ln(SNR)) = 1-exp(-|ln(SNR)|*0.726479)

The driver function (tof2snr.m) that determines the likelihood of a candidate event being physical has been modified to return a Boolean indicating if the detection delay and observed SNR ratio is within polarization averaged bounds AND the probability of the ratio being observed for polarized waves. 

Let \lambda represent the maximum likelihood estimation for the pdf and P defined above (which was 0.726479).  Using simulations coincident to those used above, the \lambda for the other detector combinations is given below:

\lambda
 LHO
 LLO
 VIRGO
 GEO
 TAMA
LHO
 N/A
 0.726479
 1.1066
 1.0602
 1.1089
LLO
 0.726479
 N/A
 1.0753
 1.1291
 1.1221



Exact Sky Locations as a Function of Detection Delay


The equation for the circle on the sky describing the source locations for a given detection delay between two detectors is difficult to separate into the two angles on the unit sphere.  The easy solution to this problem is to define a new coordinate system with the two detectors in question on the polar axis with the origin centered between the detectors (I will refer to this coordinate system as polar baseline coordinates).  The detection delay then gives a circle on the unit sphere (only a great circle when there is no detection delay) that can be completely described by a declination angle.  Points on this circle can then be rotated back into Earth coordinates.  These concepts are outlined in Saulson's book.

exactskycircle.m first retrieves the WGS-84 coordinated of the two detectors in question and calculates the distance, D, between the detectors.  The declination of the circle on the unit sphere describing the exact solution to the detection delay, \Theta, is given by (eq 15.1 in Saulson's book):

\Theta = arcsin((c*dt)/D)

where c is the speed of light and dt is the detection delay.  Given a desired angular resolution of this circle, discrete coordinates for this circle are computed.  Then the Euler angles needed to rotate this circle into Earth coordinates are computed and the rotation performed.  The exact coordinates producing a zero detection delay between LHO and LLO are plotted (in blue) on top of the seeking method described above (in black):



The results of this method solve the issue of over sampling the coordinate space near the poles when seeking solutions in the Mercator projection of Earth's coordinates since the exact solution is evenly sampled in the baseline polar coordinates and then rotated in Earth coordinates.  This is especially obvious when looking at a detection delay (here -0.004 seconds) which has a sky circle which passes near Earth's poles:

Using exact solutions also has the potential to be used in reverse for efficient and accurate determinations of possible source locations.  If an externally triggered candidate event is under consideration, just the use of these timing codes can be of value.  Say, only two detectors were taking science data when the external trigger happened.  The likelihood of a candidate event can be estimated without serious consideration of the SNR ratios observed (at least initially) since one will have prior knowledge of the location of the external trigger; if the sky circle does not coincide with the external trigger location, then it cannot be a gravitational wave associated with that event.  Code such as this can also be used in the future to externally trigger other EM observatories if the online analysis from 3+ detectors give a near real time event since 3 detectors narrow the sky locations to two areas and 4 would give a single area (area is defined to be the point and the surrounding source locations included within the timing error).

- Amber L. Stuver