Category Archives: Orbital Mechanics

Tidal Theory

Tidal Wave

Tidal Wave (Source)

We’re used to thinking of various astronomical bodies having gravity and gravitationally interacting with each other, and it’s certainly true that they do, and that these interactions can be exploited to detect extrasolar planets through their pull on their parent stars (either by astrometry or Doppler spectroscopy). So far in our study of extrasolar planetary systems, we have treated astronomical bodies as point masses — that is, the entirety of their gravity is assumed to come from a single point, and that the entirety of their matter exists at that point. But this is, of course, an oversimplification. The consideration of astronomical bodies as extended objects introduces a range of new effects that are collectively called “tides.”

We may be familiar with the concept of the tides in the ocean, where the gravitational pull of the Moon on the oceans of Earth causes them to be raised twice a day, and so this is where we will start our consideration of tides. The gravity of Earth’s moon pulls on Earth, causing it to stretch into a more elliptical shape. A diagram illustrating this is shown below.

We see then that the rise over the side of Earth opposite the Moon is caused by the distortion of the shape of Earth. It’s important to understand that this is a very general phenomenon: It is not exclusive to the Earth-moon system, but any two-body system can be represened this way. Furthermore, while it is not shown in the above image for simplicity, the primary body (Earth in the case of the Earth-Moon system, the star in the case of a hot Jupiter system) will also raise a tide on the secondary mass.

The tidal deformation of the primary body is not an instantaneous process. There is inertia to overcome and so if the primary mass is rotating, the tidal bulge will be dragged away from the line connecting the two bodies. A diagram of this is shown below.

In the above image, we assume a prograde orbit for the satellite, and r is a line joining the centre of the two bodies. Note that the rotation of the primary forces its long axis to lead ahead of the orbiting body. For a circular orbit, the long axis of the planet will lead ahead of r by some constant angle, even while the planet rotates. This is why there are two lunar tides each day. Aside from causing boats to rise and descend in the water, relative to fixed points on land, these tidal effects have some important implications.

Firstly, it affects the orbits of the two bodies. In the case of the Earth-Moon system, the oceanic tidal bulge facing the Moon has a stronger gravitational attraction on the Moon than the tidal buge on the opposite hemisphere on Earth, simply because it is closer (and gravitational attraction scales inversely with r^2). Since the nearer tidal bulge is pushed ahead of the r line, the net gravitational attraction of the Moon toward Earth is also displaced from r. Since the tidal bulge is pushed ahead of the Moon in its orbit, the gravitational attraction of the Moon toward Earth has a small component in the direction of its orbital motion. This causes the Moon to accelerate (it’s orbital velocity increases), which then causes its orbital altitude (the semi-major axis) to increase. Any readers casually aware of orbital mechanics is familiar with this effect. If you’re unfamiliar with orbital mechanics, a full description of it is beyond the scope of this writing, but I would encourage you to try Kerbal Space Program.

Conversely, if the tidal bulge lagged behind the r line, the orbit of the Moon would decay. This can be done either by having the Moon in a retrograde orbit as is the case for Neptune’s moon Triton, or an orbit that is closer to the planet than the geosynchronous orbit radius, as is the case for Mars’ moon Phobos. Triton and Phobos are thus doomed.

Let’s define two factors (given by Dobbs-Dixon et al. 2004, Ferraz-Mello, et al. 2008), \hat{s} and \hat{p}, as the strength of the stellar and planetary tides, respectively,

\displaystyle \hat{s}\equiv\frac{9}{4}\frac{k_0}{Q_0}\frac{m_1}{m_0}R_0^5,\qquad\hat{p}\equiv\frac{9}{2}\frac{k_1}{Q_1}\frac{m_0}{m_1}R_1^5

Where k_0 and k_1 is the Love number of the primary and secondary, respectively, Q_0 and Q_1 is the tidal dissipation factor of the primary and secondary (a quantity that encapsulates how a body responds to tidal deformation, but it is very hard to measure, and usually poorly known), respectively, m_0 and m_1 as the masses of the two bodies, and R as their radii. The Love number quantifies the mass concentration of the body and will not be expanded on significantly here. The separation between the two bodies will change with a rate of

\displaystyle \dot{a}=-\frac{4}{3}na^{-4}\hat{s}\left[(1+23e^2)+7e^2(\hat{p}/2\hat{s})\right]

Where a is the semi-major axis of the orbit, e is the eccentricity of the orbit, and n is the

Secondly, the tidal deformation of a body will be resisted by inertia, causing its rotation to slow (and its synchronous orbit radius to expand outward). This will continue until either the body’s spin period mathes its orbital period (as is the case for the Moon), or is some integer ratio of the orbital period other than 1:1 (as is the case for Mercury), as is often the case for bodies in eccentric orbits.

For fluid bodies (such as gas giant planets or stars) in circular orbits, the end result of tidal spindown will result in a perfectly synchronous rotation, via a process called tidal synchronisation. In eccentric orbits, the planet moves around the star at varying speeds, making a synchronous rotation impossible. For such planets, the end result of tidal spindown is “pseudosynchronisation,” with a planet in a pseudosynchronous rotation such that the rotation rate of the planet matches the angular velocity around the star at periapsis (where the planet moves fastest), but with the rotation of the planet being faster than the angular velocity elsewhere. In rigid bodies (such as the terrestrial planets), tidal synchronisation is achieved in the form of a spin-orbit resonance, such that the long axis of the body is aligned with $r$ at each periapsis. It doesn’t need to be the same hemisphere of the body, as we see in the case of Mercury, which presents alternating hemispheres toward the sun at each periapsis.

Thirdly, tidal deformation of a body in a non-synchronous, eccentric orbit will cause the orbital eccentricity to tend to zero. This process, called “tidal circularisation,” is responsible for the circular orbits we see for the hot Jupiter population. As a body swings around periapsis, the tidal deformation intensifies and inertia resists this motion. In addition to causing heating to the object, it also saps some of its orbital energy, causing it to slow down. As this happens at periapsis more strongly than anywhere else, the net result is that the eccentricity of the orbit tends to zero.

The rate of the change in the eccentricity is given by

\displaystyle \dot{e}=-\frac{2}{3}nea^{-5}\hat{s}\left[9+7(\hat{p}/2\hat{s})\right]

Tidal deformations are a significant source of heating of a form that is collectively called “tidal heating.” An extreme example of it is seen at Io, where Jupiter’s other large moons perturb Io, pumping up its eccentricity, which is subsequently damped by tidal interaction with Jupiter, causing heating, and some rather spectacular volcanism.

Understanding the influences of tides on planetary systems is a vital component of understanding their evolution, formation and habitability. Tidal interactions between moons and planets can completely melt those moons, and tidal interaction between planets can tidally heat those planets and potentially completely melt them, too. Tides can also dominate the final rotational behaviour of a planet, with effects for its climate that are still being examined.


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.

The Path of the Wanderers

The Heliocentric Model of the Solar System

While the ancients did not understand the motions of the “wandering stars,” Newton and Kepler made great strides toward understanding them.

Assuming the simple case of a single body, A, orbiting another, B, the two bodies will orbit around a common barycentre that represents the centre of mass of the system. Thus while it is acceptable to say that one body orbits another, it is more physically accurate to say that both bodies are in a mutual orbit about each other.

Assume the existence of two objects, A and B, with masses MA and MB, respectively. Let aA be the distance from A to the barycentre of a system, and let aB be the distance from B to the barycentre. The ratio of aA / aB is determined purely by inverse ratio of the two masses. For example if MA = MB, then aA = aB, and if MA = 3MB, then aA = (1/3)MB. Physically, this is because A is more massive than B, and so it has more gravity and is thus B needs to be farther from the centre of mass to balance the system.

In the latter case, the gravitational acceleration of A on B is also greater, causing B to have a higher velocity than A. Thus despite B having a wide orbit, it will complete the orbit at the same time as A, and thus, the orbital period, P, of the two will be equal. The fact that is true for all values of is important for constraining the orbital period of one object if the other is not visible for whatever reason. The orbital period of both bodies may be defined as

\displaystyle P = 2\pi \sqrt{a^3 / \mu}

Where a is the semi-major axis, roughly the average distance between A and B (and is defined as half the length of the orbits longest axis, the “major axis”), and μ is defined as μ = GM, the product of the gravitational constant, G, multiplied by the total masses of the two bodies, M. It is necessary to consider the sum of the masses because the gravitational pull of both objects on each other will contribute to the orbital period. In cases where MA > MB, the contribution of B may be considered negligible and the orbital period determined more or less accurately without taking it into account (for example, calculating the orbital period of Earth about the sun may be done reasonably well without taking Earth’s gravitational pull on the sun into account).

Not all orbits are perfect circles. Some may be eccentric. Mathematically, an orbit should never be thought of as a circle, but rather as an ellipse with the barycentre at one focus. However an orbit with an eccentricity of e = 0 is physically indistinguishable from circular.

For an eccentric orbit (any orbit for which e > 0), there exists a point in the orbit where the bodies are closest to each other, this is referred to as the periapsis, ap. The point in the orbit where the two bodies are most separated is the apoapsis aa, though these terms vary depending on which body is orbited. Aphelion and perihelion are used for bodies orbiting the sun, and apastron and periastron are used for bodies orbiting other stars.

If e = 0, the orbit may be considered circular. If 0 < e < 1, then the orbit is eccentric. If e = 1, then the orbit is parabolic, and not bound to the orbited body. If e > 1, then the orbit is hyperbolic and is not bound to the orbited body. In the case of e ≥ 1, there exists a periapsis, but no apoapsis.

Regardless of the eccentricity of an orbit, the orbital period will remain constant for a constant value of a. However the orbital velocity will change in such a way so that a line drawn from the planet to the star will cover a constant amount of area per unit time. As such, a planet will orbit faster at periastron, and slower at apastron.

Conveniently, the eccentricity of a planet around the barycenter will be equal to the eccentricity of the star’s orbit around the barycentre. Note that if a body in a circular orbit accelerates (be it by its own power like a spacecraft, or by a gravitational pull from another body), its orbit will become more eccentric.

Inclination for the planets in our solar system is defined as the angle between the plane of the planet’s orbit and the plane of Earth’s orbit, which defines the “ecliptic.” However this is invalid for other planetary systems for which Earth is not presently orbiting to provide such a reference. As such, we use a different reference plane upon which the sky is. We view the sky as if in the centre of an all-encompassing sphere. When looking at a planetary system from the viewpoint of Earth, we may consider the background sky to be a perpendicular plane, and determine the inclination of an orbit as the angle between that body’s orbital plane and the plane of the sky.

Inclination of a planetary system

Thus an edge-on orbit around some distant star would be perpendicular to the plane of the sky, and would be considered to have an inclination of 90°, and face-on orbits are parallel to the plane of the sky and have an inclination of 0° or 180°. As with the orbital period and eccentricity, both the planet and the star will have an equal inclination (and indeed their orbital planes are identical).

One may rotate an orbit about the orbited body within the orbit plane, such that the rotation axis is perpendicular to the orbit plane. This is called the longitude of perihelion, ω, because it rotates the perihelion longitudinally around the orbit, and is measured in degrees. For an orbit with e = 0, this does not make a noticeable change in the orbit, and in such cases, it is customary to set ω = 0°. The reference angle for ω is the line between the observer and the star. There is a 180° difference between ωplanet and ωstar.

The longitude of the ascending node, Ω, is the angle around the reference plane where the orbital plane of the planet in an inclined orbit crosses the reference plane. For the specific case of the ascending node, the planet rises above the reference plane (more northward). For the specific case of when the planet dips below the reference plane, this is called the longitude of the descending node. From visual observation alone, it is impossible to discern the ascending node from the descending node for a distant planetary system. The value for Ω is equal for both the planet and the star. For non-inclined orbits we let Ω = 0°.

The mean anomaly of an orbit, M, is the position of the planet around the orbit, in degrees, and is unrelated to the shape of the orbit. It is typically given with an ephemeris, or a date (or more scientifically, ‘epoch’) for which the value of M is true. There is a 180° offset between the values of M for the star and the planet for the same epoch.

The elements of an orbit