LISA Science Performance Evaluation (by Parameter Estimation)

This is a Wiki; login to edit. If you don’t have a login or you don’t have enough privileges, e-mail Curt.

How to upload files.

• Teleconferences

In an effort to cross validate the various black hole parameter estimation codes that have been developed the following exercise has been suggested:

The first step is to adopt the MLDC conventions for the LISA orbit. A brief version of the MLDC conventions document can be found here. The noise model to use is available in the form of a c-code Noise.c, which uses the header file Cnst.h. The LISA orientation parameters kappa and lambda should both be set to their MLDC default values of zero and the low frequency cut-off should be set at 1e-5 Hz. It is also very important that all the codes use the same values for the physical constants found in Cnst.h.

Four test examples are provided in the files mbh1a.txt, mbh1b.txt, mbh2a.txt, mbh2b.txt. The “a” set is for one year of observations and a time to merger of one year. The “b” set is for the same systems, but with time to merger of 1.05 years. The “b” sets should be helpful since the results will not depend on how the waveforms are truncated at the ISCO.

For the purposes of the comparison exercise both the amplitude and the phase should be computed to 2PN order.

The masses are given in solar mass units, and correspond to the physical masses. The redshifted masses are (1+z) times larger.

The sky location is given in terms of the ecliptic co-latitude θ_s and longitude Φ_s. The ecliptic latitude, lat_s = π/2 - θ_s is also quoted to remind people to be careful about this convention. The orbital orientation is given in terms of co-latitude θ_l and longitude Φ_l. The orientation is also given in terms of inclination ι and polarization angle ψ using the standard MLDC conventions. The quantity φ_c is the ORBITAL phase at coallesence. All angles are quoted in radians.

• The general spin + higher harmonic code developed by Cornish & Hughes can handle zero spin systems. Using the conventions described above it returned the following parameter errors: mbh1a_ch2.txt, mbh1b_ch2.txt, mbh2a_ch2.txt, mbh2b_ch2.txt. The Fisher matrix entries are also quoted, as are the Fisher correlations, Γ_ij/√(Γ_iiΓ_jj). The results match the independently developed Cornish & Porter code to very high accuracy.

• The waveforms are posted here as gzipped tar balls for Case 1a and Case 2a. After unzipping you will have files such as KerrA_1a.dat and KerrE_f_1a.dat. The term “Kerr” refers to the fact that the code is able to handle full spin evolution, even though these particular cases are non-spinning. The A or E refers to which LISA channel (like Curt’s I and II), and the subscript f indicates that the waveforms are in the frequency domain. The time domain waveforms are two column, with the first column giving the time of the sample. The frequency domain results are three column, and list the frequency, followed by the real and imaginary parts of the strain spectral densities.

• Waveform comparison between Cornish & Hughes (CH) and Sintes & Trias (ST). There are two things to have in mind: 1st) CH are using Cutler’s 98 convention (u₊ = 1+(L·n_DS)² and w_x = - 2 L·n_DS, where n_DS is the unit vector detector → source) while TS are using Blanchet’s 96 (u₊ = - (1+(L·n_DS)²) and w_x = 2 L·n_DS), so there is an overall sign in h(t) and h(f) to be corrected in order to compare waveforms. Notice that an overall sign does not modify SNRs and errors. 2nd) We were not using the same phase at coalescense value, in fact TS must use phic = 5.796995 (instead of 4.04657) for mbh1a and phic = 4.27908695 (instead of 2.697954) for mbh2a. After correcting that, these are some examples of the results we get: CH_ST_mbh1a_h1t.pdf, CH_ST_mbh1a_h2f.pdf, CH_ST_mbh2a_h1f.pdf and CH_ST_mbh2a_h2t.pdf. In the time domain there are small differences, just before merger, due to CH’s attenuation and in the frequency domain, we do not match exactly all details because CH has less resolution; but except for that, matching is pretty good. Note that the Cornish & Hughes power spectral density has been multiplied by the observation time in order to compare to the Sintes & Trias power spectra. Modulo the usual finite FFT ringing (see Phys. Rev. D59, 124016 (1999)), and the additional ringing caused by the filter that terminates the time domain signal, despite the constant phase difference, the agreement is pretty good.

• Update: February 6th, 2008. At the last telecon it was suggested that we look at the waveforms harmonic by harmonic. The tar ball Case 1a Harmonics contains separate files for the 0PN, 0.5PN, 1PN, 1.5PN and 2PN harmonics. In all cases the phase is computed to the full 2PN order, so the break down is by the PN order of the amplitudes. In other words, if you sum the data in the 5 files for the A or E channels you will recover the full waveforms we posted earlier. It was also suggested that we compute the SNRs using a high frequency cut-off of 1 mHz. This made no difference for Cases 1a and 1b as the signals terminate before 1 mHz. For Case 2a the 1 mHz cut-off reduced the SNR from 247.0 to 101.2. If we use a frequency cut-off of 0.1 mHz, the Case 1a SNR drops from 167.3 to 53.6.

• Update: February 14th, 2008. The parameter estimation errors for Case 1a with a 0.1 mHz cut-off are here Case 1a cut and those for for Case 2a with a 1 mHz cut-off are here Case 2a cut. The Fisher matrices and their inverses for these two cases are here: Case 1a Fisher 0.1mHz and Case 2a Fisher 1mHz

• Main results. These are our results (SNRs and errors) for the different cases: mbh1a_st.txt, mbh1b_st.txt, mbh2a_st.txt and mbh2b_st.txt. We are also posting our Fisher matrices for Detector 1: 1a_fisher_I, 1b_fisher_I, 2a_fisher_I and 2b_fisher_I, and for the two combined detectors (I+II): 1a_fisher_IpII, 1b_fisher_IpII, 2a_fisher_IpII and 2b_fisher_IpII. The ordering of the parameters is the following: phi_N, phi_L, cos(theta_N), cos(theta_L), tc, phic, logDL, mt, deltam.
• As it was suggested at the last telecon, we have computed SNRs (and also errors) using a high frequency cut-off. For mbh1a case results were not modified by using f_cut-off = 1 mHz, so we used 0.1 mHz and the new SNR we get is 53.3 (here are the new errors: mbh1a_st_1em4.txt). For the mbh2a using f_cut-off = 1 mHz the new SNR is 100.9 and the errors: mbh2a_st_1em3.txt.

• Waveforms. The waveforms we obtain for the ‘1a’ and ‘2a’ cases (here we’re using a phic value to match CH’s waveforms ; also notice that h_TS = - h_CH) are in the following tar ball files: mbh1a_st_hs.tgz and mbh2a_st_hs.tgz. After untaring them (tar -zxvf file.tar), you will find six files corresponding to the waveforms in the time and frequency domain for the detectors I and II. In the time domain, the strain is real (ampl*cos(phase)), so the 1st column corresponds to the time values and the 2nd one to the strain. In the frequency domain (computed analytically using the Stationary Phase Approximation) the 1st column corresponds to the frequency values, the 2nd for the real part and the 3rd one for the imaginary part. We have also added two files (*_allharm.dat) separating the contribution of different harmonics, so 1st col. correspond to frequency values, 2nd & 3rd to real and imaginary part of j=1 harmonic, 4th & 5th (j=2), ..., 12nd & 13rd (j=6) and 14th & 15th (total). In our original parameter estimation code we only worked with the waveform in the frequency domain (computed analyticaly using SPA); we modified a little bit that code in order to also obtain the waveform in the time domain, which is computed directly from Blanchet expressions, NOT as an invFT of the other ones. See comments and plots comparing our results with Cornish & Hughes’ ones in their section (now everything matches pretty well between CH and ST, we were wrong computing the waveform in the time domain and now we also have included all the harmonics in the frequency domain waveform).

• SNRs in the RWF case. The results we get for the SNRs working with the restricted waveform are in the following file st_SNRs.txt (for better comparison, we also have written, between brackets, Sathya’s results). We agree in the frequency bounds (Fbegin and Fend), but not in the SNR values (theirs are always higher than ours). Due to simplicity of this exercise, what we have done is to create a Mathematica notebook which is doing the same calculations as our Fortran/C code in the computation of the restricted SNR. Here you have the files for both cases: 1a_case.nb and 2a_case.nb. With this, we are making public all the details of how we are computing the restricted SNR and therefore everybody can look inside, change parameters, ... and look for the differences with their procedure. We hope this will help to start getting similar results, at least in this simplified exercise.
• Here are our SNRs going from 0PN to 2PN in amplitude (working always with 2PN in phase): SNRs_allpn.txt.

• This is the draft we are sending to CQG for the proceedings of GWDAW 12. Any suggestions or comments are welcome.

• We thought it would be useful to take a step back to the restricted waveform and first compare SNRs as well as beginning and ending frequencies. Alicia and Miquel provided us with their Mathematica notebook for the RWF and we found a mistake in the way the PSD is coded up in our notebook. After correcting this error, the SNR discrepancies between our results and those of Alicia and Miquel go away.

• After correcting our PSD, we find that our SNRs for FWF are now closer to those of Miquel and Alicia but still not quite identical (the Sintes & Trias SNRs are higher than ours, typically by about 7 percent, as seen in SRNs.txt). This difference can only come from the waveform itself. We compared SNRs going from 0PN to 2PN in amplitude in steps of 0.5PN (keeping the phasing fixed at 2PN): Orders.txt. The first non-trivial amplitude order where we can make a comparison with Alicia and Miquel is 1PN, and we already see a small difference there. This may well be related to a technical detail in the way the amplitudes are handled. When the SPA is taken, 1/sqrt(Fdot) appears in the amplitudes (where Fdot is the time derivative of the orbital frequency), and this needs to be expanded out and truncated to the desired order. However, the amplitudes also have other PN contributions. Strictly speaking, in the end one should re-expand the amplitudes as a whole and truncate them at the desired order, and this is what we do.

• Comparison between AISSV-mod and TS. We modified our code in order to do exactly the same truncation as Miquel and Alicia are doing. These are the results we get: TS-AISSV.txt, and all numbers match perfectly with TS’s ones.

• Update: February 26th, 2008. Comparison between all the groups. Coming back to our original code, we have computed both SNRs and errors for the cases proposed at the last telecon: mbh1a with high frequency cut-off of 1e-4 Hz and mbh2a with 1e-3 Hz. We also are using the new phic values that made TS’s and CH’s waveforms to match. We want to compare our results with our groups: AISSV-CH-TS.txt.

• The higher harmonic code developed by Cornish & Porter returned the following parameter errors: mbh1a_cp.txt, mbh1b_cp.txt, mbh2a_cp.txt, mbh2b_cp.txt. The Fisher matrix entries are also quoted, as are the Fisher correlations, Γ_ij/√(Γ_iiΓ_jj). The results match the independently developed Cornish & Hughes code to very high accuracy.

External Link

The MLDC version of the Barack-Cutler kludge waveforms provide a good starting point for the comparison of EMRI parameter estimation codes. The most useful data sets for this purpose are the Challenge 1.3.x and 1B.3.x training data, as the training data sets come with noise-free versions of the waveforms. These noise free waveforms are very useful for checking codes. The MLDC EMRI waveforms are fully documented here.

The EMRI search code developed by MtGWA group has been used to perform Fisher and MCMC studies of the Challenge 1.3.x training data. The following link is copied from the MLDC taskforce wiki:

* 2007/5/3: First look at EMRI posterior distributions for 1.3 training data EMRI notes

In those notes the explicit values of the Fisher matrix uncertainties were not listed. Using the key file parameter values as input, here are the Fisher uncertainties for Training 1.3.1. ( See revised results from 2008/4/8 below )

For reference, the fitting factor for the MtGWA EMRI code (which uses the low frequency approximation to the LISA response and interpolated waveforms to cut template generation times to less than a second) is 0.9993 for the 1.3.1 training data.

* 2008/4/8: It looks like the numerical central differencing offset of 1e-6 used in the Fisher computation may have been too large for some of the parameters. This quantity has been reduced to 1e-9 (consistent results have been found for 1e-10 and 1e-11), and the Fisher matrices have been recomputed. Here are the results for MLDC training 1.3.x: Training 1.3.1, Training 1.3.2, Training 1.3.3, Training 1.3.4, Training 1.3.5. It is worth noting that these values do not agree with the MCMC results (at least for 1.3.1), so more investigation is needed. The MCMC runs are currently being repeated.

* 2008/4/10: The problems with the Fisher estimates continue. I have increased the number of samples in the waveform generation and gone over to 4th order accurate numerical differencing (up from 2nd order accurate). The problems have not gone away. Over a range of numerical differencing parameters epsilon, the entries in the Fisher matrix are fairly stable - at most ten percent differences in the values. But this is enough to change the error estimates by orders of magnitude. Note that only the entries involving the mass, spin, eccentricity and lambda angle are significantly affected. The extrinsic parameters are stable. Here are some results from the training 1.3.1 case using the 4th order accurate numerical differencing epsilons of 1e-8 and 1e-9. Note that the Fisher matrices Fisher matrix 1e-8 and Fisher matrix 1e-9 differ by at most 6% in any entry, while the parameter errors Sigmas e-8 and Sigmas e-9 differ by orders of magnitude in some of the parameters. The reason for this can be traced to the large condition number of the Fisher matrices, which makes the inversions very sensitive to small changes.

* 2008/8/28: PROBLEM SOLVED. Motivated by today’s telecon, I returned to the EMRI Fisher matrix problem and found a simple solution: The eigenvalues of the Fisher matrix come in two sets, one set of 5-6 has very large eigenvalues (generally related to the intrinsic parameters), while the rest have very small eigenvalues. Some time ago I tried just deleting the small eigenvalue terms (in other words, using an SVD approach), and while that helped with the intrinsic parameters, it messed up everything else. The eventual fix is very simple: when computing the errors in the intrinsic parameters, delete the small eigenvalue contributions. When computing the errors in the extrinsic parameters, keep all the eigenvalues. This completely fixed the instability problem, and gives results in excellent agreement with the new MCMC runs. See the attached graph for a comparison of the Fisher matrix prediction in blue and the MCMC derived posterior in red for MLDC source 1B.3.1.

17/1/2008: We also have a Monte Carlo search code which we have been using for the MLDC. This employs an approximate waveform model, which is based on Barack and Cutler but with some simplifications to speed up waveform computation. As a first investigation, we are computing posteriors for the Round 1.3 training data sets. We have preliminary results that use a normalized waveform search code, i.e., with maximization over distance rather than including distance as a search parameter. Results from this code for training set 1.3.2 are here.

We are now debugging a search code that uses un-normalized waveforms, and should have results from this for all of the Round 1.3 training sets within two weeks.

13/2/2008: We now have results that can be directly compared to those of Neil Cornish, i.e., these are distributions over the full parameter space, and use parameters at t=0 rather than parameters at plunge. At present we still have results for 1.3.2 only, but are running the code for the remaining four data sets now. Results for training set 1.3.2 are here.

28/2/2008: We now have results for the other four training data sets, computed in the same way as the 1.3.2 results described above. A document summarizing results for all five training data sets is available here.

09/04/2008: We have carried out parameter estimation using the start and plunge parameter Fisher matrices for all of the 1.3 training data sets. We find that the calculation of our Fisher matrix is extremely sensitive to the parameter offset. A document detailing our parameter error and convergence tests is available :here. We need to run some MCMCs to further investigate the Fisher matrix for each challenge data set.

22/05/2008: We have tried to carry out dumb MCMCs using uniform priors from the most optimistic/pessimistic error estimates from our previous analysis. We confined ourselves to the Challenge 1.3.2 and 1.3.4 training data sets. We found that for 1.3.2 our pessimistic estimate was too large by orders of magnitude, and our most optimistic estimate was too small by a factor of 3-5. For 1.3.4 we found that the most optimistic error estimate was too small by over an order of magnitude. We have now calculated the FIM as the second derivative of the expectation value of the log likelihood. While this FIM is more convergent, it still looks to have issues of ill-conditioning. Furthermore it is expensive to compute so its use is limited. A document describing this can be found :here.

2/10/2008: We have done some investigations to test Neil’s suggestion that using the full Fisher Matrix to estimate extrinsic and phase parameter errors, while using SVD to compute intrinsic parameter errors might yield robust and accurate results. A document describing these results can be found here. We find results that are consistent with Neil’s suggestion for sources 1.3.1-1.3.3, but it is not so clear for 1.3.4-1.3.5. We are now looking at only the intrinsic-intrinsic and extrinsic-extrinsic sub Fisher Matrices to see if these yield useful estimates. As yet, we have not compared our results to the numerically computed posteriors obtained earlier.

These are preliminary notes on black hole formation models by Marta and Emanuele, including some references and topics for discussion.

Preliminary results on spin evolution due to mergers and accretion are presented here. This is a short paper summarizing the main results.

Here are corrected histograms obtained by running Scott/Neil’s MBHB PE code (including both higher harmonics and spin precession) on Marta’s model universes. (A previous version had failed to convert from dN/dV to dN/dt.) The document now shows 2 cases, corresponding to high and low spins. 3/26/08: Actually it turns out even the above results are a little off. Due to a typo in an email and one data file that was not completely read, the parameter distributions for the above results are a little bit off. By eye they looked fine. In any case, those problems were fixed before the following results re two descope options were generated.

Here are MBH param estimation results that were generated in response to a request by John Morse (head of Astrophysics at NASA) that the LISA Project investigate the cost/science loss of descope options. The U.S. LIST fixed on the following 3 noise curves–baseline, descope1, and descope2. For the MBH part of the report, I used the following 3 source distribution models generated by Marta Volonteri. Results from running Scott/Neil’s code on these 3 distributions for the 3 LISA configs are summarized in this MBHB results document. I should note that this MBH work was part of a much larger effort, especially by Neil Cornish and Jon Gair, to look at the effect of descopes on Galactic binaries, EMRIs, and MBHBs.

Update: May 19th, 2008. Marta asked about the parameter distribution of detectable sources. This document gives histograms of the physical parameter distribution for sources with total (A+E) SNR >10, for the baseline noise and the iso-eff source distribution. [Note that by mistake an earlier version of this document plotted the redshifted masses. The histograms now plot the locally-measured masses.] Also, Miquel asked about how angular resolution and distance resolution depend on source redshift. That’s shown in this accuracy versus distance,(updated June 13) which is again just for the iso-eff distribution, and 2 noise curves: baseline and descope1. I can make more plots like this, but there are hundreds of plots one COULD make, so want to know which ones people want the most.

Results document for the Proceedings, Oct.’08 Curt is putting results from running Scott/Neil’s code on Marta’s 4 models here(results). For the Proceedings, we agreed to only consider 2 noise curves: “baseline” and ‘6link’, where these designations have the same meaning as in the current working version of the LISA ScRD. The input distribtions–i.e., distributions of masses, spins, redshifts,etc.–for Marta’s four models are here(input distributions).

Here is the final draft of the proceedings. This is the source file in tgz.

Please email you acknowledgments (if missing), comments and corrections to (berti at wugrav dot wustl dot edu)

If you need the CQG style files, they are in this tgz file.