MS2 UE4: Astrophysics with the computer

Hubert Baty, Joachim Köppen, Strasbourg, 2012/13

Please note

• choose your subject before the session on 9 Nov. 2012, so that you can start working on it right away ... time runs quickly ... and the number of afternoons is limited ...
• do not hesitate to contact Hubert Baty or me (joachim.koppen -at- astro....) for more information about the subjects, and if you have any questions! Ne hesitez pas de nous contacter!
• in this TP no question is stupid
• on parle le français ou l'anglais, comme vous le preferez!

• the programming language is your choice: you can use FORTRAN, C, or JAVA
• for those who like to use JAVA, I shall give you two graphics libraries which permit you to create x-y plots in the usual scientific form. One library, MorganPlot, is very good for looking at results after they have been computed; the other, JKZoomPlot, is better for displaying data during the simulation (such as n-body simulations). Contact me for any questions.

All projects are of similar difficulty and complexity, but they may differ in numerical effort (I give a rough indication by * to ***** from 'straightforward' to 'tough', but all have their numerical delicacies!). All leave space for exploration of the physics or the numerics, and all are suitable for modifications and additions, depending on your personal interest. And all can be done within the time available.

Subjects supervised by Hubert Baty:

• Exploration numérique du problème à 3 corps restreint en interaction gravitationnelle (**)
• Développement d'un code hydrodynamique 1D: traitement des chocs (**)
• Etude numérique de l'équation de Ginzburg-Landau: instabilité convective et absolue (**)
• Prise en main et tests d'un code 1D/2D MHD (***)

Subjects supervised by Joachim Köppen:

• Onedimensional Hydrodynamics: compute the time evolution of a fluid in a tube, starting with given initial profiles of density and velocity. Showing the formation of shocks and rarefraction waves. With the addition of the energy equation, the temperature profile can be followed, giving rise to the contact discontinuity. The program can be modified for spherical geometry and with the additional force of self-gravity, the collapse of a protostar or -galaxy can be calculated. (*** - Explicit solution of a system of partial differential equations, artifical viscosity, conservation of mass and momentum). The basic project can be directed to several different applications and extensions, as described in the dossier.
  the applet

• Smoothed Particle Hydrodynamics: compute the time evolution of a fluid in a tube, starting with given initial profiles of density and velocity. This will show the formation of shocks and rarefraction waves. With the addition of the energy equation, the temperature profile can be followed, giving rise to the contact discontinuity. This problem is treated by a particle-based method: A large number of particles moves in space following certain acceleration laws in such a way that their smoothed density behaves as closely as possible the behaviour of the density in a real fluid. Because there are no fixed spatial grids on which the functions are evaluated, the program can easily be modified for two and three space dimensions, without any constraint on the geometrical symmetry. (*** - solution of a system of ordinary differential equations). The basic project can be directed to several different applications and extensions, depending on the time available.
  the applet

• N-Body Simulations: compute the evolution of the positions of one or several planets or stars due to their gravitational interaction with the central sun or the other stars. (* - Solution of a system of ordinary differential equations with constant and adaptive time steps, treating the singularity at zero separation). The basic project can be directed towards several different applications, such as mentioned in the dossier. For example, one can verify the Jeans criterion for the gravitational collapse of a cloud of stars or galaxies.
  Orbits in JavaScript
  Two Body applet
  Three Body applet
  Special cases: Three Stars of Equal Mass

• Saturn's Ring: write a code to perform full n-body simulations which compute the evolution of the satellites in orbit around a planet due to their gravitational interaction with the planet and among each other. This will check a criterion for stability of such a ring of moons, found analytically by J.C.Maxwell in 1856. (** - Solution of a system of ordinary differential equations with constant and adaptive time steps). The basic project can be extended towards several different more general cases, such as mentioned in the dossier.
  links to a JAVA applet also doing the full n-body problem

• Orbits under MOND: The postulate of mysterious dark matter - which nobody has identified yet - to explain the observed rotation curves of galaxies can be avoided if we accept that Newtonian dynamics is modified under the large dimensions in the realm of galaxies. Let us write a code that explores what type of orbits are possible if we suppose the MOND formulation, and compare them with the usual Kepler orbits of standard Newtonian gravity. We can treat the one-body problem, where small test particle moves in a fixed potential, and the two-body problem of two masses orbiting each other. More complicated cases are unfortunately beyond the simple methods we can employ! (** - Solution of a system of ordinary differential equations with constant and adaptive time steps).
  links to JAVA applets which do the one- and two-body problem in Newtonian gravity and MOND

• Circular and Escape Velocities from a Galactic Disk: derive these velocities from a given mass distribution of stars in a galaxy. (*** - Integration of the contributions to the force acting on a test star by all other mass elements. Treating the singularity of Newton's gravity law at zero distance. Integration of the force out to infinity to get the potential)
• Finding the Gravitational Potential of a Galaxy: compute the gravitational potential in and around a given mass distribution. (**** - Solution of Poisson's equation by relaxation methods (Jacobi, Gauss-Seidel, Simultaneous Over-Relaxation); also possible: by solution of system of linear equations by Gaussian elimination). This project can be split up in several projects with different methods and/or models.
• Ram Pressure Stripping of Disk Galaxies: in a simplified model we consider the fate of a point mass held in a potential which is subjected to an impulsive external force, to determine under which conditions it escapes from the potential. (* - Solution of ordinary differential equations, systematic exploration of the parameter space)
  the applet

• Monte Carlo Radiative Transfer: We follow the flight of a photon through a medium, starting at its emission at the source, going through all absorptions and re-emissions, until it escapes and reaches the observer. If we do this with a sufficiently high number of photons, we can compute the observable intensity for any arbitrary geometry of the source and the interacting medium and for any amount of scattering - and any law of scattering. Let us try out this powerful method to solve radiative transfer problems! (** - Monte Carlo simulation, transformation of uniform random numbers into those distributed with given functions)
  the applet

• Temperature Fluctuations in Ionized Nebulae: compute the temperature-sensitive intensity ratio of the forbidden lines [O III] 5007/4363 Angstrom from a nebula in which the temperature fluctuates with a certain amplitude about a central value, and derive the apparent "mean temperature". Comparison with the true mean value illustrates that in the presence of such fluctuations the interpretation of nebular spectra is systematically affected. (** - Solution of the ionic rate equations by Gaussian elimination, Monte Carlo simulation of the temperature, transformation of uniform random numbers into those distributed with a given function, fitting a formula to numerical results)
  the applet

• Time-dependent Ionization of Gaseous Nebula: HII-regions are the ionized interstellar gas around hot main sequence stars which live only a few million years. Planetary nebulae are the expanding ejecta from a star which becomes very hot before it becomes a white dwarf. The hot stars ionize the gas, and normally we assume ionization equilibrium to interpret the emission line spectrum. But what happens when such a star switches on (or off) very rapidly? Under which conditions we should expect deviations from equlibrium? Let us follow the evolution of the ionization for any given time-history of the ionizing stellar flux, including those from published stellar evolutionary tracks. (* - Solution of ordinary differential equations, adaptive time step, interpolation, systematic exploration of parameter space) The basic project can be extended towards other cases, such as mentioned in the dossier.
  the applet

• Deflection of Cosmic Rays by Earth magnetic field: compute the trajectories of cosmic ray particles in the magnetic field of the Earth, starting from various initial positions - far from the Earth - and different initial velocities, in order to find out under which conditions the particle can reach the Earth surface. With these simulations one can show which parts of the Earth are more hazardous to these rays, and one can make a world map map of the probability, i.e. the cosmic ray flux as a function of energy. These simulations can be compared to the results of the pioneering study of this subject, done in the 1930s, long before computers (*** - Solution of ordinary differential equations, automatic time stepping, systematic exploration or Monte Carlo simulation for the initial positions and velocities, transformation of uniform random numbers into those distributed with a given function)
  a basic applet

• Excitation in and Emission from a Molecular Cloud: The CO molecules in an interstellar cloud are excited by collisions with hydrogen molecules into excited states, from which they relax by producing line emission in the microwave region. The emission can be absorbed by other CO molecules, leading to their excitation. The emission from the entire cloud is the result of the balance between excitations and de-excitations of the 'average' molecule. Under certains conditions of gas density and size of the cloud the upper levels of the CO molecule can be overpopulated, and the cloud will become a MASER, which is able to amplify any radiation that comes in from the exterior, e.g. the cosmic microwave background. (** - Solution of system of linear equations; iteration between radiation field and molecular states)
  the applet

• Perturbations of planetary orbits: Newton's law of gravity predicts that the orbits of two bodies moving around each other, such as a planet around the Sun, are perfect conic sections. However, the presence of another planet will cause the first planet to deviate from the Keplerian orbit. In our Solar System, the massive planets Jupiter and Saturn have such an effect on all other planets, and also on any spaceprobe that we send away. The presence of Neptune was predicted by LeVerrier from its perturbing effects on Uranus' orbit. With our present computers it is easy to compute numerically any planetary orbits, and also include the forces from the other planets. Thus we can compute these rather small perturbations by comparing the actual trajectory with the hypothetical one, where we had neglected these additional forces. We shall do this for a simplified system of two planets moving around the sun, in order to focus on the essential effects. (** - Solution of ordinary differential equations, automatic time stepping, assuring low numerical errors)
  the applet

• Precession of orbits: The extraordinary property of Newton's law of gravity is that all orbits (in the two-body problem) are closed orbits, perfect conic sections. However, any deviation from this perfect case causes elliptical orbits to precess, that is their major axis performs a slow revolution about the central body. Newton himself had worried about deviations from the inverse-square law of the force, and shown how the precession rate is linked to such a deviation. We know that precession can occur due to perturbations by other planets, deviations of the central body's shape from a sphere, any external tidal forces and due to relativistic effects ... In more recent times, observations of galactic rotation curves suggest either that we have to accept the presence of enigmatic dark matter if we want to stick to Newton's gravity law or we have to modify Newtonian dynamics. The latter approach (MOND) implies that orbits on galactric scales would deviate from Kepler orbits and show a precession. In our exercise, we want to explore under which situations this precession would make itself detectable. Could we measure it by observing certain binary stars with high precision? Which kind of binaries and which precision would be needed? (** - Solution of ordinary differential equations, automatic time stepping, assuring low numerical errors)
  a basic applet

• Bayesian Estimation of Parameters: Measuring the Radio Flux from the Moon: The temperature on the surface of the moon can be measured by observations at radio waves: The flux is proportional to the temperature. This can be done by executing a drift scan, where the moon passes through the main lobe of the radio telescope's antenna, due to the rotation of the earth. The signal goes up, then down again, in a bell-shaped curve. Its height would give the flux ... however, the data are noisy, as observational data often is! A good and efficient method to extract as much information as possible from a data set is to compare it with a model for the physics, the observational process, and the noise. In a systematic way, we scan through all possible combinations of relevant parameters. By assuming a certain type of noise (Gaussian) we can compute the likelihood that the observed data is the result adding noise to the model with a given set of parameters, and hence we can construct maps of the probability for the model parameters. In this way, we can find the most probable height of the curve - which will give the radio flux - as well as its error bar, i.e. the range of its value for a given percentage. (* - generation of random numbers, efficient scanning of multi-dimensional parameter space, interpolation, numerical integration)
  an applet (fit Gaussian to noisy data)

• Ionization structure of an HII region: Stars hotter than about 30000 K emit plenty of photons with energies above 13.6 eV, the ionization energy of hydrogen. As a consequence, any gas around the star will be ionized up to some distance which is determined by the actual number of hard photons emitted per second. We place such a star in a large gas cloud of uniform density, and compute the state of ionization of hydrogen as a function of distance from the star. As the distance increases the photon flux decreases, and so does the ionization, until the gas becomes neutral. It is possible to compute also the ionization of helium or other heavy elements, and to compute the strength of emissions lines. (** - Solution of ordinary differential equations, automatic radial stepping, integration over stellar spectrum, fix point iteration)

• The quantum states of the Hydrogen atom: Solving Schrödinger's equation for an electron in the electrical field of a proton yields the wave function of the electron, and gives the bound states of the electron. With the wave functions one can compute the transition probabilities for the spectral lines, such as the Balmer lines. (* - Solution of ordinary differential equations, automatic radial stepping, scanning through the energies and iteration to find the eigenstates)
  a basic JAVA applet

• Internal structure of a star: stars are spherical gas clouds in hydrostatic equilibrium, held together and compressed by their own gravity, and with densities and temperatures high enough for thermonuclear reactions. With these known physical processes and using the known properties of the matter (equation of state, opacity, energy production rate) it is possible to compute the structure of stars of given mass. With guess values for the density and temperature at the centre one integrates the equations radially outwards, and changes the guess values until the conditions of low density and temperature at the surface are met. (** - Solution of ordinary differential equations, automatic radial stepping, manual or automatic adjustment of initial conditions in order to meet surface conditions)
  a simple JAVA applet

• Evolution of viscous disks: When matter falls into the gravitational well of a compact object, like a black hole, a star, or the central region of a galaxy, it forms a disk, due to its initial angular momentum. In this disk of dense gas, hydrodynamic interactions between the different layers apparently give rise to viscosity, and conversion of kinetic energy into thermal energy. In this way the gas can lose the angular momentum by emission of radiation. We shall consider the thin disk model as it is used to interpret the behaviour of X-ray binaries. It is possible to add the recipes that explain the periodic outbursts of cataclysmic variables. (*** - Solution of partial differential equations)

• Interpretation of emission lines: radial velocity and internal kinematics of a planetary nebula. Emission lines from ionized nebulae are usually optically thin, therefore their shape tells about the distribution of radial velocities on the line of sight. The wavelength of the line centre give the radial velocity of the nebula with respect to us, which is useful information of its motion within the galactic context. The line width gives information about the internal motions of the gas in the nebula, i.e. thermal and turbulent motions as well as any systematic motion of the gas with respect to the centre of the nebula, such as expansion of the nebular shell. We have some real observational data, consisting of nine spectra taken at different positions over the face of a planetary nebula. To obtain the parameters of an observed emission line, we compare it with gaussian profiles of systematically changed parameters (width, height, central wavelength etc) to compute the likelihood of this fit. This Bayesian techniques allows to determine the probability for e.g. the width of the observed line, hence its most probable value and the confidence interval. (* - Generation of random numbers, Generation of simulated noisy observational data, Scanning through a large parameter space, Integration of marginal probability distributions)
  Observational data: part 1, part 2

• Radio continuum emission from the Sun: Given a simple model for the solar atmosphere - i.e. the density and temperature as a function of height above the photosphere - one can compute the emission from the sun at various radio wavelengths. Since the atmosphere is a plasma, the frequency-dependent opacity due to the electrons need to be taken into account. The emission itself is due to free-free emission by the electrons (or "Bremsstrahlung"). The emitted intensity is obtained by integrating over various lines of sight, thus allowing to construct a radio-image of the sun. This project could be developed further by taking into account the refraction of the waves in the solar atmosphere ... (* - Computation of opacities and emissivities, integration along lines of sight through a spherical atmosphere)
  a simple JAVA applet

• Evolution of a spherical gas cloud embedded in hot gas: What happens to a gas cloud in a thin hot medium? Under certain conditions the gas will get heated and the cloud will evaporate. Under other conditions the hot gas at the surface could cool because of thermal conduction, and condense onto the cloud. There are a number of situations in the interstellar medium and in winds from stars where the fate of such a blob of gas is of interest. Here we want to solve the fluid equations for a one-dimensional system in order to see whether we can characterize the physical conditions for condensation and evaporation, and compare our findings with criteria by other authors. (*** - Solution of partial differential equations, automatic time stepping)

• Formation of absorption lines in a stellar atmosphere: The decrease of temperature with height in a stellar atmosphere gives rise to spectral lines to appear in "absorption" in the spectrum of the star. For a simple formulation of the stellar atmosphere - the height-dependency of the temperature - and a simple representation of the optical properties of the gas - continuum and line absorption - one can compute the emergent spectrum. This allows to study how lines are formed and how their strength depends on the parameters of the atmosphere and the atom. (** - Computation of absorption coefficients, optical depths, Integration along the line of sight through the atmosphere, efficient scanning of the line parameter)
  simple JAVA applet

• Dynamical evolution of an HII region: When a massive star is born and lights up in the interstellar gas that surrounds it, the emergent photons will start ionizing this gas. The gas becomes transparent and therefore the photons can ionize gas further away. This progressive ionization of its neighbourhood continues until an equilibrium is reached between the star's photon output and the recombinations taking place in the ionized region, called a Strömgren sphere. (**** - Solution of ordinary differential equations, numerical integration over stellar flux, iteration for ionization, automatic predictor/corrector method for radial integration, automatic time stepping)

• Chemical evolution of the solar neighbourhood: By thermonuclear fusion in their central parts, stars are not only generating energy that makes them shine, but they also produce heavy elements from hydrogen and helium. While massive stars produce most elements like oxygen and iron, intermediate mass stars also contribute significantly to helium, carbon, and nitrogen, which are ejected into the interstellar medium with an appreciable time delay because of the longer life times of these stars. We want to construct a code that computes the evolution of the chemical composition of the interstellar gas dues to the birth and death of stars, taking into account the details of the production of elements in stars of various mass ranges. Other ingredients, such as the inflow of gas of arbitray composition or the outflow of interstellar gas due to supernova explosions etc can be added. (*** - Solution of ordinary differential equations, numerical integration of the contributions to the elements, interpolation in tables)
  Applet: One Zone Model

  Simple Infall Model

  One Zone Model with Gas Flows

  Applet with Exponential Infall

  One Zone Model: Oxygen and Nitrogen

  Applet with fully detailed treatment (a bit complex)

• Big Bang Nucleosynthesis (underway)

| to my HomePage |

last update: 12 March 2013 J.Köppen