Computing Galactic Orbits with Astropy: A data-driven ride through the cosmos.

Recently, I've been teaching myself a bit of Physics, as well as Data Science, and realized that obviously, I could learn both faster, if I could find a way to make them overlap.

The question was, of course, where to start, and where would I get this Physics data from?

After a few days of thinking, it hit me: there's obviously an absolutely massive amount of data generated whenever we do scientific experiments, measurements, or observations, in addition to the fact that we're just passively collecting data 24/7 (in terms of telescopes and detectors for gravitational events, and particles), and considering a large swathe of the Data Science ecosystem is Python based, I figured, if I can find a Python package for programming/visualizing physical events, such as collisions, or the effect of a gravitational field on some object, then they would probably dovetail relatively agreeably, and allow me to experiment within this overlap.

The next step was to consult the Great Oracle of General Learning Expertise, a.k.a "The Google", as my Graphic Design professor of exactly one semester used to refer to it, which was, in fact, an absolutely fantastic idea (who knew?), and turned up a number of amazing libraries related to both Python & Physics.

• EinsteinPy: EinsteinPy is a Python package dedicated to working on the concepts contained in General Relativity/gravitational physics.

• Astropy: The Astropy project is a community effort to create a core Python package intended to be used for all things astronomy and astrophysics, and sports a robust ecosystem of Affiliated Packages for research, data processing, and data analysis tasks related to astronomy & astrophysics.

• Poliastro: Poliastro is a Python package for solving problems related to Astrodynamics, and orbital mechanics.

I settled on Astropy, as it seemed to have the best balance of functionality, manageable learning curve, and documentation for my current needs.

So, what goals are set before us?

• We want to query the Gaia data release 2 dataset, and retrieve data about related to nearby stars.

• With this data, we will be able to delineate the high-mass, and low-mass stellar objects, based on color-magnitude.

  • High mass stars tend to burn their fuel at a higher rate than low mass stars, and high-mass stars, are also where black holes come from. The mass of stars is based on our sun. We normalize the mass of the sun to 1, called a solar mass, and then see how many of these would fit inside other stars.

• From here, we will calculate the orbits of both high, and low mass stars, within the Milky way, to observe how younger stars, tend to have higher masses, and smaller vertical excursions.

Before we get into the data, it would probably be useful to provide some background on Gaia. Gaia is a project, attempting to chart a 3D map of our galaxy, in order to learn more about it, by surveilling over one billion stars, checking each star roughly 70 times, over the course of five years.

The data we'll be working with is from the second release, and comprises roughly 22 months of observations of stars, extra-solar planets, brown dwarfs and an expected ~500k quasars.

Setting Up Our Workstation

pip install astro-gala astroquery astropy matplotlib numpy

# from astropy

import astropy.coordinates as coord

    # Provides classes for generating celestial coordinates,

    # and their velocity components.

from astropy.table import QTable

    # Stores data in a way similar to numpy (likely matrices).

import astropy.units as u

    # Handles definition, conversion and arithmetic related to physical quantity.

from astroquery.gaia import Gaia

    # The gaia data set that is attempting to map the milky way in 3D

# 3rd party imports

import matplotlib as mpl

    # Graphing library for visualizing our data.

import matplotlib.pyplot as plt

    # Plotting package building our visualization.

import numpy as np

    # Scientific computing/matrices

%matplotlib inline

import gala.coordinates as gc

    # Special versions of astropy.coordinates.

import gala.dynamics as gd

    # Contains classes and functions for expressing gravitational dynamics.

import gala.potential as gp

    # For generating potential properties related to gravity.

from gala.units import galactic

query = '''

SELECT TOP 4069 ra, dec, parallax, pmra, pmdec, radial_velocity,

phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag

FROM gaiadr2.gaia_source

WHERE parallax_over_error > 10 AND

    parallax > 10 AND

    radial_velocity IS NOT null

ORDER BY random_index

'''

job = Gaia.launch_job(query)

# This, will launch a job in the Gaia archives tasked with retrieving the contents of our query.

gaia_data = job.get_results()

gaia_data.write('gaia_data.fits')

Also, let's analyze this SQL query, both to gain clarity about the data we're asking for, and what it actually means, with respect to the real world phenomena it describes.

• SQL stands for Structured Query Language, and is a tried and true language for describing the data you require, and selecting it from a given database.

Gaia Query Parameters:

• ra: Right Ascension: This is the term for the celestial equivalent of longitude (the east-west position on the surface of a sphere). It's a bit more complicated than terrestrial longitude, as it exists on an abstract sphere in space, called the Celestial Sphere, but for all intents and purposes, we can think about it similarly.

• dec: Declination: This is the celestial equivalent of latitude, again existing on the Celestial Sphere as opposed to a an existing object.

• radial_velocity: Generally measured in km/sec, can be considered the rate of change of distance between an object and a given point.

• parallax: The tendency of objects to appear to change position as your viewpoint changes, also causes objects in the distance to appear to be moving more slowly than those up close.

Proper Motion: The proper motion of a star is the change in position on the Celestial Sphere due to the intrinsic motion of the star, rather than the apparent motion as observed from earth.

• From this we have:

  • pmra: Proper Motion in right ascension, which is one of two values that can be used to define a stars proper motion, using the ra parameter mentioned before.

  • pmdec: Proper Motion in descension, is the other value, and uses the aforementioned dec parameter.

Stellar Classification: One of the ways stars are grouped together, is by drawing a relationship between the radiation the stars give off, and the electromagnetic spectrum. This is possible, because the nature of the radiation, and thus it's "placement" on the spectrum, is dependent upon the chemical composition, temperature and size of the star.

• In this dataset, the Morgan-Keenan system is used, wherein O, B, A, F, G, K, and M are the spectral classes, and each classification, has a scale from 0-9, to indicate an intensity (typically luminosity) within these classes, however, there are other systems for classifying stars based on the same properties.

• From this we have:

  • phot_g_mean_mag: The average magnitude of the G band of a star. In this context, magnitude refers to the apparent brightness of a celestial object. The G here, refers to one of the spectral classes mentioned above. By the way, you're already quite familiar with this class of star, as a resident of planet Earth, as our home star is in fact this type of star, often more formally stated as a G-type main-sequence star.

  • phot_bp_mean_mag: The average magnitude of the Blue Photometer, which is an instrument for measuring wavelengths of light within the visible spectrum (typically wavelengths of 330-680 nanometers).

  • phot_rp_mean_mag: The average magnitude of the Red Photometer, which is the same type of instrument as mentioned directly above, except it measures wavelengths in the 640-1050 nanometer range.

Alright, so now that we understand the contents of our data set a bit more intimately, let's get started!

gaia_data = QTable.read('gaia_data.fits')

What we've done here is convert the data returned from our query, into an Astropy Table, that contains 4096 stars, within 100 parsecs of our sun.

Q: Wait a sec, whats a parsec?

B: A parsec, is a unit of distance in Astronomy, that is equivalent to ~3.26 light years.

Given 100 of these, let's figure out what this distance is in human terms:

• We need to calculate the light year equivalent of 100 parsecs.

• Then, we need to turn that number of light years into miles.

from decimal import *

getcontext().prec = 28

parsec = 3.26

star_distance = Decimal((parsec * 100))

ly_miles = Decimal(9.46073047258 * (10 ** 12))

true_distance = Decimal((ly_miles * star_distance))

true_distance

Yep, that's right, we're measuring anything within a range that spans a little over 3 quadrillion miles.

gaia_data[:4]

So, we have our data here, with our features, and the data type as well, but we need to make sure that our data meets our standards, and that there aren't any outliers that might invalidate or dilute our analysis. We'll be creating a Distance object, out of the parallaxes of our stars, to make sure all of the observations in our dataset fall within our desired range (0-100).

dist = coord.Distance(parallax=u.Quantity(gaia_data['parallax']))

dist.min(), dist.max()

At this point, we've confirmed that all of our star observations exist within the range we desire, but now, we need to convert some of these measurements from their heliocentric/spherical values to Galactocentric or Cartesian values.

• Galactocentric simply means that we use the sun as an origin, and plot a direct line from it, to the center of our galaxy, and from this we can begin building a galactic coordinate system, that is analogous to latitude and longitude on earth.

• Cartesian implies the use of the x, y axis to plot things, and representing these locations as ordered pairs of numbers, both indicating a location on their respective axes.

In order to do this, we'll need a little assistance from the astropy.coordinates package, which we aliased as coord earlier, to turn our data into a SkyCoord object.

s_c = coord.SkyCoord(ra=gaia_data['ra'],

                   dec=gaia_data['dec'],

                   distance=dist,

                   pm_ra_cosdec=gaia_data['pmra'],

                   pm_dec=gaia_data['pmdec'],

                   radial_velocity=gaia_data['radial_velocity'])

s_c[:4]

Now we've got the ability to transform our heliocentric coordinates into galactic ones, by appending the .galactic attribute to the SkyCoordinate object.

s_c.galactic[:4]

However, if we run the code below, we see that our frame is still centered on the heliocentric coordinate system, using the sun as the center point, based on the ICRS Coordinate text in the coordinate frame data. The ICRS is the International Celestial Reference System, which takes the barycenter of our solar system as the origin point (the set of coordinate points equaling zero on all available axes, or rather, the point where all axes intersect).

• Q: What's a barycenter?

• A: The barycenter is the point around which two gravitational bodies that are orbiting each other rotate, effectively, it's the center of mass of these two objects.

s_c.frame

In order to complete the heliocentric conversion, we can use the Galactocentric class, which will allow us to apply our own conventions for distances between multiple celestial objects, or the elevation of an object over the galactic plane.

s_c.galactocentric

From here, we want to update the galcen_distance from 8.3kpc, to 8.1kpc, which is a bit more accurate, and in line with measurements gained from Gravity data, and then change our solar height to 0, which we will do with the transform_to() method.

Q: What is a kpc?

• A unit of astronomical measurement equivalent to a thousand parsecs. Earlier, we calculated the distance of one hundred parsecs, as roughly 3 three quadrillion miles, so this is equivalent to ~30q.

galcen = s_c.transform_to(coord.Galactocentric(z_sun=0*u.pc,

                                            galcen_distance=8.1*u.kpc))

galcen[:4]

Since we've completed the transformation, the placement of our stars are now available to us in x, y, z, coordinates, which we can access using .x , .y, and .z

plt.hist(galcen.z.value, bins=np.linspace(-110, 110, 32))

plt.xlabel('z [pc]')

Remember, earlier we aliased Matplotlib as plt, and here we call the histogram function from that library on our data. A histogram which is a quick way to visualize data by grouping numbers into ranges, and then seeing how many observations "stack up", within each range. This allows us to make quicker observation about our data sets as a whole. For instance, with this histogram, we can see that there are very few stars that are exactly 100pc away from our sun, and that the most dense placement of stars is seems to be within 50 parsecs of our sun, in either direction.

We also have access to the velocity measurements of each of these objects, using .v_x , .v_y, and .v_z, to create a UV plane velocity plot of the velocities of these stars mapped onto the galactic plane (you can think of it like a ~disc).

The UV plane, is how we get around the inability of an interferometer to take a picture of the entirety of the sky, in one fell swoop.

• Interferometers, are used to measure the radiation produced by celestial objects, and one can reason about the distances (and other things) between two emitters (or reflectors) of (electromagnetic) radiation. They do this by splitting the of radiation into composite wavelengths using Fourier transforms, and then recombine them in a way so as to cause interference, from which we can infer information about the objects that produced the original beam.

  • This radiation, is usually a beam of light, which, of course has a wavelength, which is why Fourier Transforms are applicable here..

• We place this object in various places, and then assume an abstract plane perpendicular to the direction in which the interferometer is pointed, upon which we can map a Fourier Transform of the distribution of brightness in the sky, and measurements on this place, have a set of coordinates {u, v}, hence the "UV" in the name of this type of plot.

fig, ax = plt.subplots(1, 1, figsize=(6,6))

ax.plot(galcen.v_x.value, galcen.v_y.value, marker='*', linestyle='none', alpha=0.2)

ax.set_xlim(-125, 125)

ax.set_ylim(200-125, 200+125)

ax.set_xlabel('$v_x$ [{0:latex_inline}]'.format((u.km/u.s)))

ax.set_ylabel('$v_y$ [{0:latex_inline}]'.format(u.km/u.s))

From this plot we can appreciate the intensity of the rate of increase of distance, between us, and celestial objects, the most extreme of which, from our sample, appear to be moving outward at almost 300 km/s^-1, and rightward at a little over 100 km/s^-1.

Moving forward, remember that Morgan-Keenan system for classifying stars?

• A quick refresher: we can classify stars by analyzing the wavelengths of the light they radiate, which contains information about the star, such as chemical composition, mass, temperature, and more.

Since we have access to that information via our dataset, we can use it as well, to group the stars according to this classification, using the G band value, with the difference of the Blue and Red photometer values (which we learned about while analyzing our SQL query) as the criteria by which they will be grouped. We'll be computing the absolute magnitude, using the dist variable we made earlier, which is an object that has a property, distmod, that provides the Distance Modulus. The Distance Modulus provides a logarithmic representation of distance at astronomical scales, and these scales are based on astronomical magnitude system, where magnitude refers the brightness of celestial objects, often within a restricted segment of the light spectrum, referred to as a passband.

M_G = gaia_data['phot_g_mean_mag'] - dist.distmod

BP_RP = gaia_data['phot_bp_mean_mag'] - gaia_data['phot_rp_mean_mag']

fig, ax = plt.subplots(1, 1, figsize=(6, 6))

ax.plot(BP_RP, M_G,

        marker='.', linestyle='none', alpha=0.3)

ax.set_xlim(0, 3)

ax.set_ylim(11, 1)

ax.set_xlabel('$G_{BP}-G_{RP}$')

ax.set_ylabel('$M_G$')

hi_mass_color = 'tab:red'

lo_mass_color = 'tab:purple'

Here, we see that the majority of stars lie in the main sequence, the term scientists use to refer to an grouping of stars that tend to appear on a plot arranging stars by color(x-axis) vs brightness(y-axis). We should take care to not get caught up on a bit of confusion: the main sequence, refers to the stars on the plot, which are also Main Sequence stars, which refers to the longest part of the lifecycle of a certain class of star, and this classification is based on how the star "obtains" or "generates" the radiation it emits. In short, the main sequence is filled with a lot of Main Sequence stars.

(Boy, there sure are a lot of ways to observe/classify stars).

We can further demarcate between the high mass, and low mass stars in the main sequence, using the calculations below to set up constraints for the values of the stars, which we will then map onto colors, which will allow stars that are outliers, in terms of mass, relative to the rest of the stars in our observation, to stand out in our visualization.

np.seterr(invalid="ignore")

hi_mass_mask = ((BP_RP > 0.5*u.mag) & (BP_RP < 0.7*u.mag) &

                (M_G > 2*u.mag) & (M_G < 3.75*u.mag) & (np.abs(galcen.v_y - 220*u.km/u.s) < 50*u.km/u.s))

lo_mass_mask = ((BP_RP > 2*u.mag) & (BP_RP < 2.4*u.mag) &

                (M_G > 8.2*u.mag) & (M_G < 9.7*u.mag) & (np.abs(galcen.v_y - 220*u.km/u.s) < 50*u.km/u.s))

A quick breakdown down of the above calculation:

For selecting high mass stars, we require the difference of the two photometers to be between 0.5 * the magnitude unit (abbreviated as mu going forward), and 0.7 * the mu. (We get this unit from our astropy.units module, which we aliased as u). Additionally, we want out the average magnitude of the G band to be between 2 * the mu and 3.75 * the mu.

For selecting low mass stars, we require the difference of the two photometers to be between between 2 * the mu, and 2.4 * the mu, while we want the average magnitude of the G band to be between 8.2 * the mu and 9.7 * the mu.

fig, ax = plt.subplots(1, 1, figsize=(6, 6))

ax.plot(BP_RP, M_G,

        marker='.', linestyle='none', alpha=0.1)

for mask, color in zip([lo_mass_mask, hi_mass_mask],

                       [lo_mass_color, hi_mass_color]):

    ax.plot(BP_RP[mask], M_G[mask],

            marker='*', linestyle='none',

            alpha=0.5, color=color)

ax.set_xlim(0, 3)

ax.set_ylim(11, 1)

ax.set_xlabel('$G_{BP}-G_{RP}$')

ax.set_ylabel('$M_G$')

We can now observe that high mass stars tend to have the lowest difference in their Blue & Red photometers (you can think of this as spanning the smallest section of the light spectrum, where a larger difference covers a larger section of it), and tend to be the brightest (the absolute magnitude scale ranges from -1 (brightest), to 10 (dimmest)), whereas low mass stars have the highest difference in their photometers, and of course, tend to be the dimmest.

Integrating & Computing Galactic Orbits with gala.

Next, taking these high and low mass stars, we'll generate visualizations of their orbits, and confirm that high mass stars tend to have lower vertical excursions than low mass stars.

For computing these orbits, we need to factor in the gravity of these objects, and the greater celestial systems they exist in, as well as objects that they exist near, from galaxies, to nebulae, to other stars. We'll use the gala python package, which allows us to bring in in galactic dynamics, to account for the effect of gravity on the orbits. gala offers a simple mass model for the Milky Way that we can use as a "world" to house our orbits, and combine them into a single visual.

milky_way = gp.MilkyWayPotential()

milky_way['bulge']

milky_way['disk']

print(milky_way['halo'])

This gives us an object representing the Milky Way, consisting of a spherical nucleus, a bulge, a disk, and a a spherical NFW dark matter halo, a model defined in Bovy, 2015.

different_disk_potential = gp.MilkyWayPotential(disk=dict(m=8e10*u.Msun))

different_disk_potential

In order to get things up and running, the first step to integrating these orbits, is to combine the mass model with a reference frame, to generate a Hamiltonian object.

• Hamiltonian refers to a re-imagination of Classical Mechanics, which is equivalent to Newton's laws in describing a classical physical system, created by William Rowan Hamilton. Hamiltonian Mechanics are a reformulation of Lagrangian Mechanics, which is a reformulation of the original Classical Mechanics. It's particularly useful for describing the energy of a system over time, and we can think of an orbit as a system oscillating between peaks of both potential, and kinetic energy- somewhat like a pendulum that completes full rotations instead of half rotations.

H_mw = gp.Hamiltonian(milky_way)

Now that we have our Hamiltonian object, we need to select the criteria by which we will select stars for orbital analysis. In gala, initial conditions are specified by creating PhaseSpacePosition objects, which can be created from a Galactocentric object, like the one we created earlier, by extracting the data with Cartesian representations, using the .cartesian attribute, and appending it to the end of our galcen frame.

• Q: What is Phase Space?

• A: When dealing with dynamical systems, a phase space is the space in which all possible states (or the set of all states) of a system are represented, and each possible state corresponds to exactly one unique point in the phase space. If the dynamical system is also mechanical, the phase space usually consists of all possible values of position and momentum variables.

w0_hi = gd.PhaseSpacePosition(galcen[hi_mass_mask].cartesian)

w0_lo = gd.PhaseSpacePosition(galcen[lo_mass_mask].cartesian)

w0_hi.shape, w0_lo.shape

Here we see 183 high-mass stars, and 574 low mass stars, that have orbits we want to analyze, which we can do, by passing initial conditions (the starting state of a system), to the integrate_orbit() method on our Hamiltonian object. We also have to specify a time-step, or a time interval, as well as the length of time, or number of these intervals, we want to walk through. Our time-step will be 1 Myr, and and the length of time will be 500 Myr (roughly two revolutions around the galaxy with an elliptical orbit), and Myr, if you haven't guessed it, is short for million years.

orbits_hi = H_mw.integrate_orbit(w0_hi, dt=1*u.Myr,

                              t1=0*u.Myr, t2=500*u.Myr)

orbits_lo = H_mw.integrate_orbit(w0_lo, dt=1*u.Myr,

                              t1=0*u.Myr, t2=500*u.Myr)

With our orbital objects, we can compare the high-mass orbits, to the low-mass orbits, we can use the plot function from these orbit objects, to visualize them in a cartesian coordinate system.

fig = orbits_hi[:, 0].plot(color=hi_mass_color)

_ = orbits_lo[:, 0].plot(axes=fig.axes, color=lo_mass_color)

fig = orbits_hi[:, 0].cylindrical.plot(['rho', 'z'],

                                       color=hi_mass_color,

                                       label='high mass')

_ = orbits_lo[:, 0].cylindrical.plot(['rho', 'z'], color=lo_mass_color,

                                     axes=fig.axes,

                                     label='low mass')

fig.axes[0].legend(loc='upper left')

fig.axes[0].set_ylim(-0.3, 0.3)

At last, having plotted the vertical excursions of a high mass, and low mass star, we can observe the range of the high mass stars vertical excursion, to be well inside the range of the vertical excursion of the low mass star. We've at least confirmed our theory, but in the name of being rigorous, we must apply this to all of our stars, and see if our hypothesis is true throughout the data set.

We can use the zmax() function on our orbit object, in order to get the maximum vertical height, of all of the high-mass stars, and low-mass stars, and plot them in a histogram, where the x-axis will be "bins", that are ranges of the height of vertical (z-axis) excursions, and the y-axis will be the masses of these stars, which will make plain the correlation between masses and certain orbital properties.

zmax_hi = orbits_hi.zmax(approximate=True)

zmax_lo = orbits_lo.zmax(approximate=True)

bins = np.linspace(0, 2, 50)

plt.hist(zmax_hi.value, bins=bins,

         alpha=0.4, density=True, label='high-mass',

         color=hi_mass_color)

plt.hist(zmax_lo.value, bins=bins,

         alpha=0.4, density=True, label='low-mass',

         color=lo_mass_color);

plt.legend(loc='best', fontsize=14)

plt.yscale('log')

plt.xlabel(r"$z_{\rm max}$" + " [{0:latex}]".format(zmax_hi.unit))

At last! We see quite plainly, that the higher a stars mass is, the smaller the range of it's vertical excursion, which is in direct accordance with our hypothesis!

If you found this interesting, head on over to Astropy to learn more about the library we just used, or you can check out Gaia, the space exploration program that provided this data set, and lastly you can check out The Olympia Academy, and keep up with me and a friend/our community, as we foray boldly into the world of Physics.