Description
For this assignment you are going to write a code that will predict the future trajectory
of M33 in the frame of M31. You will need the COM orbit files for M33 and M31 that
you computed in Homework 6 so you can plot their relative positions and velocities as a
function of time and compare them to the analytic solutions.
See Patel, Besla, Sohn 2017a
(https://arxiv.org/abs/1609.04823) for more details and background on M33’s past orbital
trajectory.
Note that G = 4.498768e−6 in units of kpc3/M /Gyr2
(define this as a global variable,
but let’s forget about storing the units and use values in this homework).
1 M33AnalyticOrbit
Create a Class called M33AnalyticOrbit.
In this class we will create a series of functions that will determine the acceleration M33
feels from M31 and integrate its current position and velocity forwards in time.
Initialize the class (def init ), taking as input a filename for the file in which you
will store the integrated orbit. At the beginning of the class, also initialize the following
quantities:
• The relative position and velocity of M33 to M31
– self.x, self.y, self.z : the relative COM position vector of M33 to M31
calculated using the disk particles of both galaxies at Snapshot 0 (location of
M33 today) from Assignment 4
– self.vx, self,vy, self.vz : the relative COM velocity vector of M33 to M31
as calculated in Assignment 4
• The scale lengths and masses for each component in M31
– self.rdisk and self.Mdisk : set self.rdisk=5 kpc and use the disk mass
computed in Assignment 3 (the value of mass should be in units of M , same for
below)
1
– self.rbulge and self.Mbulge : set self.rbulge=1 kpc and use the bulge mass
from Assignment 3
– self.rhalo and self.Mhalo : for self.rhalo, use the Hernquist scale length
(a) computed in Assignment 5 and use the halo mass from Assignment 3
2 Define the Acceleration Terms
Define functions that will compute the gravitational acceleration vectors from 3 components
of the M31 galaxy: Halo, Bulge and Disk.
2.1 Halo and Bulge Acceleration
1. The gravitational acceleration induced by a Hernquist profile is given by:
a = −∇ΦHernquist = −
GM
r(ra + r)
2
r =
GM
r(ra+r)
2 x
GM
r(ra+r)
2 y
GM
r(ra+r)
2 z
(1)
where M is the total halo or bulge mass, ra is the corresponding scale length, r is the
magnitude of the relative position vector.
2. Define a function HenquistAccel that takes 6 inputs: self, M, r a, x, y, z.
3. This function returns the acceleration vector from a Hernquist potential.
4. This function will be used for both the halo and the bulge of M31 (set by the input M
and r a).
2.2 Disk Acceleration
1. For the disk of M31, we will use an approximation that mimics the exponential disk
profile at distances far from the disk. It is called a Miyamoto-Nagai 1975 profile. The
potential for this profile is:
ΦMN(r) = −GMdisk
√
R2 + B2
, where R =
p
x
2 + y
2 and B = rd +
q
z
2 + z
2
d
, (2)
where rd is self.rdisk in our code and zd can be defined as self.rdisk/5.0 in this
homework.
The acceleration vector is thus given by:
a = −∇ΦMN =
ax
ay
az
=
−
GMdisk
(R2 + B2
)
1.5
x
−
GMdisk
(R2 + B2
)
1.5
y
−
GMdiskB
(R2 + B2
)
1.5
p
z
2 + z
2
d
z
. (3)
2
Note the z component of the acceleration vector is a bit more complicated.
2. Define a function MiyamotoNagaiAccel that takes 6 inputs: self, M, rd, x, y,
z.
3. This function returns the acceleration vector from a Miyamoto-Nagai profile.
2.3 M31Acceleration
Define a new function M31Accel that sums all acceleration vectors from each galaxy component (be careful about the vector summation; you can always take advantage of the array
operations provided by numpy). This function takes as input the 3D position vector (x,y,z)
and returns a 3D vector of the total acceleration.
3 Build an Integrator
We need to solve the orbit of M33 by integrating the following equation of motion forwards
in time
¨r = −∇Φtotal = −∇
Φhalo(r) + Φbulge(r) + Φdisk(r)
(4)
⇒ ¨r = ahalo + abulge + adisk, (5)
where Φ represents the potential for each galaxy component. To do this, we will treat M33
as a point mass and adopt a variant of the “Leap Frog” integration scheme given that a is
a pure function of r.
In your code, define a function LeapFrog that takes as input:
• a time interval for integration dt (∆t)
• a starting position vector x,y,z for the M33 COM position relative to the M31
• a starting velocity vector vx, vy, vz for the M33 relative to M31
In each integration step, your code will update the positions and velocities using standard
kinematic equations, according to the following (see also Figure 1):
1. Predict the 3D position vector of M33’s center of mass (COM) at the middle of the
timestep ∆t using the current COM velocity and position vectors according to
rn+ 1
2
= rn + vn
∆t
2
. (6)
2. Then the COM position and velocity vectors are advanced a full time step using the
acceleration at the 1/2 timestep (calculated based on the rn+ 1
2
) according to
vn+1 = vn + an+ 1
2
∆t, (7)
rn+1 = rn +
1
2
vn + vn+1
∆t, (8)
3
�! �!”# �!”$
�
!”
#
$
�
!”
%
$
Step 1: �
!”!
”
= �! + �!
∆$
%
Step 2.2: �!”& = �! + [�#”�#$!]
%
∆�
Step 2.1: �!”& = �! + �
!”!
”
∆�
Figure 1: This figure illustrates the LeapFrog integration scheme to evolve the equation of
motion of M33. The subscript (n, n + 1, etc.) denotes the steps along the time axis.
where the acceleration is determined by calling the function self.M31Accel. Note that
∆t can be positive or negative because Leap Frog integrators are symplectic, meaning they
can be used for calculations that run both forward and backward in time. For this assignment,
you’ll want positive values since we are calculating future orbits. In addition, Equation 8 is
equivalent to
rn+1 = rn+ 1
2
+ vn+1
∆t
2
. (9)
You may use either one as you see how your code fits.
4 Integrate the Orbit
Now we will loop over the LeapFrog integrator to solve the equations of motion and compute
the future orbit of M33 for 10 Gyr into the future.
• Define a function OrbitIntegrator that takes as input: self, a starting time to, a
time interval ∆t and a final time tmax.
• Supply the starting COM position and velocity of M33 relative to M31 as the initial
conditions of this orbital integration.
• Define a variable t and start the integration from to. Continue looping (i.e. a while
loop) over LeapFrog until you reach tmax, updating the time, positions and velocities
along the way. Store the results in an array that you initialize outside the loop (as you
4
did when you calculated the COM orbits from the simulations). Don’t forget to store
the initial positions and velocities before you begin the loop.
• Once the loop is complete, store the array into a file, like in Homework 6.
5 Analysis
1. Create a plot of your predicted M33 orbit from to = 0 Gyr to tmax = 10 Gyr. Start
with 0.5 Gyr intervals for ∆t and refine (for example, to 0.1 Gyr) once you know the
code is working. Overplot the solution to Assignment 6 for M33’s orbit with respect
to M31 from the simulation. Do this for both the total position and total velocity as
a function of time.
2. How do the plots compare?
3. What missing physics could make the difference?
4. The MW is missing in these calculations. How might you include its effects?
5