A calculator for Roche lobe properties
Denis A Leahy1 &
Janet C Leahy2
The Roche lobe geometry is important to understand and study the properties of the mass-losing component in a semi-detached binary system. However it is not easy to calculate accurately, and existing tables usually do not include the parameters of the binary system under study, nor do they allow for non-synchronous rotation.
A calculator for properties of the Roche lobe is presented in two formats. An easy-to-use Java version has a graphic interface, and a Fortran 90 version has a command line interface. The Fortran version allows for easy modifications by the user. Both versions have two basic output options: one provides values of a set of various quantities (such as the Lagrange points along the binary axis, and area and volume of the Roche lobe); the second provides \(R(\theta,\phi)\), the distance from the stellar center to the stellar surface for any specified polar angle. A single set of input parameters can be entered directly or a large set of input parameters can be specified in a text file. The calculator includes the options to have non-synchronous rotation of the star, or to have the star underfill its Roche lobe. It can be used to calculate Roche lobe properties for the case of elliptical orbits, with some restrictions.
We present a convenient software tool for quickly and accurately calculating Roche lobe properties for mass ratio in the range 0.01 to 100, for Roche lobe fill-out factor in the range 0.1 to 1.0, and for dimensionless rotation rate of the star in the range 0.1 to 2.0. This will allow anyone working with a binary star system to obtain the Roche lobe or stellar surface geometry for their system.
The Roche potential is the potential energy per unit test mass which is orbiting the center-of-mass of a binary star system at the same rate as the two stars. A circular orbit is assumed for the binary system. The Roche potential includes both gravitational and centripetal energy, such that its derivative gives the force (and thus acceleration) on the test mass.
The Roche potential has many uses, including finding the locations in the binary system where the acceleration is zero- the Lagrange points L1, L2, L3, L4 and L5. In the case of a binary star system where one of the stars (star 1) increases its size, star 1 will reach a limiting surface. This surface is called the Roche lobe, which is sketched in Figure 1. The Roche lobe is also defined by the surface that has the potential equal to the L1 potential. The Roche lobe has the property that any material outside this surface is unbound from star 1 and will be lost.
Roche lobe geometry. Illustration of the Roche lobe for star 1, with x, y, z axes and the various Roche radii (\(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\)).
Many works have studied and calculated properties of the Roche potential and the Roche lobe, starting with the pioneering work of Kopal (1959). Plavec and Kratochvil (1964) present tables of distances for several points on the equipotential surfaces which pass through the L1 and L2 Lagrange points. This was done for a range of mass ratios from 0.1 to 1.0. Eggleton (1983) calculated volumes, V, of the Roche lobe for a range of mass ratios to obtain the equivalent volume radius \(R_{eq}\) (defined by \(V=4 \pi R_{eq}^{3}/3\)). He then presented a simple fitting formula for \(R_{eq}\) as a function of q, accurate to better than 1%, which is often used and now known as the Eggleton formula. Mochnacki (1984) presented tables of potentials and positions of the L1 and L2 Lagrange points as a function of mass ratio. He defined a fill-out factor F in terms of the potential of an equipotential surface relative to the potential of the surfaces through the L1 and L2 points. F was defined so that \(F<1\) corresponds to surfaces inside the Roche lobe, and \(1< F<2\) corresponds to surfaces between the Roche lobe and the equipotential through the L2 point. He gave tables of volume radius, area, average gravity and average inverse gravity as a function of mass ratio and F. Pathania and Medupe (2012) present series expansions for the radius of Roche equipotential surfaces, and compare various orders of the expansion to numerically calculated values.
The purpose of the present work is to present an freely available software tool, written in two versions (Java and Fortran), which calculates radii of the Roche lobe. It does this for any specified direction, and gives some other commonly used quantities (such as the Lagrange points and values of the potential). The calculator is designed to be accurate for any mass ratio q between 0.01 and 100, for dimensionless rotation rates between 0.01 and 2, and fill-out factors between 0.1 and 1.0. The calculator may work for parameters outside these limits but has not been tested for accuracy or errors outside these limits.
The coordinates are Cartesian (x, y, z) and spherical-polar (r, θ, ϕ), centred on star 1 (with mass \(M_{1}\)), with the x-axis (\(\theta=\pi/2\), \(\phi=0\)) pointing towards star 2 (with mass \(M_{2}\)). The z-axis is perpendicular to the orbital plane. The mass ratio is defined as \(q=M_{2}/M_{1}\) and the binary separation is a.
We use the dimensionless form of the Roche potential (or potential energy per unit mass) obtained by dividing by \(G M_{1}/a\), where G is Newton's gravitational constant, and using distances in units of the binary separation a. The equation for the resulting dimensionless potential \(\Omega(r,\theta,\phi)\) in the case of synchronous rotation is:
$$ \begin{aligned}[b] \Omega(r,\theta,\phi)={}& \frac{1}{r} + q \biggl( \frac{1}{\sqrt{1-2r\sin{\theta}\cos{\phi}+r^{2}}} \\ &{}-r\sin{\theta}\cos{\phi}\biggr) +\frac{q+1}{2}r^{2}\sin{ \theta}^{2}. \end{aligned} $$
This is the same form used by Kopal (1959) and Pathania and Medupe (2012). When non-synchronous rotation is included Limber (1963), with \(p=\Omega_{star}/\Omega_{binary}\) is the rate of stellar rotation divided by the binary rotation rate, the potential becomes:
$$ \begin{aligned}[b] \Omega(r,\theta,\phi)={}& \frac{1}{r} + q \biggl( \frac{1}{\sqrt{1-2r\sin{\theta}\cos{\phi}+r^{2}}} \\ &{}-r\sin{\theta}\cos{\phi}\biggr) +\frac{q+1}{2} p^{2} r^{2}\sin{\theta}^{2}. \end{aligned} $$
This potential can be further generalized to the case of elliptical orbits, as discussed in Sepinsky et al. (2007). In the quasi-static approximation, valid for dynamical timescale of the star much less than the tidal timescale, the instantaneous shape of the star can be approximated by the instantaneous surface of constant potential, even though the instantaneous surface is changing with orbital phase. As shown by Sepinsky et al. (2007), the instantaneous potential for an elliptical orbit can be written in the same form as equation (2) above if one replaces \(p^{2}\) by the function
$$ A(p,e,\nu)=\frac{p^{2} (1+e)^{4}}{(1+e \cos(\nu))^{3}} $$
with e is the eccentricity and ν is the true anomaly for the position of the star in its elliptical orbit.
For a physical system, one knows the actual values of the masses \(M_{1}\) and \(M_{2}\). Either one specifies the binary separation a, or calculates it from the orbital period P, using Kepler's third law, \(4 \pi^{2} a^{3}= G (M_{1}+M_{2}) P^{2}\). Dimensionless distances are converted into physical units by multiplying by a.
The Lagrange points are the locations where the potential (or dimensionless potential) has a maximum (or minimum) or saddle point. There are a total of 5 Lagrange points (L1, L2, L3, L4 and L5). The L1 point is the Lagrange point located between the two stars. One finds the positions of and the values of the potential at the Lagrange points by solving
$$ \nabla \Omega(r,\theta,\phi)= 0. $$
The Roche lobe (of star 1) is the surface surrounding star 1 having the potential equal to the potential of the L1 point. The star can underfill its Roche lobe, in that case its surface will still correspond to an equipotential surface. Here we use the same definition of fill-out factor, F, as that introduced by Mochnacki (1984). The dimensionless potential used here is related to the potential C of Mochnacki (1984) by:
$$ \Omega= \frac{(1+q)C}{2}- \frac{q^{2}}{2(1+q)}. $$
F is defined by \(F=C1/C\) with C1 the value of potential for the Roche lobe and C the value of the potential for the surface defined by fill-out factor F. Thus the dimensionless potential for the Roche lobe is related to C1 by \(\Omega_{L1}= \frac{(1+q)C1}{2}- \frac{q^{2}}{2(1+q)}\). The potential for the surface specified by F is here called \(\Omega_{F}\). Thus the value of this potential for a given value of F, in terms of mass ratio q and potential \(\Omega L1\) at the L1 point, is:
$$ \Omega_{F}= \frac{\Omega_{L1}+\frac{q^{2}}{2(1+q)}}{F}- \frac{q^{2}}{2(1+q)}. $$
The purpose of the calculator is to allow one to calculate the precise size of the Roche lobe along the main coordinate axes or in any specified direction (θ, ϕ). Sizes are given along the x-axis, from the center of the star to the L1 Lagrange point, toward star 2 (positive x, called \(R_{L1}\)), and on the back side away from star 2 (negative x, called \(R_{bk}\)). They are given along the y-axis (called \(R_{y}\)) and along the z-axis (called \(R_{z}\)). The Roche lobe has size \(R_{y}\) along the positive and negative y-axes and size \(R_{z}\) along positive and negative z-axes. Additional values are given for the x-coordinate of the L2 and L3 points along the x-axis (called \(X_{2}\) and \(X_{3}\)), and the values of the potential at the L1 point (\(\Omega_{1}\)), at the L2 point (\(\Omega_{2}\)) and at the L3 point (\(\Omega_{3}\)). The area, A, and volume, V, of the Roche lobe are given (in units of \(a^{2}\) and \(a^{3}\), respectively), and the equivalent volume radius \(R_{eq}\), defined by \(V=4 \pi R_{eq}^{3}/3\) and area radius \(R_{a}\), defined by \(A=4 \pi R_{a}^{2}\). For the case of \(F<1\), i.e. star 1 inside its Roche lobe, the values of \(R_{bk}\), \(R_{y}\), \(R_{z}\), A, V, \(R_{a}\) and \(R_{eq}\) are for the stellar surface. Additionally the front radius, \(R_{fr}\) of the star (along the positive x-axis) is given because it is not equal to \(R_{L1}\) for \(F<1\).
The L1, L2 and L3 Lagrange points are along the x-axis (\(\mathrm{y}=\mathrm{z}=0\)). A root-finding routine is used to find the x-value where the derivative of \(\Omega(r,\theta,\phi)\) along the x-axis is zero. The root is found to an accuracy of 1 part in 108 (Java version) and 1 part in 109 (Fortran version) using Brent's method (Brent 1973). This method combines the bisection method, the secant method and inverse quadratic interpolation. It has been shown previously that the L4 and L5 points lie in the x-y plane, and either of these when taken with the mass centers of \(M_{1}\) and \(M_{2}\) form an equilateral triangle, so are at fixed positions independent of mass ratio. For further discussion of the Lagrange points see the textbook Rosswog and Brüggen (2007).
Values of sizes of the Roche lobe along the major axes, positions of the Lagrange points and the potentials are calculated to accuracy which is limited by the root-finding algorithm (the accuracy of the root finding algorithm can be changed by changing tolerance parameter in the root-finding subroutine of the source code). The Roche lobe volume V and area A, and the radii \(R_{eq}\) and \(R_{a}\), are obtained by using a Romberg integration method in the Fortran program and have been tested to be accurate to 8 digits (see Section 4 below). A set of results of these integrated quantities were calculated using the Fortran code for ≃150 different values of q between 0.01 and 100, nearly equally spaced in \(\log(q)\). These were used as input to the Java code, within which V, A, \(R_{eq}\) and \(R_{a}\) are calculated by interpolation. As a result the integrated quantities from the Java calculator are accurate to 4 digits or better in the range of q from 0.01 to 100.
To get the radius of the Roche surface in an arbitrary direction, \(R(\theta,\phi)\), the user enters values of q, θ and ϕ (with angles in radians or degrees). The root-finding algorithm, using Brent's method, finds the radius in the requested direction which has potential equal to the L1 potential. The user has the option of entering a text file containing a list of (q, θ, ϕ) values. An example input list is given by the first 3 columns of Table 1. The output is the input file with an added fourth column, as shown in Table 1. The output can be used, for example, to plot the Roche lobe surface in 3-D or to obtain a cross section of the Roche lobe viewed from any specific direction. To get the stellar surface for the case of a star which is under-filling its Roche lobe, one specifies an input value of \(F<1\). Then the calculator calculates \(R(\theta,\phi)\) for points with the potential equal to \(\Omega_{F}\) (equation (4)). For \(F \le 0.1\) the star is very nearly spherical.
Table 1 Example input ( q , θ , ϕ ) and output ( \(\pmb{R(\theta,\phi)}\) )
Roche lobe radii dependence on q, p and F
Figure 2 shows a plot of the radii \(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) as a function of q for the range \(0.01 \le q \le 100\) with the radii calculated using the code described here. The radii obey the relation \(R_{z}< R_{y}<R_{bk}<R_{L1}\) for all q. \(R_{eq}\) is larger than \(R_{y}\) for \(q>0.38\), and smaller than \(R_{y}\) for \(q>0.38\). Figure 3 shows the radii as a function of dimensionless rotation rate, p, for two values of q (1 and 10). The dependence of the radii on p is very similar for the two values of q: \(R_{L1}\) and \(R_{z}\) depend most strongly on p; \(R_{y}\) and \(R_{bk}\) depend most weakly on p, and \(R_{eq}\) has an intermediate dependence on p. The dependence on fill-out factor, F is considered next. For \(F<1\) the front radius of the potential surface \(R_{fr}\) is smaller than \(R_{L1}\). Figure 4 shows \(R_{fr}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) depend on F for two values of q (1 and 10) and with \(p=1\). \(R_{fr}\) depends strongly on F, especially for F near to 1. The other radii depend on F less strongly. The shape of the dependence on F is very similar for different values of q.
Roche radii vs. mass ratio. The various Roche radii (\(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) in units of a) as a function of mass ratio q.
Roche radii vs. rotation rate. The various Roche radii (\(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) in units of a) as a function of dimensionless rotation rate p, for \(q=1\) and \(q=10\).
Roche radii vs. fill-out factor. The various Roche radii (\(R_{fr}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) in units of a) as a function of fill-out factor F, for \(q=1\) and \(q=10\), with \(p=1\).
To carry out the calculations of the Roche Lobe properties, a Java version with a graphic interface (see Additional file 1), and a Fortran 90 version, with a command line interface (see Additional file 2), were created.
Fits to Roche lobe radii
The Eggleton formula has been widely used to calculate volume equivalent radius of the Roche lobe for the case of synchronous rotation \(p=1\) (here we note that the q we use is equivalent to \(1/q\) used by Eggleton (1983)). It has accuracy better than 1%. However, no similarly accurate formula for the other radii (\(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\)) has been given. Here, we check the accuracy of the Eggleton formula and provide a somewhat improved version. We also show that using a cubic spline fit yields accurate approximations of the various radii.
First, a wide variety of different analytic formulae were tested to search for a suitable approximation formula for the radii (\(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\)) as a function of q. These formula are more complex than the Eggleton formula (which limits their utility), so that they are not presented here except for the simplest generalization. The simplest generalization which gave good results for \(R_{eq}(q)\) is:Footnote 1
$$ R_{eq}= \frac{a1 q^{a4}}{a2 q^{a5} + \ln(1+a3 q^{a4+1/3})}. $$
The values of the coefficients were obtained by least-squares fitting of the formula to ≃150 numerically calculated values of \(R_{eq}\) for values of q ranging from 0.01 to 100. The resulting values of the coefficients are: \(a1=0.64334\), \(a2=0.86907\), \(a3=1.2809\), \(a4=-0.74303\), and \(a5=0.73103\). The root-mean-square (r.m.s.) of the residuals for the formula values compared to the accurate numerical values is \(2.2\times10^{-5}\). For the above \(R_{eq}\) formula, the r.m.s. increases to 4.5×10−5 with 4 digits, and to 7.1×10−4 with 3 digits. For comparison the r.m.s. for the Eggleton formula compared to the accurate numerical values is \(1.5\times10^{-3}\).
A simple alternative to using complex analytical formula is a cubic spline fit to a small table of radii vs. q. This gives good results, and is simple to implement. A cubic spline fit to the radii vs. \(\log(q)\) was found to give significantly better results than a fit to the radii vs. q, so \(\log(q)\) is used as the independent variable here. The small input table for the cubic spline fit consists of numerically calculated values of \(R_{eq}\), \(R_{L1}\), \(R_{bk}\), \(R_{y}\) and \(R_{z}\) (accurate to 8 decimal places), at 28 values of q between 0.01 and 100. The q values were spaced in ratios of square root of 2 (then rounded to 2 digits) resulting in 28 q values. The resulting cubic spline is compared to values of the various radii numerically calculated at 153 values of q between 0.01 and 100. When the cubic spline is evaluated as a function of q, then the r.m.s. of the spline fit compared to accurate numerical values ranges from \(1.2\times10^{-5}\) to \(1.9\times10^{-5}\) for the different radii. When the cubic spline is evaluated as a function of \(\log(q)\), then the r.m.s. of the spline fit compared to accurate numerical values ranges from \(1.1\times10^{-7}\) to \(1.4\times10^{-7}\) for the different radii. For \(R_{eq}\), the r.m.s. is \(1.4\times10^{-7}\), which is 104 times smaller than the Eggleton formula. The results of the residuals of the cubic spline fits from accurate numerical values are shown in Figure 5, together with the residuals of the Eggleton formula from accurate numerical values.
Accuracy of fits to Roche radii. The difference between a cubic spline fit (see text for details) for \(R_{L1}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\) and \(R_{eq}\) and the numerically calculated values, as a function of mass ratio q. Also shown is the difference between the Eggleton formula (\(R_{Egg}\)) and the accurate \(R_{eq}\), divided by 1,000. This shows that the spline fit is a factor of ∼10,000 better than the Eggleton formula for \(R_{eq}\). The vertical axis has been scaled up by a factor of 106.
Discussion and conclusions
Here we present a freely available calculator for various quantities related to the Roche lobe of a star in a binary system. The calculator is provided in two versions: a Java version which is easily installed and run on any computer system with Java, and a Fortran 90 version, which can be easily used and edited by anyone familiar with Fortran 90. The calculator provides two general types of output. The first type yields values of \(R_{L1}\), \(\Omega_{1}\), \(X_{2}\), \(\Omega_{2}\), \(X_{3}\), \(\Omega_{3}\), \(R_{bk}\), \(R_{y}\), \(R_{z}\), A, V and \(R_{eq}\), where the various symbols are defined above, given input values of mass ratio, q, stellar rotation parameter, p, and stellar fill-out factor, F. The second type gives the Roche radius in any direction (θ, ϕ) specified by the user. Both types allows a single input parameter set, or multiple sets of input parameters in the form of a text file, so that a large number of values of Roche radius can be calculated quickly and conveniently.
The present calculator can be used to calculate the instantaneous Roche lobe for elliptical orbits, for synchronous or non-synchronous stars, by using the value of \(\frac{p (1+e)^{2}}{(1+e \cos(\nu))^{3/2}}\) instead of p as an input, as explained in Section 1 above using equations (2) and (3). The calculator has only been tested for the range \(0.1 \le p \le 2\), which restricts the corresponding range of p, e and ν for the elliptical orbit case.
To check the accuracy and precision of the calculations done here we compared calculated values from the Fortran 90 version to the most recent published tables of values for the Roche radii, done by Pathania and Medupe (2012). For \(R_{L1}\) our values, in all cases they provide, agree to the 6 digits given by Pathania and Medupe (2012). For \(R_{bk}\) our values agree in all cases to 5 digits, and differ in the 6th digit in about half the cases. For \(R_{y}\) and \(R_{z}\), we agree to 4 digits, but differ in the 5th digit in many cases. Because we have calculated the radii using 3 different codes which agree, and the Fortran 90 version has good error control, we believe that the values we present are more precise than previously published values.
To check for accuracy to more than 5 digits we rely on tests that we have carried out. The Fortran 90 code was written with all calculations done in double precision. The Fortran version was run on three different Linux machines with two different flavours of Linux (RedHat Fedora and CentOS) and gave the same numerical results in all cases. Additionally, a MathCad version (Parametric Technology Corporation) was written on a Windows 7 computer system to calculate all quantities accurate to 8 digit accuracy. The Mathcad version was very slow to run for error tolerance 10−9 or less so it was not used to calculate higher accuracy results. Also the Mathcad program was written to use different integration and root-finding algorithms than the Fortran program as a check. The MathCad values agreed with the Fortran 90 values to 8 digits or better. The calculated values from the Java version and the Fortran 90 version agree to 7 digits or better for \(R_{L1}\), \(\Omega_{1}\), \(X_{2}\), \(\Omega_{2}\), \(X_{3}\), \(\Omega_{3}\), \(R_{bk}\), \(R_{y}\) and \(R_{z}\). For quantities obtained by interpolation in the Java version (A, V and \(R_{eq}\)) the accuracy is lower (but still better than 4 digits). The result is that we believe, for the above reasons, the Fortran 90 program gives results which are accurate and precise to 8 digits.
The above software is available in both Java and Fortran 90 versions. The Java version is called the Roche Radius Calculator, and is freely available as a zip file, with simple installation and usage instructions, from the Quark Nova Project website at the University of Calgary (at http://quarknova.ucalgary.ca/, follow the 'software' tab to the entry 'RoLo' as the link 'Download'). The Fortran 90 version is available as Fortran 90 source code in a single file at the same website under the 'software' tab to the entry 'RoLo' as the link RocheLobe.f90.
Dr. P. Eggleton suggested to include the large q limit that \(R_{eq}\propto q^{-1/3}\) which implies that the exponent of q in the logarithm term is \((a4+1/3)\).
Brent, RP: Some efficient algorithms for solving systems of nonlinear equations. SIAM J. Numer. Anal. 10, 327-344 (1973)
Eggleton, PP: Approximations to the radii of Roche lobes. Astrophys. J. 268, 368 (1983)
Kopal, Z: Close Binary Systems. Chapman and Hall, London (1959)
Limber, DN: Surface forms and mass loss for the components of close binaries-general case of non-synchronous rotation. Astrophys. J. 138, 1112 (1963)
Mochnacki, SW: Accurate integrations of the Roche model. Astrophys. J. Suppl. Ser. 55, 551-561 (1984)
Pathania, A, Medupe, T: Radius of the Roche equipotential surfaces. Astrophys. Space Sci. 338, 127-145 (2012)
Plavec, M, Kratochvil, P: Tables for the Roche model of close binaries. Bull. Astron. Inst. Czechoslov. 15, 165 (1964)
Rosswog, S, Brüggen, M: Introduction to High-Energy Astrophysics. Cambridge, Cambridge (2007)
Sepinsky, JF, Willems, B, Kalogera, V: Equipotential surfaces and Lagrangian points in nonsynchronous, eccentric binary and planetary systems. Astrophys. J. 660, 1624-1635 (2007)
DAL thanks Prof. Didier Barret and Dr. Natalie Webb for hospitality at Institut de Recherche en Astrophysique et Planétologie (CNRS/UPS/OMP), Toulouse, France, where part of this work was completed. The author thanks the referees for valuable suggestions which were implemented. This work was supported by the Natural Sciences and Engineering Research Council of Canada.
Department of Physics and Astronomy, University of Calgary, 2500 University Dr. NW, Calgary, Canada
Denis A Leahy
University of Calgary, 2500 University Dr. NW, Calgary, Canada
Janet C Leahy
Search for Denis A Leahy in:
Search for Janet C Leahy in:
Correspondence to Denis A Leahy.
DAL carried out the overall design, analytical work, and Fortran 90 and MathCad programming and calculations. The Java programming and Java interface design were done by JCL. All authors read and approved the final manuscript.
Below are the links to the electronic supplementary material.
Roche Radius Calculator.
Fortran 90 source code.
Leahy, D.A., Leahy, J.C. A calculator for Roche lobe properties. Comput. Astrophys. 2, 4 (2015) doi:10.1186/s40668-015-0008-8