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>

Ambiguous word, ‘of’. Do I mean the torments endured by Prometheus or the torments caused by Prometheus? In this case, both.. though mostly the second.

As I’ve mentioned before, one of the research projects Richard and I are working on (with Carl as distant overseer) is studying Saturn’s F-ring (see also here for lots of cool pictures) which is the speciality of the in-house Cassini group.

Today, for example, at the planetary lunch meeting we had a talk by Nick Cooper on various objects that may or may not have been observed in the F-ring and the difficulties involved in fitting consistent orbits to the various candidates.

For my purposes, where we’re aiming to look at the self-gravity of the ring, the most interesting result Cooper reported was that it to get a consistent fit, they needed to let the mass of Prometheus be a free parameter. They find that Prometheus is considerably lighter than the standard value of Jacobson and French 2004 (1.80*10^17 kg, or thereabouts; I’m too lazy to check). Since the JF04 value is believed to be quite accurate, this suggests that if the fit is probing physics and not noise, then Prometheus’ apparent lower mass could be due to some other material generating an opposing force: such as the F-ring, whose mass they neglected. This handwavy line of reasoning yields a ring mass roughly a third that of Prometheus. However, as Cooper emphasized, there really isn’t enough clear signal yet to make any solid claims. Hopefully improved observations will allow for a better-determined orbital fit, at which time we can have another look at what constraints it might impose on the F-ring.

For the last week or so I’ve been running simulations, with and without self-gravity, in which the moon Prometheus orbits Saturn, interior to a circular ring of material. (The real F-ring is eccentric, which plays an important role in its dynamics, but we wanted to get the easy stuff right first.) I took the orbital parameters for the moon and ring from Carl’s 2005 Nature paper, which provided an explanation for some of the features in the F-ring in terms of encounters with Prometheus.

Unfortunately, the results don’t look much like what I expected. (Warning: large animated gifs. Not for anyone on dial-up!)

(Technical note: units are scaled so that 1 time unit is 1 orbit at 140000 km = 1 distance unit. Old habit, choosing units so that typical quantities are of order unity.)

For example, here’s a look at the simulation in the x-y plane (with the x scale considerably stretched.) Prometheus is the purple-ish square, and the F-ring is in grey. The frame is corotating with Prometheus, so you see it moving in and out on its epicycle (heh! epicycle!). When it’s closest to the ring, it produces a twist, which then shears away (moves off-image due to differences in velocity) because Prometheus is closer to Saturn and is therefore orbiting faster.. so the slower-moving F-ring material seems to fall behind.

Imagine what a fast runner on the inside track sees when he looks at the slower runners outside.. er, if he were kicking them every time they got too close..

Note also that between the time that Prometheus is at the innermost part of its orbit (its pericentre) and at the outermost part (the apocentre), the features tend to change whether they’re pointing inside or out. This explains why the images of the F-ring taken by Cassini show the same thing, and was the subject of the Murray et al. (2005) Nature paper I mentioned earlier.

Another way of looking at the data is to imagine unwrapping the ring, and plotting the radius against the angle. Again, we’re corotating with (the mean motion of) Prometheus, which is at theta=0. Here, where you can see the whole ring at once, it’s even more obvious how much damage Prometheus is doing to the ring, and how the kinks it produces change direction as Prometheus orbits.

And finally, we can imagine replacing Cassini’s rather meagre flight system with a full-fledged impulse drive and fly around the ring, to see what it’d look like from above.. this one’s to true physical scale, and it’s much harder to see the effects than in the other two animations but it’s still obvious which parts of the ring have suffered damage and which haven’t.

Problem #1: the real ring doesn’t show nearly this much distortion. The kinks are a small fraction of the ring width, not several times the ring width. Problem #2: I’d expected that the synodic period (the time it’d take for the ring to go around once and “meet up” with Prometheus at the same angle again) would be about 120 orbits or so, so I expected that after 100 it’d be almost done, but it’s got a long ways to go.

It’s true that in this particular simulation Prometheus is more eccentric than it really is, but in another simulation where it’s less eccentric I still have the problem.. and these animations are more enlightening. Spent a fair bit of time today trying to sort out whether or not I’d done anything wrong, but I couldn’t come up with anything.

Finally I sent Richard an email, explaining that I hadn’t found any problems, and he consulted with Carlos Chavez (PhD student of Carl Murray’s) who’s done a lot of work on these dynamics and was the second author on the Nature letter. He was kind enough to lend us a copy of the relevant pages of his thesis, and what do I find?

Murray et al. 2005 Prometheus semi-major axis: 139700 km
Jacobson and French 2004 Prometheus semi-major axis: 139378 km

The number I used by following Murray et al. was 300 km off, which is 15 times the F-ring core width (20 km)! I’ve only realized what probably happened now: 139378 is pretty close to 139400, and 4 is an easy number to confuse with 7..

(They’re eye-somorphic!)

(Yes, I should be shot for that.)

*sigh*

Now I’ve got to rerun the sims. This is still progress, though, we’re ramping up at a pretty impressive rate considering that neither Richard nor I understand much about ring dynamics. We spent the previous month studying clumping scales, which I should get around to posting sometime.. the animations are neat, and informative, if completely unrealistic.

Anyhow, for those who were curious, this is what I do with my days, and as unlikely as this sounds, this is how I earn a living!

I spent last Friday afternoon with Richard Nelson (my primary supervisor, in the formation group) and Carl Murray (my adjunct supervisor, in the dynamics group). The pictures on their webpages don’t quite do them justice: in real life Richard looks somewhat more professorish, though still not as scholarly as he did in his July appearance on the impressive British astronomy television show The Sky at Night, and the image doesn’t capture Prof. Murray’s mountain-like presence.

This meeting was to talk about Murray’s side of our projects. He’s a major player in the Cassini imaging team which has been taking such gorgeous pictures of the Saturn system, and also a coauthor of a great reference text on planetary dynamics; that book was what got me through my master’s thesis. Learning the secular theory (the analytic model for the solar system in which you average over the orbits and study the longer-term behaviour) from the previous standard (Brouwer & Clemence) was impossible. [Aside: I'm listed in the errata for Murray & Dermott as having spotted a missing 0 in eq 7.47, which I think was my first published science result..]

Prof. Murray has the corner office on the floor below me, with a great view of the city. The room’s also enormous, with a table almost the size of my office and two gigantic flatscreen monitors, the largest I’ve ever seen in a nonindustrial setting.

Various aspects of the Cassini data are assigned to different researchers for the first analysis, and his specialization is the F ring, which is what he gave me a brief tutorial on.

It has some very strange properties, with streamers and striations and the like. Murray’s made a good deal of progress toward understanding some of these features: for example, some of the banding can be explained by the satellite Prometheus encountering the ring and causing a roughly radial perturbation which is then sheared out due to the different orbital speeds at different radii. Two of his graduate students have been working on related models. Some of the features aren’t easily associated with Prometheus, though, and so we need to figure out what’s going on.

He also mentioned cool recent discoveries I can’t talk about yet, some of which he had only made that Friday morning! High-activity, this place, even if unlike Queen’s the building empties out at five-thirty and no one’s here on the weekends..

The PPARC grant which pays my salary claims we’re going to investigate the importance of ring self-gravity and possible moonlet formation in the F-ring. We’ll probably study ring memory first: that is, for what physical parameters can objects clump and resist the shear so that Prometheus doesn’t merely find an (effectively) undisturbed stream at re-encounter, but structures can instead persist?

Should be fun! Details to follow.