Solution (source code)

= Solution

Write $x=1/a$ for the <comet>'s inverse <semi-major axis>. Its <specific orbital energy> is $-GM_\star x/2$, so escape corresponds to $x=0$. The printed quantity
$$
s_x=\frac{10}{a_p}\frac{M_p}{M_\star}
$$
has dimensions of inverse length. It cannot, literally, be a standard variance-per-time <diffusion coefficient>. Interpret it as the characteristic <root mean square> step in $x$ per periapsis passage, as in a discrete <comet energy diffusion> model.

If successive kicks are unbiased and uncorrelated, after $N$ passages the <root mean square> displacement is $s_x\sqrt N$. Starting at $a_0$ gives $x_0=1/a_0$, hence $N_{\rm ej}\sim(x_0/s_x)^2$. Taking the characteristic interval to be the initial <orbital period> $P(a_0)$ gives
$$
\boxed{t_{\rm ej}\sim\frac{P(a_0)}{100}
\left(\frac{a_p}{a_0}\right)^2
\left(\frac{M_\star}{M_p}\right)^2
=\frac{P_p}{100}\sqrt{\frac{a_p}{a_0}}
\left(\frac{M_\star}{M_p}\right)^2}.
$$
Here $P_p=2\pi\sqrt{a_p^3/(GM_\star)}$. In particular, for $a_0\sim a_p$ this is $P_p(M_\star/M_p)^2/100$. Without specifying the initial <semi-major axis> and the time per statistically independent encounter, the PDF cannot determine a unique time.

With the convention $\langle(\Delta x)^2\rangle=2D_xt$, the actual <diffusion coefficient> would be $D_x=s_x^2/[2P(a_0)]$. The estimate requires $N_{\rm ej}\gg1$, <orbit> crossing, and sufficient phase decorrelation. Near escape the <orbital period> grows, and resonant or secular correlations invalidate the simple constant-step clock. This is a characteristic diffusion scale, not an exact mean <first-passage time>: even an unbiased <Brownian motion> on an unbounded half-line has an infinite mean time to reach its absorbing endpoint.