Issue 
A&A
Volume 641, September 2020



Article Number  A86  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202038499  
Published online  17 September 2020 
Internal circulation in tidally locked massive binary stars: Consequences for double black hole formation
^{1}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
email: bhastings@astro.unibonn.de
^{2}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
^{3}
Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Ave. Universidad s/n, Cuernavaca, Morelos 62210, Mexico
Received:
26
May
2020
Accepted:
19
July
2020
Context. Steadystate currents, socalled Eddington–Sweet circulation, result in the mixing of chemical elements in rotating stars, and in extreme cases lead to a homogeneous composition. Such circulation currents are also predicted in tidally deformed binary stars, which are thought to be progenitors of double blackhole merger events.
Aims. This work aims to quantitatively characterise the steadystate circulation currents in components of a tidally locked binary system and to explore the effects of such currents on numerical models.
Methods. Previous results describing the circulation velocity in a single rotating star and a tidally and rotationally distorted binary star are used to deduce a new prescription for the internal circulation in tidally locked binaries. We explore the effect of this prescription numerically with a detailed stellar evolution code for binary systems with initial orbital periods between 0.5 and 2.0 days, primary masses between 25 and 100 M_{⊙} and initial massratios q_{i} = 0.5, 0.7, 0.9, 1.0 at metallicity Z = Z_{⊙}/50.
Results. When comparing circulation velocities in the radial direction for the cases of a single rotating star and a binary star, it is found that the average circulation velocity in the binary star may be described as an enhancement to the circulation velocity in a single rotating star. This velocity enhancement is a simple function depending on the masses of the binary components and amounts to a factor of approximately two when the components have equal masses. After applying this enhancement to stellar models, it is found that the formation of double helium stars through efficient mixing occurs for systems with higher initial orbital periods, lower primary masses and lower mass ratios, compared to the standard circulation scenario. Taking into account appropriate distributions for primary mass, initial period and mass ratio, models with enhanced mixing predict 2.4 times more double helium stars being produced in the parameter space than models without.
Conclusions. We conclude that the effects of companioninduced circulation have strong implications for the formation of close binary black holes through the chemically homogeneous evolution channel. Not only do the predicted detection rates increase but double blackhole systems with mass ratios as low as 0.8 may be formed when companioninduced circulation is taken into account.
Key words: gravitational waves / binaries: close / stars: massive / stars: rotation / stars: evolution
© ESO 2020
1. Introduction
The existence of grandscale circulations in stars was first proposed by von Zeipel (1924), who noticed that in a single rotating star the surfaces of thermal equilibrium and surfaces of hydrostatic equilibrium do not coincide, meaning that both thermal and hydrostatic equilibrium cannot coexist anywhere in the star. This results in circulation currents being driven inside the star. The first quantitative description of such currents was given by Eddington (1929), while a more rigorous model was proposed by Sweet (1950), and therefore such currents are referred to as Eddington–Sweet circulation.
The fundamental cause of Eddington–Sweet circulation is that the effective gravity of a rotating single star becomes distorted due to a centrifugal force from a spherical pattern (as found in a nonrotating star) to an ovalised shape (where gravity is weaker at the equator). This distortion changes the thermal structure of the star through the Von Zeipel Theorem (von Zeipel 1924) and leads to the previously mentioned mismatch between thermal and hydrostatic equilibrium.
It is also possible for a star’s structure to be disturbed by a binary companion, therefore, in analogy to the single rotating case, a star in a binary system should also host internal circulations. The existence of circulations in binary stars stands on just as rigorous a mathematical footing as that for single rotating stars; however, in comparison to the rotating single star case, the characterisation of this binary circulation has received little attention, although significant efforts in this area have been made (Smith & Smith 1981; Tassoul & Tassoul 1982a).
The effects of companioninduced circulation are most apparent in close binary systems, where the main effect will be a mixing of material between the stellar core and envelope, adding fresh fuel to the burning zones (de Mink et al. 2009). As shown by Maeder (1987) and Yoon et al. (2006), such mixing can have dramatic effects on the evolution of stars. In extreme cases, the mixing is strong enough to prevent a chemical gradient developing within the star, causing homogeneous evolution where the star undergoes bluewards evolution in the Hertzsprung–Russel diagram and the expansion of the stellar envelope is repressed.
In the age of gravitational wave astronomy, there is renewed interest in binaryspecific mixing mechanisms as it has been proposed by Marchant et al. (2016) and Mandel & de Mink (2016) that double black hole merger events can originate from close massive binary systems with very efficient internal mixing such that they undergo homogeneous evolution leading to double helium stars. Without internal mixing each star in the binary system would expand until one component surpasses the second Lagrangian point of the binary orbit (socalled L2 overflow), leading to a common envelope situation which will result in stellar merger for such close binary systems (Menon et al., in prep.). The models of Marchant et al. (2016) and Mandel & de Mink (2016) did not include the circulation induced by a binary companion, therefore the range of binary systems that undergo homogeneous evolution may be affected by the inclusion of these mechanisms and could have strong repercussions for double blackhole merger events.
Evidence for internal mixing in binary systems exists for a handful of observed systems, the most convincing of which is VFTS352, a system of two 30 M_{⊙} stars in the Large Magellanic Cloud with an orbital period of 1.1 days (Almeida et al. 2015; AbdulMasih et al. 2019). When comparing VFTS352 to rotating single star models, it is found that the best match to the observed effective temperatures and luminosities is given by models that are rotating significantly faster than the observed rotation rates in the system. Rotational mixing, which increases its efficiency with rotation rate, is thus judged to be stronger in the components of VFTS352 than in single star models, suggesting that mixing processes in massive close binary stars may be more efficient than those in rotating single stars. Other binary stars that show evidence of progression along a chemically homogeneous evolutionary track are R145 (Shenar et al. 2017) and HD 5980 (Koenigsberger et al. 2014).
A significant proportion (Sana et al. 2012) of Otype stars are thought to be in fact binary stars, hence our understanding of massive stars requires a detailed description of the effects of binarity, with mixing processes being one of those effects. Internal mixing affects the luminosity and surface element abundances of stars, which are in turn used to calibrate stellar models and measure quantities such as mass and age by comparison to numerical models. Therefore a solid description of mixing processes is necessary for the detailed study of massive stars in general.
This work aims to characterise quantitatively the steadystate circulation currents in a tidally locked component of a binary system and numerically investigate the effects of this circulation in massive shortperiod binary systems using a onedimensional detailed stellar evolution code. To do this we shall utilise the mathematical results of Tassoul & Tassoul (1982a,b).
In Sect. 2 we first introduce the mathematical models of internal circulation for a single rotating star and for a binary component and summarise the results of Tassoul & Tassoul (1982a,b). A quantitative comparison between the two is made in Sect. 3 which allows us to model binary circulations in a onedimensional stellar evolution code. We investigate the implications of our theoretical results in Sect. 4.1 which outlines our numerical approach and Sect. 5 which presents the results of our simulations. Other mixing mechanisms that could feasibly operate and implications for double black hole merger events are discussed in Sect. 6. The final conclusions are give in Sect. 7.
2. A mathematical model of circulation
The problem of circulation inside a star is a complex question in 3dimensional fluid mechanics, therefore a number of simplifying assumptions must be made. We follow the mathematical formulation of Tassoul & Tassoul (1982a) to describe the circulation in the radiative zone of a tidally and rotationally distorted star. Their assumptions are as follows:

The orbit of the system is circular with each star rotating about its own axis with a period equal to the period of the binary orbit, such a system is said to be synchronised or tidallylocked.

The polar axes of the binary components are parallel to each other and perpendicular to their orbital velocity vector.

Each star is acted on only by the tidal force of its companion.

The system is detached.

The companion of each star is modelled as a point mass.

The stars are not significantly distorted away from sphericity.

The stars rotate as solidbodies.

The stars are on the zeroagemainsequence such that there are no buoyancy forces arising from chemical gradients.

The stars are nonmagnetic.

The circulation propagates in the radiative zone only.

The radiative zone is inviscid.
We now discuss the justifications, if they exist, for each assumption. Many shortperiod binary systems are observed to be synchronised and to have nearly circular orbits. This is explained by dynamical tides that circularise the orbit and enforce spinorbit coupling on timescales that are relatively short (Zahn 1975, 1977). For example the predicted synchronisation time of a binary star consisting of two 10 M_{⊙} stars with a period of 2 days is approximately 1% of the nuclear timescale (Zahn 1975). For the same system the circularisation timescale is of the order of 10% of the nuclear timescale (Hurley et al. 2002). Thus massive shortperiod binary systems are expected to be circularised and synchronised, and assumption 1 is valid in many systems.
For very luminous or extreme massratio systems, a star may be affected by the radiation of its companion, which can drive a weak circulation. Tassoul & Tassoul (1982c)find that this is unlikely to have any hydrodynamical effect or result in significant mixing, thus supporting assumption 3.
In a contact binary, Stȩpień (2009) envisioned that the circulation pattern must connect both stellar cores, and thus contact binaries require a different treatment. However, the impact of circulation currents on the evolution of a star is mainly determined by the circulation strength at the coreenvelope interface, as this is where fresh fuel is mixed into the core. In most cases the coreenvelope interface would not be greatly affected by contact, thereby contact and noncontact systems are expected to behave similarly, giving justification for assumption 4.
Assumptions 5 and 6 exist to simplify the problem. Strictly, they are only justified for systems with orbital separations much larger than the stellar radii.
Assumption 7 is valid owing to the action of the TaylerSpruit dynamo (Spruit 2002), which causes angular momentum to be transported throughout the stellar interior via magnetic interactions such that near solidbody rotation is achieved. The TaylerSpruit dynamo is used to explain relatively slow rotation rates of white dwarves (Suijs et al. 2008) and young pulsars (Heger et al. 2005).
A standard assumption is that stars enter the zeroagemainsequence with a uniform chemical composition, meaning that no buoyancy forces act upon the circulation. Although during the starformation process, material may become partially segregated, no process, such as nuclear burning, occurs in stars before the zeroagemainsequence that can introduce significant chemical gradients.
Assumption 9 is a simplifying assumption only, although magnetic fields can introduce asymmetric distortions that may influence circulation currents (Roxburgh 1963; Moss 1977). The consideration of magnetic fields is outside of the scope of this work.
Typical circulation velocities are not greater than and are mostly much slower than 10 cm s^{−1} (see Fig. B.1), while convective velocities are typically of the order 10^{4} cm s^{−1} (see also diffusion coefficients plotted in Fig. 5). Thus, the circulation currents are efficiently smoothed out by convection, and it is safe to ignore circulation in convective regions as in assumption 10.
Viscosity may be taken into account at the expense of complicating the final solutions. Given that Tassoul & Tassoul (1982b) find that their “results imply that the circulation velocities do not depend much on the values of the coefficient of viscosity”, we choose to work with the inviscid solutions.
The mathematical approach is to find the mean velocity of steadystate motion by expanding the Navier–Stokes equation about hydrostatic equilibrium to first order in the dimensionless parameter ϵ, which characterises the strength of rotation and is given by
$$\begin{array}{c}\hfill \u03f5=\frac{{\mathrm{\Omega}}^{2}{R}^{3}}{\mathit{GM}}=(\frac{\mathrm{\Omega}}{{\mathrm{\Omega}}_{\mathrm{crit}}}{)}^{2},\end{array}$$(1)
where Ω is the angular velocity of a star’s rotation, G the gravitational constant, R, M the stellar radius and mass respectively and Ω_{crit} the critical angular velocity such that the star becomes gravitationally unbound at this angular velocity. It is thus clear that ϵ can never be greater than 1 and that ϵ is a measure of the strength of rotation. In a synchronised binary, Kepler’s third law can be used to express the parameter ϵ as
$$\begin{array}{c}\hfill \u03f5=\frac{G(M+{M}_{\mathrm{comp}})}{{d}^{3}},\end{array}$$(2)
where M_{comp} is the mass of the companion and d the orbital separation. The detailed derivation of circulation in a rotating single and binary star is given in Tassoul & Tassoul (1982a,b), respectively, upon which our treatment is entirely based.
2.1. Circulation in a rotating single star
The potential of a barotrope (Krogdahl 1942) is
$$\begin{array}{c}\hfill V={A}_{2}{\varphi}_{2}(r)\frac{1}{3}{\omega}_{0}^{2}{r}^{2},\end{array}$$(3)
where A_{2} is a constant which ensures that there is a continuous potential inside and outside of the star, defined as
$$\begin{array}{cc}& {A}_{2}=\frac{1}{3}{\omega}_{0}^{2}\frac{5{R}^{2}}{3{\varphi}_{2}(R)+R{\varphi}_{2}^{\prime}(R)},\hfill \end{array}$$(4)
$$\begin{array}{cc}& {\omega}_{0}^{2}=\frac{\mathit{GM}}{{R}^{3}}\xb7\hfill \end{array}$$(5)
The function ϕ_{2}(r) satisfies the differential equation
$$\begin{array}{c}\hfill {\varphi}_{k}^{\u2033}+\frac{2}{r}{\varphi}_{k}^{\prime}+{\varphi}_{k}(4\pi G\frac{\rho {\rho}^{\prime}}{{p}^{\prime}}\frac{k(k+1)}{{r}^{2}})=0,\end{array}$$(6)
with k = 2, ρ being the mass density and the boundary conditions
$$\begin{array}{c}\hfill {\varphi}_{k}(r=0)={\varphi}_{k}^{\prime}(r=0)=0.\end{array}$$(7)
A dash denotes a derivative with respect to the radial coordinate ($p\prime \equiv \frac{\mathrm{d}p}{\mathrm{d}r}$).
By employing the expressions above, expansion about hydrostatic equilibrium gives the components (to 1st order in ϵ) of the circulation velocity in the radial and colatitudinal directions as
$$\begin{array}{cc}\hfill {u}_{r}(r,\mu )& =\u03f5{U}_{\mathrm{S}}(r){P}_{2}(\mu ),\hfill \end{array}$$(8)
$$\begin{array}{cc}\hfill {u}_{\mu}(r,\mu )& =\u03f5V(r)(1{\mu}^{2})\frac{\mathrm{d}{P}_{2}(\mu )}{\mathrm{d}\mu},\hfill \end{array}$$(9)
where μ = cos(θ) and θ is the angle from the star’s pole (colatitude), P_{2} is the second order Legendre polynomial, U_{S}(r) and V(r) are velocity magnitudes in the radial and colatitudinal directions respectively. The velocity magnitudes U_{S}(r) and V(r) are related to each other through the continuity equation, giving
$$\begin{array}{c}\hfill V(r)=\frac{1}{6\rho {r}^{2}}\frac{\mathrm{d}}{\mathrm{d}r}(\rho {r}^{2}{U}_{\mathrm{S}}(r)).\end{array}$$(10)
The most important quantity for internal mixing is the circulation velocity magnitude in the radial direction, given by
$$\begin{array}{c}\hfill {U}_{\mathrm{S}}(r)=\frac{2L{r}^{4}}{{G}^{2}{m}^{3}}\frac{n+1}{n1.5}[{h}^{\prime}+(\frac{2}{r}\frac{{m}^{\prime}}{m})h],\end{array}$$(11)
where
$$\begin{array}{c}\hfill h(r)={A}_{2}{\varphi}_{2}(r).\end{array}$$(12)
In these expressions L is the total stellar luminosity, m the Lagrangian mass coordinate, r the Eulerian radial coordinate and n is the effective polytropic index throughout the model, which is the value for which the pressure and density profiles fit the polytropic relation at a given point in the star.
We will now use Eqs. (8) and (9) and the velocity magnitudes for a Cowling point source model, as tabulated in Tassoul & Tassoul (1982b), to show the 3dimensional circulation pattern for a single rotating star, in Fig. 1. This star has electron scattering opacity, mass of 3 M_{⊙} and rotation parameter of ϵ = 0.1. For simplicity only the outermost circulation cells are shown, which define the radius of the star. The other circulation cells are nested within these and have almost the same morphology. It is seen that the currents rise at the pole and retreat back into the star at the equator. The circulation is symmetric about the equatorial plane. It is remarkable that this pattern is very ordered and simple, with currents fixed in a meridional plane. The circulation pattern is largely a product of the symmetry of the system.
Fig. 1. Circulation pattern of Eddington–Sweet circulation in a single star Cowling point source model with electron scattering opacity, M = 3 M_{⊙}, R = 1.75 R_{⊙}L = 93 L_{⊙} and ϵ = 0.1. Only the outer circulation cells in the northern hemisphere are shown. The green arrows give the direction of the motion. The central convective region is marked in red, the pole is marked with a blue cross and the radiative envelope is shown in opaque grey. This figure is a reproduction of Fig. 3 of Tassoul & Tassoul (1982b). 
2.2. Circulation in a binary component
We will now give an overview of the solutions to the problem of a star which is both tidally and rotationally distorted as given by Tassoul & Tassoul (1982b). A binary system is described in standard spherical coordinates where the origin is the centre of mass of the star in question, the xaxis points towards the centre of gravity of the companion and the zaxis points towards the pole of the star in question. In these coordinates the plane described by φ = 0 is the subbinary meridian. The mathematical treatment of the binary case is identical to that of single star case, with of course a different gravitational potential. The potential, W, of a rotating star with a binary companion is then
$$\begin{array}{c}\hfill W=\frac{1}{3}{\mathrm{\Omega}}^{2}{r}^{2}[1{P}_{2}(\mu )]+\frac{G{M}_{\mathrm{comp}}}{d}{\displaystyle \sum _{k=2}^{4}}\frac{{r}^{k}}{{d}^{k}}{P}_{k}(\nu ),\end{array}$$(13)
as found by Chandrasekhar (1933) where μ = cos(θ), ν = sin(θ)cos(φ) and P_{k} represents the Legendre polynomial of order k. We note that for the unique case of pointmass companion, Eq. (13) is exact. In all other cases, the potential would contain an infinite summation.
The potential W is a combination of 4 different potentials, one arising from the rotation of the star (the first term) and 3 arising from the binary companion (the second term). The overall 3dimensional velocity field will then be a superposition of terms originating from each one of these potentials. This field is
$$\begin{array}{cc}\hfill {u}_{r}(r,\nu ,\mu )& =\u03f5[{U}_{\mathrm{S}}(r){P}_{2}(\mu )+{\displaystyle \sum _{k=2}^{4}}{U}_{k}(r){P}_{k}(\nu )],\hfill \end{array}$$(14)
$$\begin{array}{cc}\hfill {u}_{\mu}(r,\nu ,\mu )& =\u03f5[V(r)(1{\mu}^{2})\frac{\mathrm{d}{P}_{2}(\mu )}{\mathrm{d}\mu}+{\displaystyle \sum _{k=2}^{4}}{V}_{k}(r)(1{\mu}^{2})\frac{\partial {P}_{k}(\nu )}{\partial \mu}],\hfill \end{array}$$(15)
$$\begin{array}{cc}\hfill {u}_{\phi}(r,\nu ,\mu )& =\u03f5\left[{\displaystyle \sum _{k=2}^{4}}\frac{r{V}_{k}(r)}{\sqrt{(1{\mu}^{2})}}\frac{\partial {P}_{k}(\nu )}{\partial \phi}\right]\xb7\hfill \end{array}$$(16)
The velocity V_{k}(r) is related to U_{k}(r) by employing the continuity equation to give
$$\begin{array}{c}\hfill {V}_{k}=\frac{1}{k(k+1)}\frac{1}{\rho {r}^{2}}\frac{\mathrm{d}}{\mathrm{d}r}(\rho {r}^{2}{U}_{k}).\end{array}$$(17)
As in the single star case, each contribution to this velocity field consists of a dimensionless geometric term and a velocity magnitude term. The first terms of Eq. (14) and (15) are exactly those of Eqs. (8) and (9) and represent the circulation arising from the star’s rotation. The terms in the summations in the equations above are due to the potential of the binary companion. Because the potential of a binary star can be decomposed into a summation of 3 terms (see Eq. (13)), we also get 3 different contributions to the velocity field, each with its own inviscid circulation velocity magnitude, U_{k}(r), which is defined by the following equations
$$\begin{array}{cc}& {U}_{k}(r)=\frac{2L{r}^{4}}{{G}^{2}{m}^{3}}\frac{n+1}{n1.5}[{h}_{{k}^{\prime}}+(\frac{2}{r}\frac{{m}^{\prime}}{m}){h}_{k}],\hfill \end{array}$$(18)
$$\begin{array}{cc}& {h}_{k}(r)={c}_{k}{\varphi}_{k}(r),\hfill \end{array}$$(19)
$$\begin{array}{cc}& {c}_{k}={\omega}_{0}^{2}\frac{{M}_{\mathrm{comp}}}{M+{M}_{\mathrm{comp}}}{\left(\frac{R}{d}\right)}^{(k2)}\frac{(2k+1){R}^{2}}{(k+1){\varphi}_{k}(R)+R{{\varphi}^{\prime}}_{k}(R)}\xb7\hfill \end{array}$$(20)
Again ϕ_{k} can be found by solving Eq. (6) and c_{k} assures a smooth potential between the stellar interior and exterior.
Equations (14)–(16) and (18)–(20) give the full description, inline with the initial assumptions, of circulation currents in a tidallylocked binary component. The dependence on the orbital separation is expressed in Eq. (20) through the ${\left(\frac{R}{d}\right)}^{(k2)}$ term and also in Eqs. (14)–(16) through the ϵ parameter, describing the strength of rotation. In a tidallylocked system, rotation is related to the orbital separation by Kepler’s Third Law, Eq. (2). Therefore even when k = 2 the dependence on the orbital separation is retained.
Plotted in Fig. 2 is this 3dimensional velocity field for a binary component again with electron scattering opacity, M = 3 M_{⊙}, ϵ = 0.1 and assuming a companion mass of M_{comp} = 3 M_{⊙}. We ignore the k = 3, 4 terms, as they contribute little to the overall circulation pattern (see Tassoul & Tassoul 1982a and Appendix B). Figure 2 shows that this is a very different pattern to the single rotating star shown in Fig. 1. The circulation pattern is symmetric about the equatorial plane and about the plane passing through the poles and subbinary meridian whereas the Coriolis force could potentially change this, it does not enter into the linear order theory presented here. As our goal is to model the circulation in a binary star using a onedimensional stellar evolution code, we must make some kind of measurement of these two patterns to compare them.
Fig. 2. Circulation pattern of a rotationally and tidally distorted Cowling point source model with electron scattering opacity, M = 3 M_{⊙}, M_{comp} = 3 M_{⊙}, R = 1.75 R_{⊙}L = 93 L_{⊙} and ϵ = 0.1. Three different perspectives are shown. The k = 3, 4 contributions to the velocity field due to the companion have been omitted. The companion is shown as a yellow sphere at an arbitrary orbital separation. The green arrows give the direction of the motion. The central convective region is marked in red, the pole is marked by a blue cross and the radiative envelope is shown in opaque grey. 
2.3. Limitations of the mathematical model
By far the most constricting assumption in the mathematical description of internal circulation is that the star is not significantly distorted away from sphericity. As binary stars must have small orbital separations to become tidally locked early enough during their mainsequence evolution, the departures from spherical symmetry may be nonnegligible. However, for our example system that is discussed in Sect. 5.1, the tidal deformation is smaller than 3% and 5% for the primary and secondary, respectively, based on the deformation parameter given in Zahn (2008).
The mathematical model also relies on the Roche approximation, where the radii of the components are much smaller than the orbital separation so that the companion may be treated as a point mass. Our example system (Sect. 5.1) exhibits ratios of stellar radius to orbital separation of approximately 0.35 and 0.30 for the primary and secondary star respectively, so the Roche approximation may break down for massive tidally locked binaries.
The solutions presented here are theoretical constructs and despite the fact that the assumptions are mostly justified, many phenomena exist which could alter the flowpatterns in Figs. 1 and 2. One of the most likely is the baroclinic instability which can set up wavelike azimuthal motions (see Appendix B of Tassoul & Tassoul 1982b). Other processes include horizontal turbulence (Chaboyer & Zahn 1992) and shear mixing which may only be investigated by a more sophisticated mathematical procedure.
3. A comparison of circulation patterns
For convenience Table 1 gives a summary of the different functions defining the circulation velocity in the radial direction for the single star and binary star cases. See Sect. 2 for an overview of these expressions.
Summary of the functions defining the circulation in the radial direction, u_{r}, for the single and binary star cases.
Even though the circulation patterns inside a star are not onedimensional, the effect of horizontal turbulence smooths out any latitudinal or longitudinal dependent chemical gradients so that the effective transport of chemicals can be considered a onedimensional problem (Chaboyer & Zahn 1992). To analyse the circulation patterns in one dimension, the rootmeansquared (rms) velocity is calculated in the radial direction over the whole surface of the star. The rms velocity must be considered because the mean velocity in the radial direction for any closed circulation pattern is zero. The rms radial velocity, denoted v_{rms} is
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}=\frac{\sqrt{{\int}_{\phi =0}^{\phi =2\pi}{\int}_{\theta =0}^{\theta =\pi}{u}_{r}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi}}{\sqrt{{\int}_{\phi =0}^{\phi =2\pi}{\int}_{\theta =0}^{\theta =\pi}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi}}\xb7\end{array}$$(21)
3.1. A rotating single star
Using Eqs. (8) and (21) to calculate v_{rms} gives
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}=\sqrt{{\u03f5}^{2}{U}_{\mathrm{S}}^{2}(r)}\sqrt{\frac{1}{2{\pi}^{2}}}\sqrt{{\int}_{\phi =0}^{\phi =2\pi}{\int}_{\theta =0}^{\theta =\pi}{P}_{2}{(cos(\theta ))}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi},\end{array}$$(22)
which results in
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}=\u03f5{U}_{\mathrm{S}}(r)\sqrt{\frac{11}{32}}\approx 0.59\u03f5{u}_{\mathrm{S}}(r)\xb7\end{array}$$(23)
This result is used below to compare the circulation in a binary star to that in a single star.
3.2. A binary star
The inviscid radial velocity magnitudes, U_{k} and U_{S} are defined by the solution to the differential Eq. (6). It is not evident how these solutions vary with k and therefore it is in general difficult to compare the functions U_{k} and U_{S}. The exception is U_{2} and U_{S}, which both use the k = 2 solution to Eq. (6). When both the single star and binary component do not suffer significant deformation (which is a starting assumption in Sect. 2), they have the same structure and we see that U_{2} and U_{S} differ only by the constants c_{2} and A_{2} (see Table 1), which are related to each other by
$$\begin{array}{c}\hfill {c}_{2}=3\frac{{M}_{\mathrm{comp}}}{M+{M}_{\mathrm{comp}}}{A}_{2}.\end{array}$$(24)
In terms of q_{*}, the effective mass ratio, q_{*} = M_{comp}/M, this expression becomes
$$\begin{array}{c}\hfill {c}_{2}=3\frac{{q}_{\ast}}{{q}_{\ast}+1}{A}_{2}.\end{array}$$(25)
We note that q_{*} will change depending on which member of the binary is being considered, so that q_{*} > 1 for the least massive star and q_{*} < 1 for the most massive star in the system.
Using Eq. (25) we may write
$$\begin{array}{c}\hfill {U}_{2}=3\frac{{q}_{\ast}}{{q}_{\ast}+1}{U}_{\mathrm{S}}=Q{U}_{\mathrm{S}},\end{array}$$(26)
where for convenience we have set $Q=3\frac{{q}_{\ast}}{{q}_{\ast}+1}$. Armed with this relation, two extreme cases shall be considered, where U_{2} dominates the companioninduced velocity field such that U_{2} ≫ U_{3} ≫ U_{4} and the case where all terms contribute equally such that U_{2} ≈ U_{3} ≈ U_{4}.
3.2.1. Extreme case (i): U_{2} ≫ U_{3} ≫ U_{4}
When we neglect U_{3} and U_{4} and apply the relation expressed in Eqs. (26), (14) becomes
$$\begin{array}{c}\hfill {U}_{r}(r,\nu ,\mu )=\u03f5{U}_{\mathrm{S}}(r)[Q{P}_{2}(\nu )+{P}_{2}(\mu )].\end{array}$$(27)
To calculate v_{rms} we compute
$$\begin{array}{cc}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}& =\u03f5{U}_{\mathrm{S}}(r)\sqrt{\frac{1}{2{\pi}^{2}}}[{{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}{Q}^{2}({P}_{2}{(cos(\phi )sin(\theta ))}^{2}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}+{P}_{2}{(cos(\theta ))}^{2}+2Q{P}_{2}(cos(\phi )sin(\theta )){P}_{2}(cos(\theta )))\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi {]}^{\frac{1}{2}},\hfill \end{array}$$(28)
which has the analytic solution
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}=\u03f5{U}_{\mathrm{S}}(r)[\frac{49}{256}{Q}^{2}+\frac{11}{32}Q\frac{11}{32}{]}^{\frac{1}{2}}\xb7\end{array}$$(29)
3.2.2. Extreme case (ii): U_{2} ≈ U_{3} ≈ U_{4}
In the limit that U_{2} ≈ U_{3} ≈ U_{4}, Eq. (14) becomes
$$\begin{array}{c}\hfill {U}_{r}(r,\nu ,\mu )=\u03f5{U}_{\mathrm{S}}(r)[Q{\displaystyle \sum _{k=2}^{4}}{P}_{k}(\nu )+{P}_{2}(\mu )].\end{array}$$(30)
The calculation is provided in Appendix A, giving v_{rms} as
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}\approx \u03f5{u}_{\mathrm{S}}(r)\sqrt{0.3410{Q}^{2}+0.343750.2295Q}.\end{array}$$(31)
It is apparent that Eqs. (23), (29) and (31) describing the circulation velocities in one dimension for a single rotating star and a binary component are of the same functional form, that is proportional to ϵU_{S}(r) and a factor depending on the effective mass ratio, q_{*}. To compare how much circulation is enhanced in the binary case, we can divide Eqs. (29) and (31) by Eq. (23) to give the factor by which a binary companion increases the circulation velocity in a binary star compared to an equivalent single star with the same ϵ parameter. We shall name this the velocity enhancement factor, f_{ES}.
Figure 3 shows this enhancement factor, which remarkably demonstrates that there is not much difference between the two extreme cases considered. For q_{*} < 1 there is almost no difference between the two cases while for q_{*} > 1 the difference amounts to approximately 10%. For the least massive star in the binary system, q_{*} > 1 and for the most massive star q_{*} < 1, so that circulation in the least massive star is enhanced more than it is in the most massive star. This follows simply from the fact that the most massive star imparts a larger gravitational disturbance on its companion than its companion does to it.
Fig. 3. Factor by which averaged Eddington–Sweet circulation velocity is enhanced in a tidally locked binary component as a function of the effective mass ratio, q_{*} for two approximations as given in the legend. 
In the limit that the companion mass, M_{comp} tends to 0, Eqs. (29) and (31) revert to (23), and as one would expect, the velocity enhancement factor becomes 1. It is curious to note that in the limit that the companion mass becomes infinite, q_{*} tends to infinity but $Q=3\frac{{q}_{\ast}}{{q}_{\ast}+1}$ tends to −3, hence the velocity enhancement factor has a limit of approximately 3.0 and 3.4 for our two extreme cases respectively.
A decision must be made as to which approximation is the most suitable. Using polytropic models, one can show that in general the condition U_{2} ≫ U_{3} ≫ U_{4} is the most realistic (see Appendix B). Such a result could also be inferred by studying the circulation patterns arising from U_{2}, U_{3} and U_{4} plotted by Tassoul & Tassoul (1982a). Although inspection of Fig. 3 leads us to conclude that both cases give such similar results so as to make the choice of approximation nearly inconsequential.
The final result of our analysis is that for a fixed rotation parameter ϵ, in a tidally locked binary component the averaged circulation velocity, compared to a single rotating star is increased by a factor of
$$\begin{array}{c}\hfill {f}_{\phantom{\rule{0.333333em}{0ex}}\text{ES}}=\sqrt{\frac{441}{88}(\frac{{q}_{\ast}}{{q}_{\ast}+1}{)}^{2}+1+\frac{3{q}_{\ast}}{{q}_{\ast}+1}}\xb7\end{array}$$(32)
This expression is remarkable as it only depends on the masses of the stars in the binary, while the strength of the companioninduced circulation depends on the orbital separation and the component masses. As Eddington–Sweet circulation depends on a star’s rotation, for a full description of the circulation in a binary component one needs to know about the rotation rate, the orbital separation and the component masses. In a tidally locked binary, the orbital separation is set by the rotation rate, and so only the rotation rate and component masses are required. The rotation rate sets the strength of Eddington–Sweet circulation and gives us the orbital separation, therefore the velocity enhancement factor needs only to depend on the component masses.
4. Numerical method
4.1. Model prescription
To investigate the implications of companioninduced mixing as derived in Sect. 3, version r10398 of the onedimensional detailed stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) is used. The implemented physics is the same as that of Marchant et al. (2016), where massive binary systems without enhanced mixing were computed. However since that work was done, the treatment of angular momentum loss from binary systems was improved. With the improved version of MESA, the results presented here for our models without enhanced mixing will vary slightly from those of Marchant et al. (2016).
The opacities used in our calculations are given by the OPAL project (Iglesias & Rogers 1993, 1996). For the nuclear reaction rates, we employ those of JINA REACLIB (Cyburt et al. 2010).
Convection is modelled using standard mixing length theory (BöhmVitense 1958) with mixing length parameter α = 1.5 with the Ledoux criterion. We implement semiconvection as described by Langer et al. (1983) with efficiency parameter α_{SC} = 1.0. Convective overshooting is modelled with a step function with α_{OV} = 0.345 as in Brott et al. (2011).
The stellar wind prescription follows that outlined in Yoon et al. (2006). For stars with Y_{S} < 0.4 massloss rates are given by the recipe of Vink et al. (2001). Where Y_{S} > 0.7, the prescription of Hamann et al. (1995) reduced by a factor of 10 is used. Where 0.4 < Y_{S} < 0.7, the massloss rate is determined by an interpolation between the two schemes. For both schemes a metallicity scaling of (Z/Z_{⊙})^{0.85} was used. Rotationally enhanced massloss is handled by the method of Heger & Langer (2000) whereby the relationship between massloss with and without rotation is given by
$$\begin{array}{c}\hfill \dot{M}(\mathrm{\Omega})=\dot{M}(\mathrm{\Omega}=0)[\frac{1}{1\frac{\mathrm{\Omega}}{{\mathrm{\Omega}}_{\phantom{\rule{0.333333em}{0ex}}\text{crit}}(\mathrm{\Gamma})}}{]}^{0.43},\end{array}$$(33)
with the critical rotation rate being
$$\begin{array}{c}\hfill {\mathrm{\Omega}}_{\phantom{\rule{0.333333em}{0ex}}\text{crit}}(\mathrm{\Gamma})=\sqrt{\frac{\mathit{GM}}{{R}^{3}}(1\mathrm{\Gamma})}\phantom{\rule{0.166667em}{0ex}}\end{array}$$(34)
and Γ the ratio of the stellar luminosity to the Eddington luminosity, Γ = L/L_{Edd}. As the model approaches the critical rotation rate, massloss rates become divergent. To avoid this issue, once Ω/Ω_{crit} reaches 0.98, we switch to an implicit massloss rate, which is calculated such that the fraction Ω/Ω_{crit} remains 0.98.
Internal angular momentum transport is accomplished via the TaylerSpruit dynamo (Spruit 2002), which enforces near solidbody rotation through magnetic interactions. This is implemented numerically by the method of Petrovic et al. (2005). Tidal synchronisation is modelled according to Hurley et al. (2002) and Detmers et al. (2008). The criteria governing Rochelobe overflow come from considering the Rochelobe volumes and their volumeequivalent radii as per the prescription of Eggleton (1983), while masstransfer rates are determined implicitly by demanding that the massdonor stays inside its Rochelobe.
The mixing processes taken into account are Eddington–Sweet circulation, secular and dynamical shear instabilities and the GSF instability. The rotational mixing efficiency factor f_{c} is set to 1/30, as determined theoretically by Chaboyer & Zahn (1992) and the factor controlling the effects of chemical gradients on mixing f_{μ} is set to 0.1 as calibrated by Yoon et al. (2006). The suppression of mixing by chemical gradients is included as a stabilising circulation with velocity proportional to the meanmolecular weight gradient and f_{μ} (Heger & Langer 2000). Circulation is modelled as a diffusive process according to Heger & Langer (2000), where the diffusion coefficient is calculated as the product of the circulation velocity and a typical length scale for the circulation. It is thus appropriate, as determined in Sect. 3.2, to model the effects of companioninduced circulation by increasing the diffusion coefficients of Eddington–Sweet circulation by a factor of
$$\begin{array}{c}\hfill {f}_{\mathrm{ES}}=\sqrt{\frac{441}{88}(\frac{{q}_{\ast}}{{q}_{\ast}+1}{)}^{2}+1+\frac{3{q}_{\ast}}{{q}_{\ast}+1}},\end{array}$$(35)
with q_{*} = M_{comp}/M being the effective mass ratio for the star in question. For the most massive star in the binary, q_{*} < 1, while for the least massive star q_{*} > 1. The circulation velocity enhancement is achieved through the D_ES_factor control parameter in MESA. For models without companioninduced circulation we set f_{ES} = 1.
We compute grids of binary models with primary mass, M_{1}, ranging from log(M_{1}/M_{⊙}) = 1.4–2.0 in intervals of 0.1 dex, such that the models are not massive enough to undergo the pairinstability supernova phenomenon. We choose initial orbital periods ranging from 0.5 to 2.0 days in intervals of 0.1 day. Four different initial massratios are considered, q_{i} = M_{2}/M_{1} = 1.0, 0.9, 0.7, 0.5 such that the mass of the secondary star, M_{2} is the product of the primary mass and the massratio. All models have an orbital eccentricity of 0 and rotation periods synchronised with the binary orbital period.
The calculations are performed at a metallicity of Z = Z_{⊙}/50 where the initial helium abundance is the result of a linear interpolation in Z between the primordial value Y = 0.2477 at Z = 0 (Peimbert et al. 2007) and Y = 0.28 at Z = Z_{⊙} where we take Z_{⊙} = 0.017 (Grevesse et al. 1996). This relatively low metallicity is chosen because for higher metallicity systems, massloss through stellar winds causes a widening of the binary orbit which may result in double blackhole systems that are too wide to merge within a Hubble time (cf. Fig. 4 of Marchant et al. 2016). We note that as we only investigate the relative effect of enhanced circulation on the models, and that the minimum initial rotation rate needed for homogeneous evolution to occur is approximately constant with varying metallicity (Yoon et al. 2006), we do not expect our analysis to be greatly dependent on the choice of metallicity.
Termination occurs either at core helium depletion or when one component suffers L2 overflow (upon whence a merger is assumed to happen) or when the difference between central and surface helium mass fraction exceeds 0.2 ( meaning that chemically homogeneous evolution is not occurring, upon whence L2 overflow is assumed to eventually happen). When one component depletes helium, it is treated as a point mass until the other component depletes helium and the calculation is terminated.
To determine the effects of the enhanced mixing in a binary system, we calculate two grids of models, without and with the companioninduced mixing described above. Models without enhanced mixing are named the Standard Grid, while those with enhanced mixing are named the Enhanced Grid. The inlists required to run the Standard Grid can be found in the binary test suite of MESA version r10398.
4.2. Uncertainties in numerical modelling
An issue is the treatment of massloss, which is particularly important in a binary system as the binary orbit is influenced by both mass and angularmomentum loss through stellar winds. We are considering stars with initial masses up to 100 M_{⊙}, observations of which are rare and difficult to interpret. Therefore it is not certain that the massloss rates, as found by Vink et al. (2001) and Hamann et al. (1995) for hydrogen rich and poor stars respectively, give an accurate picture of massloss in very massive stars.
Furthermore the 3dimensional mass loss pattern of a star in a close binary system is a complex affair. Due to the von Zeipel theorem (von Zeipel 1924), the temperature and flux across the stellar surface will be a function of both angular coordinates, and therefore so will the massloss rates. This may lead to a significant difference between effective massloss rates of single and binary stars.
To model contact systems, the effects of energy transfer from one star to another are ignored, which for systems with initial mass ratios far from unity can play a role in the stellar evolution of both components. It is theoretically predicted and observed (Kähler 2004; Martins et al. 2017) that for contact systems with radiative envelopes, the difference in effective temperature between the stars in the system is very small, even when the massratio is significantly different from unity. This suggests that energy is indeed transferred from one component to another during the evolution of a contact binary.
5. Model results
In this section we compare the results of our model grid without enhanced mixing and our model grid with enhanced mixing. The reader is reminded that these grids are straightforwardly named the Standard Grid and Enhanced Grid respectively.
5.1. An example system
To characterise the strength of mixing, the mixing timescale, τ_{mix} is calculated as (Schatzman 1977)
$$\begin{array}{c}\hfill {\tau}_{\phantom{\rule{0.333333em}{0ex}}\text{mix}}=\frac{4}{{\pi}^{2}}[{{\displaystyle \int}}_{0}^{R}\frac{\mathrm{d}r}{\sqrt{{D}_{\phantom{\rule{0.333333em}{0ex}}\text{tot}}}}{]}^{2},\end{array}$$(36)
where D_{tot} is the sum of diffusion coefficients of all mixing processes operating within the star. The mixing timescale is interpreted as the time which is required for a particle to be mixed from the stellar centre to the photoshpere. When the fraction of the nuclear and mixing timescale, τ_{nuc}/τ_{mix} is computed, the strength of internal mixing is easily assessed. The occurrence of homogeneous evolution is ruled by whether mixing processes can erode chemical gradients produced from nuclear burning at a fast enough rate such that chemical gradients do not inhibit mixing. Therefore chemically homogeneous evolution is expected to occur for values of τ_{nuc}/τ_{mix} ≳ 1.
In the following we shall use the convention that the initially more massive star is named Star 1 and the initially less massive star is Star 2. These names will never be switched.
To understand the effects of our companioninduced circulation model we shall examine in detail a system with a primary mass, M_{1} = 100 M_{⊙}, initial mass ratio q_{i} = 0.7 and initial period 1.5 days. For this system in the Standard Grid, Star 2 does not evolve homogeneously, whereas in the Enhanced Grid a double helium star is formed. Figure 4 shows the evolution of component masses, the ratios of nuclear and mixing timescales, τ_{nuc}/τ_{mix} and the difference between central and surface helium mass fraction for this example system in both Standard and Enhanced grids.
Fig. 4. Diagnostic plots for a system with primary mass, M_{1} = 100 M_{⊙}, initial massratio, q_{i} = 0.7, and initial orbital period 1.5 days. Top panels: system in the Standard Grid (without enhanced mixing) and bottom panels: system in the Enhanced Grid (with enhanced mixing). Left panels: evolution of each component’s mass and the mass ratio, q with a dotted line. Central panels: ratios of nuclear and mixing timescales as a function of time. Right panels: difference between central and surface helium mass fraction. Star 1 is plotted in red, Star 2 in blue. The black cross and diamond represent the times for which profiles of Star 2 are plotted in Fig. 5. 
We first consider Star 2 at early times. Figure 4b shows that in the Standard Grid Star 2 is unlikely to undergo homogeneous evolution as the mixing timescale is around double the nuclear timescale. In the Enhanced Grid however, the diffusion coefficients of Eddington–Sweet circulation are increased by a factor of 2.1, which decreases the mixing timescale such that it is approximately equal to the nuclear timescale as shown in Fig. 4e.
Looking at the chemical diffusion coefficient throughout the star in Fig. 5a, we see that at early times, the Enhanced and Standard Grid models differ only by the increased diffusion coefficients of circulation in the envelope. Here we can easily distinguish between the convective core, where the diffusion coefficients are very high, and the radiative envelope where circulation operates with lower mixing efficiency. Similarly in Fig. 5b, showing the helium mass fraction profiles at an age of 0.27 Myr, not enough time has elapsed to allow either helium to be mixed to the surface, nor a strong chemical gradient to develop. In the model with enhanced mixing, even though only a negligible amount of helium has been mixed to the surface, the helium gradient between the coreenvelope boundary and surface is shallower than in the Standard Grid model, as shown in Fig. 5b.
Fig. 5. Internal profiles of Star 2 in the example system (M_{1} = 100 M_{⊙}, q_{i} = 0.7, P_{i} = 1.5 days) showing the total diffusion coefficient, D_{tot} (left panels), and helium mass fraction, Y (right panels), as a function of Lagrangianmasscoordinate. The star in the Standard Grid is marked in green and that in the Enhanced Grid in black. Top panels: profiles at an age of 0.27 Myr, as marked by a cross in Fig. 4, bottom panels: age of 0.80 Myr, as marked by a diamond in Fig. 4. The grey dashed line in the right panels represents the initial helium mass fraction, Y_{0}. 
For the Standard Grid model this steeper helium gradient causes a slowing of the internal mixing through buoyancy forces, as evidenced by decreasing values of τ_{nuc}/τ_{mix} in Fig. 4b. Looking at the internal profiles at an age of 0.8 Myr in Fig. 4d, we see that when circulation is enhanced, appreciable amounts of helium have been mixed to the surface. In the standard mixing regime, a relatively small amount of helium has been mixed to the surface, however we see that a large step in helium abundance has developed at the coreenvelope boundary, which provides a solid barrier to mixing. The effects of this barrier can be seen in Fig. 5c which shows very low diffusion coefficients at the coreenvelope boundary. After this step in the chemical gradient has developed, mixing is shut off completely and the star evolves along the standard track, expanding until it interacts with its companion, causing the companion’s surface helium abundance to increase and causing the calculation to terminate, as evidenced by Fig. 4c.
During homogeneous evolution Fig. 4e shows that τ_{nuc}/τ_{mix} evolves to larger values, meaning that mixing becomes more efficient. This occurs because a star undergoing homogeneous evolution evolves to hotter temperatures and approaches the Eddington luminosity (Maeder 1987), which increases the light to mass ratio and thus causes the convective core to grow, this growth may be seen in panel c of Fig. 5 by comparing the enhanced model to the standard model. As the diffusion coefficients associated with convection are some seven orders of magnitude larger than those describing Eddington–Sweet circulation, the mixing timescale naturally decreases when the extent of the convective core increases. This results in the stars evolving to larger values of τ_{nuc}/τ_{mix} in Fig. 4e.
Our example system in the Enhanced Grid shows the typical behaviour of a system that forms a double helium star without either component ever filling its Rochelobe. Despite the lack of masstransfer, the massratio of the system evolves towards unity as the more massive star experiences higher massloss rates, as shown in Fig. 4d, giving a final massratio somewhat closer to unity than the initial massratio. However, owing to the fact that massloss from stellar winds is strongly metallicity dependant, very low metallicity systems will remain closer to their initial massratios.
5.2. Fates of systems
The major effect of our enhanced mixing prescription is to reduce the minimum rotation rate required for homogeneous evolution to occur, as shown in Fig. 6, where the minimum initial rotation rate for homogeneously evolving models is plotted as a function of mass for equalmass systems. This plot shows stars with masses from log(M_{1}/M_{⊙}) = 1.6 because q_{i} = 1.0 systems with primary masses less than this value did not evolve homogeneously (see left panels of Fig. 7). For these systems the circulation velocity is enhanced by approximately a factor of two in each component, which causes homogeneous evolution to occur at fractions of critical angular velocity some 10–20% lower, depending on the mass, than for systems without enhanced mixing.
Fig. 6. Minimum initial fraction of critical angular velocity required for homogeneous evolution to occur as a function of primary mass for systems with q_{i} = 1.0. Models from the Standard Grid are marked in black and those from the Enhanced grid in green. 
As stellar mass increases, the contribution of radiation pressure begins to dominate over gas pressure and the dimensionless adiabatic temperature gradient $({\mathrm{\nabla}}_{\mathrm{ad}}={\left(\frac{\mathrm{d}logT}{\mathrm{d}logP}\right)}_{\mathrm{ad}})$ decreases from 0.4, the value for an ideal gas to 0.25, that of a radiative gas (see Ch. 13.2 of Kippenhahn et al. 2012). The dependence on the subadiabaticity (∇_{ad} − ∇) in many formulations of Eddington–Sweet circulation (Kippenhahn 1974; Zahn 1983; Maeder 1987), including that of the models presented here^{1}, indicates that the stability of a fluid is described by the deviation from adiabaticity, with more stable systems having larger values of ∇_{ad} − ∇. Increased radiation pressure serves to reduce the dimensionless adiabatic temperature gradient, reducing the stability of the fluid and thus leading to circulation currents that propagate more easily than in the absence of radiation pressure. Thus, as the ratio of radiation pressure to gas pressure increases with stellar mass, more massive stars are expected to have more efficient circulation currents, hence resulting in lower rotation rates required to achieve homogeneous evolution, as seen in Fig. 6.
We identify three consequential effects of our enhanced mixing:

(i)
that stars of lower mass are able to evolve homogeneously;

(ii)
that homogeneous evolution occurs in systems with longer initial periods;

(iii)
fewer systems undergo masstransfer events and contact phases, owing to larger initial orbital separations resulting from (ii).
When the minimum rotation rate required for homogeneous evolution becomes low enough to be attained by a tidallylocked binary without suffering L2 overflow at the zeroagemainsequence, stars of lower mass can evolve homogeneously in the enhanced mixing scenario.
The lowering of the minimum rotation rate required for homogeneous evolution means that systems with longer periods can undergo strong mixing and form double helium stars in comparison to the standard mixing scenario, as in tidallylocked systems the orbital period determines the stellar rotation. A side effect of this is that with enhanced mixing, binary orbits are wider and the stars are able to expand more without interacting with their companion, meaning that masstransfer episodes become rarer, as discussed in detail in Sect. 5.4
The three effects of our enhanced mixing discussed above are visible in Fig. 7, which shows the results of our binary evolution calculations. As homogeneous evolution did not occur in any model with q_{i} = 0.5, these models are subsequently ignored. We see that enhanced mixing results in systems with larger initial periods evolving homogeneously, causes systems with lower primary mass to evolve homogeneously, and causes fewer systems to evolve into a contact phase. We note that for the Enhanced Grid models with q_{i} = 1.0 and log_{10}(M_{1}/M_{⊙}) = 2.0, we have not found the upper initial orbital period limit for homogeneous evolution to occur, at this initial massratio, we strictly consider a lower limit.
Fig. 7. Final configurations of models in the Standard Grid (top row) and Enhanced Grid (bottom row) for initial mass ratios q_{i} = 1.0, 0.9, 0.7. The colours indicate where overflow at the second Lagrangian point occurred at the zeroagemainsequence or during evolution (purple and green respectively), where the difference between central and surface helium mass fraction exceeds 0.2 (red) and where homogeneous evolution occurs resulting in a double helium star (blue). Single hatched cells indicate that the system underwent contact on the mainsequence, and crosshatched that contact occurred on the zeroagemainsequence. 
In the Standard Grid, we see that no system with initial period greater than 1.6 days produces a double helium star. In comparison our Enhanced Grid models are able to evolve homogeneously at initial periods as high as 2.0 days for systems with initial massratio q_{i} = 1, 1.9 days for q_{i} = 0.9 and 1.5 days for q_{i} = 0.7. At every initial massratio, for initial primary masses greater than log_{10}(M_{1}/M_{⊙})≈1.7 we see that the enhanced mixing increases the minimum initial orbital period required for homogeneous evolution.
5.3. Fractions of homogeneously evolving systems
The enhancement of mixing in the Enhanced Grid enables systems with lower primary masses to undergo homogeneous evolution, as compared to the Standard Grid. We see from Fig. 7 that the Enhanced Grid systems with q_{i} = 0.9 and primary masses as low as log_{10}(M_{1}/M_{⊙}) = 1.5 are able to form double helium stars, while in the Standard Grid this limit is log_{10}(M_{1}/M_{⊙}) = 1.7. Similarly for systems with q_{i} = 0.7, the minimum initial primary mass that evolves homogeneously is shifted from log_{10}(M_{1}/M_{⊙}) = 1.8 to 1.6. These relatively low mass double helium stars will have a large effect on predicted numbers of such systems owing to the initial mass function strongly favouring lower primary masses.
The predicted fractions of binary systems that will evolve to form double helium stars can be estimated by giving a weight to each model using a Salpeter distribution of masses for Star 1 (dN/dm ∝ m^{−2.35}) and a flat distribution in log_{10}(P). The weight of each model is given by
$$\begin{array}{c}\hfill W({M}_{i},{P}_{i})=\frac{{\int}_{{log}_{10}({P}_{i}\mathrm{\Delta}P)}^{{log}_{10}({P}_{i}+\mathrm{\Delta}P)}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}{log}_{10}(P){\int}_{{M}_{i}\mathrm{\Delta}M}^{{M}_{i}+\mathrm{\Delta}M}{m}^{2.35}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}m}{{\int}_{{log}_{10}(0.5)}^{{log}_{10}(2.0)}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}{log}_{10}(P){\int}_{{10}^{1.4}}^{{10}^{2.0}}{m}^{2.35}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}m},\end{array}$$(37)
where ΔP is the half difference between the period spacing of the models (0.05 days in this work), and ΔM is the corresponding half difference between the mass spacing of the models in log space, such that M_{i} + ΔM = 10^{(log(M/M⊙)+0.05)}
Summing the weights of models that formed double helium stars gives an estimate of the fraction of systems forming a double helium star with initial period spanning 0.5–2.0 days, Star 1 mass spanning log(M/M_{⊙}) = 1.4 − 2.0 at each initial mass ratio. These values are tabulated in Table 2. Finally by assuming a flat initial massratio distribution, integrating over q_{i} gives the fraction of binary systems predicted to form double helium stars in our parameter space. We calculate that the models with enhanced mixing produce 2.4 times more double helium stars than those without.
Fractions of binary systems in our model grids that are predicted to form double helium stars for each initial massratio.
5.4. Occurrence of masstransfer
As systems with larger orbital periods undergo homogeneous evolution, the separation of the stars can become large enough to avoid one component evolving to fill its Rochelobe. We note that although homogeneously evolving stars do increase their radii slightly during hydrogen burning, such that stars that are detached at the zeroagemainsequence may evolve into contact, this increase is very small compared to that seen during standard evolution.
In our Standard Grid models only a few systems with primary masses near 100 M_{⊙} and initial massratio near 1 manage to avoid contact phases. However in the Enhanced Grid, owing to the larger orbital periods for homogeneously evolving systems, many systems are able to avoid binary interaction. For the Enhanced Grid models with q_{i} = 1.0, approximately half of the systems forming double helium stars do so without any contact phases occurring. While the number of systems avoiding interaction decreases as the initial massratio becomes more extreme, such systems are remarkable because they produce double helium stars with massratios close to their initial massratio. In the Enhanced Grid, double helium stars that avoid interaction are produced at initial massratios as extreme as q_{i} = 0.7, whereas in the Standard Grid no noninteracting systems exist at q_{i} = 0.7 and only two are present at q_{i} = 0.9.
Table 3 gives the proportion of double helium star systems that do not undergo masstransfer. We see that only 5% of double helium stars in the Standard Grid avoid masstransfer episodes, while the for the Enhanced Grid this figure is around 20%. This means that in the Standard Grid, almost all of the double helium star systems have a massratio of unity, whereas for the Enhanced Grid, some 20% of double helium stars remain close to their initial massratio.
Of systems that form a double helium star, the fractions of those that do not undergo masstransfer for each initial mass ratio.
To explore the final massratios of our model systems, we plot the mass of each component at core helium exhaustion in Fig. 8. Systems with q_{i} = 1 remain at their initial mass ratio and so are not plotted. If it is assumed that at the end of helium burning, the star collapses directly to a black hole without any mass loss, Fig. 8 shows the predicted massratios of binary black holes formed through the chemically homogeneous channel. Comparing the Standard with the Enhanced Grid, we see that the enhanced mixing enables binary black holes with massratios as low as 0.8 to be produced. Without enhanced mixing the final systems all have massratios very close to unity. It is also seen that the Enhanced Grid produces systems with component masses as low as 30 M_{⊙}, whereas the Standard Grid can only make carbon stars with masses not lower than 40 M_{⊙}.
Fig. 8. Masses of Star 1 (xaxis) and Star 2 (yaxis) at core helium depletion for systems with q_{i} = 0.9, 0.7 (plotted as triangles and squares respectively). Models in the Standard Grid (without enhanced mixing) are marked in black, those in the Enhanced Grid (with enhanced mixing) in green. For clarity systems with q_{i} = 1 are not plotted as they remain at their initial mass ratio. The grey dotted lines show various massratios as marked. Red, blue, yellow and purple stars represent measured progenitor masses of gravitational wave events GW170729, GW170823, GW150914 and GW170818 respectively (Abbott et al. 2016; The LIGO Scientific Collaboration & The Virgo Collaboration 2019). 
6. Discussion
6.1. Other possible mixing mechanisms
One should keep in mind that grandscale circulations are not the only proposed mixing mechanisms operating in binary stars. As our results have shown that binary evolution calculations are sensitive to even a modest increase in mixing efficiency, it of interest to discuss other processes which may have an impact on internal mixing.
The interaction of internal circulation currents with the convective core can enhance mixing at the coreenvelope boundary in a way that mimics convective overshooting (Jermyn et al. 2018). This mechanism may explain large overshoot distances inferred from observations.
Vidal et al. (2018) suggest that certain magnetically induced instabilities can arise in the radiative zones of binary stars, causing mixing and the induction of significant magnetic fields. Such a mechanism may explain the observed magnetism of stars such as Vega. Furthermore, massive stars in close proximity to each other can be imagined to interact magnetically, which may further effect or result in mixing processes. Indeed magnetic fields have been observed in both Btype components of the binary system ϵ Lupi by Shultz et al. (2015), who conclude that the two fields are interacting.
Another agent which can mix material inside a star is pulsation. Massive stars can be pulsationally unstable, where pressure and gravity waves travel throughout the star. These waves are expected to cause some degree of mixing (Alecian 2014; Rogers & McElwaine 2017). What is more, when there exists a resonance between the orbital and pulsational periods of a binary star, pulsational instabilities can be strongly excited (Kumar et al. 1995; Witte & Savonije 2001). Such behaviour is observed in socalled heartbeat stars (Fuller 2017).
6.2. Implications for black hole merger events
Although the models presented here were computed at metallicity Z = Z_{⊙}/50 and the metallicity will have an effect on the final fates of binary systems (through the stellar radii and massloss rates), it is worth comparing our models with observed gravitational wave events. The double blackhole systems produced by our models are all expected to have merger timescales much shorter than the Hubble time (Marchant et al. 2016), and so are deemed as viable progenitors for gravitational wave events.
Shown in Fig. 8 are the component masses of the currently detected gravitational wave events with secondary masses greater than 25 M_{⊙} (Abbott et al. 2016; The LIGO Scientific Collaboration & The Virgo Collaboration 2019). Although our model grid is fairly coarse, a comparison between the models and observed blackhole merger events reveals that it would be very difficult for our Standard Grid models to produce double black holes with such low component masses or massratios as low as those observed. On the other hand, although our Enhanced Grid models do not fit the observations perfectly, taking into account metallicity effects and the coarseness of the model grid, we conclude that the Enhanced Grid models do not rule out the possibility that the four merger events could have been produced by the chemically homogeneous channel.
The detailed evolutionary models of Marchant et al. (2016) investigated blackhole merger events from the chemically homogeneous channel without modelling companioninduced circulation and predicted that below the pairinstability supernova gap, every double blackhole system has a massratio of unity, in agreement with our Standard Grid models. As demonstrated by our Enhanced Grid models, the inclusion of companioninduced circulation changes this picture significantly. We predict that around 20% of double blackhole systems will have avoided masstransfer, and estimate from Tables 2 and 3 that 12% of systems have final massratios in the range 0.80 < q < 0.95 and 8% with 0.60 < q < 0.80. Although our models predict the lowest achievable massratio to be approximately 0.8, at lower metallicities (du Buisson et al. 2020 find that systems with metallicities as low as 10^{−5} Z_{⊙} contribute to observable double blackhole merger events) weaker stellar winds will allow systems to remain much closer to their initial massratios, thus making the production of double blackhole systems with massratios around 0.7 entirely possible.
As our Enhanced Grid models produce 2.4 times the number of double helium stars than those of the Standard Grid (see Sect. 5.3), the inclusion of companioninduced circulation would increase the predicted rate of double blackhole merger events by a similar factor. Although we have considered only one metallicity, it may reasonably be assumed that the inclusion of companioninduced circulation will affect other metallicities similarly. du Buisson et al. (2020) calculate the aLIGO detection rate for systems below the pairinstability supernova gap to be 245 yr^{−1}, with our results suggesting that this rate should be revised upwards, by a factor of 2.4, to some 590 events per year, equating to nearly two detections per day. It has been predicted that gravitational wave events from the homogeneous channel could well be the dominant source for future detections (Marchant et al. 2016; Mandel & de Mink 2016; du Buisson et al. 2020). The results of this work serve to reinforce such statements.
7. Conclusions
We considered the mathematical results of Tassoul & Tassoul (1982a,b) for the steadystate circulation in the radiative zones of a single rotating star and a tidally locked binary star. By making a suitably justified approximation (namely that some components of the velocity field can be ignored), it has been shown that the rootmeansquared circulation velocity in the radial direction in a binary component may be expressed as a velocity enhancement to the Eddington–Sweet circulation. This enhancement is a simple function depending only on the masses of stars in the binary system. Such a result allows one to easily include the effects of companioninduced circulation in onedimensional stellar evolution codes.
We investigated the effects of companioninduced circulation on tidally locked binary systems with masses below the pairinstability supernova gap using the detailed stellar evolution code MESA. We find that at initial massratios in the range 0.7–1.0, binary systems with larger initial periods and slightly lower initial masses are able to form double helium stars through homogeneous evolution, compared to the standard mixing scenario. The circulation velocity enhancement also results in fewer systems undergoing masstransfer phases. To calculate the relative numbers of double helium stars produced from each grid, a Salpeter primary mass distribution, a distribution of the initial orbital period that is flat in logarithmic space, and a flat initial massratio distribution are used. Under these assumptions, we calculate that, in our parameter space models with enhanced mixing produce 2.4 times more double helium stars than models with standard mixing.
Our findings have strong implications for the formation of merging double blackhole systems. Previous detailed modelling has shown that the homogeneous evolution channel is only able to form double blackhole binaries with massratios very close to one (Marchant et al. 2016). The inclusion of companioninduced circulation into the models alters this, as we find double blackhole systems with massratios as low as 0.7 ∼ 0.8 are able to be produced from tidally locked binary stars. Our models show that systems with uneven massratios are expected to make up some 20% of the total population of merger events.
Furthermore, the previously calculated gravitational wave event detection rates from the homogeneous evolution channel (Marchant et al. 2016; Mandel & de Mink 2016; du Buisson et al. 2020) should be revised upwards in light of our findings by a factor of 2.4. This would lead to a predicted aLIGO detection rate of nearly 600 events per year, making the homogeneous evolution channel a certain contender for the dominant production source of double blackhole merger events.
MESA is based on the numerical treatment of Heger & Langer (2000) which, in turn uses the formalism of Kippenhahn (1974).
Acknowledgments
We thank our anonymous referee for useful comments on an earlier version of this manuscript. GK acknowledges support from CONACYT 252499 and UNAM/PAPIIT IN103619.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. X, 6, 041014 [Google Scholar]
 AbdulMasih, M., Sana, H., Sundqvist, J., et al. 2019, ApJ, 880, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Alecian, G. 2014, in Precision Asteroseismology, eds. J. A. Guzik, W. J. Chaplin, G. Handler, & A. Pigulski, IAU Symp., 301, 185 [Google Scholar]
 Almeida, L. A., Sana, H., de Mink, S. E., et al. 2015, ApJ, 812, 102 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
 Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chaboyer, B., & Zahn, J. P. 1992, A&A, 253, 173 [NASA ADS] [Google Scholar]
 Chandrasekhar, S. 1933, MNRAS, 93, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1957, An Introduction to the Study of Stellar Structure (USA: Dover Publication, Inc.) [Google Scholar]
 Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Detmers, R. G., Langer, N., Podsiadlowski, P., & Izzard, R. G. 2008, A&A, 484, 831 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, ArXiv eprints [arXiv:2002.11630] [Google Scholar]
 Eddington, A. S. 1929, MNRAS, 90, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Eggleton, P. P. 1983, ApJ, 268, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Fuller, J. 2017, MNRAS, 472, 1538 [NASA ADS] [CrossRef] [Google Scholar]
 Grevesse, N., Noels, A., & Sauval, A. J. 1996, in Cosmic Abundances, eds. S. S. Holt, & G. Sonneborn, ASP Conf. Ser., 99, 117 [Google Scholar]
 Hamann, W. R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151 [Google Scholar]
 Heger, A., & Langer, N. 2000, ApJ, 544, 1016 [NASA ADS] [CrossRef] [Google Scholar]
 Heger, A., Woosley, S. E., & Spruit, H. C. 2005, ApJ, 626, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752 [NASA ADS] [CrossRef] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Tout, C. A., & Chitre, S. M. 2018, MNRAS, 480, 5427 [NASA ADS] [CrossRef] [Google Scholar]
 Kähler, H. 2004, A&A, 414, 317 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kippenhahn, R. 1974, in Late Stages of Stellar Evolution, eds. R. J. Tayler, & J. E. Hesser, IAU Symp., 66, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution (New York: Springer) [Google Scholar]
 Koenigsberger, G., Morrell, N., Hillier, D. J., et al. 2014, AJ, 148, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Krogdahl, W. 1942, ApJ, 96, 124 [CrossRef] [Google Scholar]
 Kumar, P., Ao, C. O., & Quataert, E. J. 1995, ApJ, 449, 294 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207 [NASA ADS] [Google Scholar]
 Maeder, A. 1987, A&A, 178, 159 [NASA ADS] [Google Scholar]
 Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Martins, F., Mahy, L., & Hervé, A. 2017, A&A, 607, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moss, D. 1977, MNRAS, 178, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Peimbert, M., Luridiana, V., & Peimbert, A. 2007, ApJ, 666, 636 [NASA ADS] [CrossRef] [Google Scholar]
 Petrovic, J., Langer, N., Yoon, S. C., & Heger, A. 2005, A&A, 435, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Roxburgh, I. W. 1963, MNRAS, 126, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
 Schatzman, E. 1977, A&A, 56, 211 [NASA ADS] [Google Scholar]
 Shenar, T., Richardson, N. D., Sablowski, D. P., et al. 2017, A&A, 598, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shultz, M., Wade, G. A., Alecian, E., & BinaMIcS Collaboration 2015, MNRAS, 454, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. C., & Smith, D. H. 1981, MNRAS, 194, 583 [CrossRef] [Google Scholar]
 Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stȩpień, K. 2009, MNRAS, 397, 857 [NASA ADS] [CrossRef] [Google Scholar]
 Suijs, M. P. L., Langer, N., Poelarends, A.J., et al. 2008, A&A, 481, L87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sweet, P. A. 1950, MNRAS, 110, 548 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Tassoul, J. L., & Tassoul, M. 1982a, ApJ, 261, 265 [CrossRef] [Google Scholar]
 Tassoul, J. L., & Tassoul, M. 1982b, ApJS, 49, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Tassoul, J. L., & Tassoul, M. 1982c, ApJ, 261, 273 [CrossRef] [Google Scholar]
 The LIGO Scientific Collaboration & The Virgo Collaboration 2019, Phys. Rev. X, 9, 031040 [Google Scholar]
 Vidal, J., Cébron, D., Schaeffer, N., & Hollerbach, R. 2018, MNRAS, 475, 4579 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Witte, M. G., & Savonije, G. J. 2001, A&A, 366, 840 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoon, S. C., Langer, N., & Norman, C. 2006, A&A, 460, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J.P. 1975, A&A, 41, 329 [Google Scholar]
 Zahn, J.P. 1977, A&A, 57, 383 [NASA ADS] [Google Scholar]
 Zahn, J. P. 1983, in SaasFee Advanced Course 13: Astrophysical Processes in Upper Main Sequence Stars, eds. A. N. Cox, S. Vauclair, & J. P. Zahn, 253 [Google Scholar]
 Zahn, J. P. 2008, in EAS Publ. Ser., eds. M. J. Goupil, & J. P. Zahn, 29, 67 [CrossRef] [Google Scholar]
Appendix A: Calculation of the rms circulation velocity
Here we give details of the calculation of the rms circulation velocity for our extreme case considered in Sect. 3.2.2, where all U_{k} terms contribute equally to the circulation velocity. In the limit that U_{2} ≈ U_{3} ≈ U_{4}, Eq. (14) becomes
$$\begin{array}{c}\hfill {U}_{r}(r,\nu ,\mu )=\u03f5{U}_{\mathrm{S}}(r)[Q{\displaystyle \sum _{k=2}^{4}}{P}_{k}(\nu )+{P}_{2}(\mu )],\end{array}$$(A.1)
and we have
$$\begin{array}{cc}& {{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}{v}_{r}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi \hfill \\ \hfill & ={\u03f5}^{2}{U}_{\mathrm{S}}{(r)}^{2}{{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}[{Q}^{2}({\displaystyle \sum _{k=2}^{4}}{P}_{k}(cos(\phi )sin(\theta )){)}^{2}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}+{P}_{2}{(cos(\theta ))}^{2}+2Q{P}_{2}(cos(\theta )){\displaystyle \sum _{k=2}^{4}}{P}_{k}(cos(\phi )sin(\theta ))]\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi .\hfill \end{array}$$(A.2)
For convenience we shall let
$$\begin{array}{c}\hfill {{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}{v}_{r}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi ={\u03f5}^{2}{U}_{\mathrm{S}}{(r)}^{2}[A+B+C].\end{array}$$(A.3)
Where A, B, C are defined below:
$$\begin{array}{cc}\hfill A& ={Q}^{2}{{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}[{P}_{2}{(\nu )}^{2}+{P}_{3}{(\nu )}^{2}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}+{P}_{4}{(\nu )}^{2}+2{P}_{2}(\nu ){P}_{4}(\nu )]\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi ,\hfill \end{array}$$(A.4)
which has the analytic solution
$$\begin{array}{c}\hfill A={Q}^{2}{\pi}^{2}\frac{682039260864257}{1000000000000000}\approx 0.68{Q}^{2}{\pi}^{2}.\end{array}$$(A.5)
Now moving onto B,
$$\begin{array}{c}\hfill B={{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}{P}_{2}{(\mu )}^{2}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi ,\end{array}$$(A.6)
giving simply
$$\begin{array}{c}\hfill B=\frac{11{\pi}^{2}}{16}=0.6875{\pi}^{2}.\end{array}$$(A.7)
Lastly we tackle C,
$$\begin{array}{cc}\hfill C& =2Q{{\displaystyle \int}}_{\phi =0}^{\phi =2\pi}{{\displaystyle \int}}_{\theta =0}^{\theta =\pi}{P}_{2}(\mu ){P}_{2}(\nu )+{P}_{2}(\mu ){P}_{3}(\nu )\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}+{P}_{2}(\mu ){P}_{4}(\nu )\phantom{\rule{0.166667em}{0ex}}\mathrm{d}\theta \phantom{\rule{0.166667em}{0ex}}\mathrm{d}\phi ,\hfill \end{array}$$(A.8)
giving
$$\begin{array}{c}\hfill C=2Q{\pi}^{2}\frac{235}{1024}\approx Q{\pi}^{2}\times 0.4590.\end{array}$$(A.9)
Combining A, B, C gives
$$\begin{array}{c}\hfill {v}_{\phantom{\rule{0.333333em}{0ex}}\text{rms}}\approx \u03f5{u}_{\mathrm{S}}(r)\sqrt{0.3410{Q}^{2}+0.343750.2295Q}.\end{array}$$(A.10)
Appendix B: An investigation of the functions U_{k}
We will use simple models to investigate how the inviscid radial circulation velocity magnitudes arising from a binary companion, U_{2}, U_{3}, U_{4} are related to each other. Specifically we aim to see which of our extreme cases considered in Sect. 3.2, U_{2} ≫ U_{3} ≫ U_{4} or U_{2} ≈ U_{3} ≈ U_{4} is most realistic. To accomplish this we employ the theory of polytropes, for a detailed description of polytropes one should refer to Chandrasekhar (1957). A polytrope is an object which obeys the following relation between pressure and mass density
$$\begin{array}{c}\hfill P={K}_{\mathrm{P}}{\rho}^{\frac{n+1}{n}}\end{array}$$(B.1)
with K_{P} being a constant and n the polytropic index. It can be shown (Chandrasekhar 1957) that the polytropic constant, K_{P}, is given by
$$\begin{array}{c}\hfill {K}_{\mathrm{P}}=\frac{G}{n+1}[\frac{4\pi}{{\xi}_{1}^{n+1}}(\frac{\mathrm{d}\theta}{\mathrm{d}\xi}{}_{{\xi}_{1}}{)}^{1n}{]}^{\frac{1}{n}}{M}^{\frac{n1}{n}}{R}^{\frac{3n}{n}},\end{array}$$(B.2)
where $\theta ={(\rho /{\rho}_{c})}^{\frac{1}{n}}$ is the solution to the Lane–Emden equation and ξ is the radial coordinate divided by normalisation factor α, ξ = r/α. The radius of the star is represented by ξ_{1}, ξ_{1} = R/α and is found by setting θ = 0.
By combining the equation of hydrostatic equilibrium,
$$\begin{array}{c}\hfill {p}^{\prime}=\frac{Gm\rho}{{r}^{2}}\end{array}$$(B.3)
with the relation from Tassoul & Tassoul (1982b),
$$\begin{array}{c}\hfill {\rho}^{\prime}=\frac{n}{n+1}\frac{Gm{\rho}^{2}}{p{r}^{2}},\end{array}$$(B.4)
one can show that
$$\begin{array}{c}\hfill \frac{{\rho}^{\prime}}{{p}^{\prime}}=\frac{n}{n+1}\frac{\rho}{p}\xb7\end{array}$$(B.5)
Thus
$$\begin{array}{c}\hfill \frac{{\rho}^{\prime}\rho}{{p}^{\prime}}=\frac{n}{n+1}\frac{{\rho}^{2}}{p}\xb7\end{array}$$(B.6)
The fraction $\frac{{\rho}^{2}}{p}$ can be rewritten using the polytropic relation Eq. (B.1) so that
$$\begin{array}{c}\hfill \frac{{\rho}^{\prime}\rho}{{p}^{\prime}}=\frac{n}{n+1}{\rho}^{\frac{n1}{n}}{K}_{\mathrm{P}}^{1}.\end{array}$$(B.7)
Now we will make use of the fact that
$$\begin{array}{c}\hfill {\rho}_{c}=\frac{M}{4\pi {R}^{3}}{\xi}_{1}[\frac{\mathrm{d}\theta}{\mathrm{d}\xi}{}_{{\xi}_{1}}{]}^{1},\end{array}$$(B.8)
(see Chandrasekhar 1957) to rewrite Eq. (B.7) in terms of $\frac{\rho}{{\rho}_{c}}={\theta}^{n}$. This gives
$$\begin{array}{c}\hfill \frac{{\rho}^{\prime}\rho}{{p}^{\prime}}=\frac{n}{n+1}[\frac{M{\xi}_{1}}{4\pi {R}^{3}}{]}^{\frac{n1}{n}}[\frac{\mathrm{d}\theta}{\mathrm{d}\xi}{}_{{\xi}_{1}}{]}^{\frac{1n}{n}}{K}_{\mathrm{P}}^{1}{\theta}^{n1}.\end{array}$$(B.9)
Lastly we use the form of the polytropic constant, K_{P} given in Eq. (B.2) to produce
$$\begin{array}{c}\hfill \frac{{\rho}^{\prime}\rho}{{p}^{\prime}}=\frac{n}{4\pi G}{\xi}_{1}^{2}{R}^{2}{\theta}^{n1}.\end{array}$$(B.10)
To determine U_{k}(r), we first need to solve the differential equation Eq. (6), which when combined with Eq. (B.10) becomes
$$\begin{array}{c}\hfill {\varphi}_{k}^{\u2033}+\frac{2}{r}{\varphi}_{k}^{\prime}+{\varphi}_{k}[n{\zeta}_{1}^{2}{R}^{2}{\theta}^{n1}\frac{k(k+1)}{{r}^{2}}]=0.\end{array}$$(B.11)
We must now convert this from a differential equation in r to one in ξ by defining F_{k}(ξ) as an equivalent of ϕ_{k}(r), thus
$$\begin{array}{c}\hfill \frac{\mathrm{d}{\varphi}_{k}}{\mathrm{d}r}=\frac{\mathrm{d}{F}_{k}}{\mathrm{d}\xi}\frac{\mathrm{d}\xi}{\mathrm{d}r}=\frac{\mathrm{d}{F}_{k}}{\mathrm{d}\xi}\frac{1}{\alpha}\xb7\end{array}$$(B.12)
We may now recast Eq. (B.11) as
$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}{F}_{k}}{\mathrm{d}{\xi}^{2}}+\frac{2}{\xi}\frac{\mathrm{d}{F}_{k}}{\mathrm{d}\xi}+{F}_{k}[n{\theta}^{n1}\frac{k(k+1)}{{\xi}^{2}}]=0.\end{array}$$(B.13)
Happily this differential equation does not depend on the model’s mass or radius, only the polytropic index, as θ is a function of n.
One last step is to redetermine the constant c_{k} in terms of ξ and F_{k}(ξ) so that it becomes
$$\begin{array}{c}\hfill {c}_{k}={\omega}_{0}^{2}\frac{{M}_{\mathrm{comp}}}{M+{M}_{\mathrm{comp}}}(\frac{R}{d}{)}^{(k2)}\frac{(2k+1){({\xi}_{1}\alpha )}^{2}}{(k+1){F}_{k}({\xi}_{1})+{\xi}_{1}{F}_{k}^{\prime}({\xi}_{1})}\xb7\end{array}$$(B.14)
Equation (B.13) must be solved numerically, not least because in general θ has no analytic solution. We shall choose to solve for an n = 3 polytrope, giving ξ_{1} = 6.8968 and $\frac{\mathrm{d}\theta}{\mathrm{d}\xi}=0.0424$. This polytrope describes the Eddington Standard Model, where energy transport is partially radiative, so is the best choice for the envelope of a massive star. In terms of the Emden variables the functions U_{k} defined in Eq. (18) become
$$\begin{array}{c}\hfill {U}_{k}(r)=\frac{2L{\xi}^{4}{\alpha}^{4}}{{G}^{2}{m}^{3}}\frac{n+1}{n1.5}{c}_{k}[\frac{\mathrm{d}{F}_{k}}{\mathrm{d}\xi}\frac{1}{\alpha}+(\frac{2}{\xi \alpha}\frac{3}{\xi \alpha}){F}_{k}].\end{array}$$(B.15)
The mass, m, is expressed as an integral of the density
$$\begin{array}{c}\hfill m(\xi )=4\pi {\alpha}^{3}{\rho}_{c}{{\displaystyle \int}}_{0}^{\xi}{\xi}^{\prime 2}{\theta}^{n}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}{\xi}^{\prime},\end{array}$$(B.16)
which when combined with Eq. (B.8) leads to
$$\begin{array}{c}\hfill m(\xi )={\alpha}^{3}{\xi}_{1}\frac{M}{{R}^{3}}[{\xi}_{1}^{2}\frac{\mathrm{d}\theta}{\mathrm{d}\xi}{}_{{\xi}_{1}}{]}^{1}{{\displaystyle \int}}_{0}^{\xi}{\xi}^{\prime 2}{\theta}^{n}\phantom{\rule{0.166667em}{0ex}}\mathrm{d}{\xi}^{\prime}.\end{array}$$(B.17)
Equations (B.15) and (B.17) can be solved numerically to provide the functions U_{k}(r). These functions are plotted in Fig. B.1 for an n = 3 polytrope with M = 3 M_{⊙}, R = 1.75 R_{⊙}, L = 93 L_{⊙} and a companion mass of M_{comp} = 3 M_{⊙}. The luminosity is assumed to be concentrated in a central point. From Eq. (B.14) we see that for small values of R/d, the constants c_{k} decrease sharply with increasing k. This means that the functions U_{k}(r) will be most similar to each other for high R/d values, therefore we have chosen the ratio of stellar radius to orbital separation to be 0.35, which represents a system that is very close to contact (for massratios of 1 a contact system has R/d = 0.379 Eggleton 1983).
Fig. B.1. Inviscid radial circulation velocity magnitudes, U_{k}(r) plotted as a function of normalised radius for an n = 3 polytrope with M = 3 M_{⊙}, R = 1.75 R_{⊙}, L = 93 L_{⊙} and a companion mass of M = 3 M_{⊙}. The ratio of stellar radius to orbital separation is set to be 0.35, representing a nearcontact system. 
A more interesting quantity to consider is the ratios of U_{3} and U_{4} with U_{2}. We note that such ratios will not depend on the specifics of our chosen star (mass, luminosity and radius) because the solutions to Eq. (B.13) only depend on the k value and polytropic index. The only other parameter not divided out is the ratio of stellar radius to orbital separation. Figure B.2 shows the ratios of the U_{k} functions for varying R/d values. It is seen that at most U_{3} is 70% of U_{2} and U_{4} is 40% of U_{2}. For smaller values of R/d, corresponding to detached systems, U_{2} becomes even greater compared to U_{3} and U_{4}. Therefore it is concluded that the approximation U_{2} ≫ U_{3} ≫ U_{4} is more suitable than U_{2} ≈ U_{3} ≈ U_{4}.
Fig. B.2. Functions U_{3} (left panel) and U_{4} (right panel) as fractions of U_{2} plotted as a function of normalised radius for the polytropic model in Fig. B.1. Various values of the stellar radius to orbital separation, R/d are plotted as coloured lines as given by the legend. R/d = 0.38 represents a contact system and is hence an upper limit. 
All Tables
Summary of the functions defining the circulation in the radial direction, u_{r}, for the single and binary star cases.
Fractions of binary systems in our model grids that are predicted to form double helium stars for each initial massratio.
Of systems that form a double helium star, the fractions of those that do not undergo masstransfer for each initial mass ratio.
All Figures
Fig. 1. Circulation pattern of Eddington–Sweet circulation in a single star Cowling point source model with electron scattering opacity, M = 3 M_{⊙}, R = 1.75 R_{⊙}L = 93 L_{⊙} and ϵ = 0.1. Only the outer circulation cells in the northern hemisphere are shown. The green arrows give the direction of the motion. The central convective region is marked in red, the pole is marked with a blue cross and the radiative envelope is shown in opaque grey. This figure is a reproduction of Fig. 3 of Tassoul & Tassoul (1982b). 

In the text 
Fig. 2. Circulation pattern of a rotationally and tidally distorted Cowling point source model with electron scattering opacity, M = 3 M_{⊙}, M_{comp} = 3 M_{⊙}, R = 1.75 R_{⊙}L = 93 L_{⊙} and ϵ = 0.1. Three different perspectives are shown. The k = 3, 4 contributions to the velocity field due to the companion have been omitted. The companion is shown as a yellow sphere at an arbitrary orbital separation. The green arrows give the direction of the motion. The central convective region is marked in red, the pole is marked by a blue cross and the radiative envelope is shown in opaque grey. 

In the text 
Fig. 3. Factor by which averaged Eddington–Sweet circulation velocity is enhanced in a tidally locked binary component as a function of the effective mass ratio, q_{*} for two approximations as given in the legend. 

In the text 
Fig. 4. Diagnostic plots for a system with primary mass, M_{1} = 100 M_{⊙}, initial massratio, q_{i} = 0.7, and initial orbital period 1.5 days. Top panels: system in the Standard Grid (without enhanced mixing) and bottom panels: system in the Enhanced Grid (with enhanced mixing). Left panels: evolution of each component’s mass and the mass ratio, q with a dotted line. Central panels: ratios of nuclear and mixing timescales as a function of time. Right panels: difference between central and surface helium mass fraction. Star 1 is plotted in red, Star 2 in blue. The black cross and diamond represent the times for which profiles of Star 2 are plotted in Fig. 5. 

In the text 
Fig. 5. Internal profiles of Star 2 in the example system (M_{1} = 100 M_{⊙}, q_{i} = 0.7, P_{i} = 1.5 days) showing the total diffusion coefficient, D_{tot} (left panels), and helium mass fraction, Y (right panels), as a function of Lagrangianmasscoordinate. The star in the Standard Grid is marked in green and that in the Enhanced Grid in black. Top panels: profiles at an age of 0.27 Myr, as marked by a cross in Fig. 4, bottom panels: age of 0.80 Myr, as marked by a diamond in Fig. 4. The grey dashed line in the right panels represents the initial helium mass fraction, Y_{0}. 

In the text 
Fig. 6. Minimum initial fraction of critical angular velocity required for homogeneous evolution to occur as a function of primary mass for systems with q_{i} = 1.0. Models from the Standard Grid are marked in black and those from the Enhanced grid in green. 

In the text 
Fig. 7. Final configurations of models in the Standard Grid (top row) and Enhanced Grid (bottom row) for initial mass ratios q_{i} = 1.0, 0.9, 0.7. The colours indicate where overflow at the second Lagrangian point occurred at the zeroagemainsequence or during evolution (purple and green respectively), where the difference between central and surface helium mass fraction exceeds 0.2 (red) and where homogeneous evolution occurs resulting in a double helium star (blue). Single hatched cells indicate that the system underwent contact on the mainsequence, and crosshatched that contact occurred on the zeroagemainsequence. 

In the text 
Fig. 8. Masses of Star 1 (xaxis) and Star 2 (yaxis) at core helium depletion for systems with q_{i} = 0.9, 0.7 (plotted as triangles and squares respectively). Models in the Standard Grid (without enhanced mixing) are marked in black, those in the Enhanced Grid (with enhanced mixing) in green. For clarity systems with q_{i} = 1 are not plotted as they remain at their initial mass ratio. The grey dotted lines show various massratios as marked. Red, blue, yellow and purple stars represent measured progenitor masses of gravitational wave events GW170729, GW170823, GW150914 and GW170818 respectively (Abbott et al. 2016; The LIGO Scientific Collaboration & The Virgo Collaboration 2019). 

In the text 
Fig. B.1. Inviscid radial circulation velocity magnitudes, U_{k}(r) plotted as a function of normalised radius for an n = 3 polytrope with M = 3 M_{⊙}, R = 1.75 R_{⊙}, L = 93 L_{⊙} and a companion mass of M = 3 M_{⊙}. The ratio of stellar radius to orbital separation is set to be 0.35, representing a nearcontact system. 

In the text 
Fig. B.2. Functions U_{3} (left panel) and U_{4} (right panel) as fractions of U_{2} plotted as a function of normalised radius for the polytropic model in Fig. B.1. Various values of the stellar radius to orbital separation, R/d are plotted as coloured lines as given by the legend. R/d = 0.38 represents a contact system and is hence an upper limit. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.