  Characterization of transiting planets with transit timing variations (TTVs)
requires understanding how to translate the observed TTVs into masses and
orbital elements of the planets. This can be challenging in multi-planet
transiting systems, but fortunately these systems tend to be nearly
plane-parallel and low eccentricity. Here we present a novel derivation of
analytic formulae for TTVs that are accurate to first order in the planet-star
mass ratios and in the orbital eccentricities. These formulae are accurate in
proximity to first order resonances, as well as away from resonance, and
compare well with more computationally expensive N-body integrations in the low
eccentricity, low mass-ratio regime when applied to simulated and to actual
multi-transiting Kepler planet systems. We make code available for implementing
these formulae.
