The project was originally developed by Jennifer Klay and modified by Brian Granger. The idea for this project came from Jorge Moreno, who visited our department to give a colloquium on his research in Winter 2014. All content is licensed under the MIT License.
One of the fundamental questions in astrophysics is how galaxies form, evolve and interact. It might seem surprising that galaxies separated by vast distances could interact with one another, but their immense masses and the gravitational forces that govern their formation and evolution nevertheless can lead to violent interactions. The Milky Way galaxy is right now on a collision course with the Andromeda galaxy and they are predicted to collide in about 4 billion years. Even though both galaxies are made up of hundreds of billions of stars, the separation between stars in the galaxies means that the chances of direct stellar collisions are extremely small. Nevertheless, if our descendents are still around to experience it, the view of our night sky will definitely change as a result of this interaction.
Here is a computer simulation of the merger generated using data from the Hubble Space Telescope and our knowledge of the gravitational interaction among the constituents.
In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('4disyKG7XtU')
Out[1]:
So how does one go about simulating the collision of two galaxies? The video above was probably generated on massive supercomputers running full blast for a very long time to generate the paths of all of the individual stars as the galaxies merged. Is it really possible that we could generate something like that using the programming skills you have learned? The answer is yes, if we make some simplifying assumptions.
In this project you will create a simulation of galaxy mergers using the methods of Toomre and Toomre, two brothers who used state of the art computers at MIT and NYU in 1972 to investigate the dynamics of massless point particles orbiting a massive galactic nucleus $M$ in a parabolic orbit about the center of mass with another massive galactic nucleus $S$. Their work was an extension of a previous paper published in German in 1963 describing the system of equations and an early attempt to investigate how the spiral structure of galaxies emerges.
In [1]:
def numpy():
print('pls')
In [4]:
numpy()
Toomre and Toomre's 1972 paper on the simulation of galaxy mergers with Newtonian mechanics.
Pfleiderer 1963 is in German, but includes the equations and the set up of the problem. Jennifer Klay translated the relevant parts of the paper to come up with a description of the equations and variables that were used by Pfeliderer and the Toomres to do their simulation.
This is a restricted 3-body problem in which particles composing the outer disk of galaxy $M$ are massless but nevertheless interact through inverse square laws with the mass centers of their galactic central mass ($M$) and the central mass of the disrupting galaxy ($S$).
The calculation is performed in the rest frame of the mass $M$ lying at the origin, with the starting positions for the massless point particles $m_i$ (stars) in the orbits around it given by $(\boldsymbol{r_0})_i$ and the position of mass $S$ relative to $M$ given by $\boldsymbol{\Re}$. The evolution of the positions of the stars $m_i$ and galaxy $S$ relative to $M$ is dictated by the set of differential equations:
where
One starts with a position and velocity vector for $S$, and the $m_i$ (remember that we are in the rest frame of $M$ so it is at rest at the origin of the coordinate system) and then solves these differential equations to get the positions of $S$, and the set of $m_i$ ($i$ = 0,...,$N$) as a function of time under the influence of pure Newtonian gravity.
After the equations are solved, one can either draw static images of the system at specific points in time, or form an animated movie of the interaction. For this project, you will do both.
In the image from Toomre and Toomre shown below, they transformed to the center-of-mass of the $M$+$S$ system for the first 4 frames and then switched back to the rest frame of $M$ for the six subsequent frames.
Section II of the paper describes exactly the algorithm used to compute the results, followed by four examples.
The original paper describes four examples:
You should pick one/two of these cases as your base question.
Pfleiderer's paper lists a variety of interactions and initial conditions, including elliptical orbits for the point masses, from which the Toomre examples are a subset. Your two additional questions should come from this list of initial conditions.
For all initial conditions studied you should have:
| Type | Position $X$ | Position $Y$ | Velocity $U$ | Velocity $V$ | Eccentricity $\epsilon_s$ | Mass ratio $S/M$ |
|---|---|---|---|---|---|---|
| S1+ | 1 | 0 | > 0 | 0 | 1 | 1 |
| S1- | 1 | 0 | < 0 | 0 | 1 | 1 |
| S2+ | 1 | 0 | > 0 | 0 | 7 | 1 |
| S2- | 1 | 0 | < 0 | 0 | 7 | 1 |
| S3+ | 1 | 0 | > 0 | 0 | 31 | 3 |
| S3- | 1 | 0 | < 0 | 0 | 31 | 3 |
| S4 | 1 | 0 | 0 | > 0 | 31 | 3 |
| S5 | 1 | 0 | 0 | > 0 | 1 | 1 |
| S6 | 0 | 1 | > 0 | 0 | 31 | 3 |
| S7 | 0 | 1 | > 0 | 0 | 1 | 1 |
The exercises related to Ordinary Differential Equations are useful guidance for how to set up your solution. Start with a single point mass orbiting the main galaxy and its interactions with the disrupting galaxy. Once you have it working for a single test mass, add in more until you have a complete solution.
You will need to tune the error tolerances of odeint to ensure that you are accurately solving the equations of motion. One good way of making sure you have the errors under control is to compute the total energy of the system (which should be constant). This will require a bit of derivation, but would be very worth your time.
One of the main challenges you will face is how to perform animations of your results. Here are some thoughts on how to do that:
npz format). You should come up with a naming scheme for your files that makes it clear what parameters were used for each simulation. You may even want to save those parameters into a JSON file with the same naming scheme.I will provide you with some additional information about how to create animations using Matplotlib.