© ESO 2020

1. Introduction

The existence of grand-scale 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 non-rotating 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 companion-induced 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 blue-wards 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 binary-specific 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 (so-called 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 black-hole 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; Abdul-Masih 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 O-type 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 steady-state circulation currents in a tidally locked component of a binary system and numerically investigate the effects of this circulation in massive short-period binary systems using a one-dimensional 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 one-dimensional 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 3-dimensional 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:

1. 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 tidally-locked.

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

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

4. The system is detached.

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

6. The stars are not significantly distorted away from sphericity.

7. The stars rotate as solid-bodies.

8. The stars are on the zero-age-main-sequence such that there are no buoyancy forces arising from chemical gradients.

9. The stars are non-magnetic.

10. The circulation propagates in the radiative zone only.

11. The radiative zone is inviscid.

We now discuss the justifications, if they exist, for each assumption. Many short-period 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 spin-orbit 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 short-period binary systems are expected to be circularised and synchronised, and assumption 1 is valid in many systems.

For very luminous or extreme mass-ratio 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 core-envelope interface, as this is where fresh fuel is mixed into the core. In most cases the core-envelope interface would not be greatly affected by contact, thereby contact and non-contact 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 Tayler-Spruit dynamo (Spruit 2002), which causes angular momentum to be transported throughout the stellar interior via magnetic interactions such that near solid-body rotation is achieved. The Tayler-Spruit 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 zero-age-main-sequence with a uniform chemical composition, meaning that no buoyancy forces act upon the circulation. Although during the star-formation process, material may become partially segregated, no process, such as nuclear burning, occurs in stars before the zero-age-main-sequence 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⁻¹ (see Fig. B.1), while convective velocities are typically of the order 10⁴ cm s⁻¹ (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 steady-state 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 ϵ=\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 ϵ=\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₂ 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{5R^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}·\hfill \end{array}$$(5)

The function ϕ₂(r) satisfies the differential equation

$$\begin{array}{c}\hfill {\varphi }_k^{″}+\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 co-ordinate ($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 )& =ϵU_{\mathrm{S}}(r)P_2(\mu ),\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill u_{\mu }(r,\mu )& =ϵV(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₂ 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}{n-1.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 co-ordinate, r the Eulerian radial co-ordinate 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 3-dimensional 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 x-axis points towards the centre of gravity of the companion and the z-axis points towards the pole of the star in question. In these coordinates the plane described by φ = 0 is the sub-binary 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 point-mass 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 3-dimensional 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 )& =ϵ[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 )& =ϵ[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 )& =ϵ\left[{\displaystyle \sum _{k=2}^4}\frac{r{V}_k(r)}{\sqrt{(1-{\mu }^2)}}\frac{\partial P_k(\nu )}{\partial \phi }\right]·\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}{n-1.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)}^{(k-2)}\frac{(2k+1)R^2}{(k+1){\varphi }_k(R)+R{{\varphi }^{\prime }}_k(R)}·\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 tidally-locked binary component. The dependence on the orbital separation is expressed in Eq. (20) through the ${\left(\frac{R}{d}\right)}^{(k-2)}$ term and also in Eqs. (14)–(16) through the ϵ parameter, describing the strength of rotation. In a tidally-locked 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 3-dimensional 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 sub-binary 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 one-dimensional 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⊙, Mcomp = 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 main-sequence evolution, the departures from spherical symmetry may be non-negligible. 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 flow-patterns in Figs. 1 and 2. One of the most likely is the baroclinic instability which can set up wave-like 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.

Even though the circulation patterns inside a star are not one-dimensional, 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 one-dimensional problem (Chaboyer & Zahn 1992). To analyse the circulation patterns in one dimension, the root-mean-squared (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_{\text{rms}}=\frac{\sqrt{{\int }_{\phi =0}^{\phi =2\pi }{\int }_{\theta =0}^{\theta =\pi }{u}_r^2\mathrm{d}\theta \mathrm{d}\phi }}{\sqrt{{\int }_{\phi =0}^{\phi =2\pi }{\int }_{\theta =0}^{\theta =\pi }\mathrm{d}\theta \mathrm{d}\phi }}·\end{array}$$(21)

3.1. A rotating single star

Using Eqs. (8) and (21) to calculate v_rms gives

$$\begin{array}{c}\hfill v_{\text{rms}}=\sqrt{{ϵ}^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\mathrm{d}\theta \mathrm{d}\phi },\end{array}$$(22)

which results in

$$\begin{array}{c}\hfill v_{\text{rms}}=ϵ|U_{\mathrm{S}}(r)|\sqrt{\frac{11}{32}}\approx 0.59ϵ|u_{\mathrm{S}}(r)|·\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₂ 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₂ and U_S differ only by the constants c₂ and A₂ (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₂ dominates the companion-induced velocity field such that U₂ ≫ U₃ ≫ U₄ and the case where all terms contribute equally such that U₂ ≈ U₃ ≈ U₄.

3.2.1. Extreme case (i): U₂ ≫ U₃ ≫ U₄

When we neglect U₃ and U₄ and apply the relation expressed in Eqs. (26), (14) becomes

$$\begin{array}{c}\hfill U_r(r,\nu ,\mu )=ϵ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_{\text{rms}}& =ϵ|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 & +P_2{(cos(\theta ))}^2+2Q{P}_2(cos(\phi )sin(\theta ))P_2(cos(\theta )))\mathrm{d}\theta \mathrm{d}\phi {]}^{\frac{1}{2}},\hfill \end{array}$$(28)

which has the analytic solution

$$\begin{array}{c}\hfill v_{\text{rms}}=ϵ|U_{\mathrm{S}}(r)|[\frac{49}{256}{Q}^2+\frac{11}{32}-Q\frac{11}{32}{]}^{\frac{1}{2}}·\end{array}$$(29)

3.2.2. Extreme case (ii): U₂ ≈ U₃ ≈ U₄

In the limit that U₂ ≈ U₃ ≈ U₄, Eq. (14) becomes

$$\begin{array}{c}\hfill U_r(r,\nu ,\mu )=ϵ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_{\text{rms}}\approx ϵ|u_{\mathrm{S}}(r)|\sqrt{0.3410Q^2+0.34375-0.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₂ ≫ U₃ ≫ U₄ is the most realistic (see Appendix B). Such a result could also be inferred by studying the circulation patterns arising from U₂, U₃ and U₄ 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_{\text{ES}}=\sqrt{\frac{441}{88}(\frac{q_{\ast }}{q_{\ast }+1}{)}^2+1+\frac{3q_{\ast }}{q_{\ast }+1}}·\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 companion-induced 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 companion-induced mixing as derived in Sect. 3, version r10398 of the one-dimensional 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öhm-Vitense 1958) with mixing length parameter α = 1.5 with the Ledoux criterion. We implement semi-convection 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 mass-loss 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 mass-loss 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 mass-loss is handled by the method of Heger & Langer (2000) whereby the relationship between mass-loss 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 }}_{\text{crit}}(\mathrm{\Gamma })}}{]}^{0.43},\end{array}$$(33)

with the critical rotation rate being

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\text{crit}}(\mathrm{\Gamma })=\sqrt{\frac{\mathit{GM}}{R^3}(1-\mathrm{\Gamma })}\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, mass-loss rates become divergent. To avoid this issue, once Ω/Ω_crit reaches 0.98, we switch to an implicit mass-loss rate, which is calculated such that the fraction Ω/Ω_crit remains 0.98.

Internal angular momentum transport is accomplished via the Tayler-Spruit dynamo (Spruit 2002), which enforces near solid-body 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 Roche-lobe overflow come from considering the Roche-lobe volumes and their volume-equivalent radii as per the prescription of Eggleton (1983), while mass-transfer rates are determined implicitly by demanding that the mass-donor stays inside its Roche-lobe.

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 mean-molecular 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 companion-induced 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{3q_{\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 companion-induced circulation we set f_ES = 1.

We compute grids of binary models with primary mass, M₁, ranging from log(M₁/M_⊙) = 1.4–2.0 in intervals of 0.1 dex, such that the models are not massive enough to undergo the pair-instability supernova phenomenon. We choose initial orbital periods ranging from 0.5 to 2.0 days in intervals of 0.1 day. Four different initial mass-ratios are considered, q_i = M₂/M₁ = 1.0, 0.9, 0.7, 0.5 such that the mass of the secondary star, M₂ is the product of the primary mass and the mass-ratio. 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, mass-loss through stellar winds causes a widening of the binary orbit which may result in double black-hole 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 companion-induced 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 mass-loss, which is particularly important in a binary system as the binary orbit is influenced by both mass and angular-momentum 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 mass-loss 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 mass-loss in very massive stars.

Furthermore the 3-dimensional 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 co-ordinates, and therefore so will the mass-loss rates. This may lead to a significant difference between effective mass-loss 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 mass-ratio 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 }_{\text{mix}}=\frac{4}{{\pi }^2}[{{\displaystyle \int }}_0^R\frac{\mathrm{d}r}{\sqrt{D_{\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 companion-induced circulation model we shall examine in detail a system with a primary mass, M₁ = 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, M1 = 100 M⊙, initial mass-ratio, qi = 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 core-envelope 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 (M1 = 100 M⊙, qi = 0.7, Pi = 1.5 days) showing the total diffusion coefficient, Dtot (left panels), and helium mass fraction, Y (right panels), as a function of Lagrangian-mass-coordinate. 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, Y0.

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 core-envelope 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 core-envelope 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 Roche-lobe. Despite the lack of mass-transfer, the mass-ratio of the system evolves towards unity as the more massive star experiences higher mass-loss rates, as shown in Fig. 4d, giving a final mass-ratio somewhat closer to unity than the initial mass-ratio. However, owing to the fact that mass-loss from stellar winds is strongly metallicity dependant, very low metallicity systems will remain closer to their initial mass-ratios.

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 equal-mass systems. This plot shows stars with masses from log(M₁/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 qi = 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¹, 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 mass-transfer 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 tidally-locked binary without suffering L2 overflow at the zero-age-main-sequence, 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 tidally-locked 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 mass-transfer 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₁₀(M₁/M_⊙) = 2.0, we have not found the upper initial orbital period limit for homogeneous evolution to occur, at this initial mass-ratio, 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 qi = 1.0, 0.9, 0.7. The colours indicate where overflow at the second Lagrangian point occurred at the zero-age-main-sequence 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 main-sequence, and cross-hatched that contact occurred on the zero-age-main-sequence.

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 mass-ratio q_i = 1, 1.9 days for q_i = 0.9 and 1.5 days for q_i = 0.7. At every initial mass-ratio, for initial primary masses greater than log₁₀(M₁/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₁₀(M₁/M_⊙) = 1.5 are able to form double helium stars, while in the Standard Grid this limit is log₁₀(M₁/M_⊙) = 1.7. Similarly for systems with q_i = 0.7, the minimum initial primary mass that evolves homogeneously is shifted from log₁₀(M₁/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₁₀(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)}\mathrm{d}{log}_{10}(P){\int }_{M_i-\mathrm{\Delta }M}^{M_i+\mathrm{\Delta }M}{m}^{-2.35}\mathrm{d}m}{{\int }_{{log}_{10}(0.5)}^{{log}_{10}(2.0)}\mathrm{d}{log}_{10}(P){\int }_{{10}^{1.4}}^{{10}^{2.0}}{m}^{-2.35}\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 mass-ratio 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.

5.4. Occurrence of mass-transfer

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 Roche-lobe. We note that although homogeneously evolving stars do increase their radii slightly during hydrogen burning, such that stars that are detached at the zero-age-main-sequence 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 mass-ratio 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 mass-ratio becomes more extreme, such systems are remarkable because they produce double helium stars with mass-ratios close to their initial mass-ratio. In the Enhanced Grid, double helium stars that avoid interaction are produced at initial mass-ratios as extreme as q_i = 0.7, whereas in the Standard Grid no non-interacting 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 mass-transfer. We see that only 5% of double helium stars in the Standard Grid avoid mass-transfer 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 mass-ratio of unity, whereas for the Enhanced Grid, some 20% of double helium stars remain close to their initial mass-ratio.

To explore the final mass-ratios 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 mass-ratios 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 mass-ratios as low as 0.8 to be produced. Without enhanced mixing the final systems all have mass-ratios 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 (x-axis) and Star 2 (y-axis) at core helium depletion for systems with qi = 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 qi = 1 are not plotted as they remain at their initial mass ratio. The grey dotted lines show various mass-ratios 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 grand-scale 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 core-envelope 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 B-type 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 so-called 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 mass-loss rates), it is worth comparing our models with observed gravitational wave events. The double black-hole 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 black-hole 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 mass-ratios 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 black-hole merger events from the chemically homogeneous channel without modelling companion-induced circulation and predicted that below the pair-instability supernova gap, every double black-hole system has a mass-ratio of unity, in agreement with our Standard Grid models. As demonstrated by our Enhanced Grid models, the inclusion of companion-induced circulation changes this picture significantly. We predict that around 20% of double black-hole systems will have avoided mass-transfer, and estimate from Tables 2 and 3 that 12% of systems have final mass-ratios in the range 0.80 <  q <  0.95 and 8% with 0.60 <  q <  0.80. Although our models predict the lowest achievable mass-ratio to be approximately 0.8, at lower metallicities (du Buisson et al. 2020 find that systems with metallicities as low as 10⁻⁵ Z_⊙ contribute to observable double black-hole merger events) weaker stellar winds will allow systems to remain much closer to their initial mass-ratios, thus making the production of double black-hole systems with mass-ratios 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 companion-induced circulation would increase the predicted rate of double black-hole merger events by a similar factor. Although we have considered only one metallicity, it may reasonably be assumed that the inclusion of companion-induced circulation will affect other metallicities similarly. du Buisson et al. (2020) calculate the aLIGO detection rate for systems below the pair-instability supernova gap to be 245 yr⁻¹, 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 steady-state 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 root-mean-squared 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 companion-induced circulation in one-dimensional stellar evolution codes.

We investigated the effects of companion-induced circulation on tidally locked binary systems with masses below the pair-instability supernova gap using the detailed stellar evolution code MESA. We find that at initial mass-ratios 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 mass-transfer 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 mass-ratio 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 black-hole systems. Previous detailed modelling has shown that the homogeneous evolution channel is only able to form double black-hole binaries with mass-ratios very close to one (Marchant et al. 2016). The inclusion of companion-induced circulation into the models alters this, as we find double black-hole systems with mass-ratios as low as 0.7 ∼ 0.8 are able to be produced from tidally locked binary stars. Our models show that systems with uneven mass-ratios 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 black-hole merger events.

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

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₂ ≈ U₃ ≈ U₄, Eq. (14) becomes

$$\begin{array}{c}\hfill U_r(r,\nu ,\mu )=ϵ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\mathrm{d}\theta \mathrm{d}\phi \hfill \\ \hfill & ={ϵ}^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 & +P_2{(cos(\theta ))}^2+2Q{P}_2(cos(\theta )){\displaystyle \sum _{k=2}^4}{P}_k(cos(\phi )sin(\theta ))]\mathrm{d}\theta \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\mathrm{d}\theta \mathrm{d}\phi ={ϵ}^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 & +P_4{(\nu )}^2+2P_2(\nu )P_4(\nu )]\mathrm{d}\theta \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.68Q^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\mathrm{d}\theta \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 & +P_2(\mu )P_4(\nu )\mathrm{d}\theta \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_{\text{rms}}\approx ϵ|u_{\mathrm{S}}(r)|\sqrt{0.3410Q^2+0.34375-0.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₂, U₃, U₄ are related to each other. Specifically we aim to see which of our extreme cases considered in Sect. 3.2, U₂ ≫ U₃ ≫ U₄ or U₂ ≈ U₃ ≈ U₄ 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}{)}^{1-n}{]}^{\frac{1}{n}}{M}^{\frac{n-1}{n}}{R}^{\frac{3-n}{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 co-ordinate divided by normalisation factor α, ξ = r/α. The radius of the star is represented by ξ₁, ξ₁ = 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}·\end{array}$$(B.5)

Thus

$$\begin{array}{c}\hfill \frac{{\rho }^{\prime }\rho }{p^{\prime }}=\frac{n}{n+1}\frac{{\rho }^2}{p}·\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{n-1}{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{n-1}{n}}[-\frac{\mathrm{d}\theta }{\mathrm{d}\xi }{|}_{{\xi }_1}{]}^{\frac{1-n}{n}}{K}_{\mathrm{P}}^{-1}{\theta }^{n-1}.\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 }^{n-1}.\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^{″}+\frac{2}{r}{\varphi }_k^{\prime }+{\varphi }_k[n{\zeta }_1^2{R}^{-2}{\theta }^{n-1}-\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 }·\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 }^{n-1}-\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}{)}^{(k-2)}\frac{(2k+1){({\xi }_1\alpha )}^2}{(k+1)F_k({\xi }_1)+{\xi }_1{F}_k^{\prime }({\xi }_1)}·\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 ξ₁ = 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}{n-1.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\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\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 mass-ratios of 1 a contact system has R/d = 0.379 Eggleton 1983).

Fig. B.1. Inviscid radial circulation velocity magnitudes, Uk(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 near-contact system.

A more interesting quantity to consider is the ratios of U₃ and U₄ with U₂. 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₃ is 70% of U₂ and U₄ is 40% of U₂. For smaller values of R/d, corresponding to detached systems, U₂ becomes even greater compared to U₃ and U₄. Therefore it is concluded that the approximation U₂ ≫ U₃ ≫ U₄ is more suitable than U₂ ≈ U₃ ≈ U₄.

Fig. B.2. Functions U3 (left panel) and U4 (right panel) as fractions of U2 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

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⊙, Mcomp = 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, M1 = 100 M⊙, initial mass-ratio, qi = 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 (M1 = 100 M⊙, qi = 0.7, Pi = 1.5 days) showing the total diffusion coefficient, Dtot (left panels), and helium mass fraction, Y (right panels), as a function of Lagrangian-mass-coordinate. 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, Y0.

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 qi = 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 qi = 1.0, 0.9, 0.7. The colours indicate where overflow at the second Lagrangian point occurred at the zero-age-main-sequence 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 main-sequence, and cross-hatched that contact occurred on the zero-age-main-sequence.

In the text

Fig. 8. Masses of Star 1 (x-axis) and Star 2 (y-axis) at core helium depletion for systems with qi = 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 qi = 1 are not plotted as they remain at their initial mass ratio. The grey dotted lines show various mass-ratios 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, Uk(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 near-contact system.
In the text

	Fig. B.2. Functions U3 (left panel) and U4 (right panel) as fractions of U2 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