Should PN code take angular velocity as input?

[vijayvarma392]

@dongzesun, @moble, @duetosymmetry:

In the PN code here (the one being put in sxs is not exactly the same, I believe, but close): https://github.com/dongzesun/HybridizationWaveforms/blob/3a854a860e6efc13d1368bfc0ccdb38b1da6b83c/PYPostNewtonian/Code/PostNewtonian.py#L14

What is this input? Is it actually angular velocity or orbital frequency? They are not the same, right? I would have thought you need the latter for initializing PN.

In any case, here's an example comparison between the spins of NR, LAL PN code, and Dongze's PN code. Here I'm initializing the PN codes using the orbital frequency in the coprecessing frame at the first time shown.

Zooming in:

There seems to be an offset between Dongze's PN and the other two. Is this because I should have used angular velocity as input?

[vijayvarma392]

GWFrames seems to ask for orbital frequency (I think): https://github.com/moble/GWFrames/blob/master/Code/PNWaveforms.cpp#L43

[dongzesun]

The omega in my code is actually v^3, where v is the PN parameter. Should we take v as the input instead of omega?

[dongzesun]

I'm using the same omega with GWFrames, which is v^3. This is just an approximation of angular velocity

[vijayvarma392]

Hmm, ok, then maybe it should just be called orbital frequency instead?

[dongzesun]

Yeah I'll make it accurate. It should be scaled by M

[vijayvarma392]

Here's a snippet of code to reproduce this:

import numpy as np
import matplotlib.pyplot as P
import quaternion

import lal
import lalsimulation as lalsim
from lalsimulation.nrfits.pn_spin_evolution_wrapper import spin_evolution

import os, sys
sys.path.append('../HybridizationWaveforms/')
from PYPostNewtonian.Code import PostNewtonian 

def PN_Dongze(q, omega0, chiA0, chiB0):

    M = 1.0 # total mass
    frame = quaternion.quaternion(1,0,0,0) # frame at reference time

    # At this point w is in corotating frame
    w, chiA, chiB = PostNewtonian.PNWaveform(q, M, omega0, 
            chiA0, chiB0, frame, return_chi=True, ForwardInTime=False) 

    t = w.t
    chiA = quaternion.as_float_array(chiA)[:,1:]
    chiB = quaternion.as_float_array(chiB)[:,1:] 

    return t, chiA, chiB

q = 1.4619394167969584 
chiA0 = [-0.1627143,  -0.08010382, -0.36227616] 
chiB0 = [0.24580962, -0.43054159, -0.03076978]
omega0 = 0.008336652403608003

dt_pn = 0.1
omega_pn, orbphase_pn, chiA_pn, chiB_pn, LNhat_pn \
        = spin_evolution(q, chiA0, chiB0, omega0, dt=dt_pn,
                         approximant='SpinTaylorT1', spinO=6, phaseO=7)
t_pn = np.arange(len(omega_pn)) * dt_pn

# This code only allows TaylorT1
t_pn_d, chiA_pn_d, chiB_pn_d = PN_Dongze(q, omega0, chiA0, chiB0)

spin_idx = 0
P.plot(t_pn, chiA_pn.T[spin_idx], label='PN LAL')
P.plot(t_pn_d, chiA_pn_d.T[spin_idx], label='PN Dongze')
P.xlabel('t')
P.ylabel('chiA'+['x','y','z'][spin_idx])
P.ylim(0.16, 0.26)
P.xlim(6500, 11500)
P.legend()
P.show() 

And it should make this plot (not showing NR):

I'd normally say, our PN is better than LAL's, but LAL's seem to be closer to NR.

[dongzesun]

Thanks, I'll try to figure out this.

[vijayvarma392]

Another thing I noticed is the biggest difference between the two PNs occurs at the turnover points for the spins. Maybe this means the LAL code has some extra spin terms we don't?

[dongzesun]

I'm comparing the source code. Are you using the same PN order?

[vijayvarma392]

Hard to say, I'm using the max PN order they have (see the options in the snippet). I don't know exactly what you have :)

[dongzesun]

Does spinO=6 means 6 PN? Then they definitely have more terms. The default in my code is 4 PN.

[vijayvarma392]

No, that's twice the PN order.

Here's the documentation: https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_inspiral_spin_taylor_8c.html#ace0a712ca64820416f1f3457880caf05 which is not super great. This function also does the same thing, and has paper links: https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_inspiral_spin_taylor__c.html#ga35cfdf3082e09cc97cda9e11ba4c2bff So, I'd look at both of these for documentation.

[vijayvarma392]

Actually, this one as well: https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/_l_a_l_sim_inspiral_spin_taylor_8c.html#a4b7ae88c1f139d87d55f10baf7a3d8b2. This is Riccardo Struani's fault, he's terrible at documenting his code.

[vijayvarma392]

What happens when you run my snippet with 3PN spin order for your code (I don't know how to do that in your code)?

[moble]

GWFrames seems to ask for orbital frequency (I think): moble/GWFrames@master/Code/PNWaveforms.cpp#L43

The documentation that should show up in python is here, where it's described as "orbital angular frequency".

Is it actually angular velocity or orbital frequency? They are not the same, right?

I'm not sure if this is what you're asking, but just to clarify: In GWFrames, I always used Omega for orbital angular velocity/frequency and omega for waveform angular frequency, which are — of course — quite different. And I always used angular frequency as opposed to "temporal" frequency, if that's what you mean. It's also true that (with G=c=1), $(M\Omega)^{1/3}$ is defined as v in the PN literature, which is what's actually used in GWFrames. (Note that the total mass is always 1 in that function.)

Maybe this means the LAL code has some extra spin terms we don't?

Quite possibly. GWFrames was very complete for its time, but I haven't kept it up at all in years. Also note that I would tend to just include as much as I could, and not worry about terms that were missing in the literature. So you could ask for 4PN, and I'd just give you as much as I could up to that order, even if — say — spin terms were only complete to 2.5PN or something.

[dongzesun]

Ok, thanks. I'll compare the terms up to 3 PN to see if they are different

[vijayvarma392]

Ok, so the code uses "angular velocity" to mean "orbital angular velocity/frequency", which I'd define as Omega=dPhi/dt, where Phi is the phase from orbital trajectories. I suppose I've been shortening "orbital angular velocity/frequency" to "orbital frequency" and you've been shortening it to "angular velocity", thus making sure we have zero overlap in the four possible words, classic physicists. Same page about this?

[vijayvarma392]

Here are Riccardo's notes about what is in the LAL code: https://dcc.ligo.org/LIGO-T1500554/public

[dongzesun]

Yeah we are on the same page

[dongzesun]

What happens when you run my snippet with 3PN spin order for your code (I don't know how to do that in your code)?

They are not the same even for same PN order. The LAL code uses equation (21) in https://dcc.ligo.org/public/0122/T1500554/023/dLdS.pdf, while GWFrames and my python code uses equation (4.5) in https://arxiv.org/pdf/1212.5520v2.pdf. I need to have a look at the derivations.

[dongzesun]

The LAL code uses terms in equation (A1) and (A2) in https://arxiv.org/pdf/1501.01529.pdf, they should be complete at 2.5 PN order. But we only have terms up to 2 PN orders.

[vijayvarma392]

Some more information: It's not just the spin evolution (maybe), the waveforms are also different between the two codes. While waveform differences could arise from the differences in spin evolution, I suspect it's coming from differences in amplitude corrections.

So, first, here's the real part of (2,2) and (2,1) in the inertial frame for one of the HybTest_rerun cases, where I compare lal_pn and Dongze's pn:

Notable thing is that (2,1) seems better for lal_pn.

Next, amplitudes:

Notable things are (2,2) amplitude is way better for Dongze (which is why I suspect differences in amplitude corrections). (2,1) amplitude is doing weird things for Dongze.

[keefemitman]

So as for why the (2,1) mode is different, this is because Dongze's code includes the 1.5 and 2 PN memory contributions to the (2,1) mode that we computed in this paper. Not sure what's up with (2,2) though. Maybe that's because LAL includes higher PN order precession terms?

[vijayvarma392]

Ah, interesting. @dongzesun can you add an option to your PN code to turn off the memory contributions Keefe mentions, but set them to be on by default? Then, I can repeat this plot and see what happens. Also, @keefemitman, seeing as Dongze's (2,2) amplitude is closer to NR, this might suggest he has more amplitude corrections than the LAL code (but is missing spin terms potentially).

[dongzesun]

Ok I can turn off the memory terms.