Properties of a good N -body solver • Good N -body solver depends on context • Galaxy dynamics: • Geometry of system often not known a priori (e.g., merger simulation —> where do the particles end up?) • Can have range of scales: dynamics in disk (<~ 100 pc) — halo dynamics (100s kpc), more extreme closer to center • Want solver that adapts to different geometries and scales • Star clusters: Geometry typically well-known (e.g., isolated cluster does not move) • Cosmological structure formation: • Cosmological box that does not change overall shape much over evolution • Rather smooth on large scales, periodic boundary conditions • Range of scales from ‘Universe’ (Gpcs) to individual clusters and galaxies (<~ Mpc)
Direct summation
Direct summation • Most naive approach: directly compute mutual gravity using • Need to determine N(N-1)/2 mutual distances —> O ( N 2 ) • Straightforward to implement, ‘easy’ to optimize implementation (incl. hardware) • Manifestly satisfies Newton’s third law —> conservation of momentum
Example implementation (note: inefficient, because we compute all N 2 mutual distances
Example implementation
Example implementation
Example implementation
Direct summation • Method of choice for collisional simulations: star clusters, particle disks around planets, … • Currently possible to go to ~ N =10 6 using GPUs • Not used for collisionless simulations: can use smoothness of the density field to motivate approximate algorithms
Tree-based gravity solvers
Tree-based gravity solvers • Introduced by Barnes & Hut (1986; still behind paywall!) • Use hierarchical tree to divide data set into cells • Approximate gravitational force by: • Group distant bodies into large cells • Group nearby bodies into small cells • Generally, leads to reduction of O ( N ) time / force calculation —> O (ln N )
Hierarchical tree • Recursively sub-divide data set into smaller and smaller cells • For gravity, oct-tree popular choice: • First level: single cell that covers data in 3D • Second level: subdivide into eight children, 2 per dimension • Third level: subdivide each child into 8 more children • Stop when cell has <=N max bodies (e.g., 1)
Hierarchical tree
Hierarchical tree
Hierarchical tree
Hierarchical tree
Hierarchical tree
Hierarchical tree: real quad-tree example
Hierarchical tree: real quad-tree example Automatically adapts resolution —> high resolution where there are many particles
Opening angle • Key to tree-based gravity calculation is approximation of distant bodies as a single body • Criterion based on an opening angle • Cell C with center-of-mass x C spans an angle of θ jC when seen from point x j : • For point x j all bodies in cell C can be approximated as single body if the cell appears smaller than a limiting opening angle • If cell C appears to be larger than the limit, we consider each of its children and whether they appear to be smaller than θ lim , if not, split
Opening angle: example
Opening angle: example “single-unit cells”
Opening angle: example “single-unit cells”
Gravity approximation • To compute the gravity of a given point x j , we consider all single-unit cells for x j and treat each as a single body (at the center of mass: red square in next slides) • Typically have O (ln N ) single-unit cells —> gravity calculation takes O (ln N ) time / body and O (Nln N ) for all bodies • This is much faster than direct summation
Gravity approximation
Gravity approximation
Gravity approximation • Final piece of the puzzle: how to approximate each cell as a single body • Simplest: single point-particle with mass = sum of all masses in cell • Gravitational potential from cell C is then • Can consider multipole expansion around the center of mass for higher accuracy (first multipole zero, bc center of mass)
Gravity approximation: example
Tree-based gravity scaling O(ln N )
Tree-based gravity scaling Tree setup takes O( N ln N ) [~N here]
Fast Multipole Method • Barnes & Hut tree-based gravity solver uses only the fact that force from all particles in a distant cell is similar • Further speed-up possible by using that force on particles that are close together is similar as well • Principle of fast-multipole method: expand gravity both around source (particle acting as source of gravity) and sink (particle on which we want the force) (Greengard & Rokhlin 1987; Dehnen 2002/2002 in astrophysics) • Leads to O ( N ) gravity calculation, but more complicated than Barnes/Hut (larger pre-factor)
Fast Multipole Method Dehnen & Read (2011)
Particle-mesh technique
Particle-mesh codes • Quite different approach solves the Poisson equation on a grid using Fourier transformation • Poisson equation becomes k 2 Φ ( ~ ~ k ) = 4 ⇡ G ⇢ ( ~ k ) • Easy to solve • Thus, algorithm becomes • Assign particles’ mass to grid points (some smoothing) —> O ( N ) • Fourier transform the density —> O ( n grid ln n grid ) • Solve Poisson above • Inverse fourier transform the potential —> O ( n grid ln n grid )
Particle-mesh codes • Advantages: • Fast because of fast-fourier-transform techniques • Ideal for cosmological simulations with periodic boundary conditions • Disadvantages: • Single grid-spacing —> not adaptive to regions of high density • Adaptive grids possible (without Fourier), but lead to energy conservation and issues in non-equilibrium situations
Scaling of collisional and collisions simulations
Numerical orbit integration
Numerical orbit integration • Motion of body under the force of gravity given by Newton’s second law or Hamilton’s equations • Cannot be solved analytically, except in very rare cases (point-mass, homogeneous sphere, …) • Thus, need to solve numerically (we have already done this a lot in the previous weeks) • Equation is an ordinary differential equation
Numerical orbit integration • Always solve Hamilton’s equations; system of first- order, linear ordinary differential equations (ODEs) • Can use any standard ODE solver or specialized solvers for Hamiltonian systems
General ODE solvers
General ODE solvers • Equation of motion in cartesian coordinates: • Or • Simplest solution: Euler’s method
Euler’s method Numerical Recipes, Chapter 17
Improving Euler’s method Numerical Recipes, Chapter 17
Improving Euler’s method • Use initial derivative to estimate mid-point • Use derivative at mid-point to advance the full time step • Example of a Runge-Kutta type method (specifically, second order Runge-Kutta; RK2) • Smaller errors than Euler, but two derivative evaluations
Fourth-order Runge Kutta • We can build higher-order Runge-Kutta type methods by using more estimated intermediate points to advance the point by 1 time step • Coefficients determined by requiring the solution to be precise up to given order • Most popular: fourth-order Runge Kutta
Fourth-order Runge Kutta Numerical Recipes, Chapter 17
Fourth-order Runge Kutta
Higher-order Runge Kutta methods • We can build even higher-order methods of the same type • Advantage of some higher orders: can produce multiple solutions from same set of intermediate points, accurate to different orders • Can use these multiple solutions to estimate error and use to adapt the time step —> adaptive time stepping • Examples: Dormand-Prince methods of order 5 and 8 • RK4 and Dormand-Prince methods are often used in galactic dynamics
Criticism of standard methods • Standard methods can give high precision at short time intervals • But long-term energy conservation is bad —> numerical errors build up in every time step • Useful for few dynamical times (up to hundreds), but not more (no solar system!)
Hamiltonian integration
Hamiltonian integration • Standard methods work by approximately solving the equations for the exact Hamiltonian • But we could also exactly solve the equations of motion for an approximate Hamiltonian • Such methods have the advantage that they exactly conserve dynamical quantities (phase-space volume, energy, …), but for the wrong system • If the approximate Hamiltonian is ~ the correct, this advantage typically leads to great long-term behavior
Discretizing the Hamiltonian • Let’s work in cartesian coordinates; Hamiltonian • Hamilton’s equations • Discretize the Hamiltonian in the following way with
Discretizing the Hamiltonian • Discretized Hamiltonian equal to • Thus, on timescales >> Δ t: approximate Hamiltonian ~ exact Hamiltonian • Hamilton’s equations for discretized Hamiltonian
Drift-kick • We can solve the equations of motion as a combination of drift and kick steps • From t= ε to t= Δ t- ε • From t= Δ t- ε to t= Δ t+ ε
Drift-kick • Limit of ε —> 0 • First step: drift • Second step: kick • Result similar to Euler, but evaluate force at final point —> modified Euler, also first-order
Leapfrog integrator • Simply by shifting by Δ t/2, we can build a second- order integrator that only uses a single force evaluation • This is the drift-kick-drift leapfrog integrator • Workhorse for collisionless simulations
Advantages of Hamiltonian integration • Typically have excellent long-term energy conservation • Connect to perturbation theory: can split Hamiltonian into largest contribution + perturbation • Important in solar-system (and exoplanet) integrations: can split Hamiltonian into point-mass potentials + small planet-planet perturbations (Wisdom-Holman integrator)
Energy conservation Bovy (2015)
Conservation of phase- space volume Bovy (2015)
Individual time steps • N -body simulations of collisionless systems have wide range of dynamical timescales • Use individual time step: particles with short dynamical times have short time steps, those with long dynamical times have long time steps • Somewhat unsolved problem of how to best choose these, standard methods far from optimal • Most efficient when particles are arranged in hierarchy of discrete time steps: block time-step scheme
Block time-step scheme BT08
and finally….
N -body modeling
N -body modeling • Now we have the tools to discuss N -body modeling • Two cases: collisional vs. collisionless • Collisional: solve evolution of system of N stars under Hamiltonian • Collisionless: solve collisionless Boltzmann equation
Collisional N -body modeling
Collisional N -body modeling • Need to follow equation of motion under the full Hamiltonian • Not easy to get around the following procedure: • Compute all forces using direct summation • Use numerical integrator to advance orbits • Typically use more advanced numerical integrators that use derivatives of the acceleration (+higher!) • Typically don’t use symplectic integrators for star clusters: need high precision on short timescales, less important on long • Need to explicitly deal with close encounters in efficient way
Collisionless N -body modeling
Recommend
More recommend
