Methods of teaching numerical methods to solve ordinary differential equations in the context of galaxy mergers were explored. The research published in a paper by Toomre and Toomre in 1972 describing the formation of galactic tails and bridges from close tidal interactions was adapted into a project targeting undergraduate physics students. Typically undergraduate physics students only take one Computational Physics class in which various techniques and algorithms are taught. Although it is important to study computational physics techniques, it is just as important to apply this knowledge to a problem that is representative of what computational physics researchers are investigating today. The model that was developed is capable of showing general trends in galactic evolution and is a good introduction for students who hope to expand on this topic in the future.
Computational physics is based upon the idea that a computer can simulate certain physical environments to produce novel scientific conclusions. This field allows scientists to conduct experiments on the computer that are just not realistic or possible to do in the real world because of the investigation of topics involving extremely large data sets and complicated physical interactions[1]. This paper investigates the components and processes of modeling galaxy mergers on an undergraduate academic level. The project outlined in this paper is derived from a paper written by Toomre and Toomre in 1972 that described the formation of galactic tails and bridges from close tidal interactions [2]. This project could be assigned to undergraduate students in an introductory computational physics class to demonstrate the capabilities and possibilities of what the field has to offer.
The galaxy mergers project teaches students numerical methods to solve ordinary differential equations. The galaxies modeled interact via ordinary differential equations so students are able to learn about these techniques through an applied project that is a simplified version of a cutting edge computational physics topic today[3]. A lot of the time undergraduate physics students only take one computational physics class in which they learn about various techniques and algorithms that can in theory be applied to model more complicated physical situations [4]. Although it is important to learn about the underlying theory and algorithms, it is just as important to apply this knowledge to a problem that is representative of what computational physics researchers are investigating today. The project in this paper does in fact accomplish this goal and teaches numerical methods through application.
Although this project uses ordinary differential equations to model galaxy mergers, currently computational physics researchers use supercomputers to solve partial differential equations using the Boltzmann approximation [5]. The cutting edge research today accounts for the evolution of black holes and dark matter, star formation, and the interactions between individual stars [6, 7, 8, 9]. These characteristics are important in modeling galaxy mergers, but supercomputers are needed to compute all of the calculations. Our model works well enough to see general trends in galactic shapes and is a good introduction for students who want to expand on this topic in the future.
Although Toomre and Toomre’s paper dates back to 1972, it is still a good model to use for undergraduate students in a computational physics class. Their research provided good and accurate results despite the assumptions and simplifications that they made. Researchers today can produce more precise results and investigate a larger range of topics than in 1972, but the foundation for producing these results is still the same. This project is an excellent introduction to the capabilities of computational physics and is a great learning opportunity for undergraduate physics students.
Galaxies are a collection of stars that are orbiting around a common center. Galaxies come in various shapes and sizes, but the scientist Edwin Hubble designed a classification system for galaxies that is still used today. He grouped galaxies into four categories including elliptical, spiral, disk, and irregular. Today we have expanded upon these four groups, but only to add subcategories. Today elliptical galaxies can be classified according to the ratio between their minor axis and major axis dimensions to differentiate between more circular as opposed to more elongated shapes. Spiral galaxies are now classified with how tight or loose the spiral arms are, and whether or not there are a cluster of stars in a bar shape in the center of the galaxy. Disk galaxies are considered to be a transition point between elliptical and spiral galaxies, and were added to the classification system by Hubble approximately 10 years after he created the original one. Although there are galaxy classification systems in place, a galaxy does not necessarily maintain its shape throughout its life time.
Galaxies can change shape due to galactic mergers and interactions. When one galaxy approaches or collides with another, the gravitational force between the two can be so strong that the general shapes of the galaxies can forever be changed. For instance gravitational interactions with other galaxies have been known to morph the arms of spiral galaxies, change or add the bar structure in spiral galaxies, add features such as tidal tails, and sometimes completely merge the two galaxies together [10]. This paper examines how the positions of individual stars in one galaxy can change with an interaction of another galaxy. This is important in the grand scheme of things, because this project allows students to predict how a galaxy’s shape can morph over time in various situations. Trends in galactic morphology are easily visible in this project despite the simplifications that were made.
Our modern understanding of cosmology indicates that a large fraction of the total matter in the universe is in fact dark matter. Dark matter does not emit or absorb light and is therefore virtually invisible. Dark matter is able to interact gravitationally with other matter, so it can be detected by observing trajectories of visible matter. For example the orbital velocity of visible matter as a function of radius from the center of spiral galaxies does not follow the standard Keplerian motion one would expect for a system with most of the mass concentrated at the center. In fact, the orbital velocities are larger which suggests that the mass increases more than expected along the galactic radius. This phenomenon suggests the presence of dark matter in galaxies [11]. Dark matter is distributed non-uniformly throughout the universe. When modeling galaxy mergers and interactions, advanced computer models use smooth particle hydrodynamics to model the diffuse gas and dark matter components of the galaxies. These simulations can convert the gas into new star populations while the dark matter interacts gravitationally through mean field approximations[3]. Our simplified model ignored these components of galaxies.
** Supermassive black holes are believed to reside in the center of galaxies due to their strong gravitational pull in relation to their size [12]. Recently a black hole was discovered in the center of a quasar that was 12 billion times the mass of our sun [13]. It is hypothesized that super massive black holes are actually formed through galaxy mergers [14]. Since the universe is expanding, the early universe was a lot more compact and more galactic interactions were able to occur. The kind of computational modeling used in this project therefore provides students a way to study galactic mergers and the interactions of black holes in the early universe. **
Galaxies are very complicated entities and there are a lot of factors that can contribute to the positions and velocities of stars. First of all, galaxies contain an uneven distribution of matter which is very difficult to model. Second of all, the stars themselves interact gravitationally with each other despite ultimately orbiting around the galactic center. These factors involve advanced computational techniques and equipment that is more serious than an average computer that an undergraduate student may own. Following the example of Toomre and Toomre [2], the simulation was simplified so that the stars in this project were treated as unit mass point objects which do not interact with each other but do interact with a large, concentrated galactic core point mass of 1e11 solar masses in the center of the galaxy. The unit-mass stars were oriented in rings around the galactic center. The two interacting galaxies were labeled "main" and "disruptor", with unit mass stars orbiting the main galaxy but none orbiting the disruptor. Although the disrupting galaxy could follow either hyperbolic, parabolic, or elliptical approach trajectories; the simulation only computed parabolic trajectories. These simplifications made it possible to conduct the computations at each time step using Newtonian gravity differential equations that describe motion.
**Figure 1: Diagram of the orbital geometry and variables**$Y_0$ | Initial "Y" position of disrupting galaxy |
---|---|
$\vec{V}$ | Initial velocity of disrupting galaxy |
$\vec{V}_{min}$ | Minimum velocity of disrupting galaxy |
$R_0$ | Initial separation between the two galaxies |
$R_{min}$ | Minimum separation between the two galaxies |
Key for Figure 1
$Rx_0$(Mpc) | $Ry_0$(Mpc) | $M$($M_{sol}$) | $S$($M_{sol}$) | |
---|---|---|---|---|
Direct Passage | -39 | -80 | 1e11 | 1e11 |
Retrograde Passage | -39 | 80 | 1e11 | 1e11 |
Light Mass Disruptor | -39 | -80 | 1e11 | 1e11/4 |
Heavy Mass Disruptor | -39 | -80 | 1e11 | 1e11*4 |
User Input Values
The table above displays the values that a user can input into the program. Using these numbers and the constraint of a parabolic passage, the initial velocity parameters can be computed.
This project consisted of two galaxies. One galaxy contained stars and remained at the origin of our reference frame, whereas the other galaxy had a perfectly parabolic trajectory (in relation to the origin) and was simply a point mass with no stars. This configuration allows one to clearly observe the evolution of star positions in a galaxy merger. Additionally the only cases that were examined in this project were direct passage, retrograde passage, light mass disruptor and heavy mass disruptor cases which will soon be explained.
After R_min was set, it was possible to calculate the initial velocity parameters of the disrupting galaxy. The disrupting galaxy was set to follow a perfect parabolic trajectory, so the initial velocity conditions were extremely important. These conditions pave the trajectory that the galaxy will follow. First, the Vis-Viva equation (Equation 1) was used to determine the minimum velocity of the disrupting galaxy that occured at R_min.
\begin{equation} \hspace{10cm}v^2 = GM(\frac{2}{r} - \frac{1}{a}) \hspace{10cm}[1] \end{equation}** where v is the relative speed between two bodies, r is the separation between two bodies, M is the mass of the central body, G is the gravitational constant, and a is the semi major axis.**
With parabolas the semi-major axis approaches infinity,so the Vis-Viva equation can be reduced to
$$ \hspace{10cm}v_{(min)}= \sqrt{\frac{2GM}{R_{(min)}^2}} \hspace{9cm}[2] $$Using v_min, the angular momentum, l_min can be calculated from
\begin{equation} \hspace{8cm}l = l_{(min)} = r_{(min)}*v_{(min)}\sin(90) = constant \hspace{6cm} [3] \end{equation}which is constant for the entire trajectory of the disrupting galaxy. This angular momentum value can be used to determine the parameters of the general equation of a parabola described by
\begin{equation} \hspace{11cm}y^2 = c^2 - 2cx\hspace{10cm} [4] \end{equation}where the paramter c is related to the masses and angular momentum through
\begin{equation} \hspace{10cm}c = \frac{l^2}{Gm_1m_2\mu}\hspace{11cm} [5] \end{equation}\begin{equation} \hspace{10cm}\mu = \frac{M_1M_2}{M_1 + M_2} \hspace{11cm}[6] \end{equation}A starting value for the initial y position, y0, was set, and the general equation of the parabola was used to determine the corresponding x0 value. These coordinates were ultimately used to determine the initial separation value between the two galaxies (R0) which can be calculated
**
$$
\hspace{10cm}R_0 = \sqrt{(x_0^2) + (y_0^2)} \hspace{10cm}[7]
$$
**
From R0, the initial velocity can be determined from a simplified version of the Vis-Viva equation
$$
\hspace{10cm}v_0= \sqrt{\frac{2GM}{R_0^2}}\hspace{11cm}[8]
$$
**
The modified Vis-Viva equation for the case of circular orbits is applied using
**
$$
\hspace{10cm}v_0= \sqrt{\frac{2GM}{R_0^2}}\hspace{11cm} [9]
$$
With the initial conditions set, the calcuation of the subsequent position and velocities for each time step under Newtonian graviational forces can be undertaken. The positions of the unit-mass stars are described by the differential equation ** $$ \hspace{10cm}\ddot{\mathbf{r}} = -\gamma \left\{ \frac{M}{r^3}\mathbf{r} -\frac{S}{\rho^3}\boldsymbol{\rho} + \frac{S}{R^3}\boldsymbol\Re \right\}\hspace{8cm}[10] $$ ** and that of the disrupting galaxy by ** $$ \hspace{10cm}\ddot{\boldsymbol\Re} = -\gamma \frac{M+S}{R^3}\boldsymbol\Re\hspace{11cm} [11] $$ ** with the variables and constants defined as follows: **
To validate the correctness of our simulation, we compared our results to the four cases studied by Toomre and Toomre [2]. The four cases were direct passage, retrograde passage, light mass disruptor and heavy mass disruptor. Direct passage occurs when the stationary galaxy and disrupting galaxy have equal mass, and the disrupting galaxy approaches the stationary one along a parabola with the same direction that the stars are orbiting. Retrograde passage is very similar except the disrupting galaxy approaches in the opposite direction that the stars are orbiting. A way to visualize these concepts is that the stars are orbiting counter clockwise and in a direct passage the galaxy approaches counter clockwise as well whereas the retrograde passage follows a clockwise trajectory. A light mass disruptor was also modeled where the disrupting galaxy approached along a direct approach, but had only a quarter of the mass of the stationary galaxy. The heavy mass disruptor on the other hand was also a direct approach, but the disruptor had four times the galactic mass of the stationary galaxy. With the model validated against previous research, it can be used to explore new configurations with confidence that the physics is valid.
Let us begin by examining the stationary galaxy containing stars which is set at the origin. This galaxy has a point mass, which means that all of its mass occurs at a point in its center. The galaxy contains 120 unit-mass stars which are initially distributed throughout 5 discrete, concentric rings around the galaxy. Each ring contains 12, 18, 24, 30, and 36 stars respectively (the inner most ring contains 12 stars). For the direct passage, retrograde passage, and light mass disruptor cases; the concentric rings are placed at exactly 20, 30, 40, 50, and 60 percent of the minimum distance between the two galaxies (R_min).
Figure 3: Stationary galaxy for direct, retrograde and light mass disruptor cases
On the other hand, in the heavy mass disruptor case the rings are placed at 12, 18, 24, 30, and 36 percent of R_min.
Figure 4: Stationary galaxy for heavy mass disruptor case
In all cases this minimum distance (Rmin) is set at 25 kpc. The initial time was set to t = 0 and the time step was set to exactly 10^8 years. In all of the featured cases, the heavier mass was set to 10^11 solar masses. Various functions inside the iPython notebook were used in the project including NumPy, Matplotlib, SciPy, odeint, and the iPython interact. These functions made it possible to conduct all of the necessary computations in this project.
Set_Init_R_Cond(Ry0, M, S):
The first step is to set the initial conditions of the disturbing galaxy. A function was written that depended on the initial “y” condition of the disturbing galaxy (Ry0), the mass of the stationary galaxy (M), and the mass of the disturbing galaxy (S). The minimum separation distance between the two galaxies was called Rmin, which was set at 25 kpc. The minimum velocity at Rmin was set as well following equation [8] with variable values of galactic masses. The angular momentum was determined using these values which allow for the c factor to be computed according to equations [5] and [6]. The c factor was used to compute the corresponding initial “X” position of the disturbing galaxy (Rx0) which in turn could be used to compute the distance of the disturbing galaxy from the origin (R0) according to [7]. The unit tangent vectors were then determined for the parabola trajectory which were used to set the initial “X” and “Y” velocities of the disturbing galaxy. The function returns initial “X” and “Y” positions and velocities of the disturbing galaxy.
Ring(particles, percent, G, M):
A function was set up where formulas and empty lists were set up to produce initial conditions of the stars around the stationary galaxy. This function took in arguments called particles, percent, G, and M. Particles described the number of stars that were desired in one ring around the center of the stationary galaxy. Percent described at what percentage of Rmin the distance between each star and the origin of the stationary galaxy should be. Lastly, G described a gravitational constant and M was the mass of the stationary galaxy. First the radius of each star was defined as a percent of the Rmin value set in the previous function. The stars can be plotted along a circular trajectory around the origin of the stationary galaxy. The arc length of each trajectory is divided among each ring of stars so that the stars are spaced equally apart. This is accomplished by creating a while loop that produces a star location with each turn of the loop. The loop runs as many times as there are stars to be created. Although the stars are initially in polar coordinates, their positions are converted to Cartesian coordinates and appended to an empty list. The initial velocity conditions are also created using equation [8]. These velocities are tangent to the circle that the stars are plotted along. This function returns an array of star positions and an array of star velocities.
Init_rings (G,M):
A function was set up where the gravitational constant (G) and the stationary galactic mass (M) were the inputs. The previous “ring” functions set up all of the formulas using variables to set the initial positions and velocities of the stars. This function upon the “ring” function so that the actual values of the number of stars desired and the percentage of Rmin are computed to create an array of real numbers that represent the positions and velocities of the stars. There are two ‘init-ring’ functions because different star numbers and percentages are used for different cases. Specifically, each ring in every case contains 12, 18, 24, 30, and 36 stars respectively (the inner most ring contains 12 stars). For the direct passage, retrograde passage, and light mass disruptor cases; the concentric rings are placed at exactly 20, 30, 40, 50, and 60 percent of the minimum distance between the two galaxies (R_min). On the other hand, in the heavy mass disruptor case the rings are placed at 12, 18, 24, 30, and 36 percent of R_min. This function returns an array of position values and an array of velocity values.
Unpack_rings_vel(rings, velocity):
This function takes in the arguments ‘rings’ and ‘velocity’ which were the outputs from the previous function. The purpose of this function is to organize the information in the rings and velocity arrays into separate lists of x positions of stars, y positions of stars, x velocities of stars, and y velocities of stars. This can be accomplished by writing a for loop that loops through the rings and velocity arrays and selects the desired information using indexing techniques and appending the proper information to an initially empty list. This function returns four separate arrays that describe the x and y positions and velocities of the stars.
Derivgalaxy(y, t, M, S):
This function inputted an array containing all of the position and velocity information of the disturbing and stationary galaxies, the period of time to be examined, and the stationary and disrupting galactic masses. Additionally all of the arrays of information were inputted into equations [10] and [11] that describe the position of each star and disrupting galaxy over time. This function outputs an array of the x and y velocity and acceleration values of each star and disrupting galaxy.
Make-Master_Array(Case = 1, Rx0 = -39, Ry0 = -80, M = 1.0e11, S = 1.0e11);
This function inputted the case which will be examined. Each of the four cases was set arbitrarily to numbers to make them simpler to call. The function also inputted the initial x and y position of the disrupting galaxy, and the galactic masses of the disrupting and stationary galaxies. This function runs the previous functions according the case number that is to be examined. This could be achieved by using if and elif statements to run certain functions if a certain case is called. The output of the derivgalaxy is solved by using an odeint function which is an ordinary differential equation solver. This function loops over each time step and appends all of the information of the positions and velocities of the stars and disrupting galaxy at each time step to a master array. The purpose of this is to run the calculation one time in its entirety and to store the results in an array for simple referencing. This function output a complete master array containing all of the information of the stars and disrupting galaxy at each time step.
Make_Plot_stars(results, M, S, t, dt):
This function input the results of the make_master_array function, the masses of the disrupting and stationary galaxies, the total time the user would like to examine, and the time steps which were integrated over. This function sets the stationary galaxy location at the origin despite the fact that it does technically move. Additionally, this function ensures that all of the stars will orbit and move in reference to the set origin at the stationary galactic center. The information from the make_master_array function is plotted using an iPython interact widget tool that allows the user to slide through all of the calculated time steps and observe the evolution of the galaxies in a movie like layout.
Make_Plots_Yellow_Star(results, M, S, t, dt, Yellowstar):
This function is almost identical to the previous Make_Plot_Stars function, except it contains an extra argument called yellowstar that allows the user to index a certain star in the master array to be a yellow star where all of the other stars are red. Essentially a for loop is written that plots and colors all of the stars red. Next, another for loop is written that only loops over a desired star and plots a yellow star on top of the existing red one. Although two stars are plotted in one location, the user only sees the yellow star that is covering the red one.
When creating the galaxy mergers model, one must be aware of some challenges that can impede one's progress. First of all, it is important to set initial position and velocity conditions according to the equations listed above in the "Model, Assumptions, and Simplifications" paragraph. If arbitrary values are set instead, the disrupting galaxy will not follow a perfect parabolic approach and the stars will not move along realistic trajectories.
Another challenge to be aware of is setting a proper time step. If one sets a time step that is too large, the positions of the stars are not calculated frequently enough and they may jump to the origin when the disrupting galaxy is close. This problem occurs because of a numerical precision problem that results in the magnification of error. The time step used in this project is 0.007 and is the recommended time step to use.
After a student gets this project working, they may build upon the existing program with the following suggestions. First, it may be interesting to examine the trajectory of a single star and follow its path. One way to do this is to color a particular star in a different color to make it stand out amongst the other stars. The method to accomplish this is described in the previous section called "How to set up General Code" under the "Make_yellow_star" description. This modification allows students to see which stars are transferred to the disturbing galaxy and which stars stay with their own galaxy.
Another modification may be to add stars to the disturbing galaxy following similar rules as modeling stars around the stationary galaxy. It would be interesting for a student to examine how the different stars mix together over time by color coding the stars based on which galaxy they came from. Also the student can manipulate the tilts of the interacting galaxies and study how the star trajectories are affected with this modification.
Finally it would be interesting to include gravitational interactions between the stars. This suggestion would take the most computing power and therefore it is reccommended for a student to start with only three stars. Once a student acccomplishes this extension with three stars, they may add more as long as their computer can process the information.
Computational Physics is a very broad field that covers many interesting topics. This project focuses on modeling Galaxies which is a sub-category of computational physics. Galaxies have many complicated factors such as dark matter, black holes, and star-star interactions that computational physicists must model on super computers. Super computers are not commonly available to undergraduate physics students, therefore this project was created to let students model galaxies using the iPython notebook on their laptops. The project was modeled after the 1972 Toomre and Toomre paper that described the formation of galactic tails and bridges. Toomre and Toomre imposed certain simplifications on their model including unit-mass stars, no dark matter or black holes, and the use of 120 stars. These assumptions were then adopted into this project. Although the model is not as precise as some of the cutting-edge research that is currently being conducted, the model is definitely good enough to see galactic trends and is therefore interesting to study.
I would like to thank the College Based Academic Fee (CBF) for funding my research. This money allowed me to spend my summer developing this project into what it has become. Additionally I would like to acknowledge the Spring 2014 Physics on a Computer class that introduced me to the topic of modeling Galaxy Mergers. This class is really where this project began! Last of all I would like to thank my wonderful advisor Jennifer Klay for her patience and hardwork in making this project a reality.