I told NASA an asteroid would hit Earth

I still think about this sentence that I wrote in 2016:

Through long-term integration, it was determined that 2003 LS3 will collide with the earth within 5,000 years.

That was the conclusion of a real scientific paper I wrote in 2016. I published the paper with Alex Lathem and Ana Tudor¹ at the Yale Summer Program in Astrophysics, a months-long research program of observatory nights and observations. We compressed our findings into five pages. Here’s the cover sheet.

A five-page paper in the two-column IEEE template. The abstract describes “an accurate model of its future orbit.”

We had five nights with the telescopes. Over the course of four weeks, we determined that the asteroid was going to hit the Earth. After drawing our conclusion, we checked our orbit against NASA JPL’s official solution. When we saw that we disagreed, we concluded that NASA might be wrong.

I still have all of it: the raw frames, the astronomical observations, and the final paper. At the end of the recap, I’ll also re-run the next decade of observations on the asteroid, and we’ll see who is actually right.

Five nights of observations

YSPA is a summer program at Yale’s Leitner Observatory where they bring in students to research real problems in astronomy.² I was assigned 2003 LS3, a near-Earth asteroid.³ Near-Earth is a technical category of asteroid whose perihelion is less than 1.3 astronomical units. NASA tracks more than 14,000 near-Earth asteroids to make sure none of them will become an issue for humanity.⁴

July 9th: landing in LaGuardia on the way to Yale. This is the last photo I took all summer. The next one is dated September 11th, nine weeks later.

Over the course of four weeks, I got five clear nights and 25 usable observations. For nights where the sky was cloudy in Connecticut, we looked up observations in other parts of the world.

Our observing routine was roughly systematic. We cooled the CCD to −5°C, pointed it at a bright star, usually Altair, and nudged the aim until it sat in the center of the field of view, so that the telescope and the sky agreed on coordinates. Then we’d move the telescope to the part of the sky where JPL’s ephemeris told us the asteroid should be and start taking 60-second exposures: about 20 luminance frames, then some red ones, then a bit more luminance at the end.

We did that from right at dark until the sky went completely gray in the morning. We’d reduce the data and do it again. So sleep deprivation became dire.

Asleep on the bus to the observatory.

A few nights of that, and we got some good photos. You can see it there. The asteroid is basically a star that’s not in any catalog.

Leitner 16-inch, July 16, 2016. The same stars, 23 minutes apart. You can see the circled dot has moved 22 arcseconds, about 16 pixels, and nothing else has moved. Magnitude 16.3, detected at 17σ over the noise.

It’s crazy that that little dot is a gigantic rock, 420 to 940 meters across. It’s moving at magnitude 16, and it’s thousands of times too dim for your naked eye to see.

We ran each frame through Astrometry.net to pattern-match the stars in the image against the catalog and figure out exactly which piece of sky we had photographed. From there, we pulled centroids in DS9 and turned every observation into a right ascension and a declination.

All of our observations resulted in 25 data points across what’s called the celestial sphere.

Some 1801 math

To take the dot observations and figure out the asteroid’s orbit, we used the Method of Gauss. Gauss invented this technique in 1801 to recover Ceres after it had disappeared behind the Sun. It worked so well that astronomers found Ceres again right where he said it would eventually appear.

The method takes three observations and adds in Newton’s law of gravitation:

which says that any object with position $\mathbf{r}$ and relative velocity $\dot{\mathbf{r}}$ moves on an ellipse with the Sun at one focus. If you feed that and the three data points into 13 steps of iteration, it returns the position and velocity vectors at the middle observation. Ours came out to

If you then rotate those into the plane of the solar system and use some other vector identities, you can get the six classical orbital elements: the size, shape, and tilt of the ellipse that is the asteroid’s orbit.

But Gauss’s method just gives you the preliminary orbit determination. It gives you three points that make up the orbit, but the points have noise in them, so you then have to take the preliminary orbit and fit it to everything you observed.

Using a genetic algorithm

You could do this part with least squares and a differential correction, but we used a genetic algorithm.

It’s basically natural selection for orbits:

1. You spawn 100 candidate asteroids, each with a position and velocity randomly scattered around the preliminary solution.
2. You propagate every candidate forward with a fourth-order Runge-Kutta integrator. From each candidate, you get where it thinks the asteroid should have appeared at each of the 25 timestamps, and then you score that guess by the root mean square error against what was actually measured.
3. You kill 99 out of 100 candidates, and the one lone survivor becomes the parent of the next generation.

As you do this iteratively, the random spread shrinks by a factor of 0.9.

The scoring function, as printed in footnote 1 of the paper, looks like this:

It’s a root mean square error with no square. Thus, residuals cancel out instead of accumulating, and a symmetric spread of errors gets scored to zero. That’s a typo. The intended function, which the code presumably computed, was

Survival of the best fit. The final orbit matched our observations with a root mean square residual of 0.0108.

Simulating the collision

Here’s where things get interesting. We plugged the best-fit vectors into Rebound, an n-body integrator. We built the solar system with the Sun, Earth, Jupiter, and Saturn, then we added the asteroid and integrated forward until either a collision happened or a million years passed, whichever happened first. We logged the distance to Earth at every close approach.

Each run took a little over an hour. We ran three and seeded each from a slightly different orbit within the error bars.

All three simulations collided with Earth.

Since Earth’s radius is 6,371 km, every single one of those close approaches is a point inside of the planet.

Thus our conclusion:

Through long-term integration, it was determined that 2003 LS3 will collide with the earth within 5,000 years.

Figure 2, printed in our paper. The Sun is at one of the foci, and the Earth is on the inner circle.

Overruling NASA

Next, we compared our results against accepted values. JPL Horizons publishes an orbit for 2003 LS3 built from years of observations.

If you look closely, you can see that we got the orientation of the orbit nearly perfect: the node to within a hundredth of a degree, and the inclination to within a degree. What we missed, by the way, were the numbers that determine how big the ellipse is and how far it stretches: half an astronomical unit on the semimajor axis. The difference was equal to nearly 200 times the distance from the Earth to the Moon, and the precision was meant to be one part in 3,000.

The models were clearly very different, and our error bars were even smaller than theirs, so we submitted the discrepancy to NASA:

Thus, we stay confident in our own generated orbit model, and submit the possibility that there is a significant enough difference between our model and the JPL model that the JPL Horizons information is not quite completely accurate.

We did not hear back from NASA.

What the error analysis says

In hindsight, there are some other reasons why we would have observed these discrepancies.

We only watched 1% of an orbit. 2003 LS3 takes 4.3 years to go around the Sun. Our first observation is JD 2457587.834, and our last is JD 2457604.699. This represents an arc of 16.9 days, which is only 1.07% of a period and only 14° of true anomaly. From that, we attempted to extrapolate 5,000 years. It’s sort of like watching a golf ball for the first inch of flight and trying to determine which blade of grass it will land on.

Here’s what it looks like when you draw both orbits to scale:

The two ellipses nearly overlap near perihelion. Our perihelion distance was 1.2765 astronomical units, and JPL’s was 1.2534 astronomical units. That’s only a 2% difference, but on the far side, which is the half of the orbit that we did not measure, the errors accumulate. Our aphelion is 5.04 astronomical units, while JPL’s is 4.05. That’s a full astronomical unit of difference, very far from the observations.

Our error bars measured the algorithm, but not the sky. This is probably the biggest one. A genetic algorithm with a shrinking search radius will always converge. The ±0.001 on our elements came from the optimizer’s final cluster, but not from the uncertainty of the measurements. The paper admits that the observations scattered by ±39 arcseconds, so we did have precision, but not strong accuracy.

Three samples can’t make a full distribution. Modern planetary defense analyses run thousands of virtual asteroids sampled from the uncertainty. They report an impact probability only. Since we only ran three because each run took a whole hour, we had to admit that we couldn’t state a probability.

Re-evaluating the conclusions

When I was writing this blog, I decided to go back and actually re-run all of the integrations.

Firstly, the orbit we published cannot hit Earth. Its perihelion, the closest the ellipse ever comes to the Sun, is 1.2765 astronomical units, but Earth’s aphelion, the furthest it ever gets from the Sun, is 1.017 astronomical units. The two orbits are a quarter of an astronomical unit apart at their nearest. I just installed Rebound again and built the same four-body solar system. I integrated the published orbit 5,000 years forward, and the closest approach was 0.271 astronomical units, about 105 lunar distances.

So the paper’s conclusion was actually at odds with the paper’s orbit.

Let’s try to figure out where the error was.

Going back to the equation from above:

This is the vis-viva equation with $\mu$, the Sun’s gravitational parameter, set to 1. It holds true in Gaussian units, where the clock is rescaled so the arithmetic is clean. One Gaussian day is about 58.1 ordinary days, so the velocity comes out in astronomical units per Gaussian day. If you plug in our calculated $\mathbf{r}$ and $\dot{\mathbf{r}}$, you get $a = 3.1583$, $e = 0.5958$, $i = 10.305°$. These are basically our published elements. The vectors, the formula, and the elements table all agree with each other in Gaussian units.

But above in the paper, it states the velocity vector is in astronomical units per year.

This is not true. The conversion between the two is

We had actually mislabeled our velocity vector by a factor of $2\pi$. Rebound was set up to pull the planets from JPL’s ephemeris, which gives values in astronomical units, years, and solar masses. We gave this simulation a number in astronomical units per year, but it was really astronomical units per Gaussian day. Because of this, the asteroid’s initial velocity was simulated at one-sixth of its actual speed.

And an asteroid moving at one-sixth of orbital velocity moves quite differently:

The second row shows the asteroid, very close to the Sun, on a seven-month orbit that crosses Earth’s path twice per rotation. If you give it 4,000 years, which is about 7,000 asteroid rotations, it has many, many chances to collide. When I integrate that, the closest approach drops from 0.271 astronomical units to 0.0079 astronomical units. That’s 34 times nearer, or three lunar distances.

To be clear, I cannot prove that this is what happened, because the original Rebound script is gone. When I ran the mislabeled orbit myself, it came close to Earth but never actually hit in any of my runs. What I can prove is that the orbit we published cannot produce a collision, and the units of our velocity vector are off by exactly $2\pi$. Misreading them is the best explanation I can come up with for how we determined that 2003 LS3 would have hit Earth.

It’s actually pretty interesting that our conclusion was a result of mislabeled units instead of bad observations.

What ended up happening to the asteroid?

2003 LS3 is now properly numbered. It’s called 468448 and has been promoted out of its provisional designation since the orbit was eventually pinned down. This was done over 621 observations and 8,264 days. Our 25 points over 17 days are in that data that was used.

It’s what’s known as an Amor, an asteroid that approaches Earth’s orbit but does not cross it. JPL classifies it as not potentially hazardous. Its orbit uncertainty parameter is 0, the best possible value. The nearest the two orbits ever come is about 95 times the distance to the Moon.

What about the correction I submitted to NASA? We can now measure whether JPL Horizons was accurate or not. In the ten years since I filed, with 621 observations in addition to my 25, JPL’s semimajor axis has moved by 0.0022 astronomical units. NASA did make a correction, but it was only about 0.5% of the correction we submitted.

Conclusions

The paper we originally wrote actually hedges more than the headline sentence suggests. It does say that a collision “is a possibility” and that the model “warrants further investigation and data collection for verification.” The fact that we observed 2003 LS3 at all was definitely the right instinct. It is genuinely a near-Earth asteroid, and watching these things is very important to do. This is the basis of the paper’s introduction, which certainly is still accurate, even if the conclusion was not.

We also, to our credit, did every hard part correctly. The astrometry was good. Our positions are within a few arcseconds of JPL’s ephemeris on the same nights. Gauss’s method definitely worked. The genetic algorithm converged. The orientation of the orbit, which is the plane it sits in and the direction it points, is correct to a hundredth of a degree using just a telescope on the roof of the observatory in New Haven.

If it weren’t for the “astronomical units per year,” we would have gotten an accurate conclusion, even if not as salacious.