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

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)$!

• Well-- if a craft like this is ever designed it absolutely must be called ICARUS. Aug 25, 2018 at 4:16
• Probably solution is not an analytical function. Aug 25, 2018 at 6:10
• @Heopps for a circular orbit or a slow spiral, I'll bet it is! Let's see what happens...
– uhoh
Aug 25, 2018 at 7:17

## 2 Answers

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.

• okay it looks like the onus is on me to check it, I'll give it a go in the morning. Thanks!
– uhoh
Aug 27, 2018 at 19:16
• @uhoh Sorry, I made a mistake. The power should be 2/3, not 2. Aug 28, 2018 at 2:49
• That looks better! It also appears to have the same shape as my simulation (I'll post something soon). Also, I think $T$ is "greater than (or equal to) the finish time" as well. Exponentiation to fractional powers returns complex numbers for negative arguments, for example $-1^{2/3}=-1/2+\sqrt{3/4}j$ So probably $T$ may end up being just the arrival time itself, since $r(t)=0$.
– uhoh
Aug 28, 2018 at 4:00

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()

• Way to check the answer, or answer the check, maybe. I appreciate you putting in the code, so I can play around with it, too. May 8, 2020 at 0:34
• I already like this community, because people share useful things. : ) Thanks! May 8, 2020 at 0:39