Similar to how Adrian has set up indexing into Orbit classes, I'd like to add a feature to slice into a Population so that you can quickly just access a particular set of bin_nums and get the subset of every table etc.

(This will make my life easy for grabbing things like "just the OBA stars", "just the NSNSs")

We'll need a way to save these to files for larger runs. Should be able to create separate tables within a single hdf5 file

It would be nice to actually create some sphinx docs, should be easy enough to copy from LEGWORK

At some point it could be interesting to consider how we emulate cluster dynamics with some simpler models in our simulations. Adrian has some experience doing this with his stellar streams and so we may be able to do something like that.

The current mask seems to make everything detectable as Zari+19 requires G < 16 - is this what we want?

Based on Mathieu's papers, it would be good to repeat the work from #39 but for variations on accretion efficiency and critical mass ratios.

Run a fairly large population through the entire code to get the recent history of the Milky Way.

• Restrict lookback times to < 200 Myr
  • This will limit us to the low-alpha disc based on Neige's model
• Remove any sort of mass cut
• Binary fraction = 1.0
  • This still bothers me a bit, since we're allowing any mass

It may be useful to actually

• estimate how large the population should be (or at least how much star forming mass came) from the last 200 Myr

It would be nice to be able to evolve single stars instead.

Perhaps the most convenient way is just catching if binfrac=0.0 and doing a different thing if that is the case.

Neige's galaxy model can give us metallicities outside of the allowable COSMIC range [1e-4, 0.03] - stop that!

Need to do #46, #47 and #48 first

Seg faults seem to be occuring whenever COSMIC produces a binary with NaNs and then Gala tries to get its orbit.

I should check if any binaries have NaNs before sending them to Gala (and also investigate why that's happening).

This'll be a question for @adrn.

I'm seeing that some low-mass binaries are managing to get ejected from the galaxy. Perhaps some sort of interaction is happening in the centre in gala?

Here are some plots of the orbit of one particular binary at 50 Myr, 200 Myr, 2Gyr and 10Gyr.

• #13
  • The big loop ends up being super slow but could be easily thrown into a pool
• #14
  • I just make a list of lists of strings currently - this is very human-readable but harder to create masks
  • I'm thinking a large dataframe of booleans might be the best way
  • We can discuss during our meeting on Wednesday

Maybe I can use tqdm here to check the progress?

This'll be a question for @katiebreivik.

I'm currently evolving every binary in COSMIC for 13.7 billion years - this is not correct given that they each have different birth times. We should account for this. Either:

1. Maybe there's a method in COSMIC to just set the max_ev_time individually 🤞🏻
2. Or, we do something sneaky like force bcm to include a timestep of an exact time for each binary perhaps?

Here's the repo! For comparing to Gaia we need to account for all sorts of things like dust that make measuring reality hard.

We currently just have met=0.02 for COSMIC. It would be easy enough to edit InitialBinaries with the drawn metallicities from whatever Galaxy model we use. Do it!

• Move events stuff to events.py
• Move evolving function into a method of Population in pop.py.
• Set initial galactic velocities in class also (probably add them to Population.initial_galaxy

Work on understanding querying of Gaia and getting out the runaways (Katie will help with some code I believe).

I'm unsure what we're going to require for this. I imagine it will be useful to know how many massive objects are passing between us and Galactic centre. Perhaps integrating the orbits forwards by a couple of decades would be useful for getting rates?

Need to think about this more

Where did merger events occur? We will need to recompute the orbits of the merging DCOs to get the exact location.

This is generally defined using effective temperature and absolute magnitude (=M_ZAMS) cuts. One option for working out what these should be could be to create a CMD for a grid of COSMIC data.

As a simpler first pass we could instead try to just do it based on the CMD directly and choose some cuts there (#19 will be useful for this of course).

E.g. if m1_cutoff is changed then they should really re-run the whole population to make it consistent. We probably don't want to erase previous data automatically so we should flash up a warning if things like this get changed.

Question for @katiebreivik - does COSMIC handle single degenerate Type IAs at all? (Could be interesting from a runaway perspective said @mathren)

If you do p[[1, 2, 2, 3]] you only actually get a population of 3 binaries - this might be useful but it seems misleading - need to fix the mass to account for this

Katie has a bunch of useful code for getting Gaia observables from our simulated binaries. This will account for dust extinction and bolometric corrections and stuff. Seems like I'll just need to find some MIST isochrones and it's all there to work.

• Integrate Katie's code into the package (observables.py? easy access functions in Population?)
• Grab those isochrones from somewhere

The result of all of this should be some magnitude of the object I think.

To compare with Gaia we need to know where these binaries are and how they are moving on the sky. Obviously from the simulation we know their Galactocentric positions and velocities. So we need some sort of coordinate transformation to get what Gaia sees...but I'm not entirely sure what that is (the concept of "sky position" for something not on Earth seems weird to me). I'll probably need help from @adrn on understanding this 🙂

• Transform Galactocentric pos/vel into what Gaia sees

Then once that works we can start making cuts based on what Gaia can see. The limiting magnitude is a function of sky position and we can cut out anything too faint

• Cut out binaries that are too faint for Gaia

How far are DCOs from Galactic centre?

To start: where is the population?

• v vs. v_z for observed vs. simulated OBA population

Currently I apply dustmaps to both the absolute and apparent magnitudes, as @katiebreivik pointed out we shouldn't be doing that for absolute mags

Use the Gaia IDs from the paper to acquire a sample to compare against.

We'll need a sample of Gaia O-star to compare to. This paper could be helpful in deciding how to do this.

For now it seems like a cut on the CMD could suffice?

E.g. test against Cygnus X-1 and some other types of XRBs to see whether results are realistic

To start let's implement the following classes for binaries. Keep track of the times at which the binary belongs to the class (e.g. maybe only XRB briefly).

• "Theory" Runaway
  • Any disrupted binary that has a post SN speed > 30 km/s
• "Observer" Runaway
  • Any disrupted binary with a final velocity relative to GSR > 30 km/s
• X-ray binaries (simple)
  • For now just any binary with a star and CO at any point
• Any binary with a BH or NS
• Any DCOs

We can do this in classify.py

Link to paper

• Create a table of scale heights and lengths
• Plot x vs z, split by "visible" vs. "underworld" populations
• Repeat everything for full single and binary populations

Now that the COSMIC PR is merged, we can hide that pesky zsun message 🙂

I need to set mass0_x when setting up the initial binaries!

I was chatting with Katie about this and she says we'll need to get an opinion from @adrn.

44% of systems in my latest simulation fall outside of the range of the bayestar dustmap. It seems that this map only covered "the three quarters of the sky north of a declination of -30°". The systems that are missed all have declinations lower than this.

What should we do about this? (Currently my default is to conservatively assume maximum dust extinction for all of these objects.)

• Is there an alternate dustmap?
• Katie suggested maybe mwdust could help?
• Could we do something smart for missing areas (maybe find an equivalent point reflected about the galactic plane??)

We know t_c and could integrate orbit (or assume it has escaped) to get the future merger locations

The mergers of two compact objects can produce a strong kick. Is this in COSMIC? Should we add it? It seems like it could be important for CO distributions (#42) as well as DCO locations (#46)

• Split into RLOF/wind-fed vs. Be-XRBs
• Plot locations of "strong" x-ray sources on the sky

In a grand reunion I'm going to be attempting to remember how healpy works (😅) and create maps for both our simulation and the Gaia observed O-star populations.

For each pixel we can create some simple statistics comparing the populations and see how this varies.

Should be able to use store_all parameter of integrate_orbit

As @katiebreivik suggested a while ago (and I am now remembering), it would be nice if Population could use ini files (since that's all the rage in COSMIC these days apparently 🙃).

I'm not entirely sure how this would work since I've not used the files before but I guess we would need to change __init__ and sample_initial_binaries to account for it.

I need to remember to ask @mathren about this. Currently I have it set up so that I only check for walk/runaways if the binary is disrupted. So for single stars I don't find any. But I figure this should still be possible with just a single supernova.

Just not sure if walk/runaways are fast moving stars, or fast moving stars that came out of binaries 🤷🏻

The for loop in identify_events should be thrown into the same multiprocessing pool that we use for the rest of the computation, there is a need for speed

Ideally we'd like to include the maps that Chiara suggested which are from these papers:

• https://arxiv.org/abs/astro-ph/0207156
• https://arxiv.org/abs/1610.09448

Will discuss these with @katiebreivik tomorrow.

There appear to be two weird cases in which not all disruption markers agree.

1. Supernova occurs during RLOF (~1/15000 in sample)
   • In this case the separation gets set to -1 and an evol_type of 11 appears...but the kick_info table does not include any lines that say disruption==1.0
2. Supernova kicks seems to lead to merger (~1/1000 in sample)
   • The separation gets set to 0.0 and no evol_type=11 appears...but the kick_info has rows that say disruption==1.0

The former seems to be a logging issue where the kick_info hasn't been updated correctly - though not sure about that.

The latter I can imagine could maybe occur if the supernova kicked one star into the other, but some of these natal kicks are ~700km/s and it seems insane to me that you get a merger out of that this often.

I can attach the initC table but I don't think that includes kicks? Will check with Katie tomorrow.