The Periodogram

When dealing with radial velocity data sets, the behaviour of the star is pretty unambiguous for a convenient sampling rate. For example, if a planet orbits a period over the course of a year, and the star is observed weekly, then there’s a clear sinusoidal curve in the data that can be easily picked out and modeled. However the sampling rate is often not convenient, and can even be larger than the period of the planet. Some data sets can essentially be great guessing games, where the periods are in the data but take some considerable effort to find. This effort takes the form of tuning the orbit parameters of the Keplerian fit until it finally lines up with enough data points.

For example, consider this radial velocity data set for 54 Piscium given by Fischer et al.

Radial velocity data set for 54 Piscium

It turns out there is a tool we can use to get a feel for what signals might be in the data. This tool, called a periodogram, was discussed earlier (here). Assume that we have a data set of velocities (or any variable, really). The radial velocity, v, is measured at a time t, resulting in a time series data with N data points, {v (tj), j=1,2,…,N0 }. A Fourier transform is performed on the data set as follows.

\displaystyle \begin{aligned} P_v(\omega)&=\frac{1}{N}|\text{FT}_v(\omega)|^2\\ &=\frac{1}{N_0}\left|\sum_{j=1}^{N_0}v(t_j)\exp(-\text{i}\omega t_j)\right|^2\\ &=\frac{1}{N_0}\left[\left(\sum_jv_j\cos(\omega t_j)\right)^2+\left(\sum_jv_j\sin(\omega t_j)\right)^2\right]\\ \end{aligned}

where P is the “power” of a given period ω, N is the number of data points, v is the radial velocity at time tj, and i is the imaginary unit.

With the above equation we can compute a basic periodogram for this data set. It is useful because if the data set contains a sinusoidal component with a frequency of ω0, then v (t) and exp(-iωt) are in phase. This causes a sinusoidal signal near ω to significantly contribute to the sums in the equation. Otherwise, values of ω that are far from a real sinusoidal signal in the data randomly fluctuate between positive and negative and therefore more or less cancel out to a small sum. Hence, if a sinusoidal signal of period ω exists in the data set, the value of Pv (ω) will be large. In a graph of periods on the x-axis, the presence of that period in the data set can be graphed as a “power” on the y-axis. Higher power values indicate that the periodicity is more strongly represented in the data.

With the periodogram, while we can technically look at the power of any of an infinite number of frequencies, it would be helpful to narrow our search to a finite set of frequencies. We must decide what frequencies those are. Since the periodogram of the data is symmetric (that is, P (-ω) = P (ω)), all of its information is contained in positive frequencies, n=0,1,2,…,N (=N0/2). Roughly half of the information in the data is discarded when calculating the periodogram because the absolute value part of the equation eliminates any distinction between positive and negative frequencies. As such, we need only to look at

$\displaystyle N=N_0/2$

evenly spaced frequencies.

It is not necessary to understand how the periodogram works, but what it means is important for this discussion. The periodogram takes your data set and tells you what periodicities are in it. For example, a periodogram of the above 54 Piscium data set would look like the following.

54 Piscium RV periodogram

A strong peak in the periodogram is found at ~62 days, corresponding to the location of 54 Piscium b. We can fold the data to the period of the planet and see a reasonable fit by a ~Saturn-mass planet in a 0.284 AU orbit with an eccentricity of 0.63, and a longitude of periapsis of 235 degrees.

54 Psc b fit

The periodogram, while useful, suffers a couple of problems. One, which is probably fairly obvious, is that it is noisy. A plethora of small peaks are seen that don’t really apply to anything physically real. The periodogram is noisy even when the actual data itself is only slightly noisy. Furthermore, the noise does not diminish with increasing sample size. This may be a bit counter-intuitive, as any fake signals would seem to be discouraged by adding new data, but it turns out that adding new data also adds new frequencies on top of the old, discouraged ones. As such, the noise is not actually averaged out by adding more data. However while the noise does not decrease with increasing sample-size, the signals in the data do. Specifically, the signal-to-noise ratio of a real signal in the data increases with the square root of the total number of observations.

Another major problem and contribution to noise is the problem of leakage. Signals tend to “leak” into other periodicities. In the 54 Psc periodogram above, you can see smaller peaks – called “sidelobes” – straddling the real signal of the planet. This is caused by other similar signals being sufficiently in phase that they end up having a higher P (ω) value. Such leakage to nearby periods is caused by the finite total interval over which the data is sampled.

Sidelobes are fairly easy to recognise and ignore, but a more dangerous sort of leakage occurs due to the finite size of the interval between measurements, a phenomena called aliasing. In particular, periods that are of integer ratios of the planet’s signal can be represented in the data. For a real signal of length S, aliased signals can exist at S/2, S/3, and so on. When you see a car driving down the road and the tires look like they are actually spinning very slowly, perhaps even in the opposite direction, the sampling rate of your eyes and brain are being fooled by an alias of the real signal – the actual rotation of the tires. Variable sampling rates are usually sufficient to eliminate the aliases that are most likely to be problematic.

As you can see, there are helpful statistical tools to take data sets and extract information out of them, but their interpretation is often unclear or ambiguous and great care must be taken when using them. As to the overreaching point that has been hinted at numerous times on this blog, detecting extrasolar planets is hard.

RV Fits: The Newtonian Solution

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

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

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

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.

Probability of Transit

Transiting Planets. Credit: NASA

Transiting planets are valuable items to explore the properties of planetary atmospheres. Planet searches like Kepler that focous on fields of sky tend to reap rewards amongst dimmer stars simply because there are many more dim stars in a given patch of the sky than bright ones. Transiting planets around bright stars are of particular value, though, as the increased brightness makes the system easier to study.

Radial velocity surveys tend to monitor brighter stars since spectroscopy is even more severely limited by stellar brightness than photometry, but it is not limited to observing patches of sky – telescopes performing Doppler spectroscopy tend to observe a single object at a time due to technical and physical limitations. Radial velocity surveys are also much less sensitive to the inclination angle of a planet orbit with respect to the plane of the sky. The planet doesn’t have to transit to be spectroscopically detectable. As such, radial velocity surveys tend to generate discoveries of planet candidates with unknown inclinations and true masses, but around much brighter stars than those planets discovered by the transit method.

As such, planet candidates discovered by radial velocity, especially planet candidates in short orbital periods are excellent targets for follow-up observations to attempt to detect transits. Transiting planets that have been discovered first through radial velocity have been of great scientific interest due to their host stellar brightness and thus ease of study. If more such systems are found, it would be of great benefit to understanding extrasolar planet atmosphere. While only a hand-full of transiting planets have been discovered first through radial velocity, they all orbit bright stars and are some of the best-characterised planets outside our solar system.

The probability that a planet will transit is, as has been discussed previously, given by
$\displaystyle P_{tr} = \frac{R_*}{a}$
where a is the semi-major axis of the planet orbit. This is the distance between the centre of the star and the centre of the planet. However, due to the inclination degeneracy – the reoccurring evil villain constantly plaguing radial velocity science – the star-planet separation is unknown. Remember that the period of the RV curve gives only the orbital period of the planet. If the orbital period is held constant, increasing the mass of the planet increases the star-planet separation. An increase in the total system mass requires greater separation between the two bodies to preserve the same orbital period.

For example, if radial velocity observations of a star reveal the presence of a mp sin i = 1 ME planet candidate, but the inclination is actually extremely low such that the true mass of the companion is in the stellar regime, then because the mutual gravitational attraction between the two stars will be much greater than the mutual gravitational attraction between the star and an Earth-mass planet at the same period, the two stars must have a wider separation, otherwise their orbital period would be smaller.

Mathematically, the true semi-major axis is given by
$\displaystyle a = \left(\frac{G[M_*+M_{\text{pl}}(i)]}{4\pi^2}\right)^{1/3}T^{2/3}$
Where G is the gravitational constant, and Mpl(i) is the mass of the planet at a given inclination i, and T is the period of the system. It is worth noting that the true semi-major axis is not significantly different from the minimum semi-major axis as long as the mass of the star is much greater than the mass of the planet – which is typically the case.

The fact that the true semi-major axis is a function of the unknown inclination makes for an interesting clarification: The probability that a planet of unknown inclination will transit is not simply given by Rstar/a, but is only approximated by it. If we assume that the distribution of planet masses is uniform (and extending through into the brown dwarf mass regime), then you would expect a planet with a minimum mass equal to Earth to have a much greater chance of being a bona-fide planet than a planet with a minimum-mass of 10 MJ, simply because there is a greater range of inclinations the former planet can be while still remaining in the planetary mass regime. Taking this a step further, even if both the Earth-mass planet candidate and the 10 Jupiter-mass planet candidate have the same orbital period, the probability that the latter planet transits ends up being less than the Earth-mass planet simply because of its high mass. Since its inclination is unknown, the probability that its mass is so high that the true semi-major axis is noticeably larger than the minimum semi-major axis is much higher, resulting in a likely lower transit probability.

Except it turns out that the mass distribution of planets and brown dwarfs isn’t constant. Earth-sized planets are significantly more common than Jupiter-sized planets, and super-Jupiters appear rare. It isn’t clear yet what the mass distribution planets actually is, with significant uncertainty in the sub-Neptune regime, but it is clear that for a highly accurate estimate of the transit probability, the inclination distribution cannot be thought of as completely random as it is fundamentally tied to the planet mass distribution.

Planet Mass Distribution given by Ida & Lin (Left) and Mordasini (Right)

Consider the case of a super-Jovian planet candidate, perhaps with a minimum mass of 7 or 8 Jupiter-masses. Because a significant fraction of physically allowable inclinations would place the true mass planet into a mass regime that is in reality sparsely populated, it is less likely that the planet candidate’s orbit is in those inclinations. It is thus more likely that the planet candidate’s orbit is edge-on than would be expected from the probability function of randomly oriented orbits. As such, the transit probability of a super-Jovian planet is actually boosted by ~20 – 50% over what you would expect from Ptr = Rstar/a. If this is the case, then we would expect to find an excess in the fraction of transiting planets in this mass regime then would be expected purely from the standard transit probability function. Indeed this is what we see.

Candidate planets with masses in the terrestrial planet regime are similarly affected, with broadened transit probabilies owing to the fact that terrestrial planets are more common than higher mass planets, arguing in favour of a higher inclination than the random inclination distribution function.

On the other hand, planet or brown dwarf candidates of minimum masses in the most sparsely populated region of the mass distribution are unlikely to truly have that mass. They are quite likely in orbits with low inclinations and with much higher true masses. The transit probability for companion candidates with minimum masses in this mass regime are actually reduced from the standard transit probability function.

Geometric and a posteriori transit probabilities

In the table above, taken from this preprint, we see that the geometric transit probability, Ptr,0, can be much less than the a posteriori transit probability, Ptr. The transit probability for 55 Cnc e, for example, jumps up from 28% to 36%. With these higher a posteriori transit probabilities, these short-period low-mass planets should be followed-up for transits. If transits are found, it would be of significant benefit to the extrasolar planet field.

In summary, there are various additional effects that can cause the a posteriori transit probability to be significantly different from the geometric transit probability. Planets with only minimum masses known can be more accurately assigned a transit probability when taking into account the uneven planetary mass distribution. Low-mass planets and super-Jupiters are more likely to transit than their geometric transit probability because a significant range of the inclination space is consumed by planets of masses that are simply rare. These planet candidates are more promising targets for transit follow-up than, for example, Jupiter-mass planets or intermediate-mass brown dwarfs.

The Real Ones

On the last post, we looked at recovering a periodic signal from a radial velocity plot and interpreting it as a planet. Now let’s look at a few of the complications involved in this.

A powerful statistical tool used to get an idea of what kind of periodic signals are in your radial velocity data set is a Lomb–Scargle periodogram (the mathematical details for the interested reader may be found here, but it is sufficiently complex to warrant skipping over in the interests of maintaining reader attention and reasonable post length). In the interests of brevity, further references to the Lomb-Scargle periodogram will be shortened to simply “periodogram.”

The purpose of this periodogram is to give an indication of how likely an arbitrary periodicity is in a data set whose data points need not be equally spaced (as is frequently the case in astronomy for a variety of reasons). Periodicities that are strongly represented in the data are assigned a higher “power,” where periodicities that are not present or only weakly present are given a lower power.

Let’s look at an example using a radial velocity data set for BD-08 2823 (source) If we calculate a periodogram for the data set, we come up with this

BD-08 2823 RV Data Periodogram

The dashed line represents a 0.1% false alarm probability (FAP). A clear, obvious peak is seen at 1 day, 230 days and ~700 days, implying that periodicities of 1 day, 230 days, and ~700 days are present in the data. Creating a one-planet model with a Saturn-mass planet at 238 days produces a nice fit. After subtracting this signal from the data, we’re left with the residuals. Now we may run up a periodogram of the residuals and see what’s left in the data.

BD-08 2823 Periodogram of Residuals

We see three noteworthy things. First and foremost is the emergence of a new peak in the periodogram that was not strongly present before at 5.6 days. We also see that the peak at 1 day remains. Lastly we see that the peak toward 700 days has weakened and moved further out. It would seem to suggest the 700-day signal is perhaps not real, or was an artifact of the 238-day signal.

Why was the 5.6-day signal not present in the first periodogram? The answer may lie in it’s mass: the planet has a mere 14 Earth-masses. It’s RV signal is completely dominated by the Saturn-mass planet. The giant planet forces the shape of the RV diagram and the signal of the second planet is just dragged along, superimposed on the larger signal.

On the radial velocity data plot, the two-planet fit we have come to looks like this:

BD-08 2823 Two-Planet Fit

It is important to realise that the obvious sine curve is not necessarily a bold line, but there is a second periodicity in there going up and down frantically, once every 5.6 days, compared to the Saturn-mass planet, at 237 days.

The fit has a reduced chi-squared of χ2 = 3.2, and a scatter of σO-C = 4.3 m s-1. There’s no obvious structure to the residuals and the scatter is not terribly bad, so any new signals will likely indicate planets of low mass. Let’s check in on the periodogram of the residuals to the two-planet fit and see what may be left in the data.

Periodogram of Residuals to 2-Planet Fit

That signal out toward a thousand days is stubbornly refusing to go away, despite a low χ2. It may either not be real, or it may be indicative of a low-amplitude signal with a rather long period.

Also noteworthy is that the periodicity at one day continues to exist, rather strongly. This periodicity is what’s known as an alias. Because the telescope observes only at night, the observations are roughly evenly spaced – there are (on average) twelve hour gaps between each data point. Therefore a sine curve with a period of 24 hours can be made to fit the data. To illustrate this, consider this (completely made up) data set:

Fake Signal

There’s no doubt that the data is well-fitted by the sine curve, but there is no real evidence that the periodicity proposed by it arises from a real, physical origin. What’s more, a sine curve with half this period could also equally well fit the data. So could a sine curve with a third of this period, and so on. There are mathematically an infinite number of aliases at ever-shortening periods that can be fit to this data.

Generally, if you observe a system with a frequency of $f_o$, and there exists a true signal with a frequency of $f_t$, then aliases will exist at frequencies $f_{t+i} * f_o$, where i is an integer.

Therefore we see that these aliases are caused by the sampling rate. If we could get data between the data points already available, if we could double our observation frequency, we could break this degeneracy. But the problem for telescopes on Earth is that the star is not actually up in the sky more than half the day, and a given portion of the time it is up could be during daylight hours. Therefore the radial velocity data sets of most stars can be plagued with short-period aliases since there is typically a small window of a few hours to observe any given star. It must be noted that as the seasons change and the stars are in different places in the sky at night, that window of availability will shift around a bit, allowing one some leverage in breaking these degeneracies. Ultimately, telescopes in multiple locations around the world (or one in space) would sufficiently break these degeneracies.

A real example of aliases exist in this example from an Alpha Arietis data set. In this case, the alias is not nearly so straightforward. Two signals of periods 0.445 days and 0.571 days can be modelled to fit the data.

Alpha Arietis RV Alias

So which of these two signals correspond to an actual planet? It turns out neither of them do: these radial velocity variations are caused by pulsations on the star – contracting and expansion of the star produces Keplerian-like signals in radial velocity data, too. That’s yet another thing to watch out for. This can be detected with simultaneous photometry of the star. If there is a photometric periodicity that is equivalent to your radial velocity periodicity, avoid claiming a planet at this period as if your academic credibility depends on it.

Additional observations could easily break this degeneracy, provided they are planned at times where the two signals do not overlap.

We see therefore that it is important to keep in mind that a low FAP speaks only to whether or not the signal is real, and not where or what it actually came from. The one-day periodicity is surely present in the data, but it is not of physical origin. It can also be extremely hard to tell whether or not a signal at a given period is actually an alias of another, more real period. There are times when the peak of an alias in the periodogram can be higher than the actual, real period. For reasons that include these, radial velocity fits must be considered fairly preliminary. New data may provide drastic revisions to the orbital periods of proposed planets if signals are exposed to be aliases.

Confusion over aliases have occurred before in literature. HD 156668 b and 55 Cnc e have both had their orbital periods considerably revised after it was realised that their published periods were, in fact, aliases. In the case of 55 Cnc e, the new, de-aliased orbital period ended up being vindicated after transits were detected). The GJ 581 data set, for example, is severely limited by sampling aliases that have spawned controversies over the possible existence of additional planets in that system.

In summary, periodograms are a useful tool to provide the user with a starting point when fitting Keplerian signals to radial velocity data, but they cannot distinguish real signals from aliases. Many observations with a diverse sampling rate are necessary to disentangle aliases from true planetary signals. Ultimately, a cautious approach to fitting signals to radial velocity data works best.

Making Waves

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

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

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

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^2$value 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^2$by 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

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

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.

Detecting Planets via Doppler Spectroscopy

We’ve established that an orbiting planet will gravitationally pull on its star in a barycentric motion around the system’s centre of mass (see here) and that we can detect this spectroscopically by watching the movement of absorption lines give away the radial velocity of the star (see here) Let us now, therefore, put the two together into some coherence.

Throughout the orbit of the planet, the observed stellar radial velocity will range between a maximum and a minimum value, half the difference between them we call the “velocity semiamplitude,” $K$. This can be calculated with

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

It might seem weird that there is no coefficient for the planetary mass in this. But realise that we aren’t observing the planet’s orbit, but rather we’re observing the star’s orbit around the star+planet barycentre. By taking the semi-major axis of the star and its orbital period into account, we’re accounting for the planet’s mass.

The velocity semiamplitude is related to the masses of the two components in a planetary system through the so-called “mass function”

$\displaystyle \frac{(m_* \sin{i})^3}{(m_*+m_p)^2} = \frac{P}{2 \pi G} K^3 (1 - e^2)^{3/2}$

Which, in the typical case of the star being much more massive than the planet, can be simplified and approximated

$\displaystyle m_p \sin{i} \approx \left(\frac{P}{2 \pi G} \right)^{1/3} Km_{*}^{2/3} (1 - e^2)^{1/2}$

Putting this approximation into Kepler’s third law gives us an expression for the semi-major axis of the planet.

$\displaystyle a \approx a_p \approx \left( \frac{G}{4 \pi^2} \right)^{1/3} m^{1/3}_* P^{2/3}$

Monitoring the radial velocity of a star over time gives us the semi-amplitude. The function of the radial velocity of the star with time may be expressed as

$V_r(t) = K[\cos{(v(t) + \omega)} + e \cos{(\omega)}] + \gamma$

Where $\gamma$ is the radial velocity drift of the star. A star not under the influence of a planet will still have a radial velocity value purely because of its galactocentric motion. Other sources of radial velocity drift include additional, more distant companions. For the case of an additional planet, the changing radial velocity in response to this outer planet will result in $\gamma$ not being constant, and will therefore be expressed as a rate, $d\gamma / dt$ or $\dot{\gamma}$, frequently in m s-2. Sometimes the RV drift is not linear, but curved, if significant progress of the outer planet’s orbit has been made over the time coverage of the observations, in which case it might be expressed as $d^2\gamma / dt^2$ or $\ddot{\gamma}$.

When plotting the radial velocity measurements as a function of time on a graph, one achieves a graph that resembles something similar to the graph below for the planet 51 Peg b.

We see, therefore, that fitting radial-velocity data with a Keplerian model gives four out of six of the orbital elements of the planet, the longitude of the ascending node and orbital inclination remain unknown. The mass of the star is estimated through other means but since the inclination is unknown, the Keplerian fit gives only a lower limit to the mass and semi-major axis of the planet.

The main take-away message is that with careful monitoring of a star with Doppler spectroscopy, we can measure the mass of the planet as multiplied by the sine of the (typically unknown) inclination, $m_p \sin{i}$, which in practice turns out to be the minimum mass of the planet (as $\sin{i} \leq 1$).

Update: Fixed exponent in $m \sin{i}$ equation. I inadvertently wrote it as negative.