In physics and astronomy, an N-body simulation is a simulation of a dynamical system of particles, usually under the influence of physical forces, such as gravity (see n-body problem for other applications). N-body simulations are widely used tools in astrophysics, from investigating the dynamics of few-body systems like the Earth-Moon-Sun system to understanding the evolution of the large-scale structure of the universe. [1] In physical cosmology, N-body simulations are used to study processes of non-linear structure formation such as galaxy filaments and galaxy halos from the influence of dark matter. Direct N-body simulations are used to study the dynamical evolution of star clusters.
The 'particles' treated by the simulation may or may not correspond to physical objects which are particulate in nature. For example, an N-body simulation of a star cluster might have a particle per star, so each particle has some physical significance. On the other hand, a simulation of a gas cloud cannot afford to have a particle for each atom or molecule of gas as this would require on the order of 1023 particles for each mole of material (see Avogadro constant), so a single 'particle' would represent some much larger quantity of gas (often implemented using Smoothed Particle Hydrodynamics). This quantity need not have any physical significance, but must be chosen as a compromise between accuracy and manageable computer requirements.
Dark matter plays an important role in the formation of galaxies. The time evolution of the density f (in phase space) of dark matter particles, can be described by the collisionless Boltzmann equation
In the equation, is the velocity, and Φ is the gravitational potential given by Poisson's Equation. These two coupled equations are solved in an expanding background Universe, which is governed by the Friedmann equations, after determining the initial conditions of dark matter particles. The conventional method employed for initializing positions and velocities of dark matter particles involves moving particles within a uniform Cartesian lattice or a glass-like particle configuration. [2] This is done by using a linear theory approximation or a low-order perturbation theory. [3]
In direct gravitational N-body simulations, the equations of motion of a system of N particles under the influence of their mutual gravitational forces are integrated numerically without any simplifying approximations. These calculations are used in situations where interactions between individual objects, such as stars or planets, are important to the evolution of the system.
The first direct gravitational N-body simulations were carried out by Erik Holmberg at the Lund Observatory in 1941, determining the forces between stars in encountering galaxies via the mathematical equivalence between light propagation and gravitational interaction: putting light bulbs at the positions of the stars and measuring the directional light fluxes at the positions of the stars by a photo cell, the equations of motion can be integrated with effort. [4] The first purely calculational simulations were then done by Sebastian von Hoerner at the Astronomisches Rechen-Institut in Heidelberg, Germany. Sverre Aarseth at the University of Cambridge (UK) has dedicated his entire scientific life to the development of a series of highly efficient N-body codes for astrophysical applications which use adaptive (hierarchical) time steps, an Ahmad-Cohen neighbour scheme and regularization of close encounters. Regularization is a mathematical trick to remove the singularity in the Newtonian law of gravitation for two particles which approach each other arbitrarily close. Sverre Aarseth's codes are used to study the dynamics of star clusters, planetary systems and galactic nuclei.[ citation needed ]
Many simulations are large enough that the effects of general relativity in establishing a Friedmann-Lemaitre-Robertson-Walker cosmology are significant. This is incorporated in the simulation as an evolving measure of distance (or scale factor) in a comoving coordinate system, which causes the particles to slow in comoving coordinates (as well as due to the redshifting of their physical energy). However, the contributions of general relativity and the finite speed of gravity can otherwise be ignored, as typical dynamical timescales are long compared to the light crossing time for the simulation, and the space-time curvature induced by the particles and the particle velocities are small. The boundary conditions of these cosmological simulations are usually periodic (or toroidal), so that one edge of the simulation volume matches up with the opposite edge.
N-body simulations are simple in principle, because they involve merely integrating the 6N ordinary differential equations defining the particle motions in Newtonian gravity. In practice, the number N of particles involved is usually very large (typical simulations include many millions, the Millennium simulation included ten billion) and the number of particle-particle interactions needing to be computed increases on the order of N2, and so direct integration of the differential equations can be prohibitively computationally expensive. Therefore, a number of refinements are commonly used.
Numerical integration is usually performed over small timesteps using a method such as leapfrog integration. However all numerical integration leads to errors. Smaller steps give lower errors but run more slowly. Leapfrog integration is roughly 2nd order on the timestep, other integrators such as Runge–Kutta methods can have 4th order accuracy or much higher.
One of the simplest refinements is that each particle carries with it its own timestep variable, so that particles with widely different dynamical times don't all have to be evolved forward at the rate of that with the shortest time.
There are two basic approximation schemes to decrease the computational time for such simulations. These can reduce the computational complexity to O(N log N) or better, at the loss of accuracy.
In tree methods, such as a Barnes–Hut simulation, an octree is usually used to divide the volume into cubic cells and only interactions between particles from nearby cells need to be treated individually; particles in distant cells can be treated collectively as a single large particle centered at the distant cell's center of mass (or as a low-order multipole expansion). This can dramatically reduce the number of particle pair interactions that must be computed. To prevent the simulation from becoming swamped by computing particle-particle interactions, the cells must be refined to smaller cells in denser parts of the simulation which contain many particles per cell. For simulations where particles are not evenly distributed, the well-separated pair decomposition methods of Callahan and Kosaraju yield optimal O(n log n) time per iteration with fixed dimension.
Another possibility is the particle mesh method in which space is discretised on a mesh and, for the purposes of computing the gravitational potential, particles are assumed to be divided between the surrounding 2x2 vertices of the mesh. The potential energy Φ can be found with the Poisson equation
where G is Newton's constant and is the density (number of particles at the mesh points). The fast Fourier transform can solve this efficiently by going to the frequency domain where the Poisson equation has the simple form
where is the comoving wavenumber and the hats denote Fourier transforms. Since , the gravitational field can now be found by multiplying by and computing the inverse Fourier transform (or computing the inverse transform and then using some other method). Since this method is limited by the mesh size, in practice a smaller mesh or some other technique (such as combining with a tree or simple particle-particle algorithm) is used to compute the small-scale forces. Sometimes an adaptive mesh is used, in which the mesh cells are much smaller in the denser regions of the simulation.
Several different gravitational perturbation algorithms are used to get fairly accurate estimates of the path of objects in the Solar System.
People often decide to put a satellite in a frozen orbit. The path of a satellite closely orbiting the Earth can be accurately modeled starting from the 2-body elliptical orbit around the center of the Earth, and adding small corrections due to the oblateness of the Earth, gravitational attraction of the Sun and Moon, atmospheric drag, etc. It is possible to find a frozen orbit without calculating the actual path of the satellite.
The path of a small planet, comet, or long-range spacecraft can often be accurately modeled starting from the 2-body elliptical orbit around the Sun, and adding small corrections from the gravitational attraction of the larger planets in their known orbits.
Some characteristics of the long-term paths of a system of particles can be calculated directly. The actual path of any particular particle does not need to be calculated as an intermediate step. Such characteristics include Lyapunov stability, Lyapunov time, various measurements from ergodic theory, etc.
Although there are millions or billions of particles in typical simulations, they typically correspond to a real particle with a very large mass, typically 109 solar masses. This can introduce problems with short-range interactions between the particles such as the formation of two-particle binary systems. As the particles are meant to represent large numbers of dark matter particles or groups of stars, these binaries are unphysical. To prevent this, a softened Newtonian force law is used, which does not diverge as the inverse-square radius at short distances. Most simulations implement this quite naturally by running the simulations on cells of finite size. It is important to implement the discretization procedure in such a way that particles always exert a vanishing force on themselves.
Softening is a numerical trick used in N-body techniques to prevent numerical divergences when a particle comes too close to another (and the force goes to infinity). This is obtained by modifying the regularized gravitational potential of each particle as
(rather than 1/r) where is the softening parameter. The value of the softening parameter should be set small enough to keep simulations realistic.
N-body simulations give findings on the large-scale dark matter distribution and the structure of dark matter halos. According to simulations of cold dark matter, the overall distribution of dark matter on a large scale is not entirely uniform. Instead, it displays a structure resembling a network, consisting of voids, walls, filaments, and halos. Also, simulations show that the relationship between the concentration of halos and factors such as mass, initial fluctuation spectrum, and cosmological parameters is linked to the actual formation time of the halos. [5] In particular, halos with lower mass tend to form earlier, and as a result, have higher concentrations due to the higher density of the Universe at the time of their formation. Shapes of halos are found to deviate from being perfectly spherical. Typically, halos are found to be elongated and become increasingly prolate towards their centers. However, interactions between dark matter and baryons would affect the internal structure of dark matter halos. Simulations that model both dark matters and baryons are needed to study small-scale structures.
Many simulations simulate only cold dark matter, and thus include only the gravitational force. Incorporating baryons, leptons and photons into the simulations dramatically increases their complexity and often radical simplifications of the underlying physics must be made. However, this is an extremely important area and many modern simulations are now trying to understand processes that occur during galaxy formation which could account for galaxy bias.
Reif and Tate [6] prove that if the n-body reachability problem is defined as follows – given n bodies satisfying a fixed electrostatic potential law, determining if a body reaches a destination ball in a given time bound where we require a poly(n) bits of accuracy and the target time is poly(n) is in PSPACE.
On the other hand, if the question is whether the body eventually reaches the destination ball, the problem is PSPACE-hard. These bounds are based on similar complexity bounds obtained for ray tracing.
The simplest implementation of N-body simulations where is a naive propagation of orbiting bodies; naive implying that the only forces acting on the orbiting bodies is the gravitational force which they exert on each other. In object-oriented programming languages, such as C++, some boilerplate code is useful for establishing the fundamental mathematical structures as well as data containers required for propagation; namely state vectors, and thus vectors, and some fundamental object containing this data, as well as the mass of an orbiting body. This method is applicable to other types of N-body simulations as well; a simulation of point masses with charges would use a similar method, however the force would be due to attraction or repulsion by interaction of electric fields. Regardless, acceleration of particle is a result of summed force vectors, divided by the mass of the particle:
An example of a programmatically stable and scalable method for containing kinematic data for a particle is the use of fixed length arrays, which in optimised code allows for easy memory allocation and prediction of consumed resources; as seen in the following C++ code:
structVector3{doublee[3]={0};Vector3(){}~Vector3(){}inlineVector3(doublee0,doublee1,doublee2){this->e[0]=e0;this->e[1]=e1;this->e[2]=e2;}};structOrbitalEntity{doublee[7]={0};OrbitalEntity(){}~OrbitalEntity(){}inlineOrbitalEntity(doublee0,doublee1,doublee2,doublee3,doublee4,doublee5,doublee6){this->e[0]=e0;this->e[1]=e1;this->e[2]=e2;this->e[3]=e3;this->e[4]=e4;this->e[5]=e5;this->e[6]=e6;}};
Note that OrbitalEntity
contains enough room for a state vector, where:
Additionally, OrbitalEntity
contains enough room for a mass value.
Commonly, N-body simulations will be systems based on some type of equations of motion; of these, most will be dependent on some initial configuration to "seed" the simulation. In systems such as those dependent on some gravitational or electric potential, the force on a simulation entity is independent on its velocity. Hence, to seed the forces of the simulation, merely initial positions are needed, but this will not allow propagation- initial velocities are required. Consider a planet orbiting a star- it has no motion, but is subject to gravitational attraction to its host star. As a time progresses, and time steps are added, it will gather velocity according to its acceleration. For a given instant in time, , the resultant acceleration of a body due to its neighbouring masses is independent of its velocity, however, for the time step , the resulting change in position is significantly different due the propagation's inherent dependency on velocity. In basic propagation mechanisms, such as the symplectic euler method to be used below, the position of an object at is only dependent on its velocity at , as the shift in position is calculated via
Without acceleration, is static, however, from the perspective of an observer seeing only position, it will take two time steps to see a change in velocity.
A solar-system-like simulation can be accomplished by taking average distances of planet equivalent point masses from a central star. To keep code simple, a non-rigorous approach based on semi-major axes and mean velocities will be used. Memory space for these bodies must be reserved before the bodies are configured; to allow for scalability, a malloc command may be used:
OrbitalEntity*orbital_entities=malloc(sizeof(OrbitalEntity)*(9+N_ASTEROIDS));orbital_entities[0]={0.0,0.0,0.0,0.0,0.0,0.0,1.989e30};// a star similar to the sunorbital_entities[1]={57.909e9,0.0,0.0,0.0,47.36e3,0.0,0.33011e24};// a planet similar to mercuryorbital_entities[2]={108.209e9,0.0,0.0,0.0,35.02e3,0.0,4.8675e24};// a planet similar to venusorbital_entities[3]={149.596e9,0.0,0.0,0.0,29.78e3,0.0,5.9724e24};// a planet similar to earthorbital_entities[4]={227.923e9,0.0,0.0,0.0,24.07e3,0.0,0.64171e24};// a planet similar to marsorbital_entities[5]={778.570e9,0.0,0.0,0.0,13e3,0.0,1898.19e24};// a planet similar to jupiterorbital_entities[6]={1433.529e9,0.0,0.0,0.0,9.68e3,0.0,568.34e24};// a planet similar to saturnorbital_entities[7]={2872.463e9,0.0,0.0,0.0,6.80e3,0.0,86.813e24};// a planet similar to uranusorbital_entities[8]={4495.060e9,0.0,0.0,0.0,5.43e3,0.0,102.413e24};// a planet similar to neptune
where N_ASTEROIDS
is a variable which will remain at 0 temporarily, but allows for future inclusion of significant numbers of asteroids, at the users discretion. A critical step for the configuration of simulations is to establish the time ranges of the simulation, to , as well as the incremental time step which will progress the simulation forward:
doublet_0=0;doublet=t_0;doubledt=86400;doublet_end=86400*365*10;// approximately a decade in secondsdoubleBIG_G=6.67e-11;// gravitational constant
The positions and velocities established above are interpreted to be correct for .
The extent of a simulation would logically be for the period where .
An entire simulation can consist of hundreds, thousands, millions, billions, or sometimes trillions of time steps. At the elementary level, each time step (for simulations with particles moving due to forces exerted on them) involves
The above can be implemented quite simply with a while loop which continues while exists in the aforementioned range:
while(t<t_end){for(size_tm1_idx=0;m1_idx<9+N_ASTEROIDS;m1_idx++){Vector3a_g={0,0,0};for(size_tm2_idx=0;m2_idx<9+N_ASTEROIDS;m2_idx++){if(m2_idx!=m1_idx){Vector3r_vector;r_vector.e[0]=orbital_entities[m1_idx].e[0]-orbital_entities[m2_idx].e[0];r_vector.e[1]=orbital_entities[m1_idx].e[1]-orbital_entities[m2_idx].e[1];r_vector.e[2]=orbital_entities[m1_idx].e[2]-orbital_entities[m2_idx].e[2];doubler_mag=sqrt(r_vector.e[0]*r_vector.e[0]+r_vector.e[1]*r_vector.e[1]+r_vector.e[2]*r_vector.e[2]);doubleacceleration=-1.0*BIG_G*(orbital_entities[m2_idx].e[6])/pow(r_mag,2.0);Vector3r_unit_vector={r_vector.e[0]/r_mag,r_vector.e[1]/r_mag,r_vector.e[2]/r_mag};a_g.e[0]+=acceleration*r_unit_vector.e[0];a_g.e[1]+=acceleration*r_unit_vector.e[1];a_g.e[2]+=acceleration*r_unit_vector.e[2];}}orbital_entities[m1_idx].e[3]+=a_g.e[0]*dt;orbital_entities[m1_idx].e[4]+=a_g.e[1]*dt;orbital_entities[m1_idx].e[5]+=a_g.e[2]*dt;}for(size_tentity_idx=0;entity_idx<9+N_ASTEROIDS;entity_idx++){orbital_entities[entity_idx].e[0]+=orbital_entities[entity_idx].e[3]*dt;orbital_entities[entity_idx].e[1]+=orbital_entities[entity_idx].e[4]*dt;orbital_entities[entity_idx].e[2]+=orbital_entities[entity_idx].e[5]*dt;}t+=dt;}
Focusing on the inner four rocky planets in the simulation, the trajectories resulting from the above propagation is shown below:
Angular momentum is the rotational analog of linear momentum. It is an important physical quantity because it is a conserved quantity – the total angular momentum of a closed system remains constant. Angular momentum has both a direction and a magnitude, and both are conserved. Bicycles and motorcycles, flying discs, rifled bullets, and gyroscopes owe their useful properties to conservation of angular momentum. Conservation of angular momentum is also why hurricanes form spirals and neutron stars have high rotational rates. In general, conservation limits the possible motion of a system, but it does not uniquely determine it.
In astronomy, Kepler's laws of planetary motion, published by Johannes Kepler in 1609, describe the orbits of planets around the Sun. These laws replaced circular orbits and epicycles in the heliocentric theory of Nicolaus Copernicus with elliptical orbits and explained how planetary velocities vary. The three laws state that:
Rotation or rotational motion is the circular movement of an object around a central line, known as an axis of rotation. A plane figure can rotate in either a clockwise or counterclockwise sense around a perpendicular axis intersecting anywhere inside or outside the figure at a center of rotation. A solid figure has an infinite number of possible axes and angles of rotation, including chaotic rotation, in contrast to rotation around a fixed axis.
In physics, equations of motion are equations that describe the behavior of a physical system in terms of its motion as a function of time. More specifically, the equations of motion describe the behavior of a physical system as a set of mathematical functions in terms of dynamic variables. These variables are usually spatial coordinates and time, but may include momentum components. The most general choice are generalized coordinates which can be any convenient variables characteristic of the physical system. The functions are defined in a Euclidean space in classical mechanics, but are replaced by curved spaces in relativity. If the dynamics of a system is known, the equations are the solutions for the differential equations describing the motion of the dynamics.
In physics, angular velocity, also known as angular frequency vector, is a pseudovector representation of how the angular position or orientation of an object changes with time, i.e. how quickly an object rotates around an axis of rotation and how fast the axis itself changes direction.
In science, work is the energy transferred to or from an object via the application of force along a displacement. In its simplest form, for a constant force aligned with the direction of motion, the work equals the product of the force strength and the distance traveled. A force is said to do positive work if it has a component in the direction of the displacement of the point of application. A force does negative work if it has a component opposite to the direction of the displacement at the point of application of the force.
A trajectory or flight path is the path that an object with mass in motion follows through space as a function of time. In classical mechanics, a trajectory is defined by Hamiltonian mechanics via canonical coordinates; hence, a complete trajectory is defined by position and momentum, simultaneously.
Newton's law of universal gravitation states that every particle attracts every other particle in the universe with a force that is proportional to the product of their masses and inversely proportional to the square of the distance between their centers. Separated objects attract and are attracted as if all their mass were concentrated at their centers. The publication of the law has become known as the "first great unification", as it marked the unification of the previously described phenomena of gravity on Earth with known astronomical behaviors.
In classical mechanics, the two-body problem is to calculate and predict the motion of two massive bodies that are orbiting each other in space. The problem assumes that the two bodies are point particles that interact only with one another; the only force affecting each object arises from the other one, and all other objects are ignored.
The classical rocket equation, or ideal rocket equation is a mathematical equation that describes the motion of vehicles that follow the basic principle of a rocket: a device that can apply acceleration to itself using thrust by expelling part of its mass with high velocity and can thereby move due to the conservation of momentum. It is credited to Konstantin Tsiolkovsky, who independently derived it and published it in 1903, although it had been independently derived and published by William Moore in 1810, and later published in a separate book in 1813. Robert Goddard also developed it independently in 1912, and Hermann Oberth derived it independently about 1920.
In astrodynamics or celestial mechanics, an elliptic orbit or elliptical orbit is a Kepler orbit with an eccentricity of less than 1; this includes the special case of a circular orbit, with eccentricity equal to 0. In a stricter sense, it is a Kepler orbit with the eccentricity greater than 0 and less than 1. In a wider sense, it is a Kepler orbit with negative energy. This includes the radial elliptic orbit, with eccentricity equal to 1. They are frequently used during various astrodynamic calculations.
A numerical model of the Solar System is a set of mathematical equations, which, when solved, give the approximate positions of the planets as a function of time. Attempts to create such a model established the more general field of celestial mechanics. The results of this simulation can be compared with past measurements to check for accuracy and then be used to predict future positions. Its main use therefore is in preparation of almanacs.
In theoretical physics, Nordström's theory of gravitation was a predecessor of general relativity. Strictly speaking, there were actually two distinct theories proposed by the Finnish theoretical physicist Gunnar Nordström, in 1912 and 1913 respectively. The first was quickly dismissed, but the second became the first known example of a metric theory of gravitation, in which the effects of gravitation are treated entirely in terms of the geometry of a curved spacetime.
The lattice Boltzmann methods (LBM), originated from the lattice gas automata (LGA) method (Hardy-Pomeau-Pazzis and Frisch-Hasslacher-Pomeau models), is a class of computational fluid dynamics (CFD) methods for fluid simulation. Instead of solving the Navier–Stokes equations directly, a fluid density on a lattice is simulated with streaming and collision (relaxation) processes. The method is versatile as the model fluid can straightforwardly be made to mimic common fluid behaviour like vapour/liquid coexistence, and so fluid systems such as liquid droplets can be simulated. Also, fluids in complex environments such as porous media can be straightforwardly simulated, whereas with complex boundaries other CFD methods can be hard to work with.
In relativity theory, proper acceleration is the physical acceleration experienced by an object. It is thus acceleration relative to a free-fall, or inertial, observer who is momentarily at rest relative to the object being measured. Gravitation therefore does not cause proper acceleration, because the same gravity acts equally on the inertial observer. As a consequence, all inertial observers always have a proper acceleration of zero.
In general relativity, Lense–Thirring precession or the Lense–Thirring effect is a relativistic correction to the precession of a gyroscope near a large rotating mass such as the Earth. It is a gravitomagnetic frame-dragging effect. It is a prediction of general relativity consisting of secular precessions of the longitude of the ascending node and the argument of pericenter of a test particle freely orbiting a central spinning mass endowed with angular momentum .
In celestial mechanics, a Kepler orbit is the motion of one body relative to another, as an ellipse, parabola, or hyperbola, which forms a two-dimensional orbital plane in three-dimensional space. A Kepler orbit can also form a straight line. It considers only the point-like gravitational attraction of two bodies, neglecting perturbations due to gravitational interactions with other objects, atmospheric drag, solar radiation pressure, a non-spherical central body, and so on. It is thus said to be a solution of a special case of the two-body problem, known as the Kepler problem. As a theory in classical mechanics, it also does not take into account the effects of general relativity. Keplerian orbits can be parametrized into six orbital elements in various ways.
The Einstein–Infeld–Hoffmann equations of motion, jointly derived by Albert Einstein, Leopold Infeld and Banesh Hoffmann, are the differential equations describing the approximate dynamics of a system of point-like masses due to their mutual gravitational interactions, including general relativistic effects. It uses a first-order post-Newtonian expansion and thus is valid in the limit where the velocities of the bodies are small compared to the speed of light and where the gravitational fields affecting them are correspondingly weak.
In physics, the n-body problem is the problem of predicting the individual motions of a group of celestial objects interacting with each other gravitationally. Solving this problem has been motivated by the desire to understand the motions of the Sun, Moon, planets, and visible stars. In the 20th century, understanding the dynamics of globular cluster star systems became an important n-body problem. The n-body problem in general relativity is considerably more difficult to solve due to additional factors like time and space distortions.
The Jeans equations are a set of partial differential equations that describe the motion of a collection of stars in a gravitational field. The Jeans equations relate the second-order velocity moments to the density and potential of a stellar system for systems without collision. They are analogous to the Euler equations for fluid flow and may be derived from the collisionless Boltzmann equation. The Jeans equations can come in a variety of different forms, depending on the structure of what is being modelled. Most utilization of these equations has been found in simulations with large number of gravitationally bound objects.