What is the functional form for r(t) for a solar-sail deorbit into the Sun?

Question

This answer describes a spacecraft in a heliocentric orbit using a solar sail to "deorbit" into the Sun by reflecting sunlight into roughly the prograde direction. In comments it was proposed that while the solar flux and therefore thrust of the sail increases as the craft gets closer to the Sun, the delta-v required for a given change in orbit size increases as well.

For a spacecraft of mass $m$ and perfectly reflecting sail area $A$, what does the plot of orbit radius versus time $r(t)$ look like? Is it a straight line, power law ($p \neq 1$), exponential, logarithmic, or something else?

Since the answer mentions that a reflector angle of 45°, resulting in exactly prograde thrust is not necessarily the optimum, you are welcome to consider either 45° or the optimum angle, fixed or dynamic (though I don't think it would change with distance).

"bonus points" for the derivation of $r(t)$!

Answer 1

This answer assumes that the spacecraft stays in an almost-circular orbit all the time. We have $$ \frac{dE}{dt} = F_\tau v, $$ where $E$ is spacecraft's energy (potential + kinetic), $v$ is its speed, and $F_\tau$ is the tangential component of the solar pressure force. For a circular orbit, $E\propto -\frac1r$ (so $\frac{dE}{dt}\propto \frac{1}{r^2}\frac{dr}{dt}$) and $v\propto\frac{1}{\sqrt{r}}$. As long as the angle between the sail and the direction to the sun stays constant (whether it's optimal or not), $F_\tau\propto -\frac{1}{r^2}$. So we have $$ \frac{1}{r^2}\frac{dr}{dt} \propto -\frac{1}{r^{5/2}}, $$ or $$ \frac{dr}{dt} \propto - \frac{1}{\sqrt{r}}. $$

A solution of a differential equation like this has the form $$ r(t) \propto (T - t)^{2/3}, $$ where $T$ (which is greater than the starting time) is determined by the starting conditions and the proportionality coefficient in the equation.

Answer 2

Looks like @Litho's answer nailed it!

$$ r(t) \propto (T - t)^{2/3} $$

I did a quick simulation based on the Planetary Society's LightSail 2's spec of 5kg and 32 m^2 sail area. I fixed it at a reflector angle of 45° so that the pressure from sunlight results in a radial force outward (momentum of the incident light) plus a a tangential force prograde (momentum of the reflected light).

The acceleration due to the momentum of light either striking or leaving a surface is just

$$\frac{AI}{mc} = \frac{AI_0}{mc} \left(\frac{\text{1 AU}}{r}\right)^2$$

where $I_0$ is the solar constant (intensity at 1 AU) of about 1361 W/m^2. See this answer for more on solar pressure and acceleration by solar sails. Remember to divide the area of the sail by $\sqrt{2}$ to get the projected area at 45°.

I started in a circular orbit at 1 AU and integrated for 15.35 years.

It turns out that $T$ is the arrival time, so in the first plot I just compare $r$, the distance to the Sun in the simulation, to the simple expression:

$$ \text{1 AU} \left(1 - \frac{t}{T}\right)^{2/3} $$

and voilà a perfect fit! The wiggles are due to the fact that I started with a heliocentric circular orbit of 1 AU and a velocity of $\sqrt{GM_{Sun}/1 AU} =$ 29783 m/s with the solar pressure effects at full strength (deceleration, slight outward force reducing gravity, and so the orbit is very slightly elliptical.

Radial accelerations due to the Sun's gravity and incident radiation pressure are given by:

$$-\frac{GM}{r^2} \ \ \text{and} \ \ +\frac{AI_0}{\sqrt{2}mc} \frac{\text{1 AU}^2}{r^2}$$

Numerically at 1 AU they are 5.930E-03 and 2.053E-05 m/s^2 respectively, and because both scale as $1/r^2$ the ratio of the two is fixed and independent of distance. In this case the ratio is about 289:1.

def deriv (X, t):

    r,  v  = X.reshape(2, -1)
    nr, nv = [thing / np.sqrt((thing**2).sum()) for thing in (r, v)] # normals
    rsqAU  = (r**2).sum() / AUsq

    acc_g     = -GMs * r * ((r**2).sum())**-1.5
    acc_solar = (Area/np.sqrt(2.) * I_zero / (m * c) / rsqAU) * (nr - nv) # radially out, and prograde

    return np.hstack((v, acc_g + acc_solar))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
degs, rads        = 180/pi, pi/180

AU     = 1.495978707E+11       # m
AUsq   = AU**2                 # m^2

GMs  = 1.327E+20               # m^3/s^2

km   = 1000.                   # meters
year = 365.2564 * 24. * 3600.  # seconds

# http://www.planetary.org/explore/projects/lightsail-solar-sailing/lightsail-faqs.html
m      = 5.                    # kg
c      = 3E+08                 # m/s
I_zero = 1361.                 # 1361 W/m^2 (at 1 AU)
Area   = 32.                   # m^2

time = np.arange(0, 15.35*year, 1E+05)  # seconds

v0    = np.sqrt(GMs/AU)

X0    = np.array([AU, 0, 0, v0])

print "X0: ", X0

answer, info = ODEint(deriv, X0, time, rtol=1E-10, full_output=True)

print answer.shape

x, v = answer.T.reshape(2, 2, -1)
r    = np.sqrt((x**2).sum(axis=0))
x, y = x

if True:
    plt.figure()
    plt.subplot(2, 1, 1)
    plt.plot(x/km, y/km)
    plt.title('heliocentric de-orbit (km)')
    plt.subplot(2, 1, 2)
    plt.plot(time/year, x/km)
    plt.plot(time/year, y/km)
    plt.plot(time/year, r/km, '-r', linewidth=2)
    plt.title('x, y and r (km) vs time (years)')
    plt.show()

if True:
    T0 = time.max()
    plt.figure()
    plt.plot(time/year, r/km)
    plt.plot(time/year, AU*(1-time/T0)**(2./3)/km)
    plt.title('r and  AU(1-t/15.35)^(2/3) (km) vs time (years)')
    plt.show()