Finding the Center of Mass of the Milky Way

Over the summer, I took Multivariable Calculus at a local college. For my final project, I was asked to explain the concept of the Center of Mass and include some contemporary examples. One of my examples was attempting to calculate the center of mass of the Milky Way using open-source data. I hope you learn as much as I did during this rewarding and informative project.

This project is relatively simple in terms of supplies. I used:

1. A Macbook
2. Matlab
3. Open-source probe data
4. A basic sense of coding

In modern mathematics, the center of mass is a critical topic in multivariable calculus. It is a widespread application of multiple integrals. To find the coordinates of the center of mass of an object, first, we need to find the moment around both the x-axis and the y-axis. A moment explains how a physical quantity, for example, mass, is spread around an object. In 2D, the moment is found by taking the product of mass and multiplying by the distance from the axis. In 3D, both x and y have their moments. We also need the mass of the object. These result in the following formula.

With a constant density, these two equations give the centroid of the object. The centroid is the geometric center of an object. In the real world, this would be an object's balancing point. However, we always assume that the density is uneven; this skews the center of mass towards one side of the object. We write this change in density as p(x,y); this function illustrates that density might change depending on where you are in the object.

The next step involves double integrals. If density is constantly changing, we cannot use area times constant density. Instead, we split the object into lots of tiny shapes and use double integrals to sum them. The equation below is used to find the total mass with an inconsistent density.

In this equation, R is the region or shape we are looking at, dA is the tiny piece of the object, and p(x,y) is the density at dA. Once we have the total mass, we can utilize the moments we found earlier. This formula will be different as we are using a non-constant density.

You can see that the total mass is multiplied by the height(y), and then the double integral is applied. This is repeated for the y-axis, with it being multiplied by length(x) instead of height. After we have found the total mass and moments, we can use the equation from earlier (Figure 1) to determine the center of mass.

A modern example of the center of mass is the attempt to map the center of mass of the Milky Way galaxy using available data. This calculation will be skewed away from the supermassive black hole located in the center of the galaxy. This is an approximation using publicly available data; unfortunately, building and launching my probe is outside the scope of this class. I used data from the Gaia probe mission. This was a probe launched in 2013 by the European Space Agency (ESA). Its mission and goal were to create a 3D map of the Milky Way. It mapped millions of stars with their position, distances, and motions. A section of this data is in the public domain, so I began with that.

The Gaia Space Observatory was launched by the ESA (European Space Agency) on December 19, 2013. Its primary purpose was to create a 3D map of the Milky Way. I obtained the Gaia data from the Gaia Data Release 3 (DR3), which is part of a series of data releases published by the ESA. The data I used specifically came from the Gaia Archive, which is ESA’s official data portal for Gaia. DR3 contains precise measurements, for example, star positions (RA/Dec), parallax (used to calculate distance), proper motions (star velocities across the sky), and magnitudes in multiple bands (brightness). In order to access certain parts of DR3, I needed to utilize an ADQL query, which was a learning experience for me.

Anyone can access the Gaia data for free.

1. Go to the Gaia Archive: https://gea.esac.esa.int/archive/
2. Search for Stars or Regions:
3. You can manually search by star name or coordinates (RA/Dec).
4. You can also upload or create ADQL (Astronomical Data Query Language) queries to download bulk data.
5. Example ADQL Query:

To get basic data from a patch of sky:

This gets 1000 stars between 250–251° Right Ascension (RA) and 36–37° Declination (Dec).

I then downloaded the data I requested as a .csv file, which I was taught is the preferred file type for MATLAB.

Next, it was onto the coding, but first, the concepts.

My MATLAB code takes this input of data and separates four important data points: Right Ascension (RA), Declination, Parallax, and G-band frequency. Right ascension and declination are the equivalents to longitude and latitude, but for stellar bodies. Right ascension is how far east of the Point of Aries a star lies. The Point of Aries is where the Sun crosses the celestial equator of the Earth during the March Equinox. It is (0,0) in ecliptic coordinates.

Right ascension is then measured in times (day, month, hour). Therefore, if a star is +30 RA, it is thirty hours east of the Point of Aries. The counterpart to RA is declination. Declination is how far north or south of the celestial equator the stars lie. The celestial equator is the entire universe on the same plane as the Earth’s equator. Declination is measured in degrees; a +30 Dec would be a star located 30 degrees north of the plane of the Earth’s equator. The next data point I gather is the Parallax of the star. Parallax is how far the “close stars” differ in apparent position. A good way to visualize this is to put your finger up to your nose and close one eye. When you open it, your finger “moves”. That difference in position is parallax. The final piece of data is the G-band frequency. This is the relative brightness of the star. It helps us determine how large a star is.

The first thing my program does after assigning those points to variables is exclude any unhelpful data. The excluded data includes any negative parallax, as distance cannot be negative, and any points with missing G-band data. Next, it converts parallax from millarcseconds to parsecs. This is a transfer from time to distance done by dividing by one thousand. Distance will help me create a 3D map. The next unit conversion is RA and declination to radians, which will allow trigonometric operations to be used. Using the parallax now in parasecs, I convert the coordinates to Cartesian form, (X, Y, Z), using the standard formulas.

X ​=d⋅cos(Dec)⋅cos(RA)

Y =d⋅cos(Dec)⋅sin(RA)

Z = d⋅sin(Dec)​

Now that I have unit conversions complete and a set of Cartesian coordinates, I can find how luminous each star is and use the center of mass formula. First, I use the magnitude of the G-Band frequency to estimate the luminosity. Each luminosity is equal to 10 to the power of the G-band times -.4 according to the University of Oregon. A brighter star will have a high luminosity and a higher mass. Finally, we apply the center of mass formula. We find mass using luminosity. The approximate mass of a star to its luminosity was first discovered by Jakob Karl Ernst Halm. The common value for most stars is 3.5. Our equation is then mass = luminosity to the power of 1 divided by 3/.5.

Now that the masses have found the center of mass can be applied. First, I sum all of the masses. Then I divided each coordinated sum of mass (X, Y, Z) by the total sum. This produces the galaxy's COM in parsec units. (1 Parsec is 1.917e+13 Miles).

Finally, I created a 3D scatter plot. As the stars get brighter, the dots get darker. The red dot shown in the plot is my calculated center of mass. I utilized AI to add nice, professional comments to my code, and I was done right??

No, no, onto velocities.

After this initial data and plotting, I decided to add more data. Going through the Gaia website, which contains all data collected by the probe, I noticed that it contains velocity data. The bulk motion of each star will help create a more accurate plot of the stars and show which direction the COM is moving. This new data also contains 500,000 stars compared to 275,000 in the previous plot.

With this new data, I needed to do more unit conversions and calculations to factor in the movement of each star. The first piece of new data provided is Proper Motion in RA and Dec. This is the angular motion of stars in the sky. Also provided is Radial Velocity (RV), which is the line of sight velocity either towards or away from us in kilometers per second. We are going to convert these numbers into a velocity vector - V = [Vx, Vy, Vz].

To do this, we must convert all of the provided data into a formula to find linear transverse velocity. This is the velocity any given object is moving perpendicular to a given point, in this case, Earth.

In the formula, 4.74 is found from converting AU/year to kilometers per second. An AU or astronomical unit is the distance from the Earth to the sun, approximately 150 million meters.

An issue I ran into is that to find the center of mass, I need to project angular motion into Cartesian Coordinates. To do this, I use cos() to represent the motion parallel to the base, and sin() to represent the motion perpendicular or upwards from the base. Using these directional cosines, I can project the angular motion into x,y, and z. The final step is combining the transverse motion with the radial motion provided. Adding these together gives the total space velocity of each star. By summing all of the directional velocities and dividing by the total mass, we can find the velocity of the COM and produce a new graph. Dividing the vector we found by its magnitude also provides the direction of the movement seen in the graph below.

I hoped you enjoyed reading this write-up. I learned a ton about astrophysics and open-source data during this project, as well as being able to apply the MATLAB I learned in the spring. Hopefully, I get an A!

Until next time,

Alex