Hunting Zebras in Outer Space

I once took a graduate course: ”Physiology of Excitable Cells I” run by an eccentric British professor who wanted to teach in a far-reaching classical sense. (In the first class on hearing my surname he commanded ”Next week you will tell the class the connection between your name and Sherlock Holmes”; perhaps a bit off-topic given the course title). One thing he said that has stayed with me is that common problems have common solutions. He was discussing how all the medical students he taught were quick to assume a patient had a rare condition, instead of something more mundane. The tag line they teach medical students is:

“When you hear the sound of hooves, think of horses not zebras”.

I use this phrase a lot when it comes to debugging issues in software. It’s far too easy to believe I am the victim of some fiendishly clever interaction in my system or (so help me) a compiler bug. The answer is (to first order) always that I have done something very dumb or just failed to see an edge case. After a long time writing code I am a lot less likely to get caught by this. However, it did happen to me with N-Body and I was so amazed at how dumb I had been I had to close the laptop and go for walk.

The issue related to the use of sprites for a collision. They seemed to work pretty well except on a level where I had a number of other masses moving around. The latent trap was that I felt guilty about the quick way I had pushed sprites into the code. I decided not to worry too much about optimizing them and allocated a sprite group on the fly. I know it’s better to pre-allocate a sprite pool and all that – but at the time I added them my main focus was elsewhere and I knew I would have to come back and clean them up.

So when I hit an issue where collisions on a busy level resulted in a serious pause in the game loop – I decided the time had come to finally put in sprite pools. That HAD to be the source of the issue. Off I went, with the added fun of feature creep in the refactoring. All in all I spent about ten hours. It made no difference.

The real issue showed itself when I was doing the feature creep part of the refactoring. The loop nesting was wrong and I was running the sprite evolution inside a loop that was running over each body in a level. When there were five masses, I evolved the sprites five times per clock tick and things got really slow. Instead of looking hard for a simple error I had jumped to the conclusion that it must be due memory inefficiency and my use of on-the-fly allocations. Doh! I was amazed I had gone on such a long Zebra hunt.

I have used the Zebra hunting metaphor for such a long time that sometimes I now use it without setting the context properly. I was working with our very bright overseas team in India but thought that this time they were on the wrong track, chasing a theory that seemed too intricate. When my manager asked me how they were making out I responded ”Oh, I think they’re off Zebra hunting”. He replied ”Really?! Do they have a holiday or something?”. I then explained my metaphor as well as the fact that India does not have Zebras…

N-body Choreographies

You would think that after three hundred years the consequences of the law of gravity for “simple” systems would have been pretty well studied and nothing surprising would emerge. Amazingly there is such richness and complexity in the physics of three masses that new discoveries are still being made.

The opening of the door was the paper “Braids in Classical Gravity” by Christopher Moore 1. This did open the door but the paper seems to have missed its audience a bit since one of its results was rediscovered a few years later 2. What these papers discovered was new families of solutions to the gravitational N-body problem not just for N=3 but for much more complex systems. They involve a mathematical assumption that the N bodies follow the same path and were termed choreographies. The gem at the center of all these solutions is a stable figure-8 solution for three equal masses.

These solutions are hypnotizing and as soon as I learned of them I had to get them into the N-body app. The solutions have very special initial conditions and when I first put them in they did not evolve as expected. The version of N-body at the time had a very simple numerical integrator, known as Leapfrog with a fixed time step 3. The choreographies required a more careful approach. Even though N-body was close to being ready to go I could just not give up on these choreographies. So I put in a 4th order Hermite integrator with some adaptive time step capabilities. This was enough to get the simple chain choreographies running. (I then found some more cool choreographies on arXiv, in fact new classes of choreographies. These will need more integrator tuning but look for them in a future update to N-body. In addition the Java gallery by Vanderbei has some solutions I’d like to include in the next release).

I’d love to understand how these solutions are found, but I’m just not there yet. Scholarpedia N-Body Choreographies  is great place to start. A useful expository paper that also finds even more solutions (including some new non-choreography solutions) is Vanderbei’s New orbits for the N-body problem. His paper provides code to calculate choreographies in the specialized language AMPL. This makes the details of the search quite clear, even if you (like me) have never seen AMPL before.

Based on what I’ve absorbed thus far it seems the general idea involves using a path optimization approach using variational principles. This comes from an idea which is in the very core of physics, the notion of finding the path along which a specific quantity (the action) is an extreme value. I learned how to do this when I was in graduate school but I have always been a bit baffled by how the mathematical machinery is searching over paths. There are such a huge number of possible paths! The key mathematics here is the calculus of variations. I’ve gotten as far a buying a copy of the classic staring point Calculus of Variations (aka Gelfand and Fomin) but it is still mocking me quietly from the shelf…

 

Notes:

  1. Phys.Rev.Lett 70:3675-3679. This was in turn anticipated by M. H´enon (1976) “A family of periodic solutions to the planar threebody problem, and their stability.” Celestial Mechanics 13 267–285.
  2.  A. Chenciner and R. Montgomery (2001) “A remarkable periodic solution of the three body problem in the case of equal masses.” Annals of Mathematics 152 881–901.
  3. This is a good place to acknowledge the places I read about N-body evolution codes.  Firstly, Moving Stars Around which is written in a quirky Socratic-dialog style. As part of the learning arc some of the code shown has deliberate errors that are then explained, so be careful. Everything you ever wanted to know about the subject is in the comprehensive Gravitational N-Body Simulations: Tools and Algorithms (Cambridge Monographs on Mathematical Physics). For a more concise summary of N-body codes and their pros and cons the section in Galactic Dynamics: (Second Edition) (Princeton Series in Astrophysics) is very to the point. 

+1 Changes Everything

There are times in physics when what seems like a simple extension to a problem causes the floor to drop out. The motion of three bodies due to gravity is one of those cases. For the two-body solution there is a modest differential equation and a solution that describes the motion as an ellipse, parabola or hyperbola. Adding a third body seems like “no big deal”. You might think that the math gets a bit messier so maybe the solution will not be quite as clean as a simple ellipse – but there should still be some kind of sensible solution. This is one of those cases where our optimism/intuition fails.

Newton’s work was published in the 1680s. It took other great mathematicians to get further. Euler found an unstable 3-body solution in the early 1700s. Lagrange found another unstable solution as well as some stable solutions in special cases for the 3-body problem. His work was published in 1772 – almost a hundred years later. [Lagrange also fundamentally changed physics with his Mechanique Analytique. Today huge parts of modern physics is still rooted in “the Lagrangian”].

Lagrange found a few special cases where three bodies could stay in a stable solution by making simplifications that result in what is called the “restricted three body problem”. In this case one body is “heavy”, the other has moderate mass and the third has no mass. As a concrete example think of the Sun, Jupiter and an asteroid in the same orbit around the sun as Jupiter. As a further simplification, assume Jupiter moves in a perfectly circular orbit. It is then convenient to think in terms of the position of the third massless body in relation to Jupiter i.e. in a co-ordinate frame that rotates around the sun at exactly the same rate as Jupiter. Lagrange found there were five points at which the third body experienced zero net force and could stay fixed in relation to Jupiter in this rotating frame.

This in itself was an interesting result. However these five Lagrange Points would only be really interesting if the massless body actually stayed in this relation to Jupiter even if it was bumped away from this zero force location i.e. if the motion around a Lagrange point is stable. It turns out that two of the five Lagrange points are stable but there’s a catch. The catch is that the Lagrange points L4 and L5 are stable only if Jupiter is not “too heavy” compared to the sun. Exactly how heavy? L4/L5 are stable as long as Jupiter is not more than 3.85 percent the mass of the Sun.

These points are of more than academic interest. This solution motivated a search for asteroids near these locations in Jupiter’s orbit and as the technology of telescopes improved the first such asteroid, Achilles, was found in 1906. There are now thought to be as many of these “Trojan” asteroids as there are asteroids in the belt between Mars and Jupiter. Other planets have such asteroids as well. Recently, more asteroids were found at the L4 point of Mars (http://arxiv.org/pdf/1303.0124.pdf).

What if the test particle is not quite at the stable Lagrange points L4/L5 or does not have quite the right velocity? Then the resulting motion is very interesting. This leads to orbital paths relative to Jupiter that get names such as tadpole or horseshoe orbits, due to their shape – you can get a sense of why these names arise by playing the levels in the Lagrange group in N-body. The image below shows oscillations about the L5 point in the N-body iPad app.

Tadpole

Oscillations about Lagrange points also have real-world implications. For example, the detection of a dust ring around Uranus puzzled astronomers since location of the ring could not be stable. It was proposed that there were “sheparding satellites” 1. Subsequent discovery of two small moons, Ophelia and Cordillia by Voyager lend credence to this idea. 2

Amazingly (and a testament to just how opaque the physics becomes when N=3) another stable three body solution for equal masses was found in the 1990s more than two hundred years after Lagrange. It was a result of considering topologies of paths. We’ll look at this fascinating solution and some of it’s unstable cousins in the choreography group of N-body levels.

One aside: In general relativity the step of going from one to two makes the floor fall out. In the case of Einstein’s general theory of relativity (GR) finding a solution for ONE body is really hard. The spherically symmetric solution was found by Schwarzschild in 1916 about a year after GR was published. It took until the 1960s before the axially symmetric one-body solution was found. There are no two body solutions corresponding to a simple orbit that can be written in terms of modest mathematical functions. This is part because in GR the idea is that gravitational interactions do not operate via a force but by defining the shape of the space they are part of. As bodies move around each other there are ripples in space-time (waves) that need to be taken into account. This makes things very messy. Numerical modelling and simplifications can be employed and things can be learned – but there is not a tidy mathematical description.

Notes:

  1. Goldreich and Tremaine (1979) Nature 277 97-99
  2. See the discussion in 10.5.3 of Solar System Dynamics

Did Gravity have to be an inverse square law?

How fundamental is the fact that gravity is proportional to 1/r^2 versus a more general 1/r^n or something more complicated like f(r) = e^{-r}? Is this solution unique? Does it have any special properties other forces would not have?

This was settled in a publication in to the Royal Society in 1684 by Newton. The idea that the force might be 1/r^2 was not new.  This had been proposed by Hooke, at the time the resident experimenter at the recently formed Royal Society. Hooke postulated an inverse square law but did not prove it, lacking the mathematical insight of Newton. He was always a bit bitter about not being recognized for his contribution in postulating the correct formula and was combative in his exchanges with Newton. It is possible that Newton’s “If I have seen further it is by standing on the shoulders of giants…” was partly to acknowledge Hooke. It’s also possible this was a bit of a dig at Hooke who was a tall, somewhat stooped gentleman. 1

To understand what was really done by Newton it helps to pause and think about what he had to work with. The most precious thing was accurate data – the positions of planets measured carefully over sustained periods of time. This data measured the shape of the planets orbits and their speed at different points of the orbit. Kepler had done a lot of data collection and in 1609 using this information he was able to see patterns that helped him determine relationships in the data. Kepler’s first law was the radical notion that planets moved in ellipses (not circles) and the second was the bodies moved such that the line from the body to the center swept out an equal area in equal times.

So what do you think? Is that in itself enough to determine what force gravity is? Enter Newton.

Newton proved a number of things about generic forces acting on a line between two bodies. In particular he established that for such a force the motion will be in a plane and there will be a constant of motion (angular momentum) that is conserved. He also demonstrated that the shape of the path of one body around the other is the same as each body around their common centre of mass. He showed that for a general force a body will sweep out equal areas in equal times. This fact alone does not distinguish between forces. It is true for any central force.

How about the path of the orbit? Does that uniquely determine the force? Halley (of Halley’s comet fame) had the same question and in 1684 he went to Cambridge to see Newton. He asked what shape a force of 1/r^2 would produce. Newton’s responded that it would be an ellipse and that he had proved it. Halley was amazed and asked for the proof, leaving Newton mumbling “Gee, I know that used to be around here somewhere…”. Newton promised to provide it and several months later he provided it and a lot more in his note The Motion of Gyrating Bodies (De Motu Corporum in Gyrum). This was read to the Royal Society on Dec. 10, 1684. This then sparked Newton into writing the Principia which was then published in 1687. (Halley had to foot the bill for the publication, since the Royal Society had “bet big” on a a previous publication The History of Fishes that turned out to be a flop and consumed all their funds. At one point the society paid Robert Hooke in copies of the book.)

Principia is simply amazing. Any one of the theorems in it would have given Newton a place in history but the book has dozens and dozens of ground breaking results. The way mathematics was written at the time and the terminology used makes a translation of this work pretty tough going, it’s an improvement on trying to read the original latin but not by much. Fortunately there is a marvellous book by Chadrasekhar Newton’s Principia for the Common Reader. It is a masterful book from a physicist whose ability to do difficult calculations is astounding. (His The Mathematical Theory of Black Holes (Oxford Classic Texts in the Physical Sciences)
is perhaps the most densely packed book of equations I ever studied.) The term “common reader” is perhaps a bit generous to the reader, since there are differential equation at about a second year undergrad level.

Newton did WAY more than just figure out the law of gravity. He looked at central forces in general, proved that a spiral would be the path due to a force of 1/r^3. (Depending on the specific type of spiral the force could be slightly different. For example the spiral r = k \theta is of the form f(r)\propto \frac{1}{r^3} + \frac{2 k^2}{r^5} 2)

Newtonian gravity does have one special characteristic. The orbits are closed. After one orbit the object returns to exactly where it started. That is the path is an ellipse and not a “spirograph” type pattern. This characteristic also occurs for a force of f(r)\propto r^2. No other forces have this characteristic, so in a way 1/r^2 is unique.

Note I was careful to say “Newtonian gravity”. Actual gravity is best described by Einstein’s general relativity. One of the drivers for revising our view of gravity was the fact the the orbit of Mercury was not closed. It’s point of closest approach to the sun shifts around very slightly over time (even after all the effects of the other planets are taken into account). As an approximation this can be modelled as a slightly different central force. The careful measurements of the orbits in the solar system (and increasingly the orbits of binaries and exoplanets) are still scrutinized for signs that our understanding of gravity is correct. (See e.g. http://arxiv.org/pdf/gr-qc/0302048.pdf) Another reason for re-thinking gravity is that a central classical force has a a kind of spooky action at a distance built in. If a mass doubles then that effect is immediately felt everywhere in the universe and that would be hard to reconcile with the idea that nothing moves faster than light. In general relativity the changes in the curvature propagate with finite speed so there is a time delay (although you can play games with the equations in some cases using negative energy density to make changes go at faster than the speed of light 3)

Notes:

  1. Hooke was an interesting figure who did a great deal and is probably a bit under-appreciated; see The Forgotten Genius: The Biography Of Robert Hooke 1635-1703
    .
  2. This comes from a book I used to bridge the gap between my engineering undergrad education and “real” physics. Introduction to Classical Mechanics by Atam P. Arya (first edition, 1990). Example 7.1 
  3. I first bumped into this in graduate school in a paper:M. Alcubierre, “The Warp Drive: Hyper-fast Travel within General Relativity”, Class. Quantum Grav., 11, pp.L73-L77, 1994. This is still an active area of investigation; see http://arxiv.org/pdf/1202.5708.pdf

Why are orbits not spirals?

As an engineering undergrad I had to take a course in probability and statistics. This course made me crazy because I could seldom get the right answer to the assigned problems. These were inevitably questions about counting certain groups of events and I would come up with a counting scheme that seemed logical but was often wrong. I could look at the provided answer which counted in some other way and it made sense BUT I could not see why my answer was wrong.

I was reminded of this when I showed N-body to a friend and his first question was “Why don’t the orbits spiral in?”. I kept trying to repeat the answers from my physics education but we kept coming back to “I’m sure your answer is correct, but why is mine wrong?”. To deal with that it finally occurred to me to ask him why he though they were spirals – instead of continuing to blather on about ellipses.

“If I throw a ball (even really far) – it hits the ground. This would happen in a vacuum too. Is this not a spiral path?” was the response. This is actually a good observation – which when linked with other observations about decaying orbits of satellites could easily leave you with the idea that orbits are spirals.

The explanation I came up with was that if the earth was replaced by a single point of the same mass and you threw the ball then the path would be an ellipse and it would come back to where it was. When the real earth is in the way, the ball hits the earth. To a person on the earth it does appear that the object is getting closer to the centre – it is. If you looked really carefully you would discover the path you thought was a spiral was actually a segment of an ellipse.

Why an ellipse? That’s something I used to be able to derive. It turns out to be not as simple a process as you might hope. After all it took someone with the stature of Newton to figure it out. The way we now do that calculation using differential equations is quite different from the geometric approach Newton used. Ironically it draws far more on the Calculus notation of Leibnitz (who independently developed Calculus, something Newton never acknowledged and worked hard to undermine).

As part of developing N-body I tried to dust off that part of my brain and sat down with a pad of paper to see how far I could get. Embarrassingly, I was a little rusty since my PhD is from the mid-nineties and since then I have not done much math. Nevertheless, I peeked only as much as needed and did finally get to the answer. I’ll post about that as a separate topic.

The details do really matter. That’s what science IS. One of my favourite illustrations of this viewpoint is in the movie Insignificancewhich has the characters of Marilyn Monroe and Albert Einstein. Marylyn explains special relativity in essentially correct lay-person terms and says:

Marilyn: I understand the results and the premise. I guess that’s the main thing, huh?

Einstein: That’s nothing.