Swinburne Astronomy Online Module Computers and Astrophysics Activity

Swinburne Astronomy Online Module: Computers and Astrophysics Activity: Astrophysical Simulations © Swinburne University of Technology

Summary In this Activity, we will learn about numerical simulations and their importance in astrophysics today. In particular, we will look at: (a) why we use numerical simulations in astrophysics (b) how to set up a numerical experiment (c) N-body codes for dynamic simulations (d) examples of N-body simulations (e) adding gas

Introduction Computer experiments - or numerical experiments - are mathematical models for simulating the behaviour of complex physical systems. There are various reasons why one may want to use a computer experiment, including: • to help design a complex device before using expensive technology to actually build it; • to gain information about physical systems that are not readily available in laboratory experiments (astronomy is a perfect example of this!); and • to help understand extremely complex system (astronomy is also a good example of this).

Dynamic systems that evolve with time, such as galaxies, have three coordinates of position (since they are 3 D objects) and also three coordinates of velocity (since they can move in all 3 directions). However, we can only observe two position coordinates and the line-of-sight velocity. observer So to determine the underlying structure of the Universe, we need to use computer simulations.

There are many situations in nature where we can view a physical system as a collections of particles that interact and obey certain laws. Astronomy offers many excellent example of this: the planets in the Solar System, the stars in the galaxy, and galaxies in the Universe obeying Newton’s law of gravitation. These particles have various physical attributes, such as mass, position and velocity, and each particle interacts with all others under the force of gravity. This is known as the N-body problem (where N is the number of particles used in the computer experiment), and we shall explore the N-body problem in this Activity.

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena chose the physical system that you wish to investigate

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena Mathematical model the physical system can be approximated by a mathematical model, which generally uses some simplifying assumptions to describe the workings of the physical system

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena Mathematical model Discretise the model the mathematical model must be converted from a continuous or differential equation into an algebraic approximation which the computer can solve. One must be careful at this step to introduce as few errors as possible. Both time and space must be discretised, which can produce huge numerical errors

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena Mathematical model Discretise the model Numerical algorithm chose the choice of discretisation is often related to the algorithm chosen to solve the discrete system. You need to be able to solve the discrete problem rapidly on your computer, else it’s all useless!

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena Mathematical model Discretise the model Numerical algorithm Computer program writing the code is what it’s all about! The code needs to be well engineered to make use of the available computing power, be easy of use and easy to modify

Setting up a Numerical Experiment The general steps involved in running your computer experiment are as follows: Physical phenomena Mathematical model Discretise the model Numerical algorithm Computer program Computer experiment finally you get to run your computer experiments, but you writingtothe code is what it’s all need know what you’re about! for, Thewhat codeyou’re needstrying to be to testing well engineered find, and how to to domake analysis useon of the data available that your computing experiment power, be easy of use produces. Pretty and pictures easy toare of modify vital at this stage! course

Discretising the Mathematical Model A major task when setting up a computer experiment is “discretising” the mathematical model. Mathematics usually involves continuous functions, whereas computers can only work with individual bits of information. For example, a computer does not really have a concept of a smooth sine wave: But it can handle a set of discrete points that represent the continuous sine wave. The computer can then deal with these discrete bits of information, adding them, multiply them, storing them etc.

Dynamical Simulations While there are many systems in astronomy which we can model using a computer code, in this Activity we will be concerned with dynamic simulations, which means the study of systems that evolve with time. We will not cover static simulations, which study a specific state of a system. In most astrophysical systems, the global properties of the system are determined by the force of gravity. While there also other forces involved - like electrical and chemical forces holding atoms and molecules together it is gravity that is the dominant long range force in the Universe.

Examples in Astronomy Dynamical systems in astronomy that are governed by gravity that we can model using numerical simulations include: • Solar System studies: – – the study of the dynamics and stability of the Solar System the motions and interactions of the asteroids in the asteroid belt modelling planetary rings the orbits of comets • Star systems: – the formation of stars through the collapse of giant molecular clouds – the dynamics of binary systems and small stellar clusters – the formation and evolution of globular clusters Globular cluster M 15

• Galaxy studies: – – the dynamics of spiral galaxies interactions between galaxies and galaxy mergers the evolution of galaxies galactic clusters and cluster dynamics • Cosmological studies: – the formation of structure in the Universe – the dark matter problem – the evolution of the Universe Merging galaxies NGC 2207 and IC 2163 These are just some of the systems that have been modelled over the past 20 years using dynamic computer simulation. So let’s have a look at how dynamic computer codes work.

The Mathematical Model There are two main parts to these dynamical codes: • the force calculation and • the time evolution. Both of these components can be described by a mathematical model. This means that we have a set of mathematical equations to solve which will tell us about the future state of our system given a set of initial conditions. Once we have our mathematical model, we need to discretise it and use a numerical algorithm to solve the discrete system on a computer. Let’s have a look at the two parts of the mathematical model that describe a dynamic gravitational system.

Gravitational Forces Newton’s universal law of gravity for the force between two bodies is given by Gm m 1 F= 2 r 2 where F is the resulting gravitational force, G is the universal gravitational constant, m 1 and m 2 are the masses of the two bodies, and r is the separation between the bodies. m 1 r F 1 F 2 m 2 Forces F 1 and F 2 are equal in size and opposite in direction.

What if we have five masses in our system? m 4 m 3 F 13 m 2 F 14 F 12 m 1 F 15 m 5 The force exerted on body 1 by the other 4 bodies would be given by: F 1 = G m 1 m 2 r 122 + G m 1 m 3 r 132 + G m 1 m 4 r 142 + G m 1 m 5 r 152

We also need to calculate the force on particle 2 due to the other 4 particles: and the force on particle 3 due to all the other particles: and the force on particle 4 due to all the other particles: m 4 m 3 m 2 m 1 m 5 and the force on particle 5 due to all the other particles: You can see that things get quite complicated very quickly and a computer could be of use!

We can write the force equation quite simply as: N Fi = j=1 G m i mj rij 2 where the (“sigma”) symbol means to sum or add over all other particles. We also use subscripts i and j to represent the particles. So for each N particle i we need to sum over all the other N-1 particles, which we assign with a j subscript. You can see that the mathematical model for gravitational force is quite easy to discretise for our N particles. There a couple of slight complications however. . .

Firstly, force is actually a vector quantity, not a scalar quantity, which means it has a magnitude and a direction. Our current equation for the force does N G mi mj rij F = - not tell us anything about the direction of i j=1 | rij 3 | the force, so that needs to be amended. Secondly, as the particles get closer together, the forces get larger. In fact, as particle i and j approach each other the denominator rij of the force equation approaches zero and so the force become infinite! Now computers do not like having to calculate infinite numbers and so we need to soften the gravity. Our equation becomes: N Fi = - j=1 G mi mj rij | rij 3 | +

This might not look like much of a change (since rij/|rij|3 ~ 1/rij 2), but it is! The funny vertical bar symbols |rij| are called “mod” or “absolute value” and they take only the positive value, for example |-3| = 3. That means that the remaining rij gives the direction of the force. If we imagine that the position rij is in 2 D space, with an x-direction y x component and a y-direction component, then rij will give y component rij the force its direction. x The softening parameter has to be carefully chosen - if it’s too large it will affect the physics of the system, and if it is too small then forces will become huge.

From the list of astrophysical examples of N-body system we gave earlier, you should be able to start seeing how each object (planet, star or galaxy) can act as a single “particle” with its own mass and position. And we’ve just seen how to calculate its force. In the case of the Solar System, we would only need the Sun and the nine planets - so N=10. To be more accurate, we could include the satellites, so let’s say N=80. And what about the asteroids? We might want to include 100, 000 of them, so now N=100, 080. You can see that the number of bodies in our system can rapidly increase!

If you want to model a small cluster of stars, you might want to use a few hundred or few thousand stars. For globular clusters you will need a few hundred thousand particles. This means that the force on each star in the globular cluster depends on the position and mass of the other hundred thousand stars in the cluster. You can see that a lot of calculating needs to be done! In fact, there will be N x (N-1) calculations, or O(N 2) (pronounced “order N squared”). Sometimes, however, you just cannot use a particle for every object in your system. You will not be able use a hundred million particles to represent all the stars in the Galaxy. So sometimes your “particles” actually represent a collection of physical bodies.

The Equation of Motion The evolution of the system is governed by the equations of motion, which tells us how the system evolves with time. Acceleration is defined as the rate of change of velocity, and velocity is the rate of change of position as follows: where a = acceleration v = velocity r = position t = time We can easily discretise this mathematical model by: v vf - vi r rf - ri a= and v = = = t t f - t i where i=initial f=final The symbol (“delta”) represents a small but finite change.

Newton’s second law relates force to motion via the equation F = m x a. We can use Newton’s law of gravitation to replace F by Fgrav, which leads to: We want to solve for both the position and velocity of the system at the next timestep. F v F t v = We can do this by solving the a = = m t m following equations: r v= r = v t t

If we now assume that we take a set timestep t between the old and new states of the system, and let v = vf - vi (where vf is the final velocity and vi the initial velocity) and do the same thing for the position, then F t the new velocity and position is vf = vi + given by: m rf = ri + v i t So once we have solved for the gravitational force F at the initial state of the system, we can solve for the system at the next timestep. There actually many fancier ways to discretise the equations of motion to give more accurate time stepping, but we won’t go into them just here.

Accuracy and Stability There are two more things that we need to be careful about: our choice of t and N. The timestep t controls the stability of the simulation. If we use a constant timestep t in the equations of motion, we will get large errors when two particles get close. So we need to use a numerical scheme with a variable timestep which will automatically decrease t if particles are too close and increase t as particles move apart. The particle number N gives the resolution of the simulation. Obviously we will need more that a few hundred particles to accurately resolve the dynamics of a galaxy, but of course there are limits as to how large N can be. Supercomputers can help us out here.

The N-body Algorithm So we’re now armed with our mathematical model for the gravitational force and equations of motion; we have a discrete algorithm for the mathematical model, and we’re ready to write our computer code to run our computer experiments. This is what our computer algorithm will look like: Set initial conditions Choose N and t, set initial particle mi, ri, vi, Fi Solve equations of motion ai = vi/ ti & vi = ri/ ti Calculate forces Fi = j Gmimj/rij 2 Update time counter tnew = told + t Output data rnew, vnew, Fnew, tnew

Examples of N-body Simulations In the next few slides, we will present some examples of N-body codes that have been used to simulate the dynamics of spiral galaxies, the formation of structure within the Universe, galaxy clustering and the formation of quasars. The number of particles in the simulation, N, provides the resolution, and increasing N allows more detail to be seen in the simulations. When N is much larger than about 30, 000, we need a supercomputer for its calculational speed, plus rapid storage of and access to the data.

Galaxy Dynamics This is a low resolution (N=7000) simulation of a spiral galaxy. The simulation began with a smooth distribution of particles, and the action of gravity sets up spiral arms within the galaxy. Colour represents density. Simulations like this have confirmed theories of spiral arm structures within galaxies and answered questions about the missing mass in the Universe.

Dark Matter in the Universe This is a low resolution (N=32, 768) simulation of the dark matter distribution in a box with a side length of 100 Mpc*, or 300 million light years. Each particle represents a “blob” of dark matter with a mass ten times that of the Milky Way. What are most interesting to note are the filamentary structures connecting regions where the particle density is much higher than the average density within the box. This box does not represent the entire Universe! * Mpc = mega parsec = 106 pc Click here for animation

Particle Number These images show the end state of a simulation of a cluster of galaxies comparable to the Virgo cluster. The colour represents density and the radius of the computational box is 4 Mpc. The image on the left is a low resolution (N=67, 000) CRAY simulation, while the image on the right is from a high resolution (N=1. 3 million!) parallel simulation run on a “cluster” supercomputer. As you can see, as more particles are used, more detail of the cluster’s structure can be seen.

Quasar Formation in the Early Universe One of the current problems in cosmology is this: simulations have shown that while structure forms very early in the Universe, galaxies tend to form quite late. Yet we know that quasars exist at very high redshifts. So how did they form? This sequence of images shows the formation of a quasar at a redshift of z=8 (which means very early in the history of the Universe) using a high resolution simulation.

Improving the Algorithms We’ve only considered N-body algorithms here, but there other types of algorithms currently used in astrophysics for solving for the gravitational force in a physical system. None of the simulations you have just seen actually use the “direct summation” technique that you have learned about in this Activity, where each of the N particles calculates its force due to the N-1 other particles. This requires N 2 calculations, and for N=5000, there are 25 million calculations to do, and for N=50, 000, 2. 5 billion calculations are needed! It just takes too long to run a simulation using this direct approach. And that’s just the gravitational part of the algorithm: we also need to integrate the equations of motion. If you have 50, 000 particles coming close together you have to decrease the timestep t. So computational astrophysicists spend a lot of their time writing very clever codes to calculate the gravity as quickly and efficiently as possible. This is usually done by approximating the forces from particles far away and accurately calculating (by direct summation) the forces of nearby particles.

Adding Gas We have looked at several examples of astrophysical N-body systems where we calculate the force on a set of gravitating particles that act as point masses. There are many cases in astronomy where we would also like to run computer experiments of the gas behaviour of systems. Examples include the collapse of molecular cores to form protostars and their surrounding disks; the modelling of gas jets emanating from young protostars; the expansion of shock waves into the interstellar medium from supernova explosions; and the study of the gas disks around cataclysmic variables and black holes. We might also want to add gas to our simulation of galaxies, so that we can see how the stars and gas in the arms of spirals interact.

The study of the motion of a fluid is called hydrodynamics. Fluids - both liquids and gases - are governed by pressure and viscous forces. Purely pressure forces are expansive forces, compared with the inward force of gravity. Viscosity, which depends on the properties of the gas, acts as a type of resistance in a fluid. Purely viscous forces are spreading forces. The total force on a particle of gas is given by: Ftotal = Fgravity + Fpressure + Fviscous This means that we need to include more physics into our computer code. We need a mathematical model for both pressure and viscosity, which we can then discretise.

Each particle in a gas simulation represents a collection of molecules or atoms, since of course we cannot represent each individual gas molecule, just like we can’t represent each individual star in a galaxy! In gas simulations, each particles represent a sphere of gas, which we can think of as a having a “radius of influence”. Any other particle in particle i’s sphere of influence is included in all force calculations with it. Those outside the sphere of influence are not. interaction radius not used to calculate the force on particle i i not used to calculate the force on particle i

Massive Protostellar Disks This is a simulation of a massive protostellar disk, which forms around a young stars as it collapses out of a molecular cloud. The disk is coloured by density and has the same mass as the central protostar. The gravity within the disk causes instabilities which create transient spiral arm features and small blobs to break off and merge. Could these blobs be forming planets? Click here for animation

Disks Around Young Binary Stars This next simulation shows the dynamics of gas in a circumbinary disk, which is a disk that surrounds a young binary star system. At the centre there are two stars in mutual orbit around each other, and the action of gravity sweeps out a hole in the disk, cutting it off from the central binary. You can see transient spiral arms occasionally feeding the young stars with gas. Click here for animation

Summary In this last Activity of the Unit we have seen how computer experiments are a fantastic tool in modern astronomy, allowing us to understand the structure and dynamics of many systems in the Universe. We have learned about how to set up an N-body computer experiment and seen several examples of their use in astronomy today. With the help of supercomputers, we can now increase the value of N into the millions, allowing the Universe to be studied in great detail. It is up the computational astrophysicists to write better algorithms and codes to make use of this great computing power!

Image Credits Sunflower galaxy M 63 - Satoshi Miyazaki, Suprime-Cam, Subaru Telescope, NOAJ http: //antwrp. gsfc. nasa. gov/apod/ap 000627. html Globular cluster M 15 - Hubble Heritage http: //heritage. stsci. edu/public/2000 aug 3/ngc 7078. html Interacting galaxies NGC 2207 and IC 2163 - Hubble Heritage http: //heritage. stsci. edu/public/99 nov 4/ngc 2207. html The planets - Welcome to the Planets, California Institute of Technology http: //pds. jpl. nasa. gov/planets/ Galaxy simulation - Sarah Maddison, Swinburne University of Technology (volume rendering by Robin Humble, CITA)

Image Credits Cosmological simulation and animation - by Chris Fluke of Swinburne University, using the Hydra N-body SPH code of H. M. P. Couchman & P. A. Thomas http: //hydra. mcmaster. ca/hydra/index. html Clustering simulation - HPCC, University of Washington http: //www-hpcc. astro. washington. edu/TSEGA/picture. html Quasar formation simulation - HPCC, University of Washington http: //www-hpcc. astro. washington. edu/TSEGA/picture. html Massive protostellar disk simulation - Sarah Maddison, Swinburne University of Technology Binary disk simulation - Sarah Maddison, Swinburne University of Technology

End of Activity and Unit