# Unified Rapid Mass Transfer

Natalia Ivanova Department of Physics, University of Alberta, Edmonton, T6G
2E7, Alberta, Canada  Surjakanta Kundu  Department of Physics, University of
Alberta, Edmonton, T6G 2E7, Alberta, Canada  Department of Physical Sciences,
Indian Institute of Science Education and Research Kolkata, Mohanpur 741246,
India  Ali Pourmand Department of Physics, University of Alberta, Edmonton,
T6G 2E7, Alberta, Canada Natalia Ivanova<EMAIL_ADDRESS>
(Accepted June 5, 2024)

###### Abstract

We present a method to obtain rapid mass loss rates in binary systems,
specifically at the onset of MT episodes. The method unifies atmospheric
(underflow) and $L_{1}$ stream (overflow) mass rates in a single continuous
procedure. The method uses averaged 3D properties of the binaries, such as
effective binary potential and effective binary acceleration, to both evolve
the donor and obtain properties of the matter at the $L_{1}$ plane. In the
case of underflow, we obtain atmospheric stratification. Our method can be
used for binaries with an extensive range of mass ratios, $0.01\leq q\leq
100$, and can also be applied to hot donors. The considered examples show that
the MT rates obtained with this revised formalism always differ from the
optically thin and optically thick MT rates widely used during the
computations of binary evolution.

Multiple star evolution — Binary stars — Roche lobe overflow

††facilities: Digital Research Alliance of Canada

##  1 Introduction

Mass transfer (MT) is a type of interaction in the lives of close binaries
when one star’s outer layers become too close to the effective binary
potential that passes through the first Lagrangian point, $L_{1}$, between the
two binary companions. The surface of this equipotential is known as a Roche
lobe. Proximity of the star to its Roche lobe allows the donor’s outer layers
to flow into the Roche lobe of its companion.

Stable MT can be a long-lived stage during which binaries appear for example
as X-ray binaries, including low-mass X-ray binaries, high-mass X-ray
binaries, cataclysmic variables, and more (see review on mass-transferring
systems in Tauris & van den Heuvel, 2023). A fast phase of stable MT can be
responsible for the formation of such exciting objects as Algols, stripped
supernovae, and more.

However, if the MT becomes dynamically unstable, the progenitor binary system
may enter a Common Envelope (CE) phase. During a CE event, one of the binary
companions (accretor) orbits inside the envelope of another binary companion
(donor). This transitional evolutionary phase was proposed to explain
cataclysmic variables (Webbink, 1975) and V471 Tau (Paczynski, 1976) (for a
more detailed history of the idea’s origin and its development throughout last
50 years, see the review in Ivanova et al., 2020). Currently, a CE phase is
considered to be a major transformational event responsible for forming the
progenitors of the already mentioned X-ray binaries, and also progenitors of
merging double neutron stars or double black holes.

A dynamical phase of a CE event starts when the donor overfills the
equipotential that passes through the outer Lagrangian point. Understanding
the initiation of a CE phase is crucial to answer CE outcomes. Specifically,
how much mass was lost during the pre-CE mass loss, i.e., before the donor
overfilled its outer Roche lobe?

There are several types of techniques to evolve the donor through the mass
loss episode. One, often called ”implicit”, is to assume that the donor has to
remain inside its Roche lobe. The MT rate is then found by taking away as much
mass as needed to keep the donor exactly inside its Roche lobe. This method is
excellent for finding smooth and continuous MT rates during stable MT.
However, this approach is intrinsically inapplicable for modeling the start of
a dynamically unstable MT leading to a CE event.

An opposite approach is when the donor is considered to not coincide precisely
with its Roche lobe, and the flow is formed from the layers that exceed the
Roche lobe. The converging-diverging behavior of the effective binary
potential in the $L_{1}$ neighborhood suggests that one may consider the
$L_{1}$ neighborhood as an analogy of the de Laval nozzle, where the $L_{1}$
plane acts as a transition cross section between subsonic and transonic flows
(Nariai, 1967). Adding the assumption that streamlines of mass flow do not
cross each other and follow the Bernoulli theorem, Nariai (1967) introduced
the way to find the rate of mass loss in the case of isentropic streams. While
the applicability of the de Laval formalism for converging-diverging
equipotentials was questioned (for example, see Lubow & Shu, 1975), it is
still a standard and widely accepted approximation to find the mass loss rate.

The next comprehensive analysis of gas dynamics during the MT, though with the
goal to understand the post $L_{1}$ ballistic stream trajectory, was developed
in Lubow & Shu (1975). The additional approximation discussed there and later
used elsewhere was to consider the MT using an isothermal stream, and the use
of a Taylor expansion for the converging-diverging behavior of the binary
potential near $L_{1}$.

Later, Eggleton (1983) made three dimensional (3D) integrations of the Roche
lobe and derived the now famous empirical formula for a volume equivalent
Roche lobe radius. Mochnacki (1984), also using 3D integrations but outside of
the Roche lobe, provided additional quantities specifically related to binary
potentials.

This set of approximations for a Roche lobe geometry led to the development of
what is considered now to be the two most commonly used prescriptions for MT
rates, both of which are functions of how much the donor overfills or
underfills its Roche lobe.

The Roche lobe “underfilling” MT occurs since a star has no rigid surface at
its photosphere, and the donor’s material extends beyond what is usually
called the star’s radius. Then, one can find how much mass can be lost from
the donor’s atmosphere before the donor’s radius formally becomes equal to its
Roche lobe volume-equivalent radius. The detailed consideration of this
prescription using the isothermal stream consideration for cold (low-effective
temperature) donors was done by Ritter (1988). Importantly, it used an
empirical equation for the potential that is based on Mochnacki (1984) 3D
integrations and is valid for mass ratios (the ratio of the donor mass to the
mass of its companion) $0.1\lesssim q\lesssim 3$. The surface temperature of
the donor should be cold enough so that there is no significant radiation
pressure. We will refer to this formalism as Original Ritter (OR) formalism.
The mass loss considered this way is the optically thin MT or atmospheric
overflow.

Roche lobe overflow (RLOF) can also be considered in a limit of an optically
thick stream formalism (Kolb & Ritter, 1990). We will refer to this formalism
as KR formalism. It uses a combination of an adiabatic stream, a second-order
Taylor expansion near $L_{1}$ for potential, and again a matter-dominated
approximation for the total pressure.

One has to choose which approach to use in their calculation, optically thin
or optically thick, based on other physics of the considered binary system. OR
can be used for both underflow and overflow, but it is problematic to justify
that a thick $L_{1}$ stream is optically thin. On the other hand, KR by design
can only find MT rates in case of an overflow, and hence using KR usually
implies that OR is used while the system is still in underflow. We note,
however, that the physical assumptions that are used to derive optically thin
and optically thick mass loss rates differ, and by design, there is no
continuity between the two regimes.

The donors that enter the CE event do not have to be as cold as those in
cataclysmic variables, for which OR and KR were derived. Further, a dynamical
instability of the MT is often the consequence when the mass ratio of the
donor to its companion is large, $q>3$ 111There is no unique value of the
critical mass ratio $q_{\rm crit}$ as it is highly dependent on the
evolutionary stage of the donor, the type of its companion, and the assumed
mode of the MT. As guidance, $q>q_{\rm crit}\approx 3.5$ would most likely
lead to dynamically unstable MT for radiative donors, $q_{\rm crit}\approx 2$
for low and intermediate-mass convective donors, and $q_{\rm crit}\approx 8$
for very massive donors. For more details, see the review in Ivanova et al.
(2020).. The dynamically unstable MT also proceeds with a significant degree
of overflow. That is beyond the applicability of the empirical approximation
integrated in OR formalism and the second-order Taylor expansion used in KR.
An additional degree of inconsistency lies in using the parametrization of the
binary effective potential for the stream condition with a star that is not
evolved in the same potential.

Recently, Pourmand & Ivanova (2023) (Paper I) have described a method to
obtain 3D averaged stellar properties for a star inside a binary potential,
from the center of the donor to the equipotential passing through the donor’s
outer Lagrangian point, as well as the properties at the $L_{1}$ plane, for an
extensive range of mass ratios, $10^{-5}\leq q\leq 10^{5}$. The obtained
tables now enable integrating a star’s evolution in a binary potential. The
next step is to obtain the mass loss rates which use the average 3D
properties.

In this paper, we first re-derive the atmospheric overflow and RLOF
formalisms, see § 2. The derived unified formalism represents the continuous
approach. In § 3, we compare the revised unified formalism with OR and KR
formalisms. In § 4, we evaluate the trends of the impact of the new formalism
on the mass transferring system.

##  2 Method description

###  2.1 Roche lobe overflow

For MT calculations when the donor overfills its Roche lobe, we use, as the
foundation, the optically thick stream formalism developed by our group,
Pavlovskii & Ivanova (2015), updated here to take into account the averaging
of the binary properties from 3D to 1D as published in Paper I.

As in Pavlovskii & Ivanova (2015), we obtain the MT rate by integrating the
mass flow over a “nozzle” cross-section that is taken on the plane
perpendicular to the line connecting the centers of the two stars and passing
through the L1 point:

$\dot{M}_{\rm RLOF}=\oiint_{L_{1}nozzle}\rho_{L1}(x,y)c_{s,L1}(x,y)dA\ .$ (1)

Here $(x,y)$ are the coordinates on the $L_{1}$ plane, $\rho_{L1}(x,y)$ is the
density of the donor’s material as it reaches the $L_{1}$ plane,
$c_{s,L1}(x,y)$ is the sonic velocity of the gas (we adopt the conventional
nozzle approximation, in which the speed with which gas passes the nozzle is
equal to its local sonic velocity), and $dA$ is the area element of the nozzle
cross-section at the $L_{1}$-plane.

Each location $(x,y)$ on the $L_{1}$ plane has a value of the effective binary
potential $\Phi_{l}(x,y)$. A group of points $(x,y)$ on the $L1$ plane that
has the same value of $\Phi_{l}(x,y)$ forms a closed loop. The closed loop
(circuit) containing all the points that have a value of the effective
potential $\Phi_{l}$ is ${\mathcal{L}}({\Phi_{l}})$. All thermodynamic
quantities, for example, density or sonic velocity, are the same along each
circuit ${\mathcal{L}}({\Phi_{l}})$. For each $\Phi_{l}$, we approximate the
circuit $\mathcal{L}({\Phi_{l}})$ with an ellipse $\mathcal{E}({\Phi_{l}})$ on
$L_{1}$ plane. Each ellipse $\mathcal{E}({\Phi_{l}})$ can be described by its
semi-major axis, $r_{l}$, and its semi-minor axis. The semi-major and semi-
minor axes of the ellipse, as a function of $\Phi$, are tabulated values from
Paper I, and $\mathcal{E}({\Phi_{l}})$ is the ellipse circumference. The
nozzle for this stream is limited by the effective potential of the donor’s
photosphere $\Phi_{\rm ph}$, which corresponds to the ellipse
$\mathcal{E}_{\rm ph}\equiv\mathcal{E}({\Phi_{\rm ph}})$ with the associated
semi-major axis $r_{\rm ph}$. We will simplify notations, stating that for a
quantity $\cal{Q}$ at the equipotential loop determined by $\Phi(r_{l})$, the
equivalency in notations is implied as follows:

${\cal{Q}}(\Phi(r_{l}))\equiv{\cal{Q}}(r_{l})\equiv{\cal{Q}}(\Phi_{l})\ .$

Equation 1 can be rewritten as

$\dot{M}_{\rm RLOF}=\int_{0}^{r_{\rm ph}}\rho_{L1}(r_{l})\ c_{s,L1}(r_{l})\
\mathcal{E}(r_{l})dr_{l}\ .$ (2)

We can find at which volume equivalent radius $r_{i}$ the donor would
experience the same effective binary potential: for each circuit on $L_{1}$
plane described by $r_{l}$, we find the volume-equivalent radius of the donor
$r_{i}$ such that $\Phi(r_{i})=\Phi(r_{l})$. To find $r_{i}$, we use tabulated
values from Paper I. The main improvement, as compared to Pavlovskii & Ivanova
(2015), is that now we can precisely relate the properties of a gas at the
$L_{1}$ plane with the same properties in the equipotential shells inside the
donor.

Following Pavlovskii & Ivanova (2015), and assuming that each streamline is
adiabatic (while the flow through the $L_{1}$-plane itself is not isentropic),
we obtain $\rho_{L1}(\Phi(r_{l}))\equiv\rho_{L1}(r_{l})$ and $c_{s,L1}(r_{l})$
from solving the Bernoulli equation for streamlines along the equipotentials:

$\frac{P(r_{i})}{\rho(r_{i})}+u(r_{i})=\frac{c_{s,L1}^{2}(r_{l})}{2}+\frac{P_{L1}(r_{l})}{\rho_{L1}(r_{l})}+u_{L1}(r_{l})\
.$ (3)

Pavlovskii & Ivanova (2015) have demonstrated that using the polytropic
equation of state, instead of the full tabulated equation of state as provided
by MESA, results in a difference for the final stream mass loss within 4
percent. In this work, we chose to use the polytropic equation of state to
speed up the calculations. In that case, the density and sonic velocity on the
$L_{1}$ plane can be found as

$\displaystyle\rho_{L1}(r_{l})$ $\displaystyle=$
$\displaystyle\left(\frac{2}{\gamma_{i}+1}\right)^{\frac{1}{\gamma_{i}-1}}\rho(r_{i})$
$\displaystyle c_{s,L1}^{2}(r_{l})$ $\displaystyle=$
$\displaystyle\gamma_{i}\left(\rho_{L1}(r_{l})\right)^{\gamma_{i}-1}\frac{P(r_{i})}{\rho_{L1}(r_{i})^{\gamma_{i}}}$
(4)

Here $\gamma_{i}$ is $\gamma\equiv\Gamma_{1}=\left(\partial\ln
P/\partial\ln\rho\right)_{\rm ad}$ in the donor at $r_{i}$.

###  2.2 Atmospheric overflow

For mass loss from the atmosphere, we use Ritter’s formalism (Ritter, 1988),
but strongly modified. To be specific, the original Ritter’s formalism assumes
that both the atmosphere and the stream are isothermal, and then applies a
number of approximations to find the mass loss rate. In what follows, we will
discuss the isothermal and adiabatic assumptions and detail which
approximations can be replaced with exact integrations.

We will start with the overall setup. The original Ritter’s formalism adopts
that both the atmosphere of the donor (to be specific, all the matter above
the photosphere) and the stream, while it travels to the $L_{1}$ plane, are
kept isothermal, and the material reaches the sonic speed at the nozzle
($L_{1}$ plane). Since the stream is isothermal due to external radiation,
only the form of the Bernoulli equation derived from the momentum conservation
can be used:

$\frac{1}{2}v^{2}+\int\frac{dP}{\rho}+\Phi={\rm const}$

To find the solution, OR evaluates the Bernoulli constant at the photosphere
and at the $L_{1}$ plane. This implies that the streamlines are meant to start
at the photosphere and then are expected to arrive at the $L_{1}$ plane. The
external energy source is mandatory here, as without an external energy
source, a flow of matter from a deep potential well at the photosphere to the
$L_{1}$ plane while accelerating to sonic velocity and keeping the same
temperature is against energy conservation. Using this assumption, the density
at $L_{1}$ in OR is

$\rho_{\rm L1}^{\rm is}(\Phi_{L1})=\frac{1}{\sqrt{e}}\rho_{\rm
ph}\exp\left(-\frac{\Phi_{L1}-\Phi_{\rm ph}}{c_{s,\rm ph}^{2}}\right)$ (5)

Here $\Phi_{L1}$ is the effective binary potential at $L_{1}$ and $\Phi_{\rm
ph}$ is the effective binary potential at the photosphere, $\rho_{\rm ph}$ is
the density of the photosphere, and $c_{s,\rm ph}$ is the sonic velocity in
the photosphere, atmosphere and $L_{1}$ plane. The equation above adopts the
use of a standard set of assumptions describing an isothermal astrophysical
flow: $P\propto\rho$, and isothermal sonic velocity $c_{\rm s}^{\rm
is}=\sqrt{P/\rho}={\rm const}$.

We note that while the original setup makes the stream originate from the
photosphere and end at the $L_{1}$ plane, the setup is effectively equivalent
to two separate assumptions: a) the atmosphere is in hydrostatic equilibrium
and is isothermal due to long-term irradiation by the donor; b) matter leaves
the atmosphere along the equipotentials and remains isothermal as it travels
to $L_{1}$.

Having the two separate assumptions adds more insight. First, the assumption
for the donor’s atmosphere to be in isothermal equilibrium is a frequent
assumption for a stellar atmosphere (although see §3.1). Second, matter
streaming along the equipotentials is more self-consistent and aligns with how
RLOF mass loss is seen222The ratio of Coriolis acceleration to the effective
binary acceleration on $L_{1}$ plane is comparable to the ratio of sonic
velocity at $L_{1}$ plane to the product of binary angular velocity and the
distance to $L_{1}$ plane from the center of the donor. The described common
assumption on the stream being aligned with equipotentials implies that
Coriolis acceleration is chosen to be ignored, for being a second-order
correction, see Lubow & Shu (1975).. However, unlike the assumption for the
hydrostatic atmosphere to be isothermal due to long-term donor irradiation,
the stream, as it starts to flow, does not have to remain isothermal,
especially if the considered MT episode is dynamically unstable. In what
follows, we will consider separately if the stream material remains isothermal
as it flows along an equipotential, or adiabatic.

Let us start from the assumption that the atmosphere is isothermal and is in
hydrostatic equilibrium. Then, pressure and density in the atmosphere at the
equipotential $\Phi_{i}$ that corresponds to the volume equivalent radius
$r_{i}$ are

$\displaystyle P_{\rm at}(\Phi_{i})$ $\displaystyle=$ $\displaystyle P_{\rm
ph}\exp\left(-\frac{\Phi_{i}-\Phi_{\rm ph}}{c_{s,\rm ph}^{2}}\right)$
$\displaystyle\rho_{\rm at}(\Phi_{i})$ $\displaystyle=$
$\displaystyle\rho_{\rm ph}\exp\left(-\frac{\Phi_{i}-\Phi_{\rm ph}}{c_{s,\rm
ph}^{2}}\right)$ (6)

A further assumption, for the stream to remain isothermal while it flows along
the equipotential, leads to the same expression for the density at $L_{1}$ as
in Equation 5. However, as we speak now about different equipotentials in the
atmosphere, we can also find the density at any location in the $L_{1}$ plane:

$\displaystyle\rho_{L1}^{\rm is}(\Phi_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{e}}\rho_{\rm
ph}\exp\left(-\frac{\Phi_{i}-\Phi_{\rm ph}}{c_{s,\rm ph}^{2}}\right)$
$\displaystyle c^{\rm is}_{s,L1}(\Phi_{i})$ $\displaystyle=$ $\displaystyle
c_{s,\rm ph}=\sqrt{\frac{P_{\rm ph}}{\rho_{\rm ph}}}$ (7)

In this scenario, the sonic speed remains constant along the stream and has,
at the $L_{1}$ plane, the same value as at the photosphere. The combination of
isothermal atmosphere and isothermal stream recovers Ritter’s values at
$L_{1}$.

If, instead, the matter is assumed to evolve adiabatically as it travels along
the equipotential, while it is still departing from the isothermal atmosphere,

$\displaystyle\rho_{L1}^{\rm ad}(\Phi_{i})$ $\displaystyle=$
$\displaystyle\left(\frac{2}{\gamma_{i}+1}\right)^{\frac{1}{\gamma_{i}-1}}\rho_{\rm
ph}\exp\left(-\frac{\Phi_{i}-\Phi_{\rm ph}}{c_{s,\rm ph}^{2}}\right)$
$\displaystyle c_{s,L1}^{\rm ad}(\Phi_{i})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2}{\gamma_{i}+1}}c_{s,\rm ph}$ (8)

Here $\gamma_{i}$ is $\gamma\equiv\Gamma_{1}$ at the corresponding
equipotential level of the atmosphere, and the sonic velocity of the
atmosphere here is the same quantity as was used to make the isothermal
atmosphere, $c_{s,\rm ph}=\sqrt{P_{\rm ph}/\rho_{\rm ph}}$.

Note that integration over the whole $L_{1}$ plane is implied, with $r_{\rm
out}$ reaching the outer Lagrangian equipotential in the $L_{1}$ plane. The
mass loss rate can be obtained for the isothermal or adiabatic scenarios using
the corresponding density and sonic velocity expressions. As at each location
in the $L_{1}$ plane, the product of $\rho_{L1}^{\rm ad}(\Phi_{i})c^{\rm
is}_{s,L1}(\Phi_{i})$ and the product of $\rho_{L1}^{\rm
ad}(\Phi_{i})c_{s,L1}^{\rm ad}(\Phi_{i})$ differ only by a constant factor of
about one, the integrated MT rates from an isothermal atmosphere with either
an isothermal or adiabatic stream will be about the same. For example, for
$\gamma=5/3$, an isothermal stream from an isothermal atmosphere will result
in MT rates that are about 1.078 times those assuming an adiabatic stream from
an isothermal atmosphere, and for $\gamma=4/3$ the ratio is 1.04.

In addition to isothermal and adiabatic atmospheres, one can construct a
simple atmosphere model and then apply the isothermal or adiabatic Bernoulli
equation for each atmospheric layer to obtain values of the density
$\rho_{L1}(r_{l})$ and sonic velocity $c_{s,L1}(r_{l})$ as functions of the
position (or the effective potential) at $L_{1}$ plane. There is no analytical
solution in this case. Further, the construction of the atmosphere can proceed
with several approximations, which we will list here:

  1. 1.

EG. This is a classic Eddington grey (EG) atmosphere in which both the opacity
of the photosphere $\kappa_{\rm ph}$ and the gravitational acceleration
$g_{\rm ph}=(d\Phi/dr)_{\rm ph})$ are kept constant through the atmosphere
(see, e.g., §4.3 in Hansen et al., 2004).

  2. 2.

EG$s$. This atmospheric model considers that the atmosphere can be extended,
and hence, its geometry is not fully plane-parallel; the radiative energy
transfer equation must be solved in spherical geometry in this case. Paczyński
(1969) noted that extended atmospheres have an intrinsically large sensitivity
to the adopted (initially arbitrary) value of the temperature at the true
surface (where density becomes zero). To resolve this numerical nuisance, a
further approximation must be made. For example, Paczyński (1969) proposed
adding an approximate slope in the equation for the temperature gradient (see
their Equation 7).

  3. 3.

EG$\Phi$. In this atmospheric model, the effective gravitational acceleration
is not constant.

  4. 4.

EG$\Phi\kappa$. In this variation, both the opacities and the effective
gravitational acceleration are functions of the radius.

Please note that only the EG atmosphere would be called an Eddington grey
atmosphere by default in the literature.

The total mass loss rate, in case of only atmospheric outflow, is

$\dot{M}_{\rm at}=\int_{0}^{r_{\rm out}}\rho_{at}(r_{l})\ c_{s,\rm at}(r_{l})\
\mathcal{E}(r_{l})dr_{l}$ (9)

###  2.3 The total mass loss rate.

Now, we can obtain the total mass loss rate when the donor overfills the Roche
lobe while keeping an atmospheric outflow. The combined mass loss rate is
important for transitioning from the Roche lobe underfilling (atmospheric)
outflow to the Roche lobe overfilling outflow.

$\displaystyle\dot{M}_{\rm tot}$ $\displaystyle=$
$\displaystyle\int_{0}^{r_{\rm ph}}\rho_{L1}(r_{l})\ c_{s,L1}(r_{l})\
\mathcal{E}(r_{l})dr_{l}$ (10) $\displaystyle+$ $\displaystyle\int_{r_{\rm
ph}}^{r_{\rm out}}\rho_{at}(r_{l})\ c_{s,\rm
at}(r_{l})\mathcal{E}(r_{l})dr_{l}$

###  2.4 Fast mass removal.

Currently, when only the total mass loss rate is provided by approximate
prescriptions, the mass loss in stellar codes can be applied only to a donor
as a whole. Essentially, during the numerical iterations for a stellar model
at the next timestep, the atmosphere is removed and reconstructed right above
the mass shell that previously was deep inside. That motivated the idea that
the fast mass loss has to be adiabatic, as the deep layers would have no time
to lose energy during fast mass loss. Accordingly, studies of the stability of
the MT are based on the adiabatic response of the envelope. The deviation from
adiabaticity so far was linked only to the reconstruction of the
superadiabatic layer in the very outer layers of the envelope (Pavlovskii &
Ivanova, 2015).

With our mass loss formalism, we can obtain the mass loss rate for each
individual mass shell inside the donor. This only plays a role during RLOF
mass loss, where the matter is streaming directly from the equipotential
inside the donor below the photosphere.

Consider a mass shell $m_{i}$ with a mass $dm_{i}$, thickness $dr_{i}$,
located at the volume-equivalent radius inside the donor $r_{i}$. With
$\Phi(r_{l})=\Phi(r_{i})$, $r_{i}$ corresponds to the semi-major axis $r_{l}$
at the $L_{1}$ plane, while the thickness $dr_{i}$ projects to a thickness
$dr_{l}$ at $L_{1}$ plane. The mass loss from the $i$ mass shell is

$\dot{m}(r_{i})=\int_{r_{l}}^{r_{l}+dr_{l}}\rho_{L1}(r_{l})\ c_{s,L1}(r_{l})\
\mathcal{E}(r_{l})dr_{l}$ (11)

For mass shells with a small thickness $dr_{i}$ and correspondingly small
$dr_{l}$,

$\dot{m}(r_{i})=-\rho_{L1}(r_{l})\ c_{s,L1}(r_{l})\ \mathcal{E}(r_{l})dr_{l}\
.$ (12)

However, from the point of computation time, removing mass directly from the
mass shells can be useful only when the mass loss rate approaches
dynamical—only then can a fraction of a mass be removed from each mass shell
during a timestep. However, another situation can occur before the mass loss
rate becomes dynamic.

Let us introduce “sonic” mass loss from the photosphere of a 1D star as the
following:

$\dot{M}_{\rm sonic}=4\pi r_{\rm ph}^{2}c_{s,\rm ph}\rho_{\rm ph}\ .$ (13)

During a significant Roche lobe overflow, the sonic speed and density near
$L_{1}$ can greatly exceed those near the photosphere. As a result, the mass
loss rate found by integrating over the $L_{1}$ neighborhood can exceed
$\dot{M}_{\rm sonic}$. Numerically, if the mass loss is applied only to the
outer layer, the code must recover the stellar model from the supersonic
removal of it. Hence, it leads numerically to adiabatic mass removal of the
surface superadiabic layer. However, in a real star, the outer layers do not
leave the star at supersonic speed. Instead, the mass loss takes place through
the range of mass shells, and it is everywhere moving at subsonic speed on
average. We know which mass shells are attached to the $L_{1}$ plane at the
moment $t$, and we can predict which mass shells will be in contact with the
L1 plane during the next timestep. We then remove the mass that has to be lost
from those mass shells. We term this as predictive mass shell removal.

Both dynamical mass shell removal and predictive mass shell removal are
computationally intensive. During the mass transfer calculation during the
binary evolution, we first use MESA standard photospheric removal. If the mass
loss removal approaches sonic, we switch to predictive shell mass removal. If
the mass loss approaches a dynamical timescale rate, we use dynamical mass
loss removal. However, the comparison with KR mass loss can only be made using
MESA standard photospheric removal.

##  3 Comparison of the formalisms.

To test our unified MT formalism, we use the Module for Experiments in Stellar
Astrophysics package, MESA (Paxton et al. (2011), Paxton et al. (2013), Paxton
et al. (2015), Paxton et al. (2018), Paxton et al. (2019), Jermyn et al.
(2023)), revision 23.05.1. MESA has integrated prescriptions for mass loss.
The option ”Ritter” uses the OR formalism, and the option ”Kolb” uses KR
formalism.

For our formalism, in addition to different mass loss rates, we evolve the
donors using the effective binary acceleration in a corotating frame, using
the subroutine published in Paper I. To provide smoother MT calculations, we
replaced the tables provided in Paper I with more refined ones: the new tables
are obtained for the mass ratios $-2\leq\lg_{10}q\leq 2$ with steps of 0.01 in
$\lg_{10}q$, and 500 points resolution between $L_{1}$ and $L_{\rm out}$
potential instead of 100 points in the original tables. Effective
accelerations and volume-equivalent radii for all equipotentials were found
numerically as in Paper I, with the precision $10^{-5}$ of its value or
better. In all the cases (testing OR, KR, or our formalism), the simulations
stopped when the donor reached $R_{\rm lout}$, which is the volume-equivalent
radius of the outer Lagrangian equipotential.

Figure 1: Relative density profiles for several atmospheric models in a
radiative (the left panel) and a convective donor (the right panel), as a
function of equipotential. The atmospheric models as described in §2.2,
specifically: EG is standard Eddington grey, EG$\Phi$ is the model where the
effective binary acceleration is non-constant, EG$\Phi\kappa$ is the model
where both the effective binary acceleration and the opacities are non-
constant, EGs is the default MESA model with a correction for non-planar
atmosphere, Isothermal is an isothermal atmosphere as in Ritter’s formalism,
and “Adiabatic” is when the entropy of all atmospheric layers is the same as
that of the photosphere. The initial mass of the donor is $30$ $M_{\odot}$,
the mass ratio is 6. The atmospheres were obtained when the donor’s radius is
$22.3$ $R_{\odot}$ (its effective temperature is $27200$ $K$) and $1340$
$R_{\odot}$ (its effective temperature is $3580$ $K$).

###  3.1 Atmospheres

Different atmospheric models (as described in §2.2) result in different
density and sonic velocity profiles at the $L_{1}$ plane. To demonstrate the
difference, we constructed different atmospheres. In each of the simulations,
the photosphere is located at the $L_{1}$ equipotential when we constructed
the profiles.

To be exact, EGs is the self-consistent atmospheric model that can be obtained
with MESA, which is termed there ”Eddington grey”. For the other models, we
build the density, pressure, and temperature profiles in the neighborhood of
the $L_{1}$ equipotential using the photospheric values of density,
temperature, and pressure provided by MESA as the initial condition. Those
models for fluid stratification do not extend to $\tau=0$ and do not provide
feedback on photospheric conditions.

The resulting density profiles, as functions of the relative effective
potential, are shown in Figure 1.

The comparison of isothermal and adiabatic atmospheres shows that the
isothermal atmosphere has a more shallow power-law density decrease at large
distances, and hence extends further. The adiabatic atmosphere has a stronger
core – i.e., its density near the $L_{1}$-neighborhood remains larger – and
then it abruptly drops. We also note that an isothermal atmosphere, while
being built using the assumption that the $P/\rho$ ratio remains constant,
does not result in an atmosphere that is truly isothermal. For example, as
density drops by two orders of magnitude in Figure 1, the temperature obtained
from the full equation of state drops from 3580K to 3190K for the convective
donor (and would drop to 2380K if density drops by three orders of magnitude).
The reason for a drop in temperature is that in the atmosphere, a significant
fraction of pressure is provided by radiation, with all pressure being
provided by radiation at the “true” surface ($\rho=0$). It serves as a
reminder that employing the equation for an ideal gas in the classic limit for
the atmosphere, as is adopted in OR formalism, can be inherently inconsistent,
and its utilization is best approached with caution.

The comparison of the EG-based models shows that, at least for the tested
models, the difference between the standard (plane-parallel) atmosphere (model
EG) and the one that introduces a correction for geometry (model EGs) is
insignificant. Varying opacities for the above models do not provide a
significant effect either (not shown in Figure 1, as it blends with the EGs
model). The most significant effect on the density profile is when the
effective acceleration is a function of the effective binary potential (model
EG$\Phi$). Finally, an additional correction to that model is provided by
using opacities that are a function of the local temperature and density
(model EG$\Phi\kappa$). EG$\Phi\kappa$ behavior has some adiabatic properties
when very close to $L_{1}$, and resembles an isothermal power law at larger
distances. Interestingly, for a convective donor, EG$\Phi\kappa$ results in
the most “isothermal” stratification among all considered models, in the sense
that for the given drop in density by 100 times from the photospheric value,
the decline in temperature is the smallest – it drops to 3290K. For a
radiative donor, the EG$\Phi\kappa$ atmosphere resembles most the EG
atmosphere.

The latter model, EG$\Phi\kappa$, is one that we adopt in what follows as the
default atmospheric model for obtaining atmospheric overflow properties, while
the isothermal and adiabatic models are used for some comparisons.

Figure 2: Atmospheric mass loss rates as a function of distance to the Roche
lobe. A radiative donor is shown in the left panel, and a convective donor is
shown in the right panel. The models are described in §3.2, specifically: II
means Isothermal atmosphere and Isothermal stream, IA means Isothermal
atmosphere and Adiabatic stream (almost identical to II in the left panel), AA
means Adiabatic atmosphere and Adiabatic stream, EGA means EG$\Phi\kappa$
atmosphere and Adiabatic stream, and OR stands for Original Ritter formalism.
The initial mass of the donor is $30$ $M_{\odot}$.

###  3.2 Atmospheric overflow

Here, we will compare the mass loss rates for atmospheric overflow with the
following assumptions:

  * •

OR: Original Ritter formalism

  * •

II: Isothermal atmosphere and Isothermal stream

  * •

IA: Isothermal atmosphere and Adiabatic stream

  * •

AA: Adiabatic atmosphere and Adiabatic stream

  * •

EGA: EG$\Phi\kappa$ atmosphere and Adiabatic stream

In the Ritter formalism, the total mass rate is obtained as

$\dot{M}_{\rm at}^{\rm Ritter}=\rho_{L1}c_{s}Q\ .$ (14)

Here, $\rho_{L1}$ is the density at the $L_{1}$ point, $c_{s}$ is the sonic
velocity at the photosphere, and $Q$ is the effective cross-section of the
flow. The sonic velocity in the original Ritter formalism uses the ideal gas
in a classic limit, with the assumption that the stream evolves isothermally.

In Figure 2, we show the comparison of the original Ritter formalism OR with
the II,AI, AA and EGA models. As expected, the difference between the $II$ and
$IA$ cases is minimal, and confirms the expected ratio (see S 2.2).

One would have expected that OR and II would produce the same MT rate, but
that is not what we find. A significant part of the difference, a factor of
few, comes from the definition of the isothermal sonic velocity in our method,
which is different from the one adopted in OR ($\sqrt{P/\rho}$ in our case and
$\sqrt{{\cal{R}}T/\mu}$ in OR case, as they have specifically considered cold
donors). The MT rates in OR are proportional to the sonic velocity to the
power of 3, so our rates are always higher for hot donors, where the different
values of isothermal sonic velocity can cause a difference in MT rates by a
factor of 3-4. Only a linear proportionality can be easily seen in Equation14.
The other two powers are hidden in the quantity $Q$. This quantity was
originally proposed in Meyer & Meyer-Hofmeister (1983). It synthesizes several
approximations, where the two most important ones are the Roche lobe
approximation and the assumption on the density drop off with the distance
from $L_{1}$ in units of pressure scale height. In the latter assumption, the
sonic velocity implicitly enters, to the power of two. The remainder of the
difference comes from the simplification of the Roche lobe geometry and the
value of the pressure scale for the isothermal atmosphere; it is hard to
disentangle the detailed role of each of the remaining approximation terms for
comparing the used assumptions with the direct integrations.

Accounting for realistic atmospheric profiles for radiative donors results in
higher mass loss rates in all atmospheric models as compared to OR when the
star is close to RLOF, with the EGA and AA models having lower MT rates when
the donor is still far from the RLOF. In convective donors, the difference is
smaller than in radiative donors for most atmospheric models. Usually, EGA
mass loss rates are slightly less than those predicted by OR for the same
fraction underflow.

###  3.3 Unified formalism in unperturbed donors

Figure 3: Mass loss rates as a function of distance to the Roche lobe. A
radiative donor is shown in the left panel, and a convective donor is shown in
the right panel. The models are described in §3, specifically: ”Unified” means
our unified scheme without taking into account the effective binary potential,
“OR” stands for OR formalism, and “KR” stands for KR formalism. KR is shown
only for the case of Roche overflow, while ”Unified” and OR are shown both for
Roche underflow (atmospheric mass loss rate) and Roche overflow cases. The
bottom branch in each case is for Roche underflow, and the upper branch is for
Roche overflow. The initial mass of the donor is $30$ $M_{\odot}$, and the
mass ratio $q=6$.

For this comparison, we select to use the EGA atmospheric model for our
unified model, although for very fast mass loss, the AA atmospheric model is
also acceptable. In Figure 3, we show the comparison of our mass loss rate
with the original Ritter formalism applied for both atmospheric overflow and
for the donor’s overflow, and optically thick mass loss rate using Kolb &
Ritter formalism (Kolb & Ritter, 1990). For this Figure, we found the MT rate
obtained assuming that the same unperturbed donor produces mass loss as a
function of the assumed underflow or overflow. No effective binary potential
was used to modify the donor. For the comparison purposes of the prescriptions
only during the binary evolution, we use default MESA options except allowing
for L-S coupling. If this option is permitted, the binaries of the same
initial orbital separation may start the MT at different orbital separations.
No direct comparison of the MT rate prescriptions would then be possible.

As before, we find that for the considered radiative donor, the MT rates are
lower in our formalism than in OR formalism when the donor is further away
from the Roche lobe but becomes larger when the distance to the Roche lobe
decreases. For convective donors, the difference is minimal. For overflow MT
rates with very small overflow, the unified MT rates are larger than predicted
by either Ritter or Kolb. For a more significant overflow, the unified MT rate
is slightly larger in the radiative donor than by the Kolb prescription and is
slightly smaller in the considered convective donor.

##  4 Impact on the binary evolution

###  4.1 Effect of modified gravity

Figure 4: Comparison of $L$ and $T_{\rm eff}$ positions for donors that have
the same radius at the start of RLOF, but evolved using modified gravity
(green stars), full synchronization (red diamonds), or with the standard MESA
binary package only (blue circles). These are shown for binaries with a mass
ratio of 10. Solid lines show the evolution of single stars of the given
masses. Grey dotted lines correspond to the constant radii.

First, we considered the difference in the evolution before RLOF, with no
atmospheric mass loss allowed, to capture only the effect of modified gravity.
In addition to the already described standard binary evolution and our
modified gravity evolution, we use MESA’s standard binary package, including
tidal evolution and orbital synchronization. We consider five initial masses
for the donor star: $M_{\rm d}=1,3,7,15$ and 30 $M_{\odot}$. For each mass, we
pick up to 7 stellar models with the star’s radius from the set of
$5,10,20,50,100,200,500$, and $1000R_{\odot}$ (stars with some initial masses
never have either some small or some large radius values from this set). The
calculations were done for several mass ratios between 1 and 10, but for the
mass ratio $q=10$, the deviation between the evolutionary paths is most
noticeable on the HR diagram, see Figure 4.

We consider the cases of standard binary evolution and fully synchronized
cases to be the two limiting cases, representing no rotation and the maximum
allowed rotation. Our modified gravity case is expected to be close to the
fully synchronized case (however, see Paper I, where we showed that the
effective acceleration in a binary could differ from the solid rotation case
by up to 20 percent). Overall, modified gravity by itself does not seem to
affect the evolution of the donor significantly, and the change it provides is
within the expected range. In rare cases, the modified gravity evolution path
deviates slightly from being located precisely between the standard evolution
and fully synchronized. It usually occurs when the donor transitions from
having a radiative envelope to a convective envelope, as modified gravity
causes a slightly different evolutionary path and envelope depths.

Figure 5: Mass loss rates as a function of the Roche lobe overflow (the left
panel) and the lost mass (the right panel) for a 30 $M_{\odot}$ donor in a
binary system with mass ratios $q=1,3,6$ and 10, and with such initial orbital
separations that the radius of the donor is $R_{1}=$350 $R_{\odot}$ at RLOF.
Solid lines show our method and dashed lines show the KR prescription. In the
left panel, we show behavior when the donor is close to RLOF. In the right
panel, we show only the evolution of MT rates during the time when the donor
loses $1~{}M_{\odot}$.

###  4.2 Unified mass loss formalism

We considered the set of donors with masses of $1,3,7,15$ and 30 $M_{\odot}$,
the set of orbital separations from $a=10,20,40,80,160,320,640,1280,1800,2400$
$R_{\odot}$, and the set of mass ratios $6$ and $10$ (and several cases with
$q=1$ and $3$ and additional donor’s radii at RLOF). We used only those
orbital separations for each donor where the donor in the considered binary
system starts RLOF. The comparison was made between standard MESA binary setup
with KR mass loss and the same MESA binary setup using our formalism while
evolving the donor using binary effective acceleration.

We find that there are two kinds of MT endings. The first is when both
formalisms (KR and unified) lead to the overflow of the outer Lagrangian
equipotential, and no sonic mass loss occurs. This situation is more typical
for extended stars, although it also occurs in 1 $M_{\odot}$ donors with all
considered separations. The comparison between the formalisms can be done
either through the evolution of the MT episode, or for the amount of mass
removed before $R>R_{\rm lout}$.

For example, in Figure 5, we show the comparison of our revised rates and the
KR rates through the MT episode as a function of the RLOF and the lost mass.
The comparison is made for the donor of the same mass and radius at RLOF,
placed in binaries with different mass ratios. We find that here atmospheric
mass loss rates are always lower with our prescription for the same overflow,
but they are higher for the same lost mass. The difference is especially
noticeable, up ten times difference, for $q>1$. In the RLOF regime, our rates
also consistently differ from those of KR. In the shown mass transfer tracks,
sonic mass loss $\dot{M}_{\rm sonic}$ was never achieved, and in all cases the
donor has lost more than $1M_{\odot}$. We note that the two ways to show
tracks are not explicitly related, as the time spent at a specific overflow
with a corresponding mass loss rate is an independent quantity, not shown.

Figure 6: Mass loss rates as a function of the Roche lobe overflow (left
panel) and of the lost mass (right panel) for a $15M_{\odot}$ donor in a
binary system with mass ratio $q=6$. The color indicates the initial orbital
separation. Solid lines show our prescription, dashed lines show the KR
prescription. In both panels, we show only the evolution of MT rates while the
donor loses the outermost $1~{}M_{\odot}$, where some binaries end their MT
early. Black-bordered light red and light blue lines show the unified
formalism with mass shell removal for $a=20~{}R_{\odot}$ and
$a=160~{}R_{\odot}$.

The second ending is when both formalisms lead (numerically) to MT exceeding
the donor’s sonic mass loss, sometimes to as high a MT as a factor of 100
times $\dot{M}_{\rm sonic}$. The sonic situation occurs for stars close to
their main sequence, and the reason for the code to stop is usually an
inability to converge the donor’s stellar model. In this case we can not
compare the final (still converged) models, but can compare the MT rates
before the sonic situation occurs.

For binaries where the MT rate exceeds the surface sonic mass loss, we do not
find much difference in $\dot{M}_{\rm tot}-$remaining mass tracks during most
of the RLOF mass transfer, although the responses to overflow ($\dot{M}_{\rm
tot}-$overflow tracks) are different between the two formalisms. In Figure 6,
we show the comparison of our revised MT rates and the KR rate through the MT
episode as a function of the lost mass, in the case of $M_{\rm d}=15M_{\odot}$
and $q=6$, for several initial orbital separations. The sonic mass loss rates
are achieved during the shown tracks for binaries with the initial $a=20$
$R_{\odot}$, at $\log_{10}(\dot{M}_{\rm tot}/M_{\odot}$ yr${}^{-1})=-3.8$, and
$a=160$ $R_{\odot}$, at $\log_{10}(\dot{M}_{\rm tot}/M_{\odot}$
yr${}^{-1})=-3.2$. It can be seen that when the outer layers are removed at a
”supersonic rate”, the donor starts quickly shrinking instead of continuing to
expand; this effect is purely a numerical artifact.

In verifying runs, where the predictive mass removal is switched on when
$\dot{M}_{\rm tot}\geq 0.01\dot{M}_{\rm sonic}$, the MT rates show the same
response for most of the RLOF episode but provide a different $\dot{M}_{\rm
tot}$ for the same remaining donor mass and RLOF when the
$\dot{M}>10\dot{M}_{\rm sonic}$, see Figure 6. The predictive mass loss allows
us to evolve to $R_{1}=R_{\rm Lout}$, where the total mass lost by the end can
be either smaller or larger than in the case of no predictive mass loss. The
$\dot{M}_{\rm tot}$-overflow tracks appear very similar to ”subsonic” tracks
in more expanded donors.

In the more extended donors, where the sonic mass loss rate is not achieved,
the difference between the MT rates found using our formalism and the KR
formalism is similar to the cases shown in Figure 5, and can be seen for
$a=640$ $R_{\odot}$ and $a=1280$ $R_{\odot}$ in Figure 6. Usually, the mass
loss rates are higher for the same removed mass and lower for the same
overflow. The exception is the case of $a=$1280 $R_{\odot}$, where the
combined effect of the evolution of the donor in a binary gravity field and
the use of the EG$\Phi\kappa$ atmosphere, instead of isothermal, led to a
lower rate of atmospheric outflow for a long period of time, resulting in a
larger mass removed at a much lower mass loss rate. The RLOF for this binary
starts when the mass loss rate is $\log_{10}(\dot{M}_{\rm tot}/M_{\odot}$
yr${}^{-1})\approx-3.2$ for our method and $\log_{10}(\dot{M}_{\rm
tot}/M_{\odot}$ yr${}^{-1})\approx-3$ for the KR formalism. Switching the
formalism leads to a difference of about 100 times in both the rate of
atmospheric outflow and the time spent in outflow, for example with
$\dot{M}>10^{-8}M_{\odot}$ yr-1, in the case of unified mass loss, vs. the
case of KR formalism.

For the case with $a=640~{}R_{\odot}$, the RLOF starts at
$\log_{10}(\dot{M}_{\rm tot}/M_{\odot}$ yr${}^{-1})\approx-4.0$ for our
formalism and $\log_{10}(\dot{M}_{\rm tot}/M_{\odot}$ yr${}^{-1})\approx-3.1$
for the KR formalism. The main difference between the two formalisms in the
stream mass loss rates (as a function of the overflow) and the outcomes for
this donor is caused not by the MT prescription, but by evolving the donor
using the effective binary potential. While only the evolution while losing
the outermost $1M_{\odot}$ is shown, the endings of the binaries are very
different: the KR prescription leads to overflowing $R_{\rm lout}$ after
losing $0.3M_{\odot}$. In comparison, in the case of our unified prescription,
almost the entire envelope is transferred ($9.1M_{\odot}$) before the donor
reaches $R_{\rm lout}$. An opposite but similarly large difference is the case
with the donor of $30M_{\odot}$, $a=640R_{\odot}$ and $q=10$: a donor evolved
using our method can only transfer $\sim 1.1M_{\odot}$, but with the KR method
the donor lost $16.6M_{\odot}$ before it reached $R_{\rm lout}$.

Overall, we find that the new method always leads to the donor losing a
different amount of mass before the onset of a CE event than if the same
binary evolution model is used with the KR method for MT. We note that the
values of the mass lost in each case should be found by considering what the
accretor is (an average star, a neutron star, or a black hole) and adopting
the angular momentum loss mode appropriate to that specific binary.

##  5 Conclusion

We present a revised method to obtain mass loss rates in 1D stellar codes. Our
method aims to improve the treatment of mass loss at the onset of the MT
episodes, especially in systems with a high mass ratio where MT evolves into
an unstable regime.

Our formalism uses the 3D average properties of the donor in a binary system
for effective binary potentials and accelerations, and can be applied in the
range of mass ratios $10^{-2}\leq q\leq 10^{2}$. The mass loss rates are found
for stream MT and continuing atmospheric overflow simultaneously, integrating
over the $L_{1}$ plane. For atmospheric overflow, we use atmospheric solutions
obtained using opacities and gravitational acceleration as a function of the
volume-equivalent radii during RLOF. To find the stream mass loss rate during
RLOF, for the density and sonic velocity profiles of the stream at the $L_{1}$
plane, we use the donor’s density and sonic velocity profiles from the same
equipotentials, using both the donor interior and its atmosphere continuously.
Our MT formalism is not restricted to cold donors. We also discussed how the
mass has to be removed from the donor’s stellar model in the case of a fast
MT.

We find that our unified mass loss rates often differ by a factor of a few
from those predicted by the Ritter (1988) and Kolb & Ritter (1990) formalisms.
Specifically, we find that our atmospheric mass loss rates, depending on the
donor, can be larger by up to a factor of 10 or smaller by a factor of 100.
This vastly different atmospheric mass loss may be important for binaries that
evolve through long episodes of atmospheric overflow, like symbiotic binaries.
The difference between our method and the KR method for overflow MT rates is
smaller, and again, our rates can either exceed or be smaller than those
predicted by the KR formalism. In some cases, the primary cause of the
difference is the MT prescription. In other cases, it is the use of the binary
effective potential. Overall, we find that our method removes a different mass
before the moment the donor exceeds the outer Lagrangian equipotential.

The influence of the revised formalism on all existing mass-transferring
systems (donor masses, mass ratios, orbital separations, donor and accretor
types) is impossible to predict. The goal of this paper is to present the new
formalism and test it against the OR and KR formalisms using the same
assumptions about binary evolution. Realistically, the adopted model of the
angular momentum loss will play a significant role in the outcome. For
example, mass losses from the accretor’s outer Lagrangian point would speed up
the system’s shrinking, and the assumptions on how conservative the mass loss
is also will change the outcomes. Simulations using appropriate models of
angular momentum loss will be the subject of future studies of specific target
groups of binaries.

N.I. acknowledges funding from NSERC under Discovery Grant No.
RGPIN-2019-04277. S.K. acknowledges funding from a MITACS Globalink Research
Internship. This research was enabled in part by support provided by Prairies
DRI and the Digital Research Alliance of Canada (alliancecan.ca). The authors
thank an anonymous referee for their helpful comments that improved the
quality of the manuscript.

## References

  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * Hansen et al. (2004) Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar interiors : physical principles, structure, and evolution
  * Ivanova et al. (2020) Ivanova, N., Justham, S., & Ricker, P. 2020, Common Envelope Evolution, doi: 10.1088/2514-3433/abb6f0
  * Jermyn et al. (2023) Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15, doi: 10.3847/1538-4365/acae8d
  * Kolb & Ritter (1990) Kolb, U., & Ritter, H. 1990, A&A, 236, 385
  * Lubow & Shu (1975) Lubow, S. H., & Shu, F. H. 1975, ApJ, 198, 383, doi: 10.1086/153614
  * Meyer & Meyer-Hofmeister (1983) Meyer, F., & Meyer-Hofmeister, E. 1983, A&A, 121, 29
  * Mochnacki (1984) Mochnacki, S. 1984, The Astrophysical Journal Supplement Series, 55, 551
  * Nariai (1967) Nariai, K. 1967, PASJ, 19, 564
  * Paczyński (1969) Paczyński, B. 1969, Acta Astron., 19, 1
  * Paczynski (1976) Paczynski, B. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 75
  * Pavlovskii & Ivanova (2015) Pavlovskii, K., & Ivanova, N. 2015, MNRAS, 449, 4415, doi: 10.1093/mnras/stv619
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pourmand & Ivanova (2023) Pourmand, A., & Ivanova, N. 2023, ApJ, 952, 126, doi: 10.3847/1538-4357/acd4c1
  * Ritter (1988) Ritter, H. 1988, A&A, 202, 93
  * Tauris & van den Heuvel (2023) Tauris, T. M., & van den Heuvel, E. P. J. 2023, Physics of Binary Star Evolution. From Stars to X-ray Binaries and Gravitational Wave Sources, doi: 10.48550/arXiv.2305.09388
  * Webbink (1975) Webbink, R. F. 1975, PhD thesis, University of Cambridge, UK


A Confidence Interval Rule sec:confidence-intervals

The QPET rule outputs a point estimate on the $t$-quantile (i.e., the number
of documents to review), and with a point estimate of recall equal to $t$.
However, we surely should feel less good about these estimates if they are
based, say, on a sample of 10 positive examples than on a sample of 1000
positive examples.

A confidence interval is an estimate that consists of a pair of scalar values
with an associated confidence level, conventionally expressed as $1-\alpha$
[Chapter 2]hahn2011statistical. A confidence interval estimator specifies a
closed interval $[X,Y]$, where $X$ and $Y$ are random variables defined in
terms of the sample. We say that such an estimator produces a $1-\alpha$
confidence interval for a population value $z$ when the probability is at
least $1-\alpha$ that, over many draws of a random sample of the specified
type, the sample-based realizations $x$ of $X$ and $y$ of $Y$ are such that
$x\leq z\leq y$.

Confidence interval estimators are used in two ways in TAR. The first is to
use a power analysis on the estimator as a guide to sample size
ryan2013sample. A review manager will draw a random sample large enough to
guarantee that any confidence interval estimate produced from it will have
some property, such as a maximum margin of error of 0.05.

The second use of confidence intervals is in reporting, i.e., in making a
statistical claim based on a labeled random sample about the effectiveness
achieved by a TAR process. If the only use of the random sample is in
reporting, this is unproblematic. But if (as is common) the same sample is
used to decide when to stop the review, the reported confidence interval
estimate will have sequential bias webber2013sequential.

Quantile Confidence Intervals

As with point estimates, the quantile perspective can rescue confidence
intervals on recall from sequential bias. A quantile confidence interval is a
confidence interval on a quantile [Chapter 4]wilcox2011introduction. To avoid
distributional assumptions about the TAR process, we can use a nonparametric
quantile confidence interval estimator [Chapter 7]david2003order. This takes
the form of a pair of order statistics, $[D_{i},D_{j}]$. The estimator
determines the values $i$ and $j$ based on the quantile level $t$, sample size
$r$, and confidence level $1-\alpha$. It provides the guarantee that, with at
least $1-\alpha$ probability over draws of the random sample, the $t$-quantile
$b_{t}$ falls within the sample-based realization $[d_{i},d_{j}]$.

If an estimator of this form is available, we can define the stopping point of
a TAR review to be the value of $d_{j}$ for our positive random sample, and
have $1-\alpha$ confidence that the $t$-quantile in B falls within
$[d_{i},d_{j}]$. By definition of the $t$-quantile, we thus have $1-\alpha$
confidence that stopping at $d_{j}$ gives a recall of at least $t$.

For many uses of confidence intervals we want estimators that make the
interval narrow (the realization $d_{j}-d_{i}$ is likely to be small) and/or
symmetric (i.e., $d_{j}-b_{t}$ and $b_{t}-d_{i}$ are likely to be similar,
where $b_{t}$ is some point estimate). For a stopping rule, however, the most
important criterion is that $d_{j}$ is likely to be small, since this reduces
the number of documents the TAR process must review before stopping.

We can minimize the likely value of $d_{j}$ by using a nonparametric one-sided
upper confidence interval (UCI) on a quantile [Chapter 5]hahn2011statistical.
Such an interval has the form $[D_{0},D_{j}]$, where $D_{0}$ is the $0$th
order statistic, i.e., the lowest logically possible value. For us this is
$D_{0}=1$ (the lowest A-rank); so the interval is $[1,D_{j}]$. We refer to the
pair as an 1-s UCI, and the upper end of the interval $D_{j}$ as a 1-s UCB
(one-sided upper confidence bound).

The estimator must choose $j$ such that the realization $d_{j}$ will be, with
$1-\alpha$ probability, a $t$-quantile or higher. This is equivalent to
requiring a probability $1-\alpha$ or higher that fewer than $j$ elements of
positive random sample D have A-rank less than the $t-$quantile. Suppose there
are $R$ positives in B, and that our sample of positives is of size $r$. Then
our estimator should choose the smallest $j$ such that:

$\sum_{k=0}^{j-1}\frac{\binom{\lceil tR\rceil-1}{k}\binom{R-\lceil
tR\rceil+1}{r-k}}{\binom{R}{r}}\geq 1-\alpha$ (1)

In a TAR setting we do not know $R$. However, if $R$ is large relative to $r$,
the binomial distribution is a good approximation to the above hypergeometric
distribution [thompson1997theory, Chapter 3]. In this condition, we want the
smallest $j$ such that

$\sum_{k=0}^{j-1}\binom{r}{k}t^{k}(1-t)^{r-k}\geq 1-\alpha.$ (2)

In fact, we can use the binomial approximation safely even when we are not
confident that $r$ is small relative to $R$. For values of $t$ greater than
0.5, the fact that the binomial has larger variance than the hypergeometric
means that the $j$ chosen using the binomial approximation will never be less
than the one chosen using the hypergeometric. Values of recall less than 0.5
are rarely of interest, but if needed we could find a similarly conservative
value of $j$ for such a $t$ by running the summation from $r$ downwards
instead of $0$ upwards. Based on the above analysis, we define the Quantile
Binomial Confidence Bound (QBCB) stopping rule. Given a sample size $r$,
recall target $t$, and confidence level $1-\alpha$, it specifies stopping a
one-phase TAR process when the $j$th positive sample document is found. Here
$j$ is smallest integer such that $[1,D_{j}]$ contains the $t$-quantile from
the unknown population of positive examples with $(1-\alpha)$ confidence,
based on the binomial approximation.

##  The QBCB Rule and Recall Estimation

We observed that the recall goal $t$ can be used as a conservative point
estimate of recall at the QPET stopping point. By the same logic, $t$ is a
conservative point estimate of recall at the QBCB stopping point. If we prefer
an interval estimate, we can use a $1-\alpha$ one-sided lower confidence
interval (1-s LCI) (or one-sided lower confidence bound, 1-s LCB) estimator
[hahn2011statistical, Chapter 2]). This defines a pair $[L,1.0]$ where, with
probability at least $1-\alpha$ over random samples of size $r$, the
realization $[\ell,1.0]$ contains a desired population value. Given the
definition of $t$-quantile, we know that $[t,1.0]$ is a $1-alpha$ 1-s LCI on
recall at the QBCB stopping point. This interval estimate on recall may seem
unsatisfying: it is identical regardless of sample size. However, this simply
reflects the task we have set for the QBCB rule: stop as soon as one has
confidence that a recall goal has been met. Larger sample sizes translate to
earlier stopping, not a tighter 1-s LCI.

##  What to Expect from the QBCB Rule

We can also compute more conventional estimates of recall at the QBCB stopping
point. As long as those estimates depend only on $j$ (which is fixed as soon
as we choose $r$ and $t$) and not on $d_{j}$ (the actual A-rank at which we
stop), these estimates are not affected by sequential bias. These estimates
give insight into the behavior of the QBCB rule. Table 1 shows the QBCB values
of $j$ for recall goal 0.8, confidence level 95% ($1-\alpha=0.95$), and
selected sample sizes $r$ from 14 to 457. (The choice of the sample sizes is
discussed in Section LABEL:sec:cost-overruns.) Sample sizes 8 to 13 are also
included. However, with these sample sizes, the only 95% 1-s UCI based on
order statistics that includes the 0.8-quantile is the trivial interval
$[D_{0},D_{r+1}]$ = $[1,N]$ (using the convention that the $r+1$’th order
statistic for a sample of size $r$ is the maximum population value). So for
$r<14$, the QBCB value of $j$ is $r+1$, and the rule does not provide
meaningful stopping behavior. For these sizes we instead show $j^{*}=j-1=r$,
the largest non-trivial stopping point. We also show both the QBCB $j$ and
$j^{*}=j-1$ for the case r=21, discussed in Section LABEL:sec:cost-overruns.
Rows with QBCB $j^{*}=j-1$ values are indicated by ”*”. We show the value of
three estimates of recall based solely on $j$ or $j^{*}$. The first is a 95%
1-s LCI, but for recall rather than for the $t$-quantile. In particular, we
use the Clopper-Pearson exact interval [brown2001interval]. Second, we show
the plug-in estimate $j/r$ discussed earlier for the PET rule. Finally, we
show a 95% 1-s UCI on recall, again computed using the Clopper-Pearson method.
For rows with the QBCB $j$ value, the lower end of the 95% 1-s LCI is always
at or above 0.80, but fluctuates and is closer to 0.80 when the sample size is
larger. This reflects the fact that the Clopper-Pearson LCI is based on the
same binomial approximation used in the QBCB rule. The only difference is that
the QBCB computation solves for integer $j$ based on fixed real $t$, while the
LCI computation solves for real $t$ based on fixed integer $j$. The QBCB
requirement that $j$ be an integer means that the $t$ at the lower end of the
LCI is typically slightly more than 0.8, with the difference decreasing as $r$
increases and more $j$ values are available. The plug-in point estimates
(which are simply $j/r$ or $j^{*}/r$ depending on the row) for small sample
sizes are much higher than 0.8. We can think of these as the estimated recall
at which the naive PET rule would need to stop to achieve the same confidence
bound as the QBCB rule, and reflects how uncertain recall estimates from small
samples are. The last column shows a 95% 1-s UCI on recall at the QBCB
stopping point. This estimate shows that, as sample sizes increase, we slowly
become more confident that the QBCB stopping point will not have very high
recall. Section LABEL:sec:cost-overruns discusses why, counterintuitively, we
should want such confidence.

Table : Stopping points within a positive sample for the QBCB rule with 0.80 recall and 95% confidence estimates, plus three conventional estimates of recall at the stopping point. |  | Recall Estimators Applied at j
---|---|---
$r$ | $j$ | 95% 1-s LCI | Plug-in | 95% 1-s UCI
8 | 8* | [0.688, 1.000] | 1.000 | [0.000, 1.000]
9 | 9* | [0.717, 1.000] | 1.000 | [0.000, 1.000]
10 | 10* | [0.741, 1.000] | 1.000 | [0.000, 1.000]
11 | 11* | [0.762, 1.000] | 1.000 | [0.000, 1.000]
12 | 12* | [0.779, 1.000] | 1.000 | [0.000, 1.000]
13 | 13* | [0.794, 1.000] | 1.000 | [0.000, 1.000]
14 | 14 | [0.807, 1.000] | 1.000 | [0.000, 1.000]
21 | 20* | [0.793, 1.000] | 0.952 | [0.000, 0.998]
21 | 21 | [0.867, 1.000] | 1.000 | [0.000, 1.000]
22 | 21 | [0.802, 1.000] | 0.955 | [0.000, 0.998]
29 | 28 | [0.847, 1.000] | 0.965 | [0.000, 0.998]
30 | 28 | [0.805, 1.000] | 0.933 | [0.000, 0.988]
31 | 29 | [0.811, 1.000] | 0.936 | [0.000, 0.988]
37 | 34 | [0.804, 1.000] | 0.912 | [0.000, 0.978]
44 | 40 | [0.804, 1.000] | 0.909 | [0.000, 0.968]
50 | 45 | [0.801, 1.000] | 0.900 | [0.000, 0.960]
63 | 56 | [0.801, 1.000] | 0.889 | [0.000, 0.947]
76 | 67 | [0.803, 1.000] | 0.882 | [0.000, 0.937]
88 | 77 | [0.802, 1.000] | 0.875 | [0.000, 0.928]
106 | 92 | [0.801, 1.000] | 0.868 | [0.000, 0.918]
129 | 111 | [0.800, 1.000] | 0.861 | [0.000, 0.908]
158 | 135 | [0.800, 1.000] | 0.854 | [0.000, 0.898]
198 | 168 | [0.800, 1.000] | 0.849 | [0.000, 0.889]
255 | 215 | [0.800, 1.000] | 0.843 | [0.000, 0.879]
332 | 278 | [0.800, 1.000] | 0.837 | [0.000, 0.870]
457 | 380 | [0.800, 1.000] | 0.832 | [0.000, 0.860]