# RV Fits: The Keplerian Solution

Johannes Kepler

We looked at deriving the radial velocity equation a little here but for the sake of this discussion, we will recover some of the material and expand on it. Furthermore, while modelling the radial velocity curve was described here, it was done so from the perspective of statistics and data, as opposed to actually modelling the stellar reflex velocity from of the gavitational influence of a planet. In hind sight, maybe it is this article which should have been titled “making waves.” I will seek to clarify the physical modelling of the radial velocity curve itself, and how the orbital orientation of the planet plays into how the Keplerian radial velocity curve is shaped.

The radial velocity equation is heavily reliant on Kepler’s Laws of Planetary Motion which govern the shape and periods of the orbits of planets and stars. As such, any fits to radial velocity data using the radial velocity equation is called a Keplerian fit. A Keplerian fit is therefore a fit to RV data that is effectively a model of the motion of the star as expected from Kepler’s Laws.

If you recall, the radial velocity equation is given by

$\displaystyle V_r(t) = K[\cos(\epsilon+\omega)+e\cos(\omega)]+\gamma+\dot{\gamma}+\ddot{\gamma}$

Where ω is the longitude of perihelion of the orbit, e is the eccentricity of the orbit, $\gamma$ is a term used to represent the radial velocity offset of the system (the radial velocity of the system arising simply from its motion through space), $\dot{\gamma}$ is any linear radial velocity trend in the data, measured in m s-1, and $\ddot{\gamma}$ is a term for a quadratic trend in the RV data, where the RV data may be curved but not in such a way that it permits the period to be readily known. This is typically observed when the aforementioned trend has been observed long enough that it begins to show some curvature in the RV data. If there is no need for any of the $\gamma$ terms, they may be set to zero. Lastly, K is the amplitude of the RV curve and directly corresponds to the mass of the planet through the semi-major axis of the star around the system barycentre, as well as being dependent on the unknown inclination of the system and the eccentricity of the orbit.

$\displaystyle K=\frac{2\pi a_*\sin{i}}{P\sqrt{1-e^2}}$

The nature of $\epsilon$ was not expanded on in the previous aforementioned page. It represents the “eccentric anomaly” and can be understood to be representative of where the planet is in its orbit at a given time. For circular orbits such that e = 0, the eccentric anomaly is simply equal to the mean anomaly, and is understood to be 0 at the start and 2π = 360o after a full orbit. For non-circular orbits such that e > 0, then it must be calculated.

Eccentric Anomaly (source)

The relationship between the mean anomaly and the eccentric anomaly is given by

$\displaystyle M = \epsilon - e \sin{\epsilon}$

It turns out that it’s impossible to solve for $\epsilon$ mathematically (in the same sense that it is impossible to solve for π) so a recursive algorithm will be needed to approximate it. This algorithm is called Newton’s Method and is given as

$\displaystyle x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

An initial guess, x0 is plugged into the algorithm and is used to determine x1, which is plugged back into the algorithm, ad infinitum. This algorithm will converge on the actual value for x regardless of the initial guess, with it returning the exact value at x, though of course after a handfull of iterations, a sufficiently close guess is arrived at to avoid the need for excessively redundant calculations. The initial guess is plugged into the original equation being solved for, divided by the initial guess plugged into the derivative of the original equation. As you might imagine from the orbit diagram above, the mean anomaly itself is a fairly decent guess at the the eccentric anomaly, so one is safe using this for x0. In the case of the expression for the eccentric anomaly, we first set the expression equal to zero: $0 = \epsilon-e\sin{\epsilon}-M$ and then with Newton’s Method

$\displaystyle \epsilon_{n+1} = \epsilon_n - \frac{\epsilon -e\sin{\epsilon}-M}{1-e\cos{\epsilon_n}}$

This may be done until $\epsilon$ is determined to the desired precision.

Graphing the RV equation is rather straightforwad with some computer scripting, though it can be somewhat challenging with a graphing calculating given the difficulty in determining $\epsilon$ to sufficient precision throughout the orbit (it turns out that astronomy uses both math and computers!). Let us examine how the shape of the RV curve varies with the orbital parameters (for a description of the various orbital parameters, see here). Obviously changing the period widens the the wavelength of the curve, and obviously varying K affects the amplitude. Changing the eccentricity will change the shape of the orbit significantly and will therefore affect the y-axis symmetry of the RV curve.

Manifestation of e on the RV curve

The “jagged” appearance of RV curves arising from eccentric planets is distinctive and difficult to replicate through other phenomena such as pulsations. As such, the RV signals from eccentric planets is less ambiguously of planetary origin than for planets in circular orbits.

Changes in the mean anomaly of the orbit simply effect a phase shift and move the RV curve of the planet along the x-axis (time axis). Changing the longitude of perihelion, ω, rotates the planetary orbit about its plane in space and therefore will cause what appears to be a phase offset to the orbit for a circular orbit, but will cause noticeable changes to the RV curve shape for eccentric planets. Changing ω for eccentric orbits causes a loss of x-axis symmetry.

Manifestation of ω on the RV curve

As before, the shape of the RV curve here is hard to replicate through processes intrinsic to the star. This “sawtooth” shape of the RV curve is much more reason to believe the RV curve is due to an orbiting companion rather than other issues such as starspots or pulsations.

Now multiple planets can be modeled by simply adding multiple RV curves over each other. For a system of n planets, the RV curve is given by
$\displaystyle V_r(t) = \gamma + \dot{\gamma} + \ddot{\gamma} + \sum_{i=1}^n K_i[\cos(\epsilon_i + \omega_i)+e_i\cos(\omega_i)]$

Keplerian modelling of RV data is simple, fast, and computationally easy. However it leaves out very important aspects which we will look at later, such as planet-planet interactions. That being said, Keplerian fitting is still the most widely used method to deriving extrasolar planet properties from RV data.