Mass changing orbits

A few weeks ago my good friend John Burk posted some intriguing questions about what happens to planetary orbits as the sun loses mass (all that heat has to come from somewhere!). I’ve been thinking about it ever since and finally got around to doing some modeling to see if I could answer any of the questions.

My first thought was to see if there were ways to simplify the differential equation solving approach. What John was doing was the full 3D version assuming the sun’s mass stays significantly above the mass of the planet. So he’s already doing some simplifications because he’s not bothering with the center of mass frame. I think he’s right to do that because after a page or so of notes, I’ve realized that the center of mass frame gets pretty ugly when one of the participants is losing mass.

So what else can be simplified? The great thing about central force problems is that you can reduce it all the way down from six variable to one (in fact that’s one of my standards when I teach Theoretical Mechanics):

1. Each of the two masses has 3 variables (x, y, and z) so you start with six.
2. The center of mass approach lets you model the problem as a fictitious mass (=m1m2/(m1+m2)) that’s the same distance from a fictitious force center as the two actual masses are from each other. Now you’re down to three.
3. If it’s a central force, the angular momentum is conserved. That means the fictitious particle has to stay in a plane. Now you’re down to two variables.
4. If the angular momentum is conserved, you can treat the rotational part of the kinetic energy as an effective potential energy, leaving only the radius variable. Now you’re down to one.

You can model the complex 6-dimensional problem as a single (reduced) mass experiencing a potential energy function given by:

where is the reduced mass and U is the potential energy that is only a function of r. To get the force this fictitious one dimensional particle feels, you just need to take a derivative (and add a negative sign).

So I gave it a try. The first thing I did was try to see how far into the future I could integrate using Mathematica. It turns out I could go quite a ways! Here’s a plot of the radius as a function of time.

As you can see, I was able to go out several billion seconds of integration time. This turned out to be around a billion “years” given the simple parameters I chose. If the radius grows, we expect the years to take longer. Here’s a plot of the instantaneous “year time” over the simulation:

as expected!

So it seems that a very slow loss of the sun’s mass would just slowly increase both the circle radius and the year time of the planet.

Another question that John asked was whether there might be an analytical solution to this. I quickly tried DSolve instead of NDSolve in Mathematica and got no joy (I wasn’t overly hopeful). I did ask a good friend of mine who’s a real expert in differential equations whether he knew of a particular decay function I could use that might have an analytical solution. He couldn’t think of one, but did point out that if you had the mass changes be discrete you could “easily” build up the solution since in between the mass changes you’ve got a simple inverse square law orbit that does have an analytical solution.

What he means is that if you start with, say, a circular orbit, you can predict exactly where the planet will be and what direction it’s traveling (and what speed) when the first mass change happens. When it does, you now have initial conditions for a slightly different inverse square law problem. Because the sun has lost mass, it doesn’t pull as hard as before so the circular orbit becomes an ellipse. That ellipse is fully analytical and you can figure out everything you need to know about the planet at the next mass drop. Repeat this to your hearts content and you’ve got a piecewise analytical solution.

This sounded intriguing, but I started to wonder what the physical differences would be between the two approaches. I figured I could check on a relatively small time scale and look for differences. So I coded up both a continuous and a discrete mass loss model, where they connect with each other after each mass loss. Here’s a plot of both mass loss functions:

Here’s the animated result (the blue dots are the points of the mass change for the discrete (blue) model:

As you can see, there’s a pretty noticeable difference in the orbits. Admittedly this is only because the mass jumps are pretty big, but it still makes me nervous.

Here’s a plot of the radius function for both:

Note how the continuous one seems to never get any close to the sun while the blue one is clearly showing more of an elliptical motion (it gets closer and further from the sun during every orbit). To see that more clearly, here’s a plot of r'[t] or the rate of change of the radius:

It sure looks like the orange line doesn’t go negative (indicating the planet never gets closer to the sun). Here’s a zoom in of the orange one to see it better:

Yep, never negative!

Your thoughts? Here are some starters for you

• This is cool, but I’d like to hear more about . . .
• This is dumb. Why didn’t you talk about this instead:  . . .
• It looks like you were doing a Hamiltonian approach in your notes. I thought you hated the Hamiltonian approach!
• You do know that the perturbations due to Jupiter alone would totally wash out these small effects, right?
• As soon as I saw that you don’t bother to label your axes, I stopped reading. Thanks for saving me some time.
• I thought you said you were going to try to do everything in python from now on? Liar!
• Why didn’t you set up your constants so that a year takes a year? Seems obvious to me.
• Can you share your code?
• How well do your students do on the 6->3->2->1 standard?
• In the piecewise analytical approach, could you look at the solution when the gap time between mass changes goes to the limit of zero?