MESAHub / MESAHub/mesa

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

Open
#507 2 comments 0 reactions 1 assignee Claimed by @matthiasfabry View on GitHub
binary bug
Dominant language
Fortran
Stars
249
Forks
82
Avg merge
1d 3h
Merged PRs (30d)
3

Description

**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](https://github.com/MESAHub/mesa/files/10942326/inlist.zip)

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.