Wolfram Computation Meets Knowledge

Exploring Exoplanet Systems with the Wolfram Language

Sun, TRAPPIST, Jupiter and Mercury with exoplanet orbit

Exoplanets are currently an active area of research in astronomy. In the past few years, the number of exoplanet discoveries has exploded, mainly as the result of the Kepler mission to survey eclipsing exoplanet systems. But Kepler isn’t the only exoplanet study mission going on. For example, the TRAnsiting Planets and PlanetesImals Small Telescope (TRAPPIST) studies its own set of targets. In fact, the media recently focused on an exoplanet system orbiting an obscure star known as TRAPPIST-1. As an introduction to exoplanet systems, we’ll explore TRAPPIST-1 and its system of exoplanets using the Wolfram Language.

To familiarize yourself with the TRAPPIST-1 system, it helps to start with the host star itself, TRAPPIST-1. Imagine placing the Sun, TRAPPIST-1 and Jupiter alongside one another on a table. How would their sizes compare? The following provides a nice piece of eye candy that lets you see how small TRAPPIST-1 is compared to our Sun. It’s actually only a bit bigger than Jupiter.

estarrad = QuantityMagnitude[["Radius"], "Kilometers"];

Although the diameter looks to be about the same as Jupiter’s, its mass is quite different—actually about 80 times the mass of Jupiter.

UnitConvert[["Mass"], "JupiterMass"]

And it has only about 8% of the Sun’s mass.

UnitConvert[["Mass"], "SolarMass"]

TRAPPIST-1 is a thus very low-mass star, at the very edge of the main sequence, but still allowing the usual hydrogen fusion in its core.

The exoplanets in this system are what actually gained all of the media attention. All of the exoplanets (blue orbits) found in the TRAPPIST-1 system so far orbit the star at distances that would be far inside the orbit of Mercury (in green), if they were in our solar system.

sats = Entity["Star", "TRAPPISTMinus1"]["Satellites"];

As a more quantitative approach to study the planets in this system, it is useful to take a look at the orbital periods of these planets, which lie very close together. Planets in such close proximity can often perturb one another, which can result in planets being ejected out of the system, unless some orbital resonances ensure that the planets are never in the wrong place at the wrong time. It’s easy to look up the orbital period of the TRAPPIST-1 planets.

pers = QuantityMagnitude[EntityValue[sats, "OrbitPeriod"], "Hours"]

Divide them all by the orbital period of the first exoplanet to look for orbital resonances, as indicated by ratios close to rational fractions.


These show near resonances with the following ratios.

{24/24., 24/15., 24/9., 24/6., 24/4., 24/3.}

TRAPPIST-1 h has an inaccurately known orbital period so it’s not clear whether it partakes in any resonances.

Similarly, nearest-neighbor orbital period ratios show resonances.

Table[pers[[i + 1]]/pers[[i]], {i, 1, Length[pers] - 1}]

Which are close to:

{8/5., 5/3., 3/2., 3/2., 4/3.}

An orbital resonance of 3/2 means that one of the planets orbits 3 times for every 2 of the other. Pluto and Neptune in our solar system exhibit a near 3:2 orbital resonance.

Rationalize[  Entity["MinorPlanet", "Pluto"]["OrbitPeriod"]/   Entity["Planet", "Neptune"]["OrbitPeriod"], .01]

This can help explain how so many planets can be packed into such a tight space without experiencing disruptive perturbations.

What about the distances of the exoplanets from their host star? If you placed TRAPPIST-1 and its planets alongside Jupiter and its four Galilean moons, how would they compare? The star and Jupiter are similar in size. The exoplanets are a bit larger than the moons (which are hard to see here) and they orbit a bit farther away, but the overall scales are of similar magnitude. In the following graphic, all distances are to scale, but we magnified the actual diameters of the planets and moons to make them easier to see.


The sizes of the planets can be compared to Jupiter’s four biggest moons for additional scale comparisons.

exoSizeGraphic[{rad_, nam_}] :=


Exploring Other Exoplanet Systems

At the time of this writing, we have curated over 3,400 confirmed exoplanets:


Most of the confirmed exoplanets have been discovered since 2014, during missions such as Kepler.

dydata = Gather[    ExoplanetData[EntityClass["Exoplanet", All], "DiscoveryYear"]];

With[{data =     Sort[Join[{#[[1]]["Year"], Length[#]} & /@        dydata, {{1990, 0}, {1991, 0}}]]},   BarChart[data[[All, 2]],    ChartLabels -> Placed[data[[All, 1]], Axis, Rotate[#, 90 Degree] &],    PlotLabel -> "number of exoplanet discoveries per year"]]

You can query for data on individual exoplanets.

emass = EntityValue[Entity["Exoplanet", "Kepler11e"], "Mass"]

emass/EntityValue[Entity["Planet", "Earth"], "Mass"]

In addition, there are various classifications of exoplanets that can be queried.


You can perform an analysis to see when the exoplanets were discovered.

dydata = Gather[    EntityValue[EntityClass["Exoplanet", "SuperEarth"],      "DiscoveryYear"]];

With[{data =     Sort[Join[      Map[Function[        yearlist, {First[yearlist]["Year"], Length[yearlist]}],        dydata], {{1990, 0}, {1991, 0}}]]},   BarChart[data[[All, 2]],    ChartLabels ->     Placed[data[[All, 1]], Axis,      Function[angle, Rotate[angle, 90 Degree]]],    PlotLabel -> "number of super\[Hyphen]Earth discoveries per year"]]

You can also do a systematic comparison of exoplanet parameters, which we limit here to the radius and density. We are only considering the entity class of super-Earths here. The red circle marks the approximate location of the TRAPPIST-1 system in this plot.

rddata = EntityValue[    EntityClass["Exoplanet", "SuperEarth"], {"Radius", "Density"}];

ListPlot[rddata, Frame -> True,   FrameLabel -> {"radius (Earth radii)",     Row[{"density (g/", Superscript["cm", 3], ")"}]},   TargetUnits -> {"EarthEquatorialRadius", "Grams"/"Centimeters"^3},   Epilog -> {Red, Circle[{.951, 4.9}, {.178, 1.6}]}]

Here is another example of systematic comparison of exoplanet parameters by discovery method, indicated by color coding. Once again, the TRAPPIST-1 system is shown, with red dots at its mean values.

pmdata = Reverse@    SortBy[GroupBy[      Select[EntityValue[        EntityClass["Exoplanet", "SuperEarth"], {"OrbitPeriod", "Mass",          "DiscoveryMethod"}],        Function[val, ! MatchQ[val[[3]], Missing["NotAvailable"]]]],       Last -> Most], Function[data, Length[data]]];

ListLogLogPlot[pmdata, Frame -> True,   FrameLabel -> {"period (days)", "mass (Earth masses)"},   PlotStyle -> PointSize[Medium],   PlotLegends ->    Placed[PointLegend[Keys[pmdata], LegendLabel -> "discovery method",      LegendFunction -> "Frame", LegendLayout -> "Row"], Bottom],   TargetUnits -> {"Days", "EarthMass"},   Epilog -> {Red, PointSize[.02], Point[{Log[24.32], Log[.9]}]}]

Atmospheric Modeling

In addition to data specific to exoplanets, the Wolfram Language also includes data on chemical compounds present in planetary atmospheres.

ThermodynamicData["Methane", {"Pressure" ->     Quantity[0.2, "Megapascals"],    "Temperature" -> Quantity[300, "Kelvins"]}]

You can use this data to show, for example, how the density of various atmospheric components changes with temperature.

temps = Table[Quantity[t, "Kelvins"], {t, 400, 900, 10}]; methanedensities =    ThermodynamicData["Methane",     "Density", {"Pressure" -> Quantity[0.2, "Megapascals"],      "Temperature" -> temps}]; co2densities =    ThermodynamicData["CarbonDioxide",     "Density", {"Pressure" -> Quantity[0.2, "Megapascals"],      "Temperature" -> temps}]; codensities =    ThermodynamicData["CarbonMonoxide",     "Density", {"Pressure" -> Quantity[0.2, "Megapascals"],      "Temperature" -> temps}]; h20densities =    ThermodynamicData["Water",     "Density", {"Pressure" -> Quantity[0.2, "Megapascals"],      "Temperature" -> temps}]; ListLinePlot[{Transpose[{temps, methanedensities}],    Transpose[{temps, co2densities}], Transpose[{temps, codensities}],    Transpose[{temps, h20densities}]}, AxesLabel -> Automatic]

Photometric Light Curves

As a more concrete application, the Wolfram Language also provides the tools needed to explore collections of raw data. For example, here we can import irregularly sampled stellar light-curve data directly from the NASA Exoplanet Archive for the HAT-P-7 exoplanet system.

rawimport =    Import["http://exoplanetarchive.ipac.caltech.edu:80/data/\ ETSS//KeplerDV/000/866/60/kplr010666592_q1_q16_tce_01_dvt_lc.tbl",     "Table"];

Then we can remove data points that are not computable.

exodata = Cases[rawimport[[36 ;;, {1, 3}]], {_?NumberQ, _?NumberQ}];

This raw data can be easily visualized, as shown here.

ListPlot[exodata, PlotStyle -> PointSize[Tiny], PlotRange -> All]

The following subsamples the data and does some additional post processing to both the times and magnitudes.

exodataDS = exodata[[1 ;; -1 ;; 100]];

times = exodataDS[[All, 1]] - First[exodataDS[[All, 1]]]; values = exodataDS[[All, 2]] - Mean[exodataDS[[All, 2]]];

Plotting the data shows evidence of eclipses, appearing as a smattering of points below the dense band of data.

ListPlot[Transpose[{times, values}], PlotRange -> All]

Deriving an Orbital Period—Method 1

A Fourier-like algorithm can identify periodicities in this data.


Plot[periodogramD[w, values, times], {w, .01, 2}, PlotPoints -> 500,  PlotRange -> All,   AxesLabel -> {"\[Omega] (radians/days)",     Row[{Subscript["P", "D"], "(\[Omega])"}]}]

Zoom into the fundamental frequency, at higher resolution.

Plot[periodogramD[w, values, times], {w, .453, .454},   PlotPoints -> 400,  PlotRange -> All,   AxesLabel -> {"\[Omega] (radians/days)",     Row[{Subscript["P", "D"], "(\[Omega])"}]}]

We find that the fundamental peak frequency occurs at .453571 radians/day, for which the reciprocal gives an estimate of the corresponding orbital period in days.

1/NArgMax[{periodogramD[w, values, times], .453 < w < .454}, w]

Deriving an Orbital Period—Method 2

With some additional processing, we can apply a minimum string length (MSL) algorithm to the raw data to look for periodicities.

stringLengths[per_?NumberQ, vals_, times_] :=   Module[{foldedtimes = FractionalPart[times/per], diffs},    diffs = Differences[vals[[Ordering[foldedtimes]]]]^2/(1 +        Differences[foldedtimes]^2);   Total[diffs]]

We can apply the MSL algorithm to a range of potential periods to try to find a value that minimizes the distance between neighboring points when the data is phase folded.

Plot[stringLengths[per, values, times], {per, 2.204, 2.2055},   PlotPoints -> 400, PlotRange -> All,   AxesLabel -> {"period (days)", "string length"}]

Clearly, the minimum string length occurs at about 2.20474 days, in close agreement with method 1.

period = NArgMin[   stringLengths[per, values, times], {per, 2.204, 2.2055}]

We can also validate this derived value with the value stored in the header of the original data.

StringCases[  StringJoin[((Function[row, StringJoin[ToString /@ row]] /@       rawimport[[1 ;; 35]]))],   ShortestMatch["\\ORBITALPERIOD=" ~~ val__ ~~ "\\"] :> val]

This orbital period corresponds with that of exoplanet HAT-P-7b, as can be seen in the Wolfram curated data collection (complete with units and precision).

UnitConvert[  EntityValue[Entity["Exoplanet", "HATP7b"], "OrbitPeriod"], "Days"]

From the known orbital period, we can phase fold the original dataset, overlapping the separate eclipses, to obtain a more complete picture of the exoplanet eclipse.

foldTimes[times_, period_] := FractionalPart[times/period]

ListPlot[Sort[   Transpose[{foldTimes[exodata[[All, 1]], period],      exodata[[All, 2]]}]], Frame -> True, Axes -> False,  FrameLabel -> {"time (days) modulo period", "amplitude fluctuation"},   PlotRange -> {All, {-.007, .001}}, PlotStyle -> PointSize[0]]

Noise can be reduced by carrying out a phase-binning technique. All data points are placed into bins of width .0005 days, and the mean of the values in each bin is determined.

bin = With[{binwidth = .0005},     Sort[Function[       binneddata, {Round[binneddata[[1, 1]], binwidth],         Mean[binneddata[[All, 2]]]}] /@       GatherBy[       Transpose[{foldTimes[exodata[[All, 1]], period],          exodata[[All, 2]]}],        Function[val, Round[val[[1]], binwidth]]]]];

lc = ListPlot[bin, Frame -> False, Axes -> True,    FrameLabel -> {"time (days) modulo period",      "amplitude fluctuation"}, PlotRange -> All, Joined -> True,    AxesOrigin -> {0, 0}]

Modeling Light Curves

This graphic, mainly for purposes of visualization, shows the host star, HAT-P-7, with its exoplanet HAT-P-7b orbiting it. All parameters, including diameter and orbital radius, are to scale. The idea is to try to reproduce the brightness variations seen in the observed light curve. For this graphic, the GrayLevel of the exoplanet is set to GrayLevel[1], which enables you to more clearly see the exoplanet go though phases as it orbits the host star.

With[{prad =     QuantityMagnitude[1.363*Entity["Planet", "Jupiter"]["Radius"],      "Kilometers"],    psemi = QuantityMagnitude[Quantity[0.0377, "AstronomicalUnit"] ,      "Kilometers"],    srad = QuantityMagnitude[Quantity[1.84, "SolarRadius"],      "Kilometers"]},   Show[Table[    Graphics3D[{{Lighting -> {{"Ambient",           GrayLevel[.3]}, {"Directional", GrayLevel[1],           ImageScaled[{0, 0, 2}]}}, GrayLevel[1],        Sphere[{0, 0, 0},         srad]}, {Lighting -> {{"Ambient", GrayLevel[0]}, {"Point",           GrayLevel[1], {0, 0, 0}}}, GrayLevel[1],        Sphere[{psemi Cos[t], psemi Sin[t], 0}, prad]}},      ViewPoint -> {0, 10000, 0}, Boxed -> False,      Background -> GrayLevel[0],      PlotRange -> {1.05 psemi {-1, 1}, 1.05 psemi {-1, 1},        srad 1.05 {-1, 1}}, ImageSize -> 600], {t, 0, 2 Pi, 2 Pi/63}]]]

Now we can do an analogous thing, generating a list of frames instead of a static graphic. In this case, the GrayLevel of the exoplanet is much reduced, as compared to the animation above. For purposes of illustration and to reduce computation time, a small set of values has been chosen around the primary eclipse.

frames = With[{prad =       QuantityMagnitude[1.363*Entity["Planet", "Jupiter"]["Radius"],        "Kilometers"],      psemi = QuantityMagnitude[Quantity[0.0377, "AstronomicalUnit"] ,        "Kilometers"],      srad = QuantityMagnitude[Quantity[1.84, "SolarRadius"],        "Kilometers"]},     Table[<|"time" -> t,       "graphic" ->        Graphics3D[{{Lighting -> {{"Ambient",              GrayLevel[0.01]}, {"Directional", GrayLevel[1],              ImageScaled[{0, 0, 2}]}}, GrayLevel[1],           Sphere[{0, 0, 0},            srad]}, {Lighting -> {{"Ambient", GrayLevel[0]}, {"Point",              GrayLevel[1], {0, 0, 0}}}, GrayLevel[.03],           Sphere[{psemi Cos[t], psemi Sin[t], 0}, prad], Opacity[.3],           Sphere[{psemi Cos[t], psemi Sin[t], 0}, 2 prad]}},         ViewPoint -> {0, 10000, 0}, Boxed -> False,         PlotRange -> {1.05 psemi {-1, 1}, 1.05 psemi {-1, 1},           prad 1.05 {-1, 1}}, ImageSize -> 1024,         Background -> GrayLevel[0]]|>, {t, {0 + 1.25,        0.0125 2 Pi + 1.25, 0.02 2 Pi + 1.25, 0.025 2 Pi + 1.25,        0.035 2 Pi + 1.25, 0.05 2 Pi + 1.25, 0.06 2 Pi + 1.25,        0.075 2 Pi + 1.25, 0.08 2 Pi + 1.25,        0.09 2 Pi + 1.25, .1 2 Pi + 1.25,        0.125 2 Pi + 1.25, .15 2 Pi + 1.25}}]];

Now, to measure how the brightness of the scene changes, we can use image processing to total all of the pixel values. It’s rasterized at a large image size so that edge artifacts are minimized (which can otherwise have measurable effects on the resulting light curve). So this code takes a minute or so to run.

pixelcounts =    Map[Function[     pair, {(pair["time"] - 1.25)/(2 Pi),       Total[Flatten[        ImageData[ImageResize[Image[pair["graphic"]], 1024]]]]}],     frames];

Next, we rescale all of the pixel counts to fit in the same vertical range as the observed light curve.

mm = MinMax[pixelcounts[[All, 2]]]

rescaledpixels =    Function[data, {data[[1]],       Rescale[data[[2]], mm, {-0.0067415, 0.000098613}]}] /@     pixelcounts;

Now compare the model data to the actual data. The red points show the model data computed at a few orbital phases around the primary eclipse.

Show[lc, ListPlot[rescaledpixels, PlotStyle -> Red], ImageSize -> 400]

A more detailed model light curve can be constructed if you increase the computation time. The version above was done for speed. Of course, additional secondary effects can be included, such as the possibility of gravity darkening and other effects that cause brightness variations across the face of the star. Such secondary effects are beyond the scope of this blog.

Other star systems can be far more complicated and provide their own unique challenges to understanding their dynamical behavior. The Wolfram Language provides a powerful tool that allows you to explore the subtleties of stellar light curve analysis as well as the periodicities in irregularly sampled data. It would be interesting to see some of these more complicated systems tackled in similar ways to what we’ve done in this blog.

Thanks to Daniel Lichtblau for his insights into exploring periodicity in irregularly sampled data.

Download this post as a Computable Document Format (CDF) file. New to CDF? Get your copy for free with this one-time download.


Join the discussion

!Please enter your comment (at least 5 characters).

!Please enter your name.

!Please enter a valid email address.