One More Epicycle

I’ve mentioned several times recently that I’ve been working on techniques to allow me to study the formation of hot Neptune systems without the simulations taking fifty times longer than they used to.  (See, e.g. here; I don’t think people looking for “astral projection” found what they wanted!)

That factor of fifty’s not an exaggeration, by the way: the orbital period at 0.05 AU is ~0.01 yrs.  Even if I only take ten steps an orbit, which is pushing things a little, that means I have to take a timestep of 0.001 yrs.. for comparison, my usual timestep for terrestrial planet formation is 0.05 yr.  For the first part of the evolution, which lasts ~5 Myr, this means I’d need five billion steps.  Assuming each step took 10 milliseconds, this simulation would take 50 million seconds, or 578 days.. completely impractical.   I want something that’ll take a few weeks at most.

Well, I have some good news to report!  My most recent attempt is behaving quite nicely, despite the small violations of momentum conservation.  (These come about because to compute the forces between the inner and outer zones I use the “current” inner particle positions — they’re updated on a short timescale — but “expired” outer particle positions, and the sampling frequency isn’t as short as it should be to resolve the inner potential completely.)  The inner objects are treated using the same code as before, so I know the close encounters are being treated correctly; likewise the outer objects; only the transition from zone to zone (which turns out to be quite smooth) and the force coupling between the inner and outer zones (which seems okay) are potential problems.  Had to play some crazy pointer games to avoid modifying what I’d already written, but now isn’t the time to go breaking things which work.

Preliminary results are very encouraging.  For example, in a test case with two 18 Earth-mass objects undergoing strong migration in a heavy (five times standard) gas disc, it’s basically impossible to determine from the semimajor axis and eccentricity evolution which code I used (the standard approach or my new version) or where the transition point is.

Tomorrow I’ll do as I promised and plot graph after graph to convince Richard that the code is working, and then hopefully we can get the low-resolution runs started soon, if not by the end of this week then the beginning of next..

Well, here I am in the office early in the morning. Why, you ask?

Had trouble getting to sleep last night, and then woke up at about five. Tossed for an hour and change before giving up and coming in. I have renewed sympathy for people with major sleep problems (hi Dad!).. even this one night of restlessness was frustrating. Not sure why I couldn’t sleep, but a research problem I’ll describe in a minute was proving very distracting.

(All-time best remains one evening when I was out with a gorgeous young woman at a family-run hole-in-the-wall restaurant in Kingston. The food was very tasty, but the coffee was do-it-yourself and I accidentally spilled an entire enormous chunk of the grounds into my cup; it took up about a third of the volume. Rather than writing off the cup as a failure, I forced myself to drink it, and only put in a little bit of extra water.

I spent the night being exhausted beyond belief and yet wired like I’ve never been in my entire life. If I’d had sleeping pills available I’d probably be dead now, because I’d have swallowed the bottle whole and not stopped to think about appropriate dosages.)

To understand my problem, you need some background. Most of my work in planet formation involves simulating the formation process using N-body integrations. The code accepts a list of objects (N objects, or bodies, hence the name), which specifies their masses, their positions, and their velocities. Then the code computes the trajectories of the objects as the result of their mutual interactions: you compute the forces, which gives you the accelerations, and so you change the velocities a bit; then you have an updated velocity, so you change the position a bit; etc. In practice there are all sorts of clever things people have dreamed up to do this efficiently and robustly, but at heart you’re just approximating the continuous path by dividing it up into small steps and advancing each part. (The relationships between the posititions, the velocities, and the accelerations involve differential equations, so solving them is an “integration”: you add up the little changes to get the final answer.)

The particular kind of integrator that I have expertise with is called “symplectic”, which is a subcategory of geometric integrators. Basically, geometric integrators build information about the structure of the equations right into the method you use to solve them, in the hopes of improving the stability and the accuracy of the integrations.

It’s sometimes easy to get the structure right: for example, if you’re modelling a spherically symmetric system, then you only have one distance variable: the radius. You *could* use Cartesian x, y, and z coordinates, but (1) it’s more work, and (2) since the method you’re using isn’t enforcing spherical symmetry, you’d expect that little errors would build up over time and eventually your spherical shells wouldn’t be spherical any more.. which tends to be what happens.

In the case of spherical symmetry, choosing the right coordinates is so obvious that we don’t even notice we’re doing it. Not every geometric structure is so straightforward to handle, and appropriate coordinate choice (though important) isn’t the only way to get your equation solver to respect the properties of the system.. which is where symplectic integrators come in.

“Symplectic” comes from the Greek for “structure”, if I remember right, and the particular kind of structure that the integration is trying to preserve here is the relationship between the positions and the velocities that Newton’s laws give us, which turns out to be roughly equivalent to preserving areas. By designing our integrators in a way which respects this, we can get much better energy conservation (the total energy stays what it should), we can take much larger timesteps (speeding up the calculation by orders of magnitude), and generally do a better job of representing the long-term dynamics, the overall shape of the evolution.

For the planetary work that I do, there’s even more structure we can build in: the Sun is by far the dominant mass in the system, so to first order you can neglect all the other planets. My hero Kepler showed that the planets orbit the Sun on near-ellipses, so if we write the equations as a Kepler orbit around the Sun — which we can solve very easily — plus a little perturbation due to the other planets, then our integration method by construction keeps us close to where we want to be. This works amazingly well, but there are a number of complications I won’t get into here which mean that getting a symplectic integrator to be tuned to the Kepler problem was very challenging. It wasn’t done until 1991 by Wisdom & Holman on one side of the Pacific and Kinoshita, Yoshida, and Nakai on the other.

There’s actually a bit of a controversy over who should be credited with the result. WH were in a way closer to what we do now, but they derived it using a method which wasn’t easily generalizable. KYN were in a way further off, but used methods which have become standard because of their simplicity. (I think in one of Tremaine’s papers he wrote something to the effect that “we follow WH’s algorithm but KYN’s analysis”, which says it pretty well.) For a while in the 90s the Kepler trick I mentioned above was known as “mixed-variable symplectic” — the mixed variables being Cartesian (for the interplanet forces) and Keplerian (for the orbits around the Sun) — which I didn’t like much because it’s not obvious what the mixing is unless you already know. These days we call it the Wisdom-Holman mapping, which I don’t like because the name doesn’t tell you anything about what you’re doing, only one of the teams who dreamed it up. I’d have called it the standard Kepler mapping, or Jacobi-Kepler mapping if we want to be technical. (Jacobi here referring to the coordinates you use.)

Recently, i.e. post-1998, a descendant of the W-H map due to my old supervisor Martin Duncan, Hal Levison, and Man Hoi Lee has become more popular because it allows you to handle close encounters between the objects, which the original map couldn’t. And, again, there’s a naming problem. Duncan et al. called their coordinate system “democratic heliocentric” (where “democratic” refers to the fact that all objects apart from the Sun are treated equally in the system); John Chambers called the same coordinate system “mixed-centre” (because the positions are measured with respect to the centre of the Sun, but the velocities are measured with respect to the barycentre); and Jack Wisdom has called it “canonical heliocentric” (because the positions are with respect to the Sun and the velocities are what they have to be to make the coordinate system canonical, which means appropriate for the underlying equations). Who knows what they’ll be calling it next year?

This Duncan et al. approach is the method I use in the code I wrote for my thesis. Unfortunately, though it’s very good at handling the case where bodies get very close to one another, it’s not so good for handling a wide range of orbital radii. The inner planets, closer to the Sun, move much faster than the outer ones, and you need to choose one base timestep to integrate all the orbits on. This means that if I have an object close to the Sun then I have to choose a timestep small enough to resolve its orbit. But this can be hundreds of times smaller than I need to handle objects further out, and if an object is close enough then it can slow my simulation down by a depressing factor.

As I’ve mentioned before, one of the projects we’re working on involves “hot” Neptunes, where the “heat” here isn’t thermal but dynamical. They’re not hot in the don’t-touch-them sense, they’re “hot” because they’re moving so fast, because they’re very close to their parent stars. For example, in the system HD69380, there are three planets each with masses less than 1 Neptune mass all inside 1 AU (the Earth’s distance from the Sun). The closest planet is 0.0785 AU, or less than a tenth the Earth’s distance from the Sun, which is just plain silly. (“Who ordered that?”)

Right before Christmas we did some quick runs studying what initial conditions we’d need to reproduce this system, and found some promising results, but only by neglecting the detailed evolution of the inner part of the system.. which is exactly what we care about.

So I’m experimenting with ways to allow more distant objects to evolve on their own timesteps. There are a number of approaches in the literature, but it’s difficult to have multiple timescales and simultaneously (1) stay symplectic, (2) allow close encounters, and (3) stay fast. I can’t give up the last two. If I give up symplecticity it’s trivial, but if I do that then I have to take far more steps an orbit than I’d like, which causes speed problems. It’s probably still a net improvement, but not by enough.

Once — just once — I’d like to study a problem which doesn’t involve endless amounts of brainwork, one where I could simply use what we’ve already done and not have to dream up new methods.

*sigh*

I’ve been brushing up on the literature about structures formed in Saturn’s rings, and pretty much everything I know from terrestrial planet formation is wrong. Robin Canup has a 1995 paper (co-authored with Larry Esposito) on accretion in the Roche zone, and it’s a crazy regime. The tidal forces due to Saturn suppress the merger of like-size bodies and runaway accretion, although self-gravity may help counteract this, and many of the Japanese teams (e.g. Daisaka and Ida, and Tanaka) have been studying the dynamics of self-gravitation-produced clumps.

One thing I’m coming to believe is that the issue of self-gravity is unrelated (at first order, anyhow) to the problem I mentioned earlier. Recall that there are these sharp features in the F-ring, and Murray et al. explain them using a simple model in which Prometheus pulls ring particles towards it when it’s nearest the ring, and then the particles oscillate back and forth, which explains why the features change direction every half-period. (Although the Cassini guys around here don’t seem to agree, it seems to me there are some F-ring features which point the wrong way for their model. I’m going to have to crunch some numbers before I declare war on the point, though, in my animations there are moments when you stop it where the differential response leads to what look like phase errors, though everything’s behaving exactly as it should.)

However, there’s a sense in which the problem they’ve solved isn’t the interesting one. As I’ve mentioned before, their model works well on a pristine ring, but doesn’t explain why the ring (after a full synodic period, when Prometheus laps the ring, like the inner car on a racetrack) doesn’t show signs of disruption on earlier passages. This is the point at which the locals tend to wave their hands and say “self-gravity” or “collisions”. I might believe the latter, but not the former: assuming the F-ring core mass is comparable to Prometheus, then the mass in the region of a feature is down by a factor of 100 or so (roughly estimating from one of my last runs). There’s no reason to believe that’s going to wash away the spikes, with or without the clumps that Daisaka et al. find.

Collisions might do it, but it’s tough to see what effect self-gravity should have, and so far in my simulations it hasn’t had any.. at least for reasonable disc masses. (You put a Jupiter mass in the F-ring and interesting things happen, but I think it’s fair to say that’s ruled out by the observations..)

To prove this, of course, I’m going to have to modify my code to handle the collisional evolution of the rings. That’s actually not going to be as difficult as it sounds: the hard parts (the close encounter resolution and the treecode) are already done and working. In principle, all I have to do is take the near-neighbour list produced by the treecode and feed that potential-encounter list into my encounter routine. In practice, I need to be careful to make sure that (1) all forces on encountering particles are being computed directly (with no tree approximations) before they’re sent off, and (2) the number of encounter candidates is kept close to the number of real encounters. My first attempt, for merely ~0.1 M particles, produced >5M encounter pairs, but almost none of those were actually going to get anywhere near each other, much less within 10 Hill radii (where the Hill radius is the size of the region in which your gravity is more important than the central body’s.)

I think I may spend time looking at the non-self-gravitating models for a while.. one thing I’ve been wondering about is whether there are any interesting dynamics with the small ring moonlets analogous to the behaviour of Earth’s quasi-satellites. That is, for the ring-crossers, how much do their orbits change qualitatively due to interactions with Prometheus, the ring, and Pandora?

Some days you feel full of energy, like you can do anything. I thought today was going to be one of those days: it was bright and blue when I walked to the office, with a cool spring-though-it’s-autumn breeze.

Since then it’s been a series of minor disasters.

Remember the Prometheus/F-ring work I’ve been doing?

I had a look at data from some sims I ran earlier this week. First it turned out I managed to get Prometheus’ eccentricity wrong: instead of being 0.002 and 0.003, it was 0.000 and 0.000. Much less interesting, although I did notice something odd that I’ll have to talk to Richard about later. Then I thought that I’d screwed things up even worse, because the self-gravitating rings had zero mass (which made the force calculations quick but utterly unnecessary).. so I deleted the useless duplicate outputs. Always a bad idea, deleting things.

Turned out that I’d used a symbolic link to point to the ring data, and the link pointed to the right place, so I’d erased data which was actually informative. Unfortunately what it’s informing me is that it doesn’t look like there’s a lot of interesting self-gravitational effects when the ring mass is comparable to Prometheus, because the mass in the region which is being strongly influenced by Prometheus is much lower. This isn’t what I was hoping for.

So I fired off yet another suite, but ran a quick low-N test-particle case to see what to expect. As predicted, after a synodic period (120 years or so), when Prometheus has managed to lap the edge of the ring, the kicks it gives the ring don’t quite match up with the kicks it gave on the earlier go-round, so the ring structure gets increasingly noisy. (Units as before.)

I don’t understand how the real ring core manages to avoid this problem and produce such clean signals. The in-house explanation works on a pristine ring over one synodic period very nicely, but I don’t see how it works on the second pass. I’m missing something obvious, I think, and I’ll ask Carl at tomorrow’s meeting.

And then a former coauthor (I have trouble saying “collaborator”.. makes it sound like you’re working with the Nazis in Vichy France or something) dropped a note asking if he could get a look at data from an earlier paper of ours for a talk he’s giving.. a talk which is next week, which means he really needs the data right now. So now to my list of things I’m behind on I have to add digging through old files and making them fit for public viewing.

Plus I ate too many of the yogurt-coated peanuts and raisins I picked up as a snack, and they’re implanting in my stomach like alien vines would on some cheap Outer Limits knockoff.

Woe, woe is me!

</self-pity>