Designing a Planetary System: Input

This month’s post is long and quite com­plex. So I will split it in two parts: Input & Out­put. I also decid­ed to write an Exten­sion post to the series, where I will dis­cuss bina­ry and mul­ti­ple star sys­tems, and things relat­ed.

Design­ing an extreme­ly real­is­tic plan­et (along with the whole sys­tem) is not a sim­ple task, so if you have no inter­est in it, I sug­gest you don’t read math & tech parts of this arti­cle (I will mark them as such). We’ll be doing some math here and there. Also, advanced lev­el of com­put­er knowl­edge is required. If you don’t know what a com­mand line is, skip those parts and enjoy the rest. And last but not least, this task will require a few days of your time (most­ly com­put­er time though). So, if you have some­thing to do, like mow grass on the foot­ball field, now is a per­fect oppor­tu­ni­ty. With all this said I’m mov­ing onto the main thing.

Count­less plan­e­tary sys­tems exist through­out the vast­ness of space. Some of them are small, con­sist­ing only of one plan­et per star; some of them are big, like our Solar Sys­tem with 8 plan­ets and lots of oth­er stuff orbit­ing around the Sun. In fact, all plan­e­tary sys­tems dis­cov­ered up until now con­sist of small­er num­ber of plan­ets than the Solar Sys­tem, with 2 to 4 plan­ets being most com­mon.

Before design­ing your own sys­tem it might be a good idea to look through descrip­tions of those sys­tems. There are numer­ous cat­a­logues with data. For exam­ple:

List of known exo­plan­ets

The Extra­so­lar Plan­ets Ency­clopae­dia

The NASA Star and Exo­plan­et Data­base

Exo­plan­et Tran­sit Data­base

Extra­so­lar Plan­et Table on Astropical.net

For more plan­et stuff you can also read the Sys­temic blog, which is writ­ten by Greg Laugh­lin, Pro­fes­sor of Astron­o­my and Astro­physics at the Uni­ver­si­ty of Cal­i­for­nia, San­ta Cruz. Sys­temic reports recent devel­op­ments in the field of extra­so­lar plan­ets, with a par­tic­u­lar focus on obser­va­tion­al and the­o­ret­i­cal astro­nom­i­cal research work sup­port­ed by the NASA Astro­bi­ol­o­gy Insti­tute, and by the Foun­da­tion­al Ques­tions Insti­tute (FQXi).

Now, there are two ways of get­ting your sys­tem up: top-to-bot­tom or bot­tom-to-top meth­ods.

The top-to-bot­tom method sug­gests using some of the accre­tion algo­rithms to gen­er­ate your sys­tem – you’ll get all stuff auto­mat­i­cal­ly, noth­ing to invent there. And I’m not talk­ing about “plan­e­tary sys­tem gen­er­a­tors” found on the web (e.g. Star­Gen), they don’t pro­duce real­is­tic sys­tems. They tend to pro­duce sys­tems “sim­i­lar to our own”; the truth is plan­e­tary sys­tems don’t fol­low this “sim­i­lar­i­ty” plan. A tru­ly alien world will be some­thing else, as you can see from worlds dis­cov­ered up to date. Per­son­al­ly, at some point I’ve been using AstroSyn­the­sis (v.3.0 will be released some­where after Sep­tem­ber 30th). It does cool things, try it. I rec­om­mend it to those who don’t want to go deep into cal­cu­la­tions.

State-of-the-art plan­e­tary accre­tion sim­u­la­tions typ­i­cal­ly use N-body inte­gra­tions, but the­se are too com­pu­ta­tion­al­ly expen­sive and time con­sum­ing. So we’ll skip that part and get to the time when all plan­ets are already in their respect­ful places. Trust your intu­ition; it’s bet­ter than some gen­er­a­tor.

The bot­tom-to-top method implies you first imag­ine your main plan­et and do all the nec­es­sary cal­cu­la­tions for its para­me­ters (I will talk about those in detail below). The rest of the sys­tem is rather sketchy: all you decide is where to put oth­er plan­ets (semi­ma­jor axes) and how many there will be in total (you can also add some of their details as well, e.g. mass, den­si­ty and orbital ele­ments). This deci­sion should be based on your star of choice and its para­me­ters.

Plan­e­tary para­me­ters

I assume you have already decid­ed on mass, den­si­ty and cal­cu­lat­ed radius. If not, check the pre­vi­ous post here.

Anoth­er two para­me­ters required for our com­pu­ta­tions are axi­al tilt and rota­tion rate. The tilt can be any­thing you want from 0 to 180°. Note that object with tilt below 90° has a ret­ro­grade rota­tion, like Venus does. Also, plan­ets usu­al­ly move on their orbits in coun­ter­clock­wise direc­tion, with clock­wise direc­tion being ret­ro­grade orbital motion.

There are three cas­es of plan­e­tary rota­tion around its axis: nor­mal (coun­ter­clock­wise), ret­ro­grade and tidal­ly locked (syn­chro­nous rota­tion or orbital res­o­nance, e.g. Mercury’s 3:2 spin-orbit res­o­nance).

When plan­ets are born, they have their pri­mor­dial spins. Dur­ing their life­times plan­e­tary spins are affect­ed by solar tides and tides inflict­ed by anoth­er bod­ies of the sys­tem. Plan­ets and moons slow down, spin up (like Mars and Sat­urn did), or become tidal­ly locked if too close to a star or any oth­er mas­sive body (e.g. Mer­cury, the Moon). Earth’s ini­tial spin rate was about 6.5 hours com­par­ing to mod­ern 24 hour day.

Math Box 1 – Ini­tial rota­tion of a plan­et

Approx­i­mate ini­tial plan­e­tary spin can be esti­mat­ed from its mass.

First, you cal­cu­late ini­tial angu­lar veloc­i­ty using Dole’s equa­tion (out­put is in rad/s). M is the mass of the plan­et; k is the moment of iner­tia ratio (I/MR^2), a val­ue of 0.33 is tak­en for a ter­res­tri­al plan­et (you can look at plan­et moments of iner­tia in the­se fact sheets); r is planet’s radius, and j is 1.46E-20 m²/s²kg.

ω = [(2*j*M) / (k*r^2)] ^ 0.5

Then, ini­tial rota­tion rate in sec­onds is cal­cu­lat­ed as t = (2*PI)/ω

To get it in hours th = t/3600

Stud­ies sug­gest that pri­mor­dial spins are ran­dom­ly deter­mined by a few giant impacts dur­ing accre­tion, e.g. Agnor et al. (1999). So pri­mor­dial spin might be what­ev­er you want (the bound­ary is the veloc­i­ty at which the plan­et will fall apart).

Math Box 2 – Plan­e­tary rota­tion at a cer­tain age

Now, from this point onward things start to slow down, because a tidal decel­er­a­tion torque from the pri­ma­ry is applied. (On some occa­sions things might actu­al­ly spin up.)

» dω/dt = (dωe/dt)*(ke/kp)*(rp/re)*(me/mp)*((Mstar/Msol)^2)*((Rsol/Rstar)^6)

Here (dωe/dt) = (-1.3E-6 /109) rad s-1/yrs; k is as men­tioned above for ter­res­tri­als; m is a plan­e­tary mass; index let­ter e stands for Earth and p for plan­et.

You can esti­mate plan­e­tary angu­lar veloc­i­ty at any age of the sys­tem. In the fol­low­ing for­mu­la sys­tem age is in sec­onds. The new angu­lar veloc­i­ty then becomes

ωnew = Ini­tial angu­lar veloc­i­ty ω + (sys­tem age (s) * tidal decel­er­a­tion torque dω/dt)

Rota­tion rate in sec­onds is cal­cu­lat­ed » t = (2*PI)/ωnew

In hours » th = t/3600

Thus you can deter­mine how your plan­et day becomes longer. How­ev­er, the pri­ma­ry isn’t the only one influ­enc­ing rota­tion rate of a plan­et. There’s also atmos­pher­ic drag, and a small influ­ence of the system’s bod­ies. (Some of the for­mu­las above are tak­en from Fogg, 1992.)

Note: If you have a moon or two, things get more com­pli­cat­ed. You’ll have to take them into account as well.

Note 2: I sug­gest you cal­cu­late rota­tion rate this way instead of inte­grat­ing the whole sys­tem for con­sid­er­able amount of time; that could take days on an aver­age machine, espe­cial­ly if pre­ci­sion and long times­pan (more than 100×108 to X × 109 years) are cho­sen.

Angu­lar veloc­i­ty also gives us a clue to plan­e­tary flat­ten­ing and dif­fer­ence between polar and equa­to­ri­al radii.

Math Box 3 – Plan­et flat­ten­ing (dynam­i­cal ellip­tic­i­ty) and radii

The radius you deter­mined from mass and den­si­ty is the mean radius, R. If not, the time is now:

R=0.5*[(mass/density)*(6/PI)]^(1/3)

or, if given mass and log g, through rela­tion
Mstar/Msun=((Rstar/Rsun)^2)*10^(loggstar-loggsun).

Ellip­tic­i­ty e = SQRT(f*(2-f)) is relat­ed to flat­ten­ing, which can then be expressed as f = 1-SQRT(1-e^2) or f = (1/SQRT(1-e^2))-1; or, if you have only radii, flat­ten­ing becomes f=(Re-Rp)/R.

Flat­ten­ing f is also cal­cu­lat­ed as f* = (5/4)*(((ω^2)*(R^3))/(M*G)); ω — from cal­cu­la­tions above (math box 1 & 2); R – mean radius; M — plan­e­tary mass; G = 6.67259E-11. Equa­to­ri­al radius thus becomes Re = R*(1+(f/3)), and polar Rp = R*(1-((2*f)/3)). The dif­fer­ence between radii is expressed as Re-Rp, or sim­ply as mean R*f.

*Note: We’ll get some­what over­es­ti­mat­ed degree of rota­tion­al flat­ten­ing of the plan­et because the for­mu­la takes planet’s inte­ri­or as incom­press­ible flu­id of uni­form den­si­ty. In real­i­ty ter­res­tri­al planet’s core is much denser than its crust. Or, in case of gaseous plan­ets, an over­es­ti­mate will rise from the fact that they have a mass dis­tri­b­u­tion which is strong­ly con­cen­trat­ed at their cores.

Now, when we have our rota­tion rates, we can move onto the next step.

Orbits and orbital ele­ments

The tra­di­tion­al orbital ele­ments are the Kep­le­ri­an ele­ments. They describe the orbit of the body in space.

Kep­le­ri­an ele­ments can be obtained from and con­vert­ed to orbital state vec­tors (x-y-z coor­di­nates for posi­tion and veloc­i­ty) by man­u­al trans­for­ma­tions or with com­put­er soft­ware. For our inte­gra­tion Kep­le­ri­an ele­ments should be enough. If you want some­thing else, the body-cen­tered, iner­tial rec­tan­gu­lar com­po­nents of the radius vec­tor can be deter­mined from the clas­si­cal orbital ele­ments (see math box 4).

Math Box 4 – Orbital State Vec­tors

The mean lon­gi­tude λ = M + ω + Ω where M is the mean anom­aly, ω is the argu­ment of perigee, Ωis the lon­gi­tude of the ascend­ing node.

The true anom­aly ν is found from the expres­sion ν=λ-Ω-ω

The body-cen­tered, iner­tial rec­tan­gu­lar com­po­nents of the radius vec­tor can be deter­mined from the clas­si­cal orbital ele­ments as fol­lows:

Here a is the semi­ma­jor axis, e is the eccen­tric­i­ty, i is the incli­na­tion; in the­se equa­tions p is called the semi­pa­ra­me­ter of the orbit and is cal­cu­lat­ed from p=a*(1-e^2). η is the grav­i­ta­tion­al con­stant of the pri­ma­ry or cen­tral body.

Math Box 5 — The spin angu­lar momen­tum*

*this will be nec­es­sary for the main inte­gra­tion part to cal­cu­late obliq­ui­ties and spin rates. If you won’t do inte­gra­tion, skip this part.

SAM is a vec­tor, it has 3 com­po­nents. It is cal­cu­lat­ed as H = H * Ĥ

H = 0.4 * M * (R^2) * v

Here, R is planet’s equa­to­ri­al radius, M is plan­e­tary mass, and v is a spin rate (angu­lar veloc­i­ty).

The ini­tial coor­di­nates of unit vec­tor Ĥ in the ref­er­ence frame are:

X = sin ψ * sin ε
Y = cos ψ * sin ε
Z = cos ε

Where ψ is pre­ces­sion in lon­gi­tude and ε is axi­al tilt of a plan­et.

The pre­ces­sion in lon­gi­tude ψ is defined by ψ = L — Ω, where Ω is the lon­gi­tude of the ascend­ing node, and L is the incli­na­tion of eclip­tic. The angle between the planet’s equato­rial plane and incli­na­tion of eclip­tic is the obliq­ui­ty ε.

Your both pre­ces­sion angle ψ and obliq­ui­ty ε should be in radi­ans.

**You can cal­cu­late your pre­ces­sion rate as well. Note that pre­ces­sion for­mu­las use semi­ma­jor axis (a) and eccen­tric­i­ty (e). You can start inte­gra­tion at the begin­ning of the sys­tem age (if you have time and patience) to get all things neat­ly; or you can skip to the nec­es­sary age and use the “invent­ed” data as ini­tial point. What­ev­er works for you.

The final step is to mul­ti­ply each com­po­nent of Ĥ by the quan­ti­ty H.

Hx = H * X
Hy = H * Y
Hz = H * Z

Inte­gra­tion will require the­se quan­ti­ties in AU^2/day, so take mass as Mass/Mass of the star, R of a plan­et as R/149597870.7, and angu­lar veloc­i­ty in radians/day.

TECH TOOLBOX: Envi­ron­ment and Inte­gra­tor

Note: Numer­ous pro­grams exist to per­form inte­gra­tions of n-body sys­tems. If you want to learn more, google some­thing like “n-body inte­gra­tor”.

Cur­rent­ly I’m using Mer­cury, a hybrid sym­plec­tic inte­gra­tor for orbital dynam­ics devel­oped by John Cham­bers. A cor­rect­ed ver­sion of the main file can be down­load­ed here.

Mer­cury is writ­ten in For­tran77 in a slight­ly non-stan­dard code and before you can use it, you have to com­pile it. So, first of all, you need a com­pil­er. I used GNU gfor­tran for this pur­pose. Depend­ing on your oper­at­ing sys­tem there are some options:

Installing Cyg­win with com­pil­er packs (Win­dows)

Assum­ing you installed all nec­es­sary things, we now can move onto our main goal – Mer­cury pack­age.

The pack­age has the man­u­al, the file called mercury6.man and can be opened with any text edi­tor of your choice. Read it CAREFULLY and fol­low all the instruc­tions.

If using gfor­tran, com­pile the dri­vers typ­ing lines
gfor­tran -o mer­cury6 mercury6_2.for
gfor­tran -o ele­men­t6 element6.for
gfor­tran -o close6 close6.for

To run what you have com­piled, just type
./mercury6

and hit enter.

Pret­ty sim­ple.

If you have ques­tions regard­ing instal­la­tion and usage, feel free to ask. But please note that I’m run­ning all this under Win­dows XP and might not know about any issues with oth­er oper­at­ing sys­tems.

So far I had no issues with any soft­ware men­tioned here.

Inte­gra­tion of orbits

When inte­grat­ing the equa­tions of motion in clas­si­cal or rel­a­tivis­tic celes­tial mechan­ics, one has to know the val­ues of the para­me­ters, in par­tic­u­lar the mass­es and the ini­tial con­di­tions (orbital para­me­ters or state vec­tors). I won’t be going into the n-body prob­lem details here. If you are inter­est­ed I sug­gest read­ing some good lit­er­a­ture on the sub­ject, e.g. Grav­i­ta­tion­al N-body Sim­u­la­tions: Tools and Algo­rithms by Sver­re J. Aarseth (a cheap copy can be found on abebooks.com or on any of their inter­na­tion­al sites).

So far all nec­es­sary info on required para­me­ters is in the Mer­cury man­u­al.

Use­ful Julian date con­vert­er can be found here.

Plan­e­tary sys­tem sta­bil­i­ty

Are plan­e­tary sys­tems sta­ble? Prob­a­bly. You can’t tell just by look­ing at them; the invis­i­ble force out there guides all the inter­ac­tions between the bod­ies. That’s why sci­en­tists per­form numer­i­cal inte­gra­tions of plan­e­tary orbits.

#For the next part I’ll do a small sam­ple sys­tem inte­gra­tion.

Jeno Marz
JENO MARZ is a science fiction writer from Latvia, Northern Europe, with background in electronics engineering and computer science. She is the author of two serial novels, Falaha’s Journey: A Spacegirl’s Account in Three Movements and Falaha’s Journey into Pleasure. Marz is current at work on a new SF trilogy. All her fiction is aimed at an adult audience.