  First-order variational equations are widely used in N-body simulations to
study how nearby trajectories diverge from one another. These allow for
efficient and reliable determinations of chaos indicators such as the Maximal
Lyapunov characteristic Exponent (MLE) and the Mean Exponential Growth factor
of Nearby Orbits (MEGNO).
  In this paper we lay out the theoretical framework to extend the idea of
variational equations to higher order. We explicitly derive the differential
equations that govern the evolution of second-order variations in the N-body
problem. Going to second order opens the door to new applications, including
optimization algorithms that require the first and second derivatives of the
solution, like the classical Newton's method. Typically, these methods have
faster convergence rates than derivative-free methods. Derivatives are also
required for Riemann manifold Langevin and Hamiltonian Monte Carlo methods
which provide significantly shorter correlation times than standard methods.
Such improved optimization methods can be applied to anything from
radial-velocity/transit-timing-variation fitting to spacecraft trajectory
optimization to asteroid deflection.
  We provide an implementation of first and second-order variational equations
for the publicly available REBOUND integrator package. Our implementation
allows the simultaneous integration of any number of first and second-order
variational equations with the high-accuracy IAS15 integrator. We also provide
routines to generate consistent and accurate initial conditions without the
need for finite differencing.
