N-body (gala.dynamics.nbody)

Introduction

With the Hamiltonian and potential classes (potential), Gala contains functionality for integrating test particle orbits in background gravitational fields. To supplement this, Gala also now contains some limited functionality for performing N-body orbit integrations through direct N-body force calculations between particles. With the gala.dynamics.nbody subpackage, gravitational fields (i.e., any potential class from gala.potential) can be sourced by particles that interact, and optionally feel a background/external potential. To use this functionality, the core class is DirectNBody. Below, we’ll go through a few examples of using this class to perform orbit integrations

For code blocks below, we’ll assume the following imports have already been executed:

>>> import astropy.units as u
>>> import numpy as np
>>> import gala.potential as gp
>>> import gala.dynamics as gd
>>> from gala.dynamics.nbody import DirectNBody
>>> from gala.units import galactic, UnitSystem

Getting started

The DirectNBody, at minimum, must be instantiated with a set of particle orbital initial conditions along with a specification of the gravitational fields sourced by each particle — that is, the number of initial conditions must match the input list of gravitational potential objects that specify the particle mass distributions. Other optional arguments to DirectNBody allow you to set the unit system (i.e., to improve numerical precision when time-stepping the orbit integration), or to specify a background gravitational potential. Let’s now go through a few examples of using this class in practice.

Example: Mixed test particle and massive particle orbit integration

Like with Hamiltonian orbit integration, orbital initial conditions are passed in to DirectNBody by passing in a single PhaseSpacePosition object. Let’s create two initial conditions by specifying the position and velocity of two particles, then combine them into a single PhaseSpacePosition object:

>>> w0_1 = gd.PhaseSpacePosition(pos=[0, 0, 0] * u.pc,
...                              vel=[0, 1.5, 0] * u.km/u.s)
>>> w0_2 = gd.PhaseSpacePosition(pos=w0_1.xyz + [100., 0, 0] * u.pc,
...                              vel=w0_1.v_xyz + [0, 5, 0] * u.km/u.s)
>>> w0 = gd.combine((w0_1, w0_2))
>>> w0.shape
(2,)

We’ll then treat particle 1 as a massive object by sourcing a HernquistPotential at the location of the particle, and particle 2 as a test particle: To treat some particles as test particles, you can pass None or a NullPotential instance for the corresponding particle potential:

>>> pot1 = gp.HernquistPotential(m=1e7*u.Msun, c=0.5*u.kpc, units=galactic)
>>> particle_pot = [pot1, None]

Let’s now create an N-body instance and try integrating the orbits of the two particles. Here, there is no external potential, so particle 1 (the massive particle) will move off in a straight line. We’ve created the initial conditions for particle 2 so that it will remain bound to the potential sourced by particle 1, and so will orbit it as it moves. Let’s create the object and integrate the orbits:

>>> nbody = DirectNBody(w0, particle_pot)
>>> orbits = nbody.integrate_orbit(dt=1e-2*u.Myr, t1=0, t2=1*u.Gyr)
>>> fig, ax = plt.subplots(1, 1, figsize=(5, 5)) 
>>> _ = orbits[:, 0].plot(['x', 'y'], axes=[ax]) 
>>> _ = orbits[:, 1].plot(['x', 'y'], axes=[ax]) 

(Source code, png, pdf)

../_images/nbody-1.png

Example: N-body integration with a background potential

With DirectNBody, we can also specify a background or external potential to integrate all orbits in. To do this, you can optionally pass in an external potential as a potential object to DirectNBody. Here, as an example, we’ll repeat a similar integration as above, but (1) add a positional offset of the initial conditions from the origin, and (2) specify an external potential using the MilkyWayPotential class as an external potential:

>>> external_pot = gp.MilkyWayPotential()
>>> w0_1 = gd.PhaseSpacePosition(pos=[10, 0, 0] * u.kpc,
...                              vel=[0, 200, 0] * u.km/u.s)
>>> w0_2 = gd.PhaseSpacePosition(pos=w0_1.xyz + [10., 0, 0] * u.pc,
...                              vel=w0_1.v_xyz + [0, 5, 0] * u.km/u.s)
>>> w0 = gd.combine((w0_1, w0_2))
>>> pot1 = gp.HernquistPotential(m=1e7*u.Msun, c=0.5*u.kpc, units=galactic)
>>> particle_pot = [pot1, None]
>>> nbody = DirectNBody(w0, particle_pot, external_potential=external_pot)
>>> orbits = nbody.integrate_orbit(dt=1e-2*u.Myr, t1=0, t2=1*u.Gyr)
>>> fig, ax = plt.subplots(1, 1, figsize=(5, 5)) 
>>> _ = orbits[:, 0].plot(['x', 'y'], axes=[ax]) 
>>> _ = orbits[:, 1].plot(['x', 'y'], axes=[ax]) 

(Source code, png, pdf)

../_images/nbody-2.png

In this case, the orbits are so similar that it is hard to tell whether the test particle is actually bound to the secondary mass. Let’s instead now plot the position in the x-y plane of particle 2 relative to particle 1. This will look strange because we have not transformed to the frame of particle 1, but it should give us a sense of whether particle 2 is bound or unbound to this mass:

>>> dxyz = orbits[:, 0].xyz - orbits[:, 1].xyz
>>> fig, ax = plt.subplots(1, 1, figsize=(5, 5)) 
>>> ax.plot(dxyz[0], dxyz[1]) 

(Source code, png, pdf)

../_images/nbody-3.png

From this, it looks like particle 2 is indeed still bound to particle 1 as they both orbit within the external potential.

gala.dynamics.nbody Package

Classes

DirectNBody(w0, particle_potentials[, …])

Perform orbit integration using direct N-body forces between particles, optionally in an external background potential.

Class Inheritance Diagram

Inheritance diagram of gala.dynamics.nbody.core.DirectNBody