Low-thrust spiraling to escape, is the flight path angle (gamma) at C3=0 always 39 degrees?

Question

@MarkAdler's comment led me to ask Why would a slow spiral from a C3 of zero take about 2.4 times as much ΔV as an impulsive maneuver? which resulted in this tidy and efficient @MarkAdler answer which points to another thoughtful answer about slowly spiraling out from a circular orbit to escape in the limit of very weak prograde propulsion, which (at first counterintuitively) slows you down while raising your orbit.

Below that answer is yet another easter-egg-like comment gem.

Always aligned with the velocity vector. That is the most efficient use of the thrust in order to increase the specific energy. The final γ is 31°.

In this answer @Julio provides a diagram showing definitions for both $\beta$ and $\gamma$ angles which measure the angle between the instantaneous velocity vector and the radial and the tangential directions, respectively.

In this answer @TomSpilker elaborates on these angles, and in this answer I give a little more information on how to calculate them.

Now I've gone back and calculated an outwardly spiraling orbit under low thrust using various conditions. Invariably I end up with a final angle $\gamma$ (gamma) of about 39 degrees when checking the moment where C3 = 0, not 31 degrees.

I'm doing a unitless calculation where GM = 1.0 and the period of an r=1.0 orbit is $2 \pi$. In this case C3 = v^2 - 2/r.

note: For this calculation, thrust is always in the same direction as velocity $\mathbf{v}$, rather than in the tangential direction (perpendicular to $\mathbf{r}$) and I'm beginning to wonder if herein lies the difference between 31 and 39 degrees.

Question: Is this ~39 degrees at C3=0 correct, and is it expected to be invariant like this?

      starting conditions                              at C3 = 0
-------------------------------     ------------------------------------------
rstart  vstart    C3    thrust      time   delta-v  gamma(deg)    r       v        C3
 1.0     1.0    -1.0    0.01        74.5    0.745     38.9       8.78    0.477   0.000
 1.0     1.0    -1.0    0.001       856.3   0.856     39.2      27.80    0.268   0.000
 1.0     1.0    -1.0    0.0001      9192.1  0.919     39.2      87.91    0.151   0.000
 4.0     0.5    -0.25   0.0001      4192.1  0.419     39.1      87.90    0.151   0.000

def deriv(X, t):
    x, v  = X.reshape(2, -1)
    vnorm = v / np.sqrt((v**2).sum())
    acc_g = -x * ((x**2).sum())**-1.5
    acc_t = thrust * vnorm
    return np.hstack((v, acc_g + acc_t))

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

T    = 16 * twopi        # or 160, 1600

ntot = 20001
time = np.linspace(0, T, ntot)

rstart = 1.0             # or 4.0
vstart = np.sqrt(1./rstart)

X0     = np.array([rstart, 0, 0, vstart])

thrust = 0.01            # or 0.001, 0.0001

answer, info = ODEint(deriv, X0, time, full_output= True)

xx, vv = answer.T.reshape(2, 2, -1)

r   = np.sqrt((xx**2).sum(axis=0))
vsq =         (vv**2).sum(axis=0)
C3 = vsq - 2./r

nstop = np.argmax(C3>0) + 1

dotted     = (xx*vv).sum(axis=0)
rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
gamma      = np.arcsin(dotted/(rabs*vabs))   # Per Tom Spilker's answer Eq. 3

print 'C3 min, max: ', C3.min(), C3.max()
print 'nstop, ntot: ', nstop, ntot
if True:
    plt.figure()

    plt.subplot(1, 2, 1)
    plt.plot(xx[0, :nstop], xx[1, :nstop])

    plt.subplot(3, 2, 2)
    plt.plot(time[:nstop], r[:nstop])
    plt.ylabel('r')

    plt.subplot(3, 2, 4)
    plt.plot(time[:nstop], C3[:nstop])
    plt.plot(time[:nstop], np.zeros_like(C3)[:nstop], '-k')
    plt.ylabel('C3')

    plt.subplot(3, 2, 6)
    plt.plot(time[:nstop], degs*gamma[:nstop])
    plt.ylabel('gamma (deg)')

    plt.suptitle('thrust = 0.0001, start at r=4, time=4192.1, gamma=39.12 deg, r=87.90', fontsize=16)

    plt.show()

Answer 1

Sorry, must have been a typo in the comment. I went back to the original notebook in which I made the plots, and indeed the final $\gamma$ for the 0.001 acceleration case was 39.2°

It is not always 39.2°, but it does go to that asymptotically as the acceleration becomes smaller. Here is a plot of the $\gamma$ in degrees at $C_3=0$ as a function of the relative acceleration:

I am not aware of a way to determine that $\gamma$ analytically.

Below is the same plot for when accelerating tangentially, as opposed to in the velocity direction. It looks identical, except for the y-axis, where here it converges to 32.3°.

Though you wouldn't do that, since accelerating in the velocity direction is more efficient.