IMRPhenomD giving NaNs for equal mass spinless BHs

[ThibeauWouters]

Hi there!

I have been working on adding tidal contributions from NRTidalv2 to ripple, focusing on the IMRPhenomD waveforms as BBH baseline for the moment. While testing my implementation and sampling masses, I noticed that the IMRPhenomD waveform generator, gen_IMRPhenomD_hphc, returns NaN values for some mass values if enforcing the BHs to be of equal mass. The problem seems to arise from the equal mass constraint, since the problem does not arise if I nudge the smaller mass with a very small number, e.g. 1e-8. It seems that this issue arises in around 20% of cases that I checked by randomly sampling a mass from a uniform distribution between 1 and 100 solar masses (i.e. the range reported in the ripple paper0.

This issue arises if I install ripple from scratch in a fresh conda env (Python 3.10.13). I have not checked this observation for other Python versions yet. For convenience, I have written my observations in a small notebook, which you can find here to reproduce the result.

[tedwards2412]

Hi! That's great that you're working on NRTidalv2, this was actually going to be my next target waveform. Let me know if there is anything I can do to help.

My suspicion from looking at your notebook is that when solve for eta, sometimes (due to small numerical inaccuracies) it will give eta > 0.25 which will cause NaNs since eta is strictly bounded. It would probably be useful to put a check on this directly in the PhenomD evaluation so that users know this is happening.

[jacobgolomb]

Coming across this issue now; I am able to produce this with the following parameters:

{'psi': 3.14159075186615,
  'iota': 2.340374948484553,
  's2_z': -0.3226819680450276,
  's1_z': 0.9999654752321001,
  'M_c': 13.102260506719983,
  'ra': 4.648127372804028,
  'dec': -0.8857181758352013,
  't_c': 0.049999999999999975,
  'phase_c': 4.210371005335187,
  'd_L': 1937.002164088187,
  'eta': 0.24999999999999997,
}

Changing M_c and eta slightly makes it return a finite value. I tracked the issue down to this line https://github.com/tedwards2412/ripple/blob/main/src/ripplegw/waveforms/IMRPhenomD_utils.py#L41
where small numerical errors can cause the argument of sqrt to be negative (like -1e-16). This may need to be protected with a jnp.clip or something.