Why has Vanguard 1 not decayed significantly since its orbital injection in 1958?

Question

Vanguard 1 became the first satellite to use solar panels as a power source in 1958, when it was launched into an orbit with the following parameters:

654 by 3,969 kilometres (406 mi × 2,466 mi), 134.2 minute elliptical orbit inclined at 34.25 degrees on March 17, 1958

Source, sadly without citation, from Wikipedia.

As you can see, its periapsis is within the tenuous grasp of the exosphere. It was estimated at the time of launch that its orbit would decay over 2000 years, but at the time, it was not known that the exosphere "puffed up" during periods of high solar activity:

[this] caused a significant decrease in its expected lifetime to only about 240 years.

So, if this is correct (and it will only be a rough estimate as orbital decay predictions are notoriously inaccurate), it has already spent nearly a quarter of its total orbital lifetime in orbit.

Yet, the satellite tracking website N2YO (which receives regular TLE updates from current orbital objects) has the satellite currently in a 653km x 3833.5km orbit, inclined 34.3 degrees.

Would we not expect to have a seen a greater decay than what has occurred so far? If these details are correct, its apoapsis has only decayed by 136km since launch, and its periapsis has not shrunk at all, at least not appreciably. This does not appear to match its predicted orbital reentry timeframe of 240 years. What is going on here?

Answer 1

Orbital decay due to aerobraking/drag typically acts to first circularise an orbit, then slowly spiral the satellite into reentry.

This is for one main reason:you get more drag at perigee. Where a force acts on an orbit collinear with the velocity vector the satellites current altitude remains unchanged, all other parts of the orbit are changed. So the natural progression of the satellite in question is:circularisation of orbit;slow reduction of semimajor axis (under practically constant 0 eccentricity);final rentry.

Also to note, due to the exponential nature of the density of atmosphere the it takes alot longer to move from 500km to 400km than it would from 400km to reentry.

If you're interested you can use DAS to calculate the expected deorbit time yourself using JPL's Orbital Debris tool.

Answer 2

In this supplementary answer I've crudely processed all TLEs for Vanguard 1 and plotted the trends.

I used the mean motion (revs/day) to get a period, divided by 0.9975 (estimating from this answer) to undo the effects of $J_2$, then used $a^3=GM (T/2 \pi)^2$ to estimate a semimajor axis, periapsis and apoapsis based on the TLE's eccentricity value. This is not particularly accurate but it's good enough to see the trend and easier than propagating 13,000+ TLEs.

We can see that the periapsis hasn't budged since 1958, and the apoapsis has dropped only about 110 kilometers in over 60 years!

As pointed out in this answer:

1. it's the apoapsis that's dropping because the drag occurs around periapsis, so periapsis stays put and semi-major axis, period, and eccentricity all drop.
2. drops in orbit occur around periods of solar activity maxima and sunspots
3. it definitely has something like hundreds of years before it will reenter the atmosphere.

Source

import numpy as np
import matplotlib.pyplot as plt

# https://space.stackexchange.com/questions/8434/why-has-vanguard-1-not-decayed-significantly-since-its-orbital-injection-in-1958?rq=1

fnames = ('Vanguard 1 TLEs 1950 to 1980.txt', 
          'Vanguard 1 TLEs 1980 to 2000.txt', 
          'Vanguard 1 TLEs 2000 to 2020.txt', 
          'Vanguard 1 TLEs 2020 to 2040.txt')

lines = []
for fname in fnames:
    with open(fname, 'r') as infile:
        lines += infile.readlines()

pairs = zip(lines[0::2], lines[1::2])

goodies = []
for L1, L2 in pairs:
    year = int(L1[18:20])
    if year < 57:
        year += 2000
    else:
        year += 1900
    daynum = float(L1[20:32])
    n = float(L2[52:63])
    ecc = float('.'+ L2[26:33])
    inc = float(L2[8:16])
    year += daynum/365.25

    goodies.append([year, n, ecc, inc])

GM = 3.986E+14
twopi = 2*np.pi

year, n, ecc, inc = np.array(goodies).T
T = 24*3600/n
T /= 0.9975  # fudge factor for J2's contribution https://space.stackexchange.com/a/25906/12102
a = (GM * (T/twopi)**2)**(1./3)
peri, apo = a*(1-ecc), a*(1+ecc)

peri = (peri-6378137.)/1000.
apo  = (apo -6378137.)/1000.

if True:
    plt.figure()
    names = ('mean motion', 'eccentricity', 'inclination',
             'periapsis', 'apoapsis')
    limitz = (10.7, 10.88), (0.18, 0.195), (34.1, 34.4), (560, 760), (3800, 4000)
    for i, (thing, name, (L0, L1)) in enumerate(zip((n, ecc, inc, peri, apo), names, limitz)): #
        plt.subplot(5, 1, i+1)
        plt.plot(year, thing)
        plt.ylim(L0, L1)
        plt.title(name, fontsize=14)
    plt.xlim(1955, 2025)
    plt.show()

Answer 3

This additional answer is just to show the effect of the solar activity on the decay rate.

The graphs were obtained from 15279 TLEs downloaded from https://celestrak.org/NORAD/archives/request.php processed with the CSpOC's SGP4 library freely downloadable from www.space-track.org.

The following graph shows the mean radius vector and the mean air density (see at the end of the post for the meaning of the terms):

We see that the decay rate increases dramatically at peaks in air density.

The next graph shows that the apoapsis has decreased by around 125 km since launch:

The linear regression for the perigee shows an average decay rate of about 140 mm/day (about 50 m/year):

The last graph shows the minimum, mean, maximum radius vector and the actual eccentricity:

---

Definitions

T: orbital period.

1-orbit mean radius vector: numerically integrated radius vector against the eccentric anomaly (it’s the semi-major axis). The integration starts from –T/2 and ends to T/2 before and after the TLE epoch (1 orbit). Not to be confused with the osculating semi-major axis.

1-orbit minimum radius vector: the smallest radius vector from –T/2 to T/2 before and after the TLE epoch (1 orbit). Not to be confused with the osculating perigee.

1-orbit maximum radius vector: the biggest radius vector from –T/2 to T/2 before and after the TLE epoch (1 orbit). Not to be confused with the osculating apogee.

actual eccentricity = (Ra - Rp) / (Ra + Rp), where Ra is the 1-orbit maximum radius vector and Rp is the 1-orbit minimum radius vector. Not to be confused with the osculating eccentricity.

1-orbit mean air density: numerically integrated air density against the time divided by T. The integration starts from –T/2 and ends to T/2 before and after the TLE epoch. The air density is calculated at the satellite position with the NRLMSISE-00 atmosphere model updated with the solar and geomagnetic indices in the file www.celestrak.com/spacedata/SW-All.txt.

For the graphs that show the radius vector, the vertical axis is scaled to a sphere with a radius of 6371 km, just to show an approximate altitude.

Answer 4

Here's a simulation to answer your question:

Would we not expect to have a seen a greater decay than what has occurred so far?

Answer: no.

The simulation includes the Newtonian and the relativistic accelerations of all the planets, Sun and Moon.
The Earth's gravity field is modeled with the SGG-UGM-1 gravity model (computed using EGM2008 derived gravity anomaly and GOCE observation data) truncated to the degree and order 15 (to save the running time, while retaining good accuracy when compared to the full model).
For the calculation of the air density, I use the NRLMSISE-00 model along with an updated data file for the solar and geomagnetic indices: www.celestrak.com/spacedata/SW-All.txt.

The first step involves determining the best ballistic coefficient to minimize a particular simulation parameter. After 27 minutes, the program finds a ballistic coefficient of about $16\,kg/m^2$ (it’s not fixed, because the drag coefficient varies with the air composition).

Now the simulation can be started:
1) with one TLE and the CSpOC’s SGP4 propagator calculate the initial state (position and velocity) of the satellite for the TLE epoch;
2) propagate that initial state with a specially crafted propagator (my propagator is based on the 8(5,3) Dormand-Prince integrator);
3) stop the propagation at the last available TLE (20024.79812439).

Here’s the result obtained with the TLE 92001.73795324 (I start the simulation from 1992 because before that year the TLEs are a bit sparse):

The graph shows both the "actual" (TLE+SGP4) and the expected (integrated) mean radius vector (or semi-major axis) from 1992 to 2020-01-24.
We see that the integrated trajectory follows very closely the one obtained from TLE+SGP4 (the amplitude of the two plots are different, but on average the decay rate is the same).