Plotting the Roche Potential Surface of a Cataclysmic Binary Star System

by Jeff Bryant

Introduction

The shape of any star is determined by the gravitational forces within its environment.  A star that is unaffected by gravitational forces other than its own should be spherical, assuming that the star is not oblate due to rotational effects.  The spherical nature of a lone star is due to the fact that all gravitational forces acting on the star are centrally directed and symmetric about all three axes.  If a star, however, is in orbit around a companion star, the gravitational effects of that companion star play a role in the determination of the shape  of the star in question.
     Throughout this notebook, we will be dealing with a certain type of binary star system, namely, a cataclysmic binary star system.  These systems contain a white-dwarf primary star of approximately one solar-mass and about the size of the earth and a small red-dwarf secondary star about half the size of our sun.  These two stars orbit each other in a volume smaller than the volume of our sun, thus resulting in orbital periods on the order of hours.

Initializing parameters and Defining the Roche Potential

Most of the interesting gravitational distortions take place on the secondary star.  In order to mathematically calculate the surface of the secondary star, the mass ratio must first be defined.  This quantity is defined as \(TraditionalForm\`q = \(Mass\ of\ the\ secondary\ star\)\/\(Mass\ of\ the\ primary\ star\)\).  Keep in mind that the secondary star is usually less massive than the primary star, thus giving values of q less than one.  The next cell defines several parameters of the system.

q =  .4; 
μ = ((M1 q)/(M1 + q M1)); 
a = 1; 
M1 = 1; 
G = 1; 
Porb = (Sqrt[((4 (π^2) (a^3))/(G (1 + q)))]); 
Ωorb = ((2 π)/Porb); 
M2 = q M1;

With the mass ratio defined, the gravitational potential per unit mass of the two stars can be calculated.  All scaling factors have been removed such as the mass of the primary star, the gravitational constant, G, and the distance between the stars.  This enables us to constrain the resulting numbers into reasonable values.  The equation below has three terms: the first two treat each of the stars as a point mass.  This may at first glance seem an unreasonable assumption but, because most of the mass of a star lies near its core where fusion occurs, the outer layers contribute little to the overall mass.  The first term represents the potential of the primary star and the second term represents the potential of the secondary star (notice the difference in the two terms).  The third term is a result of the orbital motion of the two stars.

Φ[x_, y_, 
        z_] := -((G M1)/(Sqrt[(((x + (q/(1 + q))))^2) + (y^2) + (z^2)])) - ((G M1 q)/(Sqrt[(((x - a + (q/(1 + q))))^2) + (y^2) + (z^2)])) - (1/2) Ω(orb^2) ((((x - μ + ((q a)/(1 + q))))^2) + 
              (y^2))

Plotting the Roche Potential

Plot[Φ[x, 0, 0], {x, -1, 2}, 
      PlotRange -> {0, -15}, 
      FrameLabel -> {"x/a", "Potential"}, Frame -> True, 
      TextStyle -> {FontFamily -> "Arial", 
          FontWeight -> "Bold"}, ImageSize -> 450]

Plotting Φ(x, y, 0) allows us to view the wells as a gravitational distortion on a surface.  This view makes the ice cube bin analogy seem more appropriate.  In this view, the potential well of the secondary star is in the foreground.

Plot3D[Φ[x, y, 0], {x, -1, 2}, {y, -1.5, 1.5}, 
      PlotRange -> {0, -5}, Boxed -> False, Axes -> False, 
      ViewPoint -> {2.841,  0.144,  1.833}, PlotPoints -> 50, 
      ImageSize -> 450]

By looking at the first plot, it can be seen that the critical point at which matter transfer occurs, the L1 point (first Lagrange point), is somewhere near x = 0.6.   A second Lagrange point, L2, is located just beyond the secondary star near x = 1.5.  Any star filling its Roche lobe beyond L2 will result in matter loss from the system.  Similarly, L3 is located behind the primary star's Roche lobe near x = -1.  The following compound input statement numerically solves for these critical points.  Although L4 and L5 are not solved for, they are seen as two circular features 90 degrees on either side of the x-axis in the next plot.

deriv = D[Φ[x, 0, 0], x];
posAtL1 = FindRoot[deriv == 0, {x,  .3}];
ΦAtL1 = Φ[x /. posAtL1, 0, 0];
posAtL2 = FindRoot[deriv == 0, {x, 1.5}];
ΦAtL2 = Φ[x /. posAtL2, 0, 0];
posAtL3 = FindRoot[deriv == 0, {x, -1}];
ΦAtL3 = Φ[x /. posAtL3, 0, 
          0];
posAtL4 = 
        FindRoot[
          D[Φ[x /. posAtL1, y, 0], y] == 0, {y, 1}];
ΦAtL4 = Φ[x /. posAtL1, y /. posAtL4, 
          0];

If we now look at the x-y plane and do a contour plot and a density plot, the shape of the secondary star becomes more evident.  The secondary star assumes the teardrop shape on the right side of the plot.  However, due to the small size of the white-dwarf primary star, it sits at the bottom of its well and assumes the shape of the well at that depth, a circle (or a sphere in three dimensions).  If the mass of the primary star were kept constant but its size were increased then the primary star would begin to deviate from a circular/spherical shape.  In other types of star systems, both stars may assume teardrop shapes or may actually come in contact with each other, thus becoming contact binary stars.

\(ContourPlot[Φ[x, y, 0], {x, \(-1.5\), 2}, {y, \(-1.75\), 
        1.75}, Contours -> {ΦAtL1, \(-3.6\), ΦAtL2, ΦAtL3, \(-1.85\), \(-1.9\), \(-1.96\), \(-2\)}, ContourShading -> True, PlotPoints -> 100, 
      AspectRatio -> 1, FrameLabel -> {"x/a", "y/a"}, 
      ImageSize -> 450, ColorFunction -> \((Hue[ .8\ #] &)\), 
      TextStyle -> {FontFamily -> "Arial", 
          FontWeight -> "Bold"}]\)

Plotting the Individual Stars

Now that the potential surfaces are known, the radius of the secondary star, in units of the separation of the two stars, can be found.  This is because the secondary fills one of the potential surfaces: the one that passes through the L1 point.  This surface is not spherical, so the radius will be different depending on the direction you are measuring, especially in the x-direction.

Solve Φ(x, y, z) for the radius of the secondary star(in units of a) in different directions

This finds the x radius towards L1

RL1 = 1/(1 + q) - (x /. posAtL1)

This finds the x radius away from L1

RBack = 
      Min[Cases[
            x /. Solve[Φ[x, 0, 0] == ΦAtL1, x], 
            n_ /; Re[n] > 1/(1 + q)]] - 1/(1 + q)

This finds the radius in the y-direction

Ry = 
      Min[Cases[
          y /. Solve[Φ[1/(1 + q), y, 
                  0] == ΦAtL1, y], n_ /; n > 0]]

This finds the radius in the z-direction

Rz = 
      Min[Cases[
          z /. Solve[Φ[1/(1 + q), 0, 
                  z] == ΦAtL1, z], n_ /; n > 0]]

The next and last step is to render the entire 3-dimensional surface of the secondary star and place a small star nearby at the location of the white dwarf.  This is done by using the ContourPlot3D package, which must first be loaded.  We also need to generate transparent wireframe surfaces.  Both of these are loaded as follows.  (This actually loads more than we need, but is fine for simplicity.)

<< Graphics`

With this done, the 3-dimensional Roche potential surfaces can be plotted.  Because we only want to view the surfaces passing through the L1 and L2 points, the range on the plotted variables is adjusted so that only these surfaces are viewed.  In addition, a red color is added to the L1 surface.  The surface passing through the L2 point is drawn in wireframe so as to allow the viewing of the inner L1 surface as well.

gr = 
        Show[Insert[
            ContourPlot3D[Φ[x, y, z], {x, x /. posAtL3, 
                x /. posAtL2}, {y, - .6,  .6}, {z, - .6,  .6}, 
              Contours -> {ΦAtL1}, 
              PlotPoints -> {7, 7}, Boxed -> False, 
              ContourStyle -> SurfaceColor[RGBColor[1, 0, 0]], 
              AmbientLight -> RGBColor[0, 0, 0], 
              ViewPoint -> {-0.011, -3.199, 1.104}, 
              DisplayFunction -> Identity], 
            EdgeForm[], {1, 1}]];
gr2 = 
        Show[WireFrame[
            ContourPlot3D[Φ[x, y, z], {x, x /. posAtL3, 
                x /. posAtL2}, {y, - .6,  .6}, {z, - .6,  .6}, 
              PlotPoints -> {8, 8}, Boxed -> False, 
              Contours -> {ΦAtL2}, 
              ViewPoint -> {-0.011, -3.199, 1.104}, 
              DisplayFunction -> Identity]]];
Show[{gr, gr2}, DisplayFunction -> $DisplayFunction, 
      ImageSize -> 400, Lighting -> True]

One can focus in on the shape of the secondary star by doing another ContourPlot3D, but in this case changing the plot range and lighting effects.  The secondary star fills the equipotential surface that passes through the L1 point.

gr3 = 
      Show[Insert[
          ContourPlot3D[Φ[x, y, z], {x, x /. posAtL1, 
              1/(1 + q) + RBack}, {y, - .6,  .6}, {z, - .6,  .6}, 
            Contours -> {ΦAtL1}, PlotPoints -> {7, 7}, 
            Boxed -> False, 
            ContourStyle -> 
              SurfaceColor[RGBColor[ .6, 0, 0], RGBColor[1, 1, 1], 20], 
            ViewPoint -> {-0.011, -3.199, 1.104}, 
            LightSources -> {{{-1.5, 0, 0}, RGBColor[1, 1, 1]}}, 
            DisplayFunction -> Identity], EdgeForm[], {1, 1}], 
        DisplayFunction -> $DisplayFunction, Background -> Black, 
        AmbientLight -> RGBColor[ .5,  .5,  .5]]

The VectorAnalysis package is loaded next in order to more easily control the viewpoint.

<< Calculus`VectorAnalysis`

This next step generates a spherical star to be used as the white dwarf.

gr4 = 
      Show[Graphics3D[{EdgeForm[], SurfaceColor[White], 
            TranslateShape[
              Sphere[ .01, 20, 20], {-q/(1 + q), 0, 0}]}]]

Finally, we combine the two stars into a single frame.

Show[gr3, gr4, 
      ViewPoint -> 
        CoordinatesToCartesian[{4, N[80 Degree], N[270 Degree]}, 
          Spherical], 
      PlotRange -> {{-q/(1 + q) -  .2, 
            x /. posAtL2}, {- .3,  .3}, {- .3,  .3}}, 
      Boxed -> False, Axes -> False, Axes -> True, 
      Background -> RGBColor[0, 0, 0], 
      AmbientLight -> RGBColor[ .5,  .5,  .5], ImageSize -> 400, 
      LightSources -> {{{-q/(1 + q), 0, 0}, White}}]

We can specify an inclination angle for the system and create an animation of the orbit.

With[{inclination = 70 Degree}, 
Table[
        Show[{gr3, gr4, Graphics3D[Text["+", {0, 0, 0}]]}, 
          ViewPoint -> 
            CoordinatesToCartesian[{4, N[inclination], t}, Spherical], 
          PlotRange -> {{-(x /. posAtL2), 
                x /. posAtL2}, {- .3,  .3}, {- .3,  .3}}, 
          Boxed -> False, Axes -> False, Axes -> True, 
          Background -> RGBColor[0, 0, 0], 
          AmbientLight -> RGBColor[ .5,  .5,  .5], 
          ImageSize -> 400, 
          LightSources -> {{{q/(1 + q) Sin[t],  
                  q/(1 + q) Cos[t] Cos[inclination], -q/(1 + q) 
                    Cos[t] Sin[inclination]}, White}}, 
          SphericalRegion -> True, ViewVertical -> {0, 0, 1}], {t, 
          2 π, 0, -2 π/20}]
]