Roche potential

Let there be two masses \(m_1\) and \(m_2\) with \(m_1 \ge m_2\) and \(m_1\) centered at \((x,y,z) = (0,0,0)\) and \((x,y,z) = (a,0,0)\), where \(a\) is the separation between the masses. Both objects are located in the x,y plane with x-axis counting distance along the connecting line and z measuring distance perpendicular to the orbital plane.

For synchronous rotation and a circular orbit, the Roche potential reads (e.g., Hilditch 2001 “An introduction to close binary stars”)

\[\begin{split}\Phi &= -\frac{G m_1}{r_1} - \frac{G m_2}{r_2} - \frac{\omega^2}{2}\left[\left(x - a\frac{m_2}{m_1+m_2}\right)^2 + y^2 \right] \\\end{split}\]

with

\[\begin{split}r_1 &= \sqrt{x^2 + y^2 + z^2} \\ r_2 &= \sqrt{(x-a)^2 + y^2 + z^2} \\ q &= \frac{m_2}{m_1} \\ M &= m_1 + m_2 \\ m_1 &= \frac{M}{1+q} \\ m_2 &= \frac{Mq}{1+q} \\ \omega^2 &= \frac{GM}{a^3} \;\;\; \mbox{(Kepler's third law)} \\ x', y', z', r_1', r_2' &= \frac{x}{a}, \frac{y}{a}, \frac{z}{a}, \frac{r_1}{a}, \frac{r_2}{a}\end{split}\]

we write

\[\begin{split}\Phi &= -\frac{G M}{(1+q) r_1} - \frac{G M q}{(1+q) r_2} - \frac{GM}{2 a^3}\left[\left(x - \frac{a\,q}{1+q}\right)^2 + y^2 \right] \\ \Phi &= -\frac{G M}{(1+q) r_1'\,a} - \frac{G M q}{(1+q) r_2'\,a} - \frac{GM}{2 a^3}\left[\left(x'\,a - \frac{a\,q}{1+q}\right)^2 + y'^2\,a^2 \right] \\ \Phi &= -\frac{GM}{2 a} \left( \frac{2}{(1+q) r_1'} + \frac{2 q}{(1+q) r_2'} + \left(x' - \frac{q}{1+q}\right)^2 + y'^2 \right) \\ \Phi &= -\frac{GM}{2 a} \; \Phi_n(x', y', z', q)\end{split}\]

where \(\Phi_n\) is the dimensionless Roche potential.

Note that moving material in the rotating frame is also subject to Coriolis forces, which are non-conservative and, therefore, cannot be incorporated into the scalar Roche potential.

The module provides functions to obtain the derivative of the dimensionless Roche potential with respect to the dimensionless variables representing the three directions. To turn the result into real, physical numbers, the following conversion is required

\[\frac{\partial \Phi}{\partial x} = -\frac{G M}{2 a^2} \frac{\partial \Phi_n}{\partial x'}\]

Example of usage

from __future__ import print_function, division
from PyAstronomy import pyasl
import numpy as np
import matplotlib.pylab as plt

x, y = np.linspace(-1.5, 2, 300), np.linspace(-1.6, 1.6, 300)
xx, yy = np.meshgrid(x, y)
# Coordinates in orbital plain
z = 0

# Mass ratio
q = 0.2

# Get dimensional values of Roche potential
p = pyasl.rochepot_dl(xx, yy, z, q)

# Positions (and potentials) of Lagrange points
l1, l1pot = pyasl.get_lagrange_1(q)
l2, l2pot = pyasl.get_lagrange_2(q)
l3, l3pot = pyasl.get_lagrange_3(q)
l4, l5 = pyasl.get_lagrange_4(), pyasl.get_lagrange_5()
l4pot = pyasl.rochepot_dl(l4[0], l4[1], l4[2], q)
l5pot = pyasl.rochepot_dl(l5[0], l5[1], l5[2], q)


print("Effective (dimensionless) radii of first and second mass")
print("According to the approximation of Eggleton 1983:")
r1eff = pyasl.roche_lobe_radius_eggleton(q, 1)
r2eff = pyasl.roche_lobe_radius_eggleton(q, 2)
print("    Reff1: %5.3f" % r1eff)
print("    Reff2: %5.3f" % r2eff)
print()
print("Roche volume and effective radius from Monte Carlo integration:")
mcvol1 = pyasl.roche_vol_MC(q, 1)
mcvol2 = pyasl.roche_vol_MC(q, 2)
print("    MC Roche lobe volume 1: %6.4f +/- %6.4f" % (mcvol1[0:2]))
print("    MC Roche lobe volume 2: %6.4f +/- %6.4f" % (mcvol2[0:2]))
print("    MC effective radius 1: %6.4f +/- %6.4f" % (mcvol1[2:]))
print("    MC effective radius 2: %6.4f +/- %6.4f" % (mcvol2[2:]))

plt.contour(p, [l5pot*1.02, l3pot, l2pot, l1pot], colors=['g', 'c', 'b', 'r'], extent=[-1.5, 2, -1.6, 1.6])
plt.text(l1, 0, 'L1', horizontalalignment='center')
plt.text(l2, 0, 'L2', horizontalalignment='center')
plt.text(l3, 0, 'L3', horizontalalignment='center')
plt.text(l4[0], l4[1], 'L4', horizontalalignment='center')
plt.text(l5[0], l5[1], 'L5', horizontalalignment='center')
plt.show()

Example: Substellar, polar, and side radii of Earth and hot Jupiter planets

from PyAstronomy import pyasl
from PyAstronomy import constants as PC

# Shape of the Earth and its Roche lobe (no rotation of Earth)
pc = PC.PyAConstants()
pc.setSystem("SI")

# Mass ratio (m2/m1)
q = pc.MEarth/pc.MSun
print(f"Earth/Sun mass ratio = {q}")

# Radii of Roche lobe along Earth--Sun connecting line, polar (out-of-plane),
# and side (in plane)
rads = pyasl.get_epradius_ss_polar_side(q, pot=None)
# Convert into km
rads = [r*pc.AU/1e3 for r in rads]
print(f"Roche lobe radii (substellar, polar, side) [1e6 km] = " + \
    ", ".join(["%g"%(r/1e6) for r in rads]))

reff_earth = pc.REarth/pc.RJ
# Small epsilon because rocky Earth is a small body on a large
# orbit (by Roche standards)
erads = pyasl.get_radius_ss_polar_side(q, 1, reff_earth, eps=1e-10)

# Convert into km
erads = [r*pc.RJ/1e3 for r in erads]
print(f"Earth radii (substellar, polar, side) [km] = " + \
    ", ".join(["%g"%(r) for r in erads[0:3]]))

print()
print("Roche lobe and planetary shape of typical hot Jupiter")
# Typical hot Jupiter
q = 1e-3
# sma in au
sma = 0.02
# Effective (circular disk) transit radius [RJ]
reff = 1.4

rads = pyasl.get_epradius_ss_polar_side(q, pot=None)
# Convert into km
rads = [r*pc.AU/1e3 for r in rads]
print(f"Roche lobe radii (substellar, polar, side) [1e6 km] = " + \
    ", ".join(["%g"%(r/1e6) for r in rads]))

jrads = pyasl.get_radius_ss_polar_side(q, sma, reff, eps=1e-10)
print(f"Hot Jupiter radii (substellar, polar, side) [RJ] = " + \
    ", ".join(["%g"%(r) for r in jrads[0:3]]))

With output

Earth/Sun mass ratio = 3.0032983882201428e-06
Roche lobe radii (substellar, polar, side) [1e6 km] = 1.49152, 0.952662, 0.995462
Earth radii (substellar, polar, side) [km] = 6378.14, 6378.14, 6378.14

Roche lobe and planetary shape of typical hot Jupiter
Roche lobe radii (substellar, polar, side) [1e6 km] = 10.1264, 6.52393, 6.80657
Hot Jupiter radii (substellar, polar, side) [RJ] = 1.51623, 1.38718, 1.41294

Functionality and API

Roche lobe potential

PyAstronomy.pyasl.rochepot_dl(x, y, z, q)

Dimensionless Roche potential (\(\Phi_n\), synchronous rotation)

More massive component (\(m_1\)) is centered at (x,y,z) = (0,0,0). Less massive component (\(m_2\)) is at (1,0,0). The unit of length is the distance between the objects. Both objects are in the x,y plane (x-axis along the connecting line and z perpendicular to the orbital plane).

Parameters
x, y, zfloat or array

Location(s) at which to calculate the potential. Unit of length is the distance between the masses m1 and m2.

qfloat

Mass ratio (0 <= m2/m1 <= 1)

Returns
Potentialfloat or array

The potential at the specified location(s)

PyAstronomy.pyasl.rochepot(x, y, z, m1, m2, a)

Roche potential (\(\Phi\), synchronous rotation)

More massive component (\(m_1\)) is at (x,y,z) = (0,0,0). Less massive component (\(m_2\)) is at (a,0,0).

Parameters
x, y, zfloat or array

Position [m]

afloat

Separation of masses [m]

m1, m2float

Masses of objects [kg] (m1 >= m2)

Returns
Roche potentialfloat or array

At specified locations [m**2/s**2]

Partial derivatives

PyAstronomy.pyasl.ddx_rochepot_dl(x, q, y=0, z=0)

Derivative of dimensionless Roche potential along x-axis

Parameters
x, y, zfloat or array

Position (default for y and z is zero)

qfloat

Mass ratio m2/m1

Returns
Derivative: float or array

d/dx of dimensionless Roche potential

PyAstronomy.pyasl.ddy_rochepot_dl(y, q, x=1, z=0)

Derivative of dimensionless Roche potential along y-axis

Parameters
x, y, zfloat or array

Position (default for x is one and zero for z)

qfloat

Mass ratio m2/m1

Returns
Derivative: float or array

d/dy of dimensionless Roche potential

PyAstronomy.pyasl.ddz_rochepot_dl(z, q, x=1, y=0)

Derivative of dimensionless Roche potential along z-axis

Parameters
x, y, zfloat or array

Position (default for x is one and zero for z)

qfloat

Mass ratio m2/m1

Returns
Derivative: float or array

d/dz of dimensionless Roche potential

Roche lobe volume and radius

PyAstronomy.pyasl.roche_lobe_radius_eggleton(q, m)

Approximate effective dimensionless Roche lobe radius according to Eggelton 1983 (ApJ 268, 368).

The effective Roche lobe radius is the radius of a sphere with the same volume as the actual equipotential surface defining the Roche lobe.

\[r_L \approx \frac{0.49\,q^{2/3}}{0.6\,q^{2/3} + \ln(1+q^{1/3})}\]
Parameters
qfloat or array

Mass ratio (m2/m1)

mint, {1,2}

Calculate radius for mass 1 or 2

PyAstronomy.pyasl.roche_vol_MC(q, m=2, n=100000, pl=None, eps=0.0001, fullout=True)

Calculate (dimensionless) volume of equipotential surface such as the Roche lobe

Uses Monte Carlo (MC) integration

Parameters
qfloat

Mass ratio (m2/m1)

mint, {1,2}

Whether to calculate volume for m1 or m2

nint, optional

Number of samples to be used in the Monte Carlo integration. Default is 100000. Increase to improve accuracy.

plfloat, optional

The (dimensionless) potential level bounding the volume. If None (default), the Roche lobe potential (L1) is used.

epsfloat, optional

Margin used to avoid singularities of Roche potential. Default is 1e-4.

fulloutboolean, optional

If True (default), provides more output than volume.

Returns
V, Verrfloats

Dimensionless volume of Roche lobe and estimate of uncertainty

Reff, Refferr: floats, optional

Radius of a sphere with the same radius (dimensionless) and estimate of uncertainty. Provided if fullout is True (default).

PyAstronomy.pyasl.roche_yz_extent(q, m=2, pl=None, eps=0.0001)

Extent of equipotential surface in y and z direction

Returns
dy, dzfloat

Distance to equipotential level along y and z direction.

Roche limit

PyAstronomy.pyasl.roche_limit_fluid(q, R2)

Calculate the Roche limit for a fluid body on circular orbit

An expression for the Roche limit, \(a_R\), was given by the eponymous Édouard Albert Roche in 1849 (see, e.g., Rappaport et al. 2013, ApJL 773 for a modern reference)

\[a_R = 2.44 R_1 \sqrt[3]{\frac{\rho_1}{\rho_2}} = 2.44 R_2 q^{-1/3}\]

where \(R_{1,2}\) are the radii of the primary and secondary bodies, \(\rho_{1,2}\) the respective densities, and \(q\) the mass ratio; the second equality relies on spherical shapes.

Parameters
qfloat

Mass ratio

R2float

Radius of secondary body (assuming spherical shape). Units are arbitrary and same as result.

Returns
aRfloat

The Roche distance in the units of R2.

Lagrange points

PyAstronomy.pyasl.get_lagrange_1(q, getdlrp=True, eps=0.0001)

Get location of first Lagrange point

Parameters
qfloat

Mass ratio

getdlrpboolean, optional

If True (default), also the dimensionless Roche potential at the first Lagrange point is returned.

epsfloat, optional

The potential diverges at x=0 and x=1. The search for a root of the derivative is performed in the interval [eps, 1-eps].

Returns
xL1float

The location of the first Lagrange point (xL1, 0, 0).

Potentialfloat, optional

The dimensionless potential at the first Lagrange point (default is True).

PyAstronomy.pyasl.get_lagrange_2(q, getdlrp=True, eps=0.0001)

Get location of second Lagrange point

Parameters
qfloat

Mass ratio

getdlrpboolean, optional

If True (default), also the dimensionless Roche potential at the first Lagrange point is returned.

epsfloat, optional

The potential diverges at x=0 and x=1. The search for a root of the derivative is performed in the interval [1+eps, 1+10*rL], where rL is the estimated (effective) Roche lobe radius according to Eggleton 1983.

Returns
xL2float

The location of the second Lagrange point (xL2, 0, 0).

Potentialfloat, optional

The dimensionless potential at the second Lagrange point (default is True).

PyAstronomy.pyasl.get_lagrange_3(q, getdlrp=True, eps=0.0001)

Get location of third Lagrange point

Parameters
qfloat

Mass ratio

getdlrpboolean, optional

If True (default), also the dimensionless Roche potential at the third Lagrange point is returned.

epsfloat, optional

The potential diverges at x=0. The search for a root of the derivative is performed in the interval [-10, -eps].

Returns
xL3float

The location of the first Lagrange point (xL3, 0, 0).

Potentialfloat, optional

The dimensionless potential at the third Lagrange point (default is True).

PyAstronomy.pyasl.get_lagrange_4()

Get location of forth Lagrange point

Orbital angular momentum is supposed to point into +z direction. L4 is in direction of motion of m2.

Returns
x,y,zfloat

The location of the fourth Lagrange point (xL3, 0, 0).

PyAstronomy.pyasl.get_lagrange_5(getdlrp=True)

Get location of fifth Lagrange point

Orbital angular momentum is supposed to point into +z direction. L5 is behind motion of m2.

Returns
x,y,zfloat

The location of the fifth Lagrange point (xL3, 0, 0).

Coriolis force

PyAstronomy.pyasl.coriolisAcceleration(omega, v)

Get Coriolis acceleration

In the rotating frame, moving material is subject to the Coriolis force, determined by the circular frequency (vector) of the rotating frame and the velocity vector within it.

\[a_C = -2 \vec{\omega} \times \vec{v}\]
Parameters
omegafloat or 3-element 1d array

Circular frequency. If float is given, vector assumed to be omega*(0,0,1).

varray

Either 3-element 1d array with velocity vector in corotating frame or array with shape (nv, 3) with nv velocity vectors (i.e., [::,0] gives all x components of velocity)

Returns
Coriolis accelerationarray

Either 3-element array or [vn,3] array, depending in input shape of v (units determined by input).

Radii in the Roche geometry

PyAstronomy.pyasl.get_epradius_ss_polar_side(q, pot=None, eps=1e-06)

Get dimensionless radii of equipotential surface for secondary mass in three directions

Assumes that the equipotential surface is enclosed within the Roche lobe

Parameters
qfloat

Mass ratio m2/m1

potfloat, optional

Dimensionless potential for equipotential surface. If None (default), the L1 potential will be adopted.

epsfloat, optional

Lowest dimensionless radius considered in the search.

brlboolean, optional

Search below Roche lobe (default is True).

Returns
rxfloat

Substellar radius, height along connecting line (dimensionless)

rzfloat

Polar radius (along z-axis)

ryfloat

Side radius (along y-axis)

PyAstronomy.pyasl.get_radius_ss_polar_side(q, sma, reff, eps=1e-06)

Calculate substellar, polar, and side radii of secondary body

The function first finds the dimensionless Roche potential corresponding to the equipotential surface corresponding to the effective radius and then determines the radii of that surface along the x, y, and z directions.

The effective radius is the radius of a sphere with the same area as the ellipse spanned by the polar and side radii of the equipotential area.

Parameters
qfloat

Mass ratio m2/m1

smafloat

Semi-major axis [AU]

refffloat

Effective radius of the body [RJ]

Returns
rss, rpol, rsidefloat

Radius of equipotential surface at substellar point (along the connecting x-axis), polar radius (along the z-axis), and radius toward the side (y-axis).

pequifloat

Dimensionless potential of equipotential surface corresponding to effective radius.