Category Archives: Modelling

Orbit Projections in the Sky


Projection (Source)

Whether you are determining the barycentric motion of a star with µas precision to determine the orbit and mass of an unseen planetary companion, or directly observing planetary companions orbiting dozens of AU from their stars over the course of some decent fraction of a decade to pin down their orbits, you are basically trying to determine the motion of an object in the sky. While the orbit of the object is three-dimensional, the path it takes in the sky is not (at least not for monitoring nearby objects, especially in the solar system). An integral part of binary star astrometry is therefore projecting a 3D orbit onto a 2D plane, and fitting data to that observed 2D orbit. Since a star plus planet system is physically analogous to a binary star with a high mass ratio, the mathematics ends up being the same.

With our focus thus far in this blog on Doppler spectroscopy and transit photometry, the orbital parameters we have concerned ourselves with are the semi-major axis, a, the eccentricity of the orbit e, the inclination of the orbit, i, and the longitude of periapsis ω, which define the size of the orbit, the deviation from circular, the angle between the orbit plane and the plane of the sky, and the angle of the periapsis of the orbit from the observer, respectively. For the longitude of periapsis, ω=0° is defined in such a way that the line connecting the star to the periapsis of the orbit is perpendicular to the line of sight — the Earth-star-planet angle is a right angle. ω=90° after a rotation of the orbit 90° in the orbit plane in such a way to where if the inclination of the orbit is 90°, the periapsis would be between Earth and the star, and the transit midpoint would occur at periapsis.

Longitude of Periapsis

Changing the longitude of periapsis rotates the orbit in the orbit plane

We have virtually ignored the ascending node, Ω, which defines the rotation of the orbit plane around the line of sight. For some inclination near 90°, when Ω=0°, the planet orbits “up-and-down” and when Ω=90°, the planet orbits “left-and-right.” Note that this does not necessarily mean a polar or equatorial orbit around the star, as in this discussion, we are agnostic about the orientation of the stellar rotation axis. The ascending node of a planetary orbit has been mostly ignored in this blog because it does not actually affect either the Doppler behaviour of the star, or for the most part the shape of the transit light curve. Technically, with all else held constant, varying Ω for a transiting planet will affect the projected angle between the star spin axis and the planet orbital axis, λ, which is detectable as the Rossiter-McLaughlin effect with Doppler spectroscopy (See here). However, since we do not know the orientation of the stellar spin axis in space, we aren’t able to fit Ω as an adjustable parameter to the data.

Ascending Node

Changing the ascending node rotates the orbit around the line of sight

Plotting 3D orbits as they appear on the (locally) two-dimensional plane of the sky to determine it’s orbit is the astrometric equivalent of calculating the radial velocity behaviour of a star. We input orbital parameters into a model and try to fit the data to that model to assess how well it approximates reality.

First, we will need to define a couple of functions that address the orbital motion of the object in the x and y dimensions in the sky, let’s call these functions f_x(t) and f_y(t).

\displaystyle f_x(t)=\cos E - e \\ f_y(t)=\sqrt{1-e^2} \sin E

where e is the eccentricity of the orbit, t is the time of calculation, and E is the eccentric anomaly as derived and calculated in the same way as in the Doppler spectroscopy method discussed here.

Now we will need to bring in transformations to describe the deprojection of the 3D ellipse into a 2D plane. These take the form of the Theile-Innes constants:

\displaystyle A=a(\cos\omega\cos\Omega-\sin\omega\sin\Omega\cos i)\\ B=a(\cos\omega\sin\Omega+\sin\omega\cos\Omega\cos i)\\ F=a(-\sin\omega\cos\Omega-\cos\omega\sin\Omega\cos i)\\ G=a(-\sin\omega\sin\Omega+\cos\omega\cos\Omega\cos i)

Where a is the semi-major axis of the orbit, ω is the longitude of periapsis, Ω is the ascending node and i is the inclination of the orbit. We now combine the two to get the position of the object in the x- and y- axes as a function of time:

\displaystyle x(t)=B f_x(t)+Gf_y(t)\\y(t)=Af_x(t)+Ff_y(t)

We are now in a position to evaluate the goodness-of-fit of a given model to data using essentially the same statistical tools and techniques described in this post. The observed-minus-calculated for an astrometric position is in this case of course determined with Pythagorean Theorem, O-C=\sqrt{\Delta x^2 + \Delta y^2}. As an example, here are the positions of a star in a 15 year orbit around the supermassive black hole in the centre of the Galaxy.

Star orbiting SgrA*

Star orbiting SgrA* (source)


RV Fits: The Newtonian Solution

Isaac Newton

Isaac Newton

In the last post we looked at modelling the radial velocity of a star as the reflex motion due to an orbiting planet’s gravitational influence. This is the most prominent way used to model extrasolar planets to fit RV data. The Keplerian solution offers a lot in its simplicity, but it has a major shortcoming in that it is not particularly dynamic, and really only approximates any system with more than two bodies.

A planetary system may be thought of as an n-body problem, where you have n gravitating objects all influencing each other. This is different from the Keplerian models that take into account only the star-planet gravitational interaction. This will be a rather rudimentary look at simulating a gravitational n-body model.

First let’s look at vectors. In this case, we’ll consider a velocity vector \vec{v} defined simply with an x and a y component.

\displaystyle \vec{v} = [v_x, v_y]

This simply allows us to get some sense of where our velocity is toward in a two-dimensional cartesian plane. We’re used to saying that an object has a velocity of, for example, 5 m s-1, but this alone tells you little useful in the context of modelling a system – 5 m s-1 in what direction? By defining a vector with an x and y component, we can specify the direction that the velocity is in. For example, \vec{v} = [1, 0], then the velocity of the object is to the right, whereas \vec{v} = [0, -1] implies moving down. If \vec{v} = [1, 1], then the object is moving both to the right and up, with the actual total velocity being determined by Pythagorean Theorem.

\displaystyle v = \sqrt{v_x^2 + v_y^2}

This is called the magnitude of the velocity vector.

Now each particle will have its own gravitational field. To determine the attraction to a particle from a given coordinate, we can use Newton’s Law of Universal Gravitation.

\displaystyle \vec{F_g} = -G\frac{m}{r^2}

Where m is the mass of the particle we’re considering, let’s call it p_0, and r is the distance between it and the other particle, let’s call it p_1. Since we’re considering in this context only the gravitational attraction of p_1 on p_0 (not both particles on each other), we need not take into account the mass of both particles at this time.

Once you have determined the gravitational attraction of p_0 toward p_1, you will need to determine which direction the gravitational attraction is toward. The angle between two points is simply given by

\displaystyle \theta = \tan^{-1}{\frac{y_0-y_1}{x_0-x_1}}

The x and y components of the gravitational attraction may be determined with F_g \cos{\theta} and F_g \sin{\theta}, respectively. These x and y components may be added to the x and y components of the velocity vector to determine the new \vec{v}.

This will have to be done for every particle. You will have to calculate the gravitational influence of all particles on each other particle. For i particles, the new velocity of any given particle j, \vec{v}'_j, is given by

\displaystyle \vec{v}_j' = \vec{v}_j -  \sum_{i\neq j}^i G\frac{m_i}{r^2_{ij}}

Where r is the separaton between the two particles. For most chosen values for \vec{v}, the system will have a systemic velocity (corresponding to the γ term in the Keplerian radial velocity equation). If you want to determine \vec{\gamma}, you can simply sum up all of the initial velocity vectors and divide it by the total system mass M.

\displaystyle \vec{\gamma} = M^{-1} \sum_i m_i\vec{v}_i

N-body simulations are calculation-intensive and are typically run on computers, where gravitational attractions and their modifications of velocity vectors may be calculated over and over for each new position. Particles are typically assigned a mass, a starting coordinate and a starting velocity vector, and a computer will calculate new positions as dictated by Newtonian physics for as long as the user wants.

Let’s look at a specific example. For simplicity, I set the gravitational constant to unity and use no units. I simulate three particles, a star and two planets, with positions of (0,0), (-30, 0), (60, 0), starting velocity vectors of (0.1, 0.1), (0, 0.5), (0, -0.22), and masses of 10, 0.1 and 0.1, respectively. Letting the simulation run for 3,000 time steps and plotting the y-axis component of the star’s velocity vector provides us with the stellar radial velocity.

Simulated RV data for chaotic system

Simulated RV data for chaotic system

What is going on? The RV behaviour is decidedly non-Keplerian – there are no Keplerian models that would be able to fit this data. If we look at a diagram of the system, it becomes clear what is happening here.

Planet Scattering

Planet Scattering

Our radial velocity plot is showing us a nightmare in planetary system architectures: planet-planet scattering. The outer planet (whose positions are shown in green) is clearly interacting very strongly with the inner planet (whose positions are marked in red).

Not all dynamically interacting multi-planet systems are so chaotic. Systems often have planets in close but stable orbits where they gravitationally tug on each other, causing the stellar radial velocity profile to deviate from a simple Keplerian model.

While Keplerian models are not capable of fitting to planetary systems with interacting planets, this Newtonian fit can. It is dynamical in that it simulates the influences of the planets on each other, however it can be particularly computer-intensive. Newtonian fits are often run of systems where Keplerian fits suffice simply to probe the dynamical behaviour of a system, such as its long-term stability or possible resonances between the orbits. Knowing whether to use Keplerian or Newtonian models for planets will depend on the planetary system architecture, but it is clear that Newtonian models are powerful tools for predicting the behaviour of planets in multi-planet systems.

RV Fits: The Keplerian Solution

Johannes Kepler

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

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

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

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.

Modelling the Transit Light Curve

The Sun (Credit: NASA)

The Sun (Credit: NASA)

We looked at a simplistic understanding of planetary transits and transit light curves (here). In the interests of accuracy, it would be worth investigating a more formal understanding of the transit light curve. Specifically, this post will be aimed at actually modelling transit light curve while taking into account some more complicated dynamics such as limb darkening.

The actual shape of the transit light curve can be mathematically described by a function of the sky-projected centre-to-centre distance between the star and the planet, z, where z = 0 is understood to be the mid-transit of a centrally-crossing transit and the impact parameter, b is understood to be the minimum value of z throughout the transit, and the transit starts at about z ~ 1. More accurately, first and fourth contact start at z = R_* + R_p, second and third contacts start at z = R_* - R_p. Furthermore, let’s define p as the ratio of the radii of the planet and the star, p = R_p/R_*.

Let us now define a function of the total fractional area of the star that is not obscured by the planet. Normally you might be inclined to feel that (R_p/R_*)^2 suffices, however it describes only the area ratio between second and third contacts, and is thus inadequate to describe the entire transit. We therefore will construct the amount of obscured stellar area as a piecewise function of p and z whose intervals describe different parts of the transit. It will be given as the ratio between the total area and the unobscured area subtracted away from the whole total area. Specifically, A(p,z) = 1 - \lambda(p,z).

\displaystyle  \lambda(p,z) =  \begin{cases}  0, & 1 + p < z\\  \frac{1}{\pi} \left[ p^2\kappa_0+\kappa_1 - \sqrt{\frac{4z^2-(1+z^2-p^2)^2}{4}} \right], & \lvert 1 - p\lvert < z \le 1 + p\\  p^2, & z \le 1 - p\\  1, & z \le p - 1  \end{cases}

Where \kappa_1 = cos^{-1}[(1 - p^2 + z^2)/2z] and \kappa_0 = cos^{-1}[(p^2 + z^2 - 1)/2pz]. There’s a lot going on here but it isn’t too complicated to understand. The first case, 1 + p < z, describes when the planet is completely off the stellar disc, before first contact hence there is no (zero) obscured area. The second case is the most complicated and describes the time between first and second (third and fourth) contacts. The third case is what we covered earlier, where the amount of obscured area is simply the square of the ratio of their radii, and therefore simply p^2. The final case is only satisfied in the (unlikely) condition that the planet-to-star radius ratio > 1, and the entire area of the stellar disc is obscured.

Calculating z as a function of time may be done with knowledge of the impact parameter and use of Pythagorean theorem.

For a transiting planet with an impact parameter of zero and a radius ratio of R_p/R_* = 0.1 (producing a 1% transit depth), the above model produces the following light curve plot.

Simulated Light Curve

Simulated Light Curve

Real light curves are actually curved, however, whereas this model is clearly much more “rigid” in shape. The reason for this is an effect known as limb darkening, where the outer layers of the star are not able to isometrically scatter light from underneath, causing less light from the limb of the star to reach the observer than from the centre of the stellar disc. For this reason, stars generally appear dimmer at their limbs and brighter in their centres. The magnitude of this difference will depend on a number of factors intrinsic to each star.

Proper handling of limb darkening will be required to obtain high-fidelity fits to high-accuracy data. A particularly popular and effective model for limb darkening was presented by Claret (2000). We can consider a four-parameter limb-darkening law describing the relative brightness of a point on the star expressed as a function of the angle between the observer, the centre of the star, and a line between the stellar centre and the given point on the star.

\displaystyle I(r) = 1 - \sum^4_{n=1} c_n(1 - \mu^{n/2})

Where \mu is given as \cos \theta = \sqrt{1 - r^2}, and each value of c is called a “limb darkening coefficient.” In the case of the centre of the stellar disc, I(0) = 1. Each star will have a unique set of limb darkening coefficients, and they will change depending on what wavelength the transit is being observed in (due to how the properties of the stellar photosphere’s light scattering varies with wavelength). One could choose to have as many limb darkening parameters (and thus coefficients) as desired, but a high number is rarely required to adequately fit transit light curve data.

Limb darkening must be taken into account to get accurate measurements of radius ratios of transiting objects, and thus planetary radii. Because limb darkening will concentrate the brightness of the star into the centre of the disk, a planet that would cause a 1% transit depth in a uniformly bright star would cause a deeper transit depth in a central transit around limb-darkened star, and a shallower transit depth in a grazing transit around a limb-darkened star. Proper care in modelling the limb-darkening behaviour of the star is therefore required to ensure good measurement of the planetary radius.

We see therefore that the brightness difference is not a pure function of the ratio of the stellar and planetary radii. \Delta F = (R_p/R_*)^2 is only an approximation, and only correct for the unrealistic case of a uniformly bright stellar disc.

So what can we do to get a more comprehensive representation of a light curve? Let’s modify the equation for \Delta F we’re familiar with and consider brightnesses instead of just purely area ratios. Let’s take our limb-darkening equation as the radius of the star and determine the area, since it describes the brightness of the star from centre to limb. With this equation as the stellar “radius,” we use some simple calculus to expand it to the “brightness area” of the circle that is the stellar disc. This calculates the total brightness of the star while taking into consideration the changing brightness behaviour from centre to limb. The actual transit depth of a planet at any given distance to the centre of the star (represented as z again) is simply the ratio of the total stellar brightness in the area of the star under the planet to the total brightness of the stellar disc. Therefore, we draft a new function based off our original understanding of \Delta F, while taking into effect the above model for how much of the star’s area is covered at any given time in the transit, as well as bringing in our limb-darkening model to describe the brightness distribution of the stellar disc:

\displaystyle F(p,z) = \left[ \int^1_0 dr 2rI(r) \right]^{-1} \int^1_0 dr I(r)\frac{d[A(p/r,z/r)r^2]}{dr}

And it makes quite a difference, too. By incorporating limb darkening coefficients, we’re able to much more accurately simulate the transit of a planet across a limb-darkened star. The diagram below shows the same planet (with p = 0.1) transiting stars of different limb-darkening parameters. The shape and modification to the transit depth is readily apparent.

Effects of Limb-Darkening coefficients

Effects of Limb-Darkening coefficients

Making Waves

Photo by Mila Zinkova, 2007

Photo by Mila Zinkova, 2007

We looked at how the Doppler effect can tell you about the radial velocity of a star (here), and how you can use this information to detect an orbiting companion (here). Let’s now look a bit more at some of the later parts of this process. You’ve already finished your observations, you’re done with the telescope time and you have a collection of data points. What do you do with them?

The first thing you will want to do is to create a model that fits the data. This is done using the equations in the post referred to in the above paragraph. The model will simply manifest itself as a plot of the radial velocity behaviour you would expect for a planet with a given mass and orbit. The closer your model fits the data, the more likely your model is a good representation of reality. For example, let’s look at radial velocity data of the giant star Iota Draconis (Credit: ESO).

Iota Draconis (HR 810) Radial Velocity Data

Iota Draconis (HR 810) Radial Velocity Data

The individual points are the actual data, whereas the sine wave is the model of a 2.26 MJ planet with a 320-day period. We see that the proposed model fits the data well (being very close to the data), and is thus likely representative of reality. There are a few stray points and the model is not a perfect match to the data, but the data points can be prone to intrinsic noise, such as flares or granulation on the star’s surface, or instrumental noise with the spectrometer.

We see that the radial velocity data set covers several orbits of the planet around the star. As the planet continues to orbit the star (as it will do so for aeons) and the radial velocity data continues to come in (as we would hope would occur, more data is more opportunities to discover something new), the graph could become extremely long and complicated. For example, consider a set of radial velocity data for HD 192263 (source).

Unfolded HD 192263 RV data

Unfolded HD 192263 RV data

What we could do would be to fold the entire data set by a certain period of time, specifically the orbital period (taking the time of the measurement modulo the orbital period). Since the model is periodic, that is, repeating each orbital period, the data points (provided they do have the same periodicity that you’re folding by) will fold up into a nice, neat, singular sine curve.

HD 192263 Folded RV data

HD 192263 Folded RV data

Here it is quite convincing that the model well fits the data. But how do we quantify how good of a fit to the data our model is? A value exists called \chi^2(read “chi squared”).

\displaystyle \chi^2 = \sum \frac{(O-C)^2}{\sigma^2}

That is, calculate the difference between each observed data point and the model (or “observed minus calculated”), divided by \sigma^2, which represents the variance of the data. Then total up each of these. The higher the \chi^2, the poorer the model fits the data. The variance, \sigma^2,can be determined with

\displaystyle \sigma^2 = \frac{1}{N} \sum_{i=1}^{N}(v_i - \bar{v})^2

Which is simply an average of the differences of each point, v, from the average values of the points \bar{v}. To put this in other terms (and this is more of a discussion about general statistics than it is specific to radial velocity data), all of the radial velocity measurements can be averaged together. The average of the differences of each point from that average radial velocity value is the variance of the data, \sigma^2, and is a representation of how scattered the data is. Our \chi^2value is simply the total of each point’s difference from our model, while dividing it by the variance to compensate for highly scattered data. Models that are not rigidly affixed to the data points (high total OC) are given more forgiveness of the data itself has a lot of variance (\sigma^2) anyway. If there is low variance, then only models that are closer to the datapoints (low OC) are treated well, or given a low \chi^2 value.

We can obtain another value called the “reduced chi-squared” which takes into account the number of freely adjustable variables (or, the “degrees of freedom”). This is obtained simply by dividing \chi^2by the degrees of freedom, d.

\displaystyle \chi^2_{red} = \frac{1}{Nd} \sum \frac{(O-C)^2}{\sigma^2} = \frac{\chi^2}{d}

The need to do this stems from the fact that a quantitative analysis of data must be undertaken to seriously distinguish between two models whose goodness of fit to the data may not be as easily discriminated between in our “chi-by-eye” technique of just looking at them.

Large reduced chi-squared values indicate a poor model fit, and a value of 1 is typically sought but can be extremely difficult or impossible to achieve due to instrumental or intrinsic noise that are out of one’s ability to control. Values of the reduced chi squared that are less than 1 are typically a sign that you are “over-fitting” the model to the data. This can result from trying too hard to fit the model to noise, or overestimating the variance.

HD192263 RV Residuals

HD192263 RV Residuals

Above are the residuals to the HD 192263 dataset. This is keeping the x-axis (the time of the measurement) and plotting the O-C on the y-axis, instead of the radial velocity value for that data point. In this example, there is a trend in the data that has been modelled out which could potentially be due to a second orbiting companion that is far enough away that only a very small section of its orbit has been observed.

Other than the trend, there is no apparent structure in the residuals, and so there is no reason to suspect the existence of additional orbiting bodies.

One can calculate the standard deviation, or variance, of the residuals to get a measure of how scattered the residuals are. This is typically represented with \sigma or written as ‘rms’ in literature, and is an algebraically obvious transformation of one of the above equations,

\displaystyle \sigma = \sqrt{\frac{1}{N} \sum_{i=1}^{N}(v_i - \bar{v})^2}

The radial velocity jitter of a star scales with both its mass and chromospheric activity. A quiet M dwarf may have a jitter as low as 2 to 3 metres per second, whereas a sun-like star may have a jitter on the order of 3 – 5 metres per second. F-type stars have jitters close to 10 metres per second, and the jitter sky-rockets from there, where A-type stars are nearly unusable for using radial velocity to detect extrasolar planets. As of the time of this writing, no planets have ever been discovered around an A-type star through the use of the radial velocity method. So there is an amount of variance in the residuals we should expect to have given the star we’re studying. If the variance is much higher than expected for the star, then it could be evidence of either unusually high stellar chromospheric activity, or as-yet undetected planets with small amplitude signals.

The systems we have looked at thus far throughout this blog in the context of radial velocity have been systems whose radial velocity data graphs are dominated by a single periodicity. Indeed, our own solar system is like this, with Jupiter dominating the sun’s radial velocity profile. But not all planetary systems come so conveniently organised.

I leave you with an image of how daunting this can be and why we are all thankful that computers are here to help us: An image of the 2008 three-planet fit for HD 40307 (since discovered to have three more planets).

HD 40307 Three-Planet Fit

HD 40307 Three-Planet Fit

Update (15 Apr 2013): I have corrected a math error in the equation for variance and standard deviation. In both cases the (v_i-\bar{v}) term should have been squared.