Question about evolving the binaries

[AstrophysicsAndPython]

While evolving the binaries using legwork.evol.evol_ecc the parameter t_before is used to stop the evolution sometime before the merger happens e.g, 1 million years (by default). My questions are,

1. Why take the value of t_before as 1Myr when we will only have LISA for 4 to 10 years.
2. If the evolution stops a million years before the LISA mission is deployed how are the million years compensated in LEGWORK.
3. The legwork.source.Source class doesn't include any time boundary, does that mean it will evolve the binary till the merger happens? If yes, how is the SNR calculated for those that merge way before LISA will be deployed?
4. What is the correct path to evolve binary, directly from source class using COMPAS DCO parameters or going through evolution -> source class?

Thank you.

[TomWagg]

Thanks for the questions 🙂

1. Even though LISA is only 4-10 years, users may also want to evolve sources before the LISA missions (e.g. from birth in the Milky Way until present day) and so they could be evolving for billions of years. The choice of 1 Myr presents numerical integration warnings for binaries with eccentricities below 0.95.
2. It stops 1 million years before the merger, not necessarily 1 million years before the LISA mission
3. Yes you are correct. Any source that merges before the LISA mission is given an SNR of exactly 0. But good point, I hadn't realised we have a different default choice of avoid_merger in Source and evol.evol_ecc. I will discuss this with @katiebreivik and make them match up.
4. It's really up to you how to evolve sources. You can instantiate a full Source class and then evolve them from there. But if you're only interested in evolution then Source will contain extra things that you may not need and in this case you may prefer to use evol.evol_ecc directly.

Let me know if you have further questions!

[AstrophysicsAndPython]

Hi, thanks for your quick reply and sorry for my late response. It does clarify everything I asked, but I don't understand the integration warning problem, sorry. The other thing I also noticed is that if you increase the n_step variable in evol.evol_ecc function, the SNR increases which results in an increased number of detections.

For example, consider the binary

m_1 = 15.48664*asu.M_sun
m_2 = 1.72879*asu.M_sun
e_DCO = 0.21275
a_DCO = 0.02773*asu.AU
lookback_time = 11.99873*asu.Gyr
t_evol = 7.22095*asu.Myr
dist = 9.49449*asu.kpc

e_today, a_today = evol.evol_ecc(ecc_i=e_DCO, t_evol=lookback_time - t_evol, m_1=m_1, m_2=m_2,
                                 a_i=a_DCO, output_vars=['ecc', 'a'])
e_today, a_today = e_today[-1], a_today[-1]

snr_Source = sources.Source(m_1, m_2, ecc=e_today, dist=dist, a=a_today).get_snr_evolving(t_obs=4*asu.yr)[0]

Which results in SNR ~ 0.065647. However,

e_today, a_today = evol.evol_ecc(ecc_i=e_DCO, t_evol=lookback_time - t_evol, m_1=m_1, m_2=m_2,
                                 a_i=a_DCO, output_vars=['ecc', 'a'], n_step=int(1e4))
e_today, a_today = e_today[-1], a_today[-1]

snr_Source = sources.Source(m_1, m_2, ecc=e_today, dist=dist, a=a_today).get_snr_evolving(t_obs=4*asu.yr)[0]

This results in SNR ~ 7.198972

[TomWagg]

Integration warnings First I can explain the warnings. These warnings were cropping up a lot previously (see this issue). Essentially every time you evolved a binary near a merger you'd get ~100 warnings that crowded the terminal. These warnings were essentially the adaptive timestepping code saying that the derivatives near to mergers are so large that the computer can't take small enough timesteps to change and ended up repeating lots of timesteps which is wasteful.

In order to fix this we added the avoid_merger feature which prevents the integration from getting too close to the merger and hence stops these warnings from appearing.

Your timestepping issue This issue is a little subtler but is essentially because you're evolving this source for far longer than its merger time.

Note that the merger time of your binary is ~0.3 Gyr. But you're evolving it for around 12 Gyr. This means that LEGWORK will choose 100 timesteps equally spaced from 0 to 12 Gyr. However after 3 timesteps it will have merged and so the "final" eccentricity and separation that you get from the evolution is just the values from the second timesteps. In contrast, if you tell LEGWORK to use 10,000 timesteps then you can get much closer to the merger before the timesteps become useless.

The results of this is you end up with a much tighter binary (a_today is much smaller) when using more timesteps. But in reality this binary has already merged so both SNRs are wrong, they should both be 0.

You highlight an important point which is that evol_ecc doesn't do a great job of telling you whether you evolve a binary past its merger. I'll add a warning to prevent this in future, thanks!!

Does that help? p.s. thanks for including the code in your question, was very helpful for testing!

[TomWagg]

Oh and also I think that if you had evolved these sources using the Source class then this issue wouldn't have occurred because it handles whether binaries have merged. Though to be fair, evolve_sources doesn't let you choose a number of timesteps so you can't test this sorry haha 🤷🏻

[AstrophysicsAndPython]

The results of this is you end up with a much tighter binary (a_today is much smaller) when using more timesteps. But in reality this binary has already merged so both SNRs are wrong, they should both be 0.

This was the primary motivation for the investigation :-)

I'm also curious about the code where the time selection choice is made in snr.snr_ecc_evolving

if t_merge is None:
    t_merge = evol.get_t_merge_ecc(m_1=m_1, m_2=m_2, f_orb_i=f_orb_i, ecc_i=ecc)
t_before = 0.1 * u.yr
t_evol = np.minimum(t_merge - t_before, t_obs).to(u.s)
e_evol, f_orb_evol = evol.evol_ecc(ecc_i=ecc, t_evol=t_evol, n_step=n_step, m_1=m_1, m_2=m_2,
                                   f_orb_i=f_orb_i, n_proc=n_proc, t_before=t_before, t_merge=t_merge)

Considering the above example again,

m_1 = 15.48664*asu.M_sun
m_2 = 1.72879*asu.M_sun
e_DCO = 0.21275
a_DCO = 0.02773*asu.AU
lookback_time = 11.99873*asu.Gyr
t_evol = 7.22095*asu.Myr
dist = 9.49449*asu.kpc

The binary, when at the DCO stage, has merge_time = ~0.3Gyr. According to the code, the t_merge is calculated first (0.3Gyr) then it is compared with t_obs (4 yr) and the minimum of both is selected, which is 4 yr. Later, the binary is evolved for t_evol e.g, 4 yr. What is the reason to take this minimum of the two at the DCO stage? Would it not essentially evolve every binary for 4 years?

[TomWagg]

This was the primary motivation for the investigation :-)

Haha I see, nice test in that case!!

The reason for taking the minimum is that the evolution in the SNR function only handles the evolution during the detector mission. For LISA this is 4 years, so you only want to integrate the SNR over the 4 years of the LISA mission (when observations are taking place), not over the hundreds of millions of years as the binary inspirals. So yes it does essentially evolve every binary for 4 years, but this is actually what we intend in this case. If you want to also evolve the source until present day you'll need to do that before calculating the SNR.

[TomWagg]

@AstrophysicsAndPython I'm going to close this for now but feel free to re-open/open a new issue if you have further questions! 🙂

[AstrophysicsAndPython]

Hi, thank you for your previous answers. I do have some other queries regarding LEGWORK as well (I'll probably open separate issues for them) but as of right now I'm getting a little bit confused by the lookback time and the distances from galaxy.py. It might be an extremely basic question, but here goes,

For a galaxy I generated, these are the particular values for a single point

lookback = 8.42769448*u.Gyr
distance = 8.81452481*u.kpc
Z = 0.01784518

Lookback time The time it takes for the light from the object to reach the Earth.

Another definition from here is that lookback time is the difference between the age of the universe for today and the age of the universe when the object emitted the light.

So we're seeing the object as it was ~8.5 Gyr according to definition 1, or 5.3 Gyr according to definition 2 considering the age of the universe to be 13.8 Gyr?

distance Distance of the object from Earth.

So the object is at ~8.8 kpc, but wouldn't that mean that we'll be seeing the object as it was 0.028 Myr ago?

Am I missing something?

[TomWagg]

I assume this is galaxy.py from my other repository?

Definition 2 is correct: lookback_time = universe_age - birth_time and so yeah a lookback time of 8.5 Gyr means that the system was formed ~5.3 Gyr ago. Note you can set lookback=False in the simulate_mw function to get birth times instead.

I don't quite understand your last question, the distance doesn't really affect at what time we see it (for the Milky Way at least, it has an effect at large distances).

[AstrophysicsAndPython]

Hi, sorry for the confusion.

The last question was related to the time it takes for the GW to reach us from the source.

As the source is at ~8.8 kpc and considering the finite speed of Gw propagation, c, it was just related to that time. I got it confused with look back time, but it is clear to me now.

Cheers!

[TomWagg]

Gotcha, glad that's it clear now.

Just wondering, @AstrophysicsAndPython would you find it useful to add a module to LEGWORK that dealt with placeing the sources in a galaxy model? That's something I would consider if people were interested.

[AstrophysicsAndPython]

Hi,

Can you explain it a little bit more? Is it related to the placement of the DCO in the galaxy model related to some parameter of the binary?

[TomWagg]

Sure! I was imagining setting it up so that set up the Source class as usual, but can then apply a galaxy model. I.e. it would give each source a random position/birth time/metallicity and proceed like that.

Additionally, I'd make it so that users could call the galaxy model stuff without the source class. E.g. you'd run something like simulate_mw(model="some model") and then it would hand back randomly generated positions, birth times and metallicities.

[AstrophysicsAndPython]

1. If I'm understanding this correctly, you want the input binaries in the Source class to get their random position/birth time/metallicities from your galaxy.py.
   • However, as (in my case) we're taking the binaries from COMPAS and they may already have associated metallicities with them. So reproducing the exact metallicities that were generated randomly by COMPAS, to me, sounds a bit difficult task.
2. I didn't understand this point, is this about another model for MW, or is it about keeping the galaxy.py module separate from the Source class? I mean, it can be called from within the Source class (point # 1), but if the user wants it can be called separately as well.

Coming back to your point # 1, I did something like this for a paper I'm working on currently by comparing the metallicities from the COMPAS with those present in the simulated milky way. It did work out, but it was a bit slow for a million points in the galaxy :laughing: And make it run over several galaxies, it was stressful. So, yes, that would be very helpful indeed and at least we're interested in this :-)

[TomWagg]

Gotcha good point, I can also add an option to add the metallicities manually to hopefully make that easier! Yeah it's like another model of the MW, I imagine there are lots of good options other htan mine that we could add.

I'll add an issue to the main list to gauge interest from others as well, thanks for your thoughts! :)

[AstrophysicsAndPython]

Hi, I'm curious about a particular mask you have created for your data,

https://github.com/TomWagg/detecting-DCOs-in-LISA/blob/fa7da9b60c63805ea54fb88e5fa58a61687a20dc/simulation/src/simulate_DCO_detections.py#L247-L248

Here you're taking the inspiralling binaries to be those that have a merger time greater than the difference of their lookback time and evolution time. But isn't lookback time the opposite of what we want?

For a binary to have a lookback time of 8 Gyr, the binary is formed only ~6 Gyr ago.

Definition 2 is correct: lookback_time = universe_age - birth_time and so yeah a lookback time of 8.5 Gyr means that the system was formed ~5.3 Gyr ago. Note you can set lookback=False in the simulate_mw function to get birth times instead.

If the merger time is greater than 8 Gyr, aren't we over evolving it?

[TomWagg]

Ah, I messed up with what I wrote previously. Definition 2 is correct but I then misspoke. The correct version is: a lookback time of 8.5 Gyr means that the system was formed ~5.3Gyr after the start of the universe., or equivalently, it formed 8.5 Gyr ago.

Hopefully that makes the condition make sense?

[AstrophysicsAndPython]

Hi,

Sorry for taking too long for coming back. Regarding the mask I mentioned in my previous comment, why evolve the binary for n_step = 2. Isn't it essentially breaking the evolution time into two parts, placing the binaries midway through evolution and after that calculating the SNR for only 4 years with the source class?

The reason for taking the minimum is that the evolution in the SNR function only handles the evolution during the detector mission. For LISA this is 4 years, so you only want to integrate the SNR over the 4 years of the LISA mission (when observations are taking place), not over the hundreds of millions of years as the binary inspirals. So yes it does essentially evolve every binary for 4 years, but this is actually what we intend in this case. If you want to also evolve the source until present day you'll need to do that before calculating the SNR.

[TomWagg]

The n_step = 2 only sets the output timesteps, internally it will still take the necessary steps to resolve the evolution. But by taking n_step=2 we get only the initial and final state in the output (and we only care about the final state in this case, hence the mask afterwards)