Binary implicit mass transfer scheme convergence fails if rotation and bondi-hoyle are enabled

[ruggero-valli]

Description of the bug When evoling a binary system that undergoes a stable mass transfer, under a certain set of conditions (see below), the implicit mass transfer scheme may not be able to converge due to a problem in the interaction between Bondi-Hoyle accretion, rotation, and the implicit mass transfer convergence algorithm.

To Reproduce To reproduce the bug, evolve a binary system with

1. an implicit mass transfer scheme (I used Kolb)
2. Rotation enabled.
3. Bondi-Hoyle accretion enabled.
4. To make the issue more evident, set (in binary_controls) max_tries_to_achieve to a high number (e.g. 100) and implicit_scheme_tiny_factor=1d-99 so that the code will not assume convergence just because the bisection interval is small. It can also be useful to set report_rlo_solver_progress=.true. as a diagnostic.

When mass transfer starts, the accretor will be quickly spun up to critical rotation. As soon as the accretor reaches critical rotation (and enters in implicit_wind_limit), the implicit mass transfer scheme will be unable to converge and will redo the step until max_tries_to_achieve is reached.

Expected behavior The implicit mass transfer scheme should be able to converge.

A possible explanation Short version: when the accretor reaches critical rotation, it will shed most of the mass that is being thrown onto it. Some of this mass will be flung back to the donor through the Bondi-Hoyle accretion mechanism. This will change the mass (and radius) of the donor while the implicit scheme is trying to converge to a solution. The result is that the root of the implicit equation that the implicit mass transfer scheme is trying to solve will be shifted outside of the interval where it is supposed to be found. As a consequence, the implicit mass transfer scheme will never find it.

Long version: In binary_mdot.f90:

function_to_solve = (explicit_mdot - new_mdot) / new_mdot

The implicit mass transfer scheme (e.g. Kolb) tries to find the root of the function $f(\dot M{tr}) = (\dot M{tr}^{explicit}- \dot M{tr})/\dot M{tr}$, where $\dot M_{tr}^{explicit}$ is calculated at the end of the step and thus depends on

1. the state of the system at the beginning of the step (e.g. mass and structure of the stars, roche lobe radius, etc..)
2. the mass transfer rate $\dot M_{tr}$ that is being guessed at the current iteration of the step loop.

Crucially, the state of the system at the beginning of the step must be the same at each iteration of the step loop, in other words, the calculated value of $f(\dot M{tr})$ should only depend on $\dot M{tr}$. This is usually the case, in fact, at the end of a step loop iteration, if result == redo, the system is reset to the initial state. But if the accretor is rotating at critical velocity, the following condition is satisfied (found in run_binary_support.f90)

! Avoid repeting the accretor when using the implicit scheme plus
! rotation and implicit winds. When this happens the accretor won't
! usually care about the result of the evolution of the donor.
if (j == b% a_i .and. b% num_tries >0 .and. s% was_in_implicit_wind_limit) &
    cycle

and the accretor is not evolved nor reset. As stated in the code comments, it is true that the accretor won't care about the result of the evolution of the donor, but -if Bondi-Hoyle accretion is turned on- the donor might instead care about the state of the accretor.

At the beginning of the next iteration, the $\dot M$ of the donor due to Bondi-Hoyle is calculated as (in binary_mdot.f90)

b% mdot_wind_transfer(b% a_i) = b% s_accretor% mstar_dot * &
             b% wind_xfer_fraction(b% a_i)
...
b% s_donor% mstar_dot = b% s_donor% mstar_dot + b% mtransfer_rate - &
            b% mdot_wind_transfer(b% a_i)

At this point of the flow, the variable b% s_accretor% mstar_dot is supposed to only contain the contribution of the stellar winds, calculated shortly before the above lines with the subroutine do_set_mdot(s, L_phot, M_phot, T_phot, ierr). But if the accretor is at critical rotation (s% was_in_implicit_wind_limit == .true.) do_set_mdot is never called and b% s_accretor% mstar_dot contains the whole mass lost at the end of the previous iteration of the step loop, including all the mass that the accretor has been shedding due to being at critical rotation. This new mass that is added to the donor through Bondi-Hoyle can sometimes be sufficient to change $\dot M{tr}^{explicit}$ enough to to move the root of $f(\dot M{tr})$ outside of the interval where the implicit mass transfer scheme is seeking it, preventing the convergence of the algorithm.

Suggested solution The simplest solution would be to insure that the accretor is always reset before a redo, even if at critical rotation (s% was_in_implicit_wind_limit == .true.). But evolving the accretor every time can introduce a significant increase in computational time. Furthermore, I am not sure whether this suggested solution could lead to something else breaking.

Further information I am attaching the inlists that I have used to reproduce the behavior. I am running mesa version 22.05.1 inlist.zip

[matthiasfabry]

Hello Ruggero, thank you for this extensive report. As a first (and obvious) curiosity, do you recover your intended behavior if you simply comment out the avoid repeating the accretor conditional?

[ruggero-valli]

Hi Matthias, Yes, if I comment out that conditional the implicit mass transfer scheme is able to converge. I don't know if other parts of the coded depended on the assumption that the accretor is not being evolved, though.

Moreover, evolving the accretor every redo takes a lot of time, because it needs to do an additional loop to determine the amount of mass to lose that brings it below critical rotation (the function select_mdot_action in the implicit_mdot_loop in evolve.f90), and the result of this calculation is always more or less the same. So, if we choose to evolve the accretor we commit to repeat every time a cumbersome calculation that gives always the same result.

Maybe a better way to solve the problem could be

1. the first time the step loop is run, save somewhere the value of b% s_accretor% mstar_dot that is calculated with the subroutine do_set_mdot(s, L_phot, M_phot, T_phot, ierr).
2. in the following iterations of the step loop, if the accretor is in implicit_wind_limit, use that saved value to calculate Bondi-Hoyle instead of the current b% s_accretor% mstar_dot.

An issue that I see with this approach is that we are still not resetting the accretor to the initial state at the beginning of the step loop. Thus the mass of the accretor might be slightly different and this might reflect on a different roche lobe geometry and in turn affect the mass transfer rate. Although in the models that I have run this was never enough to hinder the convergence of the mass transfer scheme, maybe there are cases in which it does.