Mass (or other parameter) vs time plots like COSMIC

[Akarsh Sahay]

Hey,

I am an undergraduate student learning to use POSYDON for my summer project on Binary Population Synthesis. I wanted to know how can we produce graphs like Mass vs Time on finer timescales. I realise that using the pop.history, I can plot mass for discrete steps (ZAMS->ROFL->CC etc) of the evolution but is there no way to get a finer graph like in COSMIC here

The closest thing I found was this. But the input .h files to generate these graphs are not the same ones outputted when I simulate binaries using the way described here . The PSYGRID class doesn't create a grid for me when I give it a file like '' 1e+00_Zsun_population.h ". From what I can understand is that PSYGRID wants MESA grids, but population_runner.evolve does not output that for the simulated binaries. So is there no way to get a fine mass vs time plot for binaries like in COSMIC.

This would be useful in understanding how mass loss/transfer is occurring in my binaries.

Thanks

Akarsh Sahay

University of Oxford

[Jeff Andrews]

Akarsh,

Thank you for your interest in POSYDON. 

Your question is a good one. By design, POSYDON works fundamentally differently from COSMIC, which means there are a couple ways to answer your question depending on exactly what you would like to do. First, the background: POSYDON works by pre-computing binary evolution grids then interpolating between them. However, that interpolation currently only works from the initial to the final binary parameters, not the time series in between (although we are working on this). If you would like to know the time evolution of a particular parameter (say, the mass) then you will be limited to one of our pre-computed binaries. Alternatively, if you have a specific binary you would like to evolve, you could re-evolve it yourself. I am happy to help you with that if you'd like, but for now, let's assume you want the time series values for one of our pre-computed binaries.

You have two options. First, you can load up the specific PSYGRID .h5 file you are interested in and grab the specific binary you want time series data for:

from posydon.grids.psygrid import PSyGrid

path_to_lite_grid = os.path.join(PATH_TO_GRIDS, 'grid.h5')
grid = PSyGrid(path_to_lite_grid) 

grid.load()

 For each binary, we have a binary_history as well as a history1 (and a history2 for the HMS-HMS grid) which gives the mesa history files for that binary. See our documentation on the PSyGrid object for a complete description.

As a second option, you can run a population using the "nearest_neighbour" approach. This method will evolve binaries from start to finish, going through all the steps, but rather than using initial-final interpolation to calculate its evolution, it will map a binary onto its closest grid point in our binary grids as the binary evolves through its different stages of evolution. I do not recommend this method for accuracy, as our grid spacing is only so small, leading to different binaries mapping to the same trajectory. However, if you want to explore the evolution of binaries in more detail, this may be a good option. To evoke this option, make sure the interpolation_method parameter is set to nearest_neighbour in the MESA_STEP option:

# STEP CUSTOMISATION
MESA_STEP = dict(
    interpolation_method = 'nearest_neighbour'

)

This will require a somewhat custom evolution of the binary, which can be followed using the documentation here.

Does this answer your question? I am happy to help further. I recommend continuing this conversation by making a post on our GitHub Discussions page or by raising a GitHub Issue. Our developers are constantly working to improve the code. 

-Jeff

--
You received this message because you are subscribed to the Google Groups "posydon-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to posydon-user...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/posydon-users/f6e8b327-452c-42a6-af4a-baaa84cf40d1n%40googlegroups.com.

[Akarsh Sahay]

Hey Jeff,

This does answer my question and I did suspect that would be the case. Anyways that's fine because I have the option to simulate my own grids in case I want a particular evolution in detail. I had another, more important question which is the reason why I wanted the mass vs time plots. I am facing issue in simulating binaries with roche lobe overflow before the 1st supernova, as in between 2 main sequence stars or a postMS and MS star. The only RLO events I can see in my binaries happen after the 1st supernova. This is very strange considering I have simulated binaries across a vast parameter space. 

I dug deep and found some possible cues - 

1) I was initially working with the 1st Dataset (link) which contained the grids CO-HeMS, CO-HeMS_RLO, CO-HMS_RLO, HMS-HMS, single_HMS, single_HeMS. As you can see there are grids for RLO after a CO formation but not one for before it. I thought this was the issue and downloaded the newer dataset(link) which had the additional dataset HMS-HMS_RLO. But even after running simulation with this grid downloaded I could not get my pre-CO-formation RLO.

2) I then looked at the POSYDON flowchart on the github repo (link) and upon inspecting I found that there was no step corresponding to the HMS-HMS_RLO grid in the default POSYDON flowchart. But there was this function at the very end -

def initial_eccentricity_flow_chart(FLOW_CHART=None, CHANGE_FLOW_CHART=None):

"""Modify POSYDON's default flow to:

ZAMS binaries -> detached

oRLO1/oRLO2 -> HMS-HMS RLO grid

Parameters

----------

FLOW_CHART : dict or None

CHANGE_FLOW_CHART : dict or None

Returns

-------

dict

Modified flow chart.

"""

# get POSYDON's default flow chart

if FLOW_CHART is None:

if CHANGE_FLOW_CHART is None:

MY_FLOW_CHART = flow_chart()

else:

MY_FLOW_CHART = flow_chart(CHANGE_FLOW_CHART=CHANGE_FLOW_CHART)

else:

if CHANGE_FLOW_CHART is None:

MY_FLOW_CHART = flow_chart(FLOW_CHART=FLOW_CHART)

else:

MY_FLOW_CHART = flow_chart(FLOW_CHART=FLOW_CHART,

CHANGE_FLOW_CHART=CHANGE_FLOW_CHART)

# modify the default flow chart

for key in MY_FLOW_CHART.keys():

s1_state, s2_state, state, event = key

# always take ZAMS binaries to step_detached

if event == 'ZAMS' and state in BINARY_STATES_ZAMS: # check for event

MY_FLOW_CHART[key] = 'step_detached'

# Add two stars initating RLO (now coming from detached) into flow chart

# Add stripped stars (from winds) initating RLO

for s1 in STAR_STATES_H_RICH_EVOLVABLE + STAR_STATES_HE_RICH_EVOLVABLE:

for s2 in STAR_STATES_H_RICH_EVOLVABLE + STAR_STATES_HE_RICH_EVOLVABLE:

MY_FLOW_CHART[(s1, s2, 'RLO1', 'oRLO1')] = 'step_HMS_HMS_RLO'

MY_FLOW_CHART[(s1, s2, 'RLO2', 'oRLO2')] = 'step_HMS_HMS_RLO'

return MY_FLOW_CHART

Which I believe is modifying the flowchart from

ZAMS ->HMS-HMS (unless very wide binaries or initial RLO) in the default flowchart and no step to go to the HMS_HMS-RLO grid

to

ZAMS ->Detached ->(And if conditions for RLO) HMS_HMS-RLO and no step to go to the HMS-HMS grid

I have several questions regarding this choice- 

a) Why can't we have ZAMS -> HMS-HMS (unless very wide binaries in which it becomes detached or initial RLO) -> (And if conditions for RLO) HMS_HMS-RLO for the default flowchart.

This seems like the most physically intuitive default evolutionary pathway for a binary. I realise we can make custom steps, but I don't get the rationale for excluding 1 of the 2 very normal steps in the provided flowcharts. Please let me know if I am missing something very obvious.

b) Even when I used the 'initial_eccentricity_flow_chart' which has the HMS-HMS_RLO step, I wasn't able to simulate a nice conservative mass transfer between the 2 stars before the 1st CO formation. For eg here - 

This is the mass vs time ( just from the data in pop.history for a simulation). The primary undergoes what I am guessing is wind mass loss initially and then RLO (between the first 2 dotted lines) which is confirmed by the event oRLO1 there. The issue is the mass loss isn't conservative and the secondary is not gaining all the mass lost by the primary. Meanwhile in your 1st paper (link) on page 9 - "For binaries with a nondegenerate accretor (those in our grid of two H-rich stars), initially all of the mass lost by the donor through RLO is accepted by the accretor." If that's the case why can't I have conservative mass loss through the Roche Lobe in my case. To investigate this further I plotted some of your HMS-HMS_RLO grids which also do not show conservative RLO. For eg - 

Even in the grid mass loss is not conservative.

c) RLO is a very common thing to happen, I don't realise what I am missing if I am finding it so hard to simulate. Even in both your papers, on page 14 for the first paper and page 26 for the second paper, the majority of the grids undergo some form of stable RLO before the first supernova. Even when simulating the exact initial conditions from these graphs, I am not able to replicate the results.

Please let me know what I am missing in regards to simulating stable RLO before first CO formation as that's critical in forming the kind of binaries I want.

Thanks,

Akarsh

[Max Briel]

Hi Akarsh,

Thank you for the follow-up questions and the investigation you've done.

I think they boil down to:

1. Why do I not see RLO1 before supernova 1 (CC1) in the history?

2. Why is the mass transfer non-conservative?

I answer them in detail below and hope it clarifies things. I've also added some clarification on the use of RLO grids.

If it's okay with you, I will put this email chain on our GitHub Discussions page to keep track of it more easily.

Let me know if that's okay and if you still have any more questions.

Cheers,

Max Briel

1. Why do I not see RLO1 before supernova 1 (CC1) in the history?

Roche lobe overflow is happening, but it is not shown within the "history" because it is occurring within the step-HMS-HMS (within the MESA binary models). If no unstable mass transfer took place, the complete evolution of the binary until carbon depletion occurs within step-HMS-HMS, thus you reach CC1 immediately after that step. On the other hand, models that experienced unstable mass transfer will exit step-HMS-HMS with "RLO1" as a state and "oCE1/2" as an event. These go into step_CE.

The history does not show you RLO1 for stable mass transfer phases, because we do not exit step-HMS-HMS and use the MESA binary evolution until carbon depletion. I want to note that multiple mass transfer phases can take place in the binary, which are considered if all are stable.

I hope this explains why the RLO1 does not show up in the history. There are the following ways to get more detailed information on the evolutionary history:

A. For your population, you can use

pop_file_path = "your/file"

pop = Population(pop_file_path)

pop.calculate_formation_channels(mt_history=True)

This will add the "pop.formation_channels" to your population, including information about the step-HMS-HMS evolution. See this tutorial section for more details: Population.formation_channels

These should give you the channels, such as:  "ZAMS_oRLO1_CC1_oRLO2_oCE2_CC2_END"

This includes the mass transfer from step-HMS-HMS and should be what you need.

B. Manual inspection

The formation_channels use the information in the "oneline" to build the channel; specifically, it uses the columns:

"interp_class_HMS_HMS" and "mt_history_HMS_HMS". 

"interp_class" provides for each grid what the mass transfer type was and is used to determine what interpolator to use.

"mt_history", on the other hand, tells you what the last interaction was or if contact occurred. 

Example:

pop.oneline['interp_class_HMS_HMS']

pop.oneline['mt_history_HMS_HMS']

C. Grid history

The grid contains a bit more information on the HMS-HMS mass transfer history, specifically "termination_flag_2". It contains a cumulative string of all the mass transfer phases that have occurred in the specific binary model or that it has undergone a contact phase.

This can be simple: "case_A1", or complicated with multiple phases: "case_B1/C1/BB1".

Currently, only when using the nearest neighbour population synthesis can these be made available in a population. Otherwise, you will have to inspect the grid manually.

Example:

HMS_HMS_file = os.path.join(PATH_TO_POSYDON_DATA, "HMS-HMS", "1e+00_Zsun.h5")

grid = PSyGrid(HMS_HMS_file)

grid.final_values['termination_flag_2']

These 3 methods should give enough insight into what specific binaries are doing.

2. Why is the mass transfer non-conservative?

The mass transfer efficiency is defined by the response of the accretor star, as modelled within the MESA model. You mention the following line from the POSYDON v1 paper: "For binaries with a nondegenerate accretor (those in our grid of two H-rich stars), initially all of the mass lost by the donor through RLO is accepted by the accretor.".

This is correct. However, the accreted material will rapidly spin up the accretor star to near critical rotation. A bit further down in the same paragraph, the following is mentioned: "The accreted angular momentum spins up the accretor, and mass accretion is restricted when the accretor reaches critical rotation.". When it reaches this stage, the mass will be ejected with the angular momentum of the accretor. Since only a small amount of mass is required to spin up the surface of the accretor to near critical, most mass transfers are generally non-conservative. It should be noted that the efficiency depends on the initial period, mass, and mass ratio. Tighter binaries allow for more efficient accretion due to tides spinning down the accretor.

At the moment, changing the mass transfer efficiency with POSYDON can only be achieved by rerunning grids of MESA binary models.

Usage of RLO grids

I think it is worth clarifying the difference between RLO grids and non-RLO grids, and the addition of the HMS-HMS_RLO grid for others reading this later as well. The POSYDON DR1 and DR2 binary grids are evolved with the assumption of circular orbits. After receiving a natal kick during the first supernova, this assumption is no longer true. Thus, after SN1 the binaries are evolved in step_detached, where the eccentric orbit is considered until Roche lobe overflow occurs. Then the binary is circularised and goes into CO-HMS_RLO, where the grid starts at Roche lobe overflow. This is why you are seeing RLO2 occur after the first SN: The exit condition of step_detached is RLO2.

Because in the default population synthesis, no eccentricity at ZAMS is added. The binaries are already circular. So instead of first evolving the binaries detached and then through RLO, we immediately go into the HMS-HMS grid. This has the benefit that the binaries are self-consistently modelled from ZAMS until they leave the HMS-HMS grid (either at carbon-depletion or due to unstable MT), and includes all the mass transfer occurring in the HMS-HMS grid.

If eccentricity at ZAMS is needed, the alternative flow can be used with sampling eccentricity at ZAMS. This first evolves the eccentric binaries in step-detached and then goes into HMS-HMS_RLO when a Roche lobe is filled, foregoing the fully self-consistent modelling from ZAMS for including eccentricity at ZAMS. For example, see Kruckow+25, which includes it for non-interacting binaries in the context of Gaia black hole systems. The HMS-HMS_RLO grid is a processed HMS-HMS grid, where non-interacting binaries are removed, and the evolution of interacting binaries only starts at Roche lobe overflow. Thus, the HMS-HMS grid already contains the RLO phases, as described under question 1.

To view this discussion visit https://groups.google.com/d/msgid/posydon-users/6655521d-2d17-4b90-99a1-81fe9d839c78n%40googlegroups.com.