  Mercury's orbit can destabilize, generally resulting in a collision with
either Venus or the Sun. Chaotic evolution can cause g1 to decrease to the
approximately constant value of g5 and create a resonance. Previous work has
approximated the variation in g1 as stochastic diffusion, which leads to a
phenomological model that can reproduce the Mercury instability statistics of
secular and N-body models on timescales longer than 10 Gyr. Here we show that
the diffusive model underpredicts the Mercury instability probability by a
factor of 3-10,000 on timescales less than 5 Gyr, the remaining lifespan of the
Solar System. This is because g1 exhibits larger variations on short timescales
than the diffusive model would suggest. To better model the variations on short
timescales, we build a new subdiffusive phenomological model for g1.
Subdiffusion is similar to diffusion but exhibits larger displacements on short
timescales and smaller displacements on long timescales. We choose model
parameters based on the behavior of the g1 trajectories in the N-body
simulations, leading to a tuned model that can reproduce Mercury instability
statistics from 1-40 Gyr. This work motivates fundamental questions in Solar
System dynamics: Why does subdiffusion better approximate the variation in g1
than standard diffusion? Why is there an upper bound on g1, but not a lower
bound that would prevent it from reaching g5?
