Sunday, June 12, 2016

Formaline: The Provenance of Computational Astrophysics Simulations and their Results (Part I)

Key to all scientific research is reproducibility. In areas of science that rely heavily on computer simulations, the open-source approach is essential for allowing people to check the details of the simulation code and for re-running simulations for reproducing results.

Open-source is essential, but it's not enough. Anybody involved in simulations knows that the code used for running the actual simulations often differs in non-trivial ways from the actually released source code. Even when we use version control software like git, too often do we run simulations with a version of the code that is either not tagged or has uncommitted changes. And this does not take into account that complex simulation codes typically have a myriad of parameters that are set individually for each simulation and parameter files typically don't make it into publications or git repositories. A similar issue exists also with the initial conditions (i.e. the input data) from which a simulation is started. Rarely are these available for others to reproduce results.

Provenance is a word used in computational science (and in some other research areas). It describes the full history of a (computational) scientific result. So the provenance of a simulation does not only include the released open-source code, but also the full information on the exact source code that was used, perhaps even the way it was compiled and where it was run, the parameters that were set for the simulation, and all initial conditions and other ingredients that were used in whatever form (and also their respective provenances).

We have to tackle provenance and reproducibility step by step. In this two-part series of posts, my focus is on preserving the actual source code that was used to run a simulation and produce specific results or to create some other (astro)physics result (for example, a neutrino-matter opacity table used in core-collapse supernova simulations). This can be accomplished by what I will generally call Formaline. Formaline is a chemical compound, a 10% solution of formaldehyde in water; used as a disinfectant or to preserve biological specimens. Formaline for code/simulations is to preserve code and simulation/physics metadata. Formaline is an idea by Erik Schnetter, Research Technologies Group Lead at Perimeter Institute for Theoretical Physics. Erik implemented Formaline for the first time for the Cactus Computational Toolkit some time in the early 2000s.

Using Formaline is super easy! (Python example)

Formaline comes in two forms and in both forms it preserves source code: (1) It can include source code in the binary of the simulation code and, upon execution, the code writes the source code (usually in a tar ball) into the output directory. (2) It can include source code and parameter files in HDF5 data files that are now very broadly used for simulation and physics inputs to simulations (such as a neutrino interaction table).

I'll talk about Formaline (1) in my next post. Here, I'll focus on how to implement Formaline (2) and include (and later extract) a source tar ball in/from an HDF5 data file. It turns out that this is extremely easy if you have the h5py and numpy Python packages installed.

Let's consider the case in which you have some standalone code that produces an HDF5 table called table.h5 and that is generated by some Fortran 90 code in and below directory src and that takes a top-level parameter file called parameters, a Makefile, a for machine-specific build flags, and a README file.

Here is the code that puts all source code into the HDF5 file:

So, after running this script, your HDF5 table table.h5 contains a dataset that holds the tar.gz ball of our code source, input, & parameter files. You can view the table of contents of any HDF5 file with h5ls [filename].

Now, let's say 5 years later, we find table.h5 and want to know how it was created. We can now go and extract the source/input/parameter tar.gz ball! Here is the script that does just that. It first creates an output directory called saved_code to make sure nothing is overwritten.

So that's a super easy way to make results reproducible -- just store the source code and all input/parameter files with the results! I think it's a no-brainer and we all should be doing it. HDF5 is a standard, well documented format that will be readable decades from now -- so your work is safe!

I wrote this implementation of Formaline for the NuLib neutrino interaction library and and a variant of it is bundled with NuLib and available on github at

In my next post, I'll talk about Formaline (1), which includes a source tar ball inside the executable. This is the original incarnation of Formaline that Erik Schnetter came up with.

Friday, April 15, 2016

A Beginner's Super Simple Tutorial for GNU Make

(GNU) Make has been around for a very long time, has been superseded (as some may argue) by CMake, but it is still an incredibly useful tool for... for what? For developing code on the command line. Yes, that's right. We still use editors and terminals to develop code. Integrated "windowy" development environments (IDEs) are great, but often overkill. And it's really hard to run them when you are developing remotely on a supercomputer.

But even a hard-core terminal shell & editor developer needs some help for productivity. Make is a key productivity tool for developing code.

Back in the old days (say just after punch cards were no more; way before my time), most people would write code in one big source file and then would just invoke the compiler to build the executable. So your entire code would be in a single file. Not very convenient for editing. Certainly not convenient for working collaboratively. When codes got bigger, it just became impractical. Compile times became huge, since every time a single line changed in one routine, the entire code needed to be compiled.

So people started breaking up their code into pieces, each of them living in its own source file. That's natural, since code is usually broken down into some kind of logical parts (e.g., subroutines and functions in procedural programming, classes/objects in object-oriented programming). So it makes sense to put logical units of the code into their own source files. But how do you get your executable built from an array of source files? How do you prevent that not everything has to be recompiled all the time?

This is where Make comes in!

Let's assume we are working with Fortran 90+ here and compile things using gfortran (if you like C better, just replace everything Fortran specific with C; or with whateva language you like). Let's say you have three source files.


Note that I am using upper case F90, instead the more common f90. The upper case tells the compiler to pre-process the source files before compiling them. In this way, we can use preprocessor directives, see here if you are interested: We won't use these directives here, but perhaps I'll write a post about them in the future.

So that you can follow along easily, let's actually write these routines:


program myprogram

  implicit none

  call subroutine_one
  call subroutine_two

end program myprogram


subroutine subroutine_one
  implicit none
  write(*,*) "Hi, I am subroutine one"
end subroutine_one


subroutine subroutine_two
  implicit none
  write(*,*) "Hi, I am subroutine two"
end subroutine_two

Naively, you could build your program this way:
gfortran -g -O3 -o myprogram main_program.F90 subroutine_one.F90 subroutine_two.F90
So every single time you compile, everything gets compiled. This produces an executable called myprogram. By the way, the options I am using in the compile line are: "-g" -- add debug symbols and "-O3" -- optimize aggressively (not really needed here, but I do it out of habit (I want fast code!)).

Recompiling is of course no problem if the routines don't do anything of substance as in our example. But let's for a moment assume that each of these routines does something incredibly complicated that takes minutes or longer to compile. You don't want to have to recompile everything if you made a change in only one of the source files. The trick to avoid this is to first build object files (file suffix ".o"), then to link them together and only remake those whose source files changed.

Let Make take care of figuring out what changed. You tell Make what to do by writing a makefile in which you specify what you want to make and what it depends on. Here is the makefile for our example. Call it Makefile, then Make will automatically find it when you type make on the command line. Here is what will be in your makefile:

myprogram: main_program.o subroutine_one.o subroutine_two.o
     $(FC) $(FCFLAGS) -o myprogram main_program.o subroutine_one.o subroutine_two.o

main_program.o: main_program.F90
     $(FC) $(FCFLAGS) -c main_program.F90

subroutine_one.o: subroutine_one.F90
     $(FC) $(FCFLAGS) -c subroutine_one.F90

subroutine_two.o: subroutine_two.F90
     $(FC) $(FCFLAGS) -c subroutine_two.F90

     rm -f *.o

There is a very clear structure to this. The first two lines are definitions of variables you are using later -- the compiler and the compiler flags (you could call these variables anything you want). Then you have lines that start with the target on the left, followed by a colon. After the colon, you list the things that the target depends on (the dependencies). This tells Make what needs to be done before working on the current target. So in our case, myprogram depends on a bunch of object files. The rules for making these sub-targets are specified subsequently. One important thing to remember: In the line where you tell Make what to do for a given target (once the dependencies are met and up-to-date), the first character on that line has to be a tabulator character (i.e., you have to hit your "tab" key). Only this gives the actual command the correct indentation (that's something historic, I guess). But clearly this is not a show stopper!

Note the "clean" target -- I added this so that you can easily clean things up.

So, for example, let's say you want to compile your code, but you want to compile it freshly from scratch. The first thing you do is say (on the command line) make clean. This will get rid of all previously compiled object files. Next you just type make and Make will compile and link everything for you! Now you go and work on one of the source files. You save it and then you want to recompile. Just type make. This time only the source file that you changed is recompiled and then linked together with the other object files to make the executable. Voilá!

I hope you see the benefits of using Make now! Make is of course much more powerful than this simple example can tell. You can read more about more advanced uses of Make here and here and here and at many other places that your favorite search engine will find for you.

Before I close, there is one more thing I want to mention: Make is also great for compiling LaTeX documents, in particular if you are using cross-referencing and a bibliography software like BibTeX.
Here is an example makefile for compiling a LaTeX document called paper.tex to paper.pdf with BibTeX. I assume that you have your BibTeX bibliography entries in a file called paper.bib:
paper.pdf: paper.tex paper.bib
    pdflatex paper.tex
    bibtex paper
    pdflatex paper.tex
    pdflatex paper.tex
This should work. Enjoy!

Saturday, January 30, 2016

Simulation vs Modeling

For a long time, I have been thinking about how a clear distinction could be made between "modeling" and "simulation." The answer to this question likely depends on the field of study. I will limit myself to the field of (computational) astrophysics. Many colleagues use "modeling" and "simulation" interchangeably, but I believe that a semantic distinction needs to be made.

For starters, "modeling" sounds less sophisticated/complicated/detailed/involved than "simulation." But this is just purely subjective. Let's go a bit beyond the sound of these words.

A good example that can help us appreciate the differences between simulation and modeling comes from gravitational wave physics: Advanced LIGO needs templates (predictions) of gravitational waves to find the same or similar waves in its data stream. Such predictions can come from detailed, approximation-free simulations that implement Einstein's equations on a computer ("numerical relativity simulations", see Alternatively, they can come from so-called post-Newtonian "models" or from phenomenological "models" that try to approximate the numerical relativity simulations' results. These models are a lot simpler and computationally cheaper than full numerical relativity simulations. So that's why they find frequent use. But their results are not quite as good and reliable as the results of numerical relativity simulations.

Another useful example comes from supernova simulations and modeling of observables: It is not currently possible to simulate a core-collapse supernova explosion of a massive star end-to-end from first principles. The process involves the very last stages of core and shell burning, core collapse, core bounce, the postbounce phase during which the stalled supernova shock must be revived (all occurring within seconds of the onset of collapse), and the long-term propagation (up to days of physical time!) of the re-invigorated shock through the stellar envelope. So, on the one hand, detailed simulations are used to study the mechanism that revives the shock. But these simulations are too computationally intensive to carry out for more than a few hundred milliseconds, perhaps a second in 3D. On the other hand, simpler (often spherically-symmetric [1D]) modeling is applied to predict the propagation, breakout, and expansion of the ejecta and the resulting light curves and spectra. These explosion lightcurve/spectral models (most of the time) start with fake ad-hoc explosions put in at some mass coordinate in various ways.

Here is a slide from a talk on simulation and modeling of gravitational wave sources that I gave at a recent Instituto de Cosmologia y Fisica de las Americas (COFI) workshop in San Juan, Puerto Rico:

I think that the items listed on this slide are broadly applicable to many problems/phenomena in astrophysics that are currently tackled computationally. Ultimately, what we want is to simulate and have fully self-consistent, reliable, and predictive descriptions of astrophysical events/phenomena. This will require another generation of supercomputers and simulation codes. For now, it is a mix of simulation and modeling.

Saturday, March 28, 2015

A super simple Introduction to version control with Git for busy Astrophysicists

Welcome, collaborator! You may be reading this, because you, I, and others are working on some kind of joint document. It could be a manuscript to be submitted to a journal, some kind of vision document, or (as is often the case) a proposal.

Sending around LaTeX files via email is so 1990s. Let's not do this. Let's also not use Dropbox. Dropbox is great for syncing files between computers and serves as a great backup solution. But it's terrible for collaborating on text documents -- it can't handle conflicts (two people editing the same file at the same time creates a second copy of the file in the dropbox; Dropbox won't merge text for you).

Let's move on and use version control.  It will make working together on the same document infinitely easier and will minimize the time we have to deal with merging different versions of the file we are working on together. There are multiple version control packages. For our project we will use git. You can learn more about git and version control in general on If you are unconvinced, read this: Why use a Version Control System?

Here is how it works:

(1) Preparations
  1. Make sure you have git installed on your computer (I presume you have either a Mac or a PC that runs a flavor of Linux). Just open a terminal, type git, hit enter, and see if the command is found. If not, you need to install git.

  2. Make and send me a password-protected public ssh key as an email attachment (not cut & paste!). The standard path and file name for this are ~/.ssh/ or ~/.ssh/ If you don't have such a key or you don't remember if it's password protected (or you dont remember your password...), type on the command line:
    ssh-keygen -t rsa
    and follow the instructions. By all means, give your key a strong password. Then send me the file ~/.ssh/ as an email attachment (do not cut & paste it into the email body) 

(2) Cloning the repository from the server
By now I have told you details about the repository and sent you the command needed to clone it from the server. Go into your preferred working directory (that's where you keep projects, papers, proposals etc.). Type:
This will clone the repo -- in other words, it will create a directory with the repository name in your working directory and will download the contents of the repo into it. If this does not happen and instead you are asked to enter a password, then your ssh key is not loaded. In this case, load your ssh key like this:
ssh-add ~/.ssh/id_rsa
Then try again. If this still does not do the trick or if you get an error message about no "ssh-agent" running, then try
This will spit out a few lines of text. Copy them and paste them back into the prompt, hit enter, then repeat
ssh-add ~/.ssh/id_rsa
This should do the trick!

(3) Pulling in changes.
Before starting to work on a section of the document, you need to make sure to get the most recent changes. You must absolutely always do this before starting to work! Go into the repository directory.
git pull -r 
(the "-r" is to make the pull a "rebase" -- don't worry about the meaning of this, just do it).  

(4) Committing and pushing changes.
Whenever you have completed a significant addition or change to the document -- this may be a completed or edited paragraph, something that takes no longer than about 30 minutes of work, you need to commit and push your changes so that others can pull them onto their computers.
git commit -m "Your commit message goes here"
git pull -r           # this is to pull in any changes by others!
git push
Please give sensible commit messages so that we others know what you changed!

(5) Looking at the history of changes and checking the status of your repo.
That's very easy! First pull in all recent changes via git pull -r. Then use git log to browse through the repository history. Try git status to see what the status of your repo is. This will tell you which files have changed and should be committed. It will also tell you which files you have not yet added (i.e. "untracked" files).

(6) Golden Rules
Working with version controlled documents is easy and painless as long as you obey two simple rules:
  1. Always git pull -r before starting to edit.
  2.  Extremely frequently:
    git add [FILE]
    git commit -m "[COMMIT MESSAGE]"
    git pull -r
    git push
If you stick to these rules, working collaboratively on the same document will be easy and conflicts will be rare.

(7) Adding additional files to the repo
That's easy!
git add [FILENAME]
git commit -m "Added file [FILENAME] to repo"
git pull -r
git push
(8) Resolving merge conflicts
Oh boy. Somebody was editing the same portion of text you were editing (did you perhaps not sufficiently frequently commit and push?!?). You now need to resolve the conflict. Let's do this by example. Let's say you were editing a file called You added, committed, then pulled and ended up with the following git merge error:
remote: Counting objects: 5, done.
remote: Total 3 (delta 0), reused 3 (delta 0) Unpacking objects: 100% (3/3), done.
23ae277..b13d53f master -> origin/master Auto-merging
CONFLICT (content): Merge conflict in
Automatic merge failed; fix conflicts and then commit the result.
Good luck! You now need to resolve the conflict. Let's look at the (hypothetical) file in question.
$ cat 
<<<<<<< HEAD 
>>>>>>> b13d53f9cf6a8fe5b58e3c9c103f1dab84026161
git has conveniently inserted information on what the conflict is. The stuff right beneath <<<<<< HEAD is what is currently in the remote repository and the stuff just below the ======= is what we have locally. Let’s say you know your local change is correct. Then you edit the file and remove the remote stuff and all conflict stuff around it.
$ cat 
Next you git add and git commit the file, then you execute git push to update the remote with your local change. This resolves the conflict and all is good!

Wednesday, May 14, 2014

Be Scientifically Limitless

Don't you misunderstand the title of this post!

I don't mean to say there should be no limits on how you, the reader of this post, do science. There are the very clear boundaries of the scientific method and of professional conduct that you must respect.

What I mean to say is: Do not limit yourself in what you can achieveDo not let yourself be limited by circumstances or your environment.

More often than I would like, I encounter the following self-limiting behavior:

A person is confronted with a problem (perhaps a scientific problem, or a coding problem, or something else). They realize that it is hard, perhaps harder than what they have done before. They also have no experience in solving such a problem. Perhaps they have seen others fail trying to solve similar problems. Some give up at this point.

Others might try for a bit, fail to make progress, get frustrated and quit trying. They often don't ask for help, for a variety of reasons, the two major ones being: (1) They might think they are just stupid/inexperienced and should learn to solve problems on their own. (2) They are afraid of what others might think about them if they ask for help and let "others solve their problems too frequently".

Always ask for help. Those that ask are the ones that do things and go places. If you don't ask for help, you limit what you can achieve.

Another thing I have seen too often: People tend to put people into comfortable little virtual mind boxes based on what they think about them and, perhaps, what their history has been. This then limits what the person in the box thinks they can and should do now and in the future. They remain limited in what they achieve to the boundaries of the box they have been placed in. You need to realize that these boundary conditions are made up. They are made up by people no smarter than you (despite of what you may think!). You can and should and must break through them. Liberate your thoughts and dreams!

Here are two very inspirational YouTube videos that are excerpts from a longer 1994 interview with Steve Jobs. You can think about this man what you want. The advice he gives is brilliant and invaluable.

On asking for help:

On not limiting yourself:

Sunday, April 27, 2014

Nature has it in for us: Supernovae from Massive Stars in 3D

[A version of this post has appeared in the Huffington Post. It is a more generally accessible discussion of our recent Astrophysical Journal Letters paper: Mösta, Richers, Ott et al., The Astrophysical Journal, 785, L29 (2014). ADS]

We don't know what precisely happens when a massive star -- about ten times the mass of our Sun or more -- first collapses and then goes supernova and leaves behind a neutron star or a black hole. The explosion expels the products of stellar evolution into the interstellar medium from which, ultimately, new stars and planets are made. These products, carbon, oxygen, silicon, calcium, magnesium, sodium, sulfur, and iron (among other elements), are the elements of life. It's arguably quite important to understand how this works, but we don't. At least we don't in detail.

The supernova problem is so complex and rich that computer simulations are crucial to even begin to formulate answers. This sure hasn't kept people from trying.  Los Alamos National Laboratory maverick Stirling Colgate (1925-2013; "I like things that explode: nukes, supernovae, and orgasms.") was one of the first to simulate a supernova in the 1960s (think: punch cards).  Over the following decades, computers got much faster, and simulations got better and better.

Today's supercomputers are more than 100 million times more powerful than the computers of the 1960s. Yet we are still struggling with this problem.

We are witnessing a revolution in dimensionality -- We are finally able to simulate stars in three dimensions:

For decades, it was possible to cram in all the complex physics of supernovae only into simulations that were spherically symmetric. Assuming that stars are spherical, the computer codes had to deal only with one spatial dimension, described by the radial coordinate inside the star. Turns out that this is a very bad approximation. Stars are not spherical cows. They rotate, the have regions that are convective (buoyant hot bubbles moving up, cold ones moving down) and turbulent, and they have magnetic fields, which are fundamentally nonspherical.

Two-dimensional (2D) simulations were the next step up from spherical symmetry. They became possible in the early to mid 1990s. In such a simulation, a 2D slice (think of the x-y plane in a graph) of the star is simulated and assumed to be symmetric under rotation about one of the axes. If a star's core is rotating rapidly and has a strong magnetic field, then 2D simulations show that such stars can explode in a "magnetorotational" explosion. Such an explosion is mediated by a combination of rapid rotation and a strong magnetic field that pushes out matter symmetrically along the rotation axis in jets that bore through the stellar envelope in the north and south directions.  This is called a bipolar explosion.  Such an explosion can be very energetic and could possibly explain some extremely energetic supernova explosions that make up about 1% of all supernovae and are referred to as hypernovae.

But not so fast. Nature tends to be 3D. And Nature has it in for us.

The current generation of supercomputers is the first to allow us to simulate supernovae in all three dimensions without the limiting assumption that there is some kind of symmetry. Studying supernovae in 2D was already quite exciting and we thought we'd already understood how things like rotation, magnetic fields, or convection affect a supernova.  So when we set out last year to simulate a magnetorotational supernova in 3D (read about the results in Mösta et al. 2014, ApJ 785, L29), we had an expectation bias: such explosions worked beautifully and robustly in 2D. We expected them to work quite similarly in 3D, perhaps with some slight, hopefully interesting variations about the general theme of a jet-driven explosion.

Colormap of a slice through the meridional plane of a rapidly rotating magnetized stellar core. The left panel shows the axisymmetric (2D) simulation, while the three slices to the right show the full 3D simulation without symmetry constraints at different times. The color coding is specific entropy. High values (indicated by red and yellow regions) of this quantity can be interpreted as "hot", "low-density", and "disordered". 2D and 3D yield fundamentally different results.

We were wrong. Nothing is perfect in nature. Even a quickly spinning star is not perfectly symmetric about its spin axis. There will always be small local variations and if there is some amplifying process around (an "instability"), they can grow. And that is precisely what we found in our 3D simulations. An instability grows from small variations from rotational symmetry. It distorts, twists, crumples, and ultimately destroys the jets that quickly and energetically blow up 2D stars. In 3D, we are left with something fundamentally different that we are only just beginning to understand. It's certainly not the runaway explosion we were looking for.

Volume rendering of the specific entropy distribution in our 3D magnetorotational supernova simulation. Red and yellow regions indicate high entropy ("high disorder", "high temperature", "low density"). The flow structure is fundamentally aspherical. Two large-scale polar lobes have formed that slowly move outward, but are not yet running away in an explosion. This snapshot is from about 180 milliseconds after the proto-neutron star is made.
Here is a YouTube movie from our SXS collaboration YouTube Channel. It shows the dynamical 3D dynamics that drive the supernova towards the state shown in the above picture. The color coding is again specific entropy. Blue and green regions are "cold/ordered/high-density", yellow and red regions are "hot/disordered/low-density". (Viewing tip: Switch to HD and look at the movie full screen!)

Here is another movie, this time showing what is called the plasma β parameter encoding the ratio of gas pressure to the effective pressure exerted by the magnetic field. Small values of β mean that the magnetic field dominates the pressure (and thus drives the dynamics). Regions in which this is the case are color-coded in yellow in the below movie. Dark colors (black/blue) indicate dominance of fluid pressure and in red regions, the magnetic field plays a role, but does not dominate.

And we should have known!

When a spinning magnetized star collapses, it's magnetic field lines get wound up tightly about the star's spin axis, sort of like a tight coil or a phone cord. Plasma physics laboratory experiments have shown long ago that such magnetic fields are unstable to a "kink" instability. If one introduces a small kink that pushes the field lines further apart on one side, those on the other side are compressed, creating a stronger magnetic field there. This increases the force that is excerted on the stellar material, pushing it to the other side. As a result, the small kink is amplified into a bigger kink. In this way, a small microscopic deviation from symmetry will become macroscopic and globally change what is happening in the supernova. This instabilty is fundamentally 3D -- in 2D, there is no "other side," since everything is forced to be rotationally symmetric about the spin axis.

Now that we know what happens in 3D, we feel like we should have known.  We should have anticipated the magnetic field in our star to go kink unstable -- it's really textbook physics. In fact, it's the same problem that plasma physicists struggle with when trying to get thermonuclear fusion to work in Tokamak reactors!

The final word about what ultimately happens with our 3D magnetorotational supernova is not yet spoken. It could be that the explosion takes off eventually, blows up the entire star, leaving behind the central neutron star.  It's also possible that the explosion never gains traction and the stellar envelope falls onto the neutron star, which will then collapse to a black hole. We'll see.  We are pushing our simulations further and are ready for more surprises.

Meet the Team:

Science is a team sport; and this is true in particular for the kind of large-scale, massively-parallel simulations that our group at Caltech is pushing. The full author list of the Mösta et al. 2014, The Astrophysical Journal, 785, L29 (2014) includes Philipp Mösta, Sherwood Richers, Christian Ott, Roland Haas, Tony Piro, Kristen Boydstun, Ernazar Abdikamalov, Christian Reisswig, Erik Schnetter. 

Everybody on this team made important contributions to the paper, but I would like to highlight the roles of the first two people in the author list:

Philipp Mösta
Philipp Mösta is a postdoc in TAPIR at Caltech and part of our Simulating eXtreme Spacetimes program. He is funded primarily by a National Science Foundation Astronomy & Astrophysics research grant (NSF grant no. AST-1212170). Philipp received his PhD in 2012 from the Albert Einstein Institute (the Max Planck Institute for Gravitational Physics) in Potsdam, Germany. Philipp spent the past two years adding magnetic fields to our code and making the entire 3D simulation machinery work for these extremely demanding magnetorotational supernova simulations. His previous training was primarily in numerical relativity, but he picked up on supernova physics with an impressive pace. Philipp carried out all simulations that went into our new paper and he worked closely with grad student Sherwood Richers on analyzing them.

Sherwood Richers
Sherwood Richers is currently a second-year graduate student in physics at Caltech. Sherwood is independently funded by a Department of Energy Computational Science Graduate Fellowship (CSGF) and we are delighted that he has chosen to collaborate with us. Sherwood received his undergraduate degree from the University of Virginia, where he carried out research on magnetohydrodynamics (MHD) with John Hawley. Sherwood is a real expert in all things MHD and he and postdoc Tony Piro are the ones who pointed out that what we are seeing in our 3D magnetorotational supernova simulations is most likely an MHD kink instability. Sherwood also participated in visualizing our simulation output and he is the one responsible for the pretty pictures and movies that we were able to produce from our simulation data!

Saturday, November 16, 2013

Collapsing Stellar Cores, Gravitational Waves, and finally a Gender-balanced Authorlist

Once in a while, I like to blog about new papers coming out of my group. I do this only for papers that are special to me and that have an interesting story behind them. Sure, it's a bit of shameless self-promotion, but, you know, hell, why not? It's also about promoting the great young people I get to work with every day and that's a great thing. I think we all can agree on that. The following post is written for a physics/astrophysics audience, but I have tried to keep it in as simple terms as possible and link to wikipedia articles for background info wherever possible.

So here it is,  our newest paper:

It's so new that we don't even have an number at the time I write this blog post. It'll come out on arXiv on Sunday and we'll submit the paper to Phys. Rev. D in a few days. You can get an exclusive sneak peak now.

This paper is special to me for a couple of reasons.  

First, in case you haven't noticed: There are two women and two men on the author list! 50% women, 50% men. That's how physics should be done. I am very happy about this. This is the first time in my life that I am a co-author on a paper with gender-balanced authorship. Hopefully, this will happen more in the future as we work to fight the gender-imbalance in physics/astrophysics. There is a lot of work that's left to do, unfortunately.

Second, the paper is about rotating stellar collapse to neutron stars and connects back to the very first paper that I wrote back in the summer of 2003 (TEN years ago; I am getting old!). It's great to get back to my roots once in a while. Back then, I was a student at Heidelberg and writing my Diplom thesis about rotating core collapse simulations that I conducted with Adam Burrows at Arizona.
Ernazar Abdikamalov

Third, the paper wasn't even my idea! Ernazar Abdikamalov is a postdoc in TAPIR working with me on a variety of projects in core-collapse supernova theory and radiation transport. Ernazar suggested the project to me and identified it as a great project in which we could involve an undergraduate researcher. Some background on Ernazar: He received his PhD in 2009 from SISSA in Trieste and was a postdoctoral scholar at the Center for Computation and Technology before he joined my group at Caltech.

Alex DeMaio

Alex DeMaio from Rutgers joined us this past summer on a LIGO REU Summer Undergraduate Research Fellowship (SURF). Alex worked with Ernazar on carrying out and analyzing the simulations. She did a great job and contributed tremendously to this paper. Alex is now a junior at Rutgers. She has continued to work with us and helped complete the paper.

Sarah Gossan
Fourth, this is the first paper in which we take a combined approach: We present simulation results for rotating core collapse and, in the same paper, use methods of gravitational wave astronomy to study how our simulation results can be used to make statements about actual observations. The person who made this possible is TAPIR graduate student Sarah Gossan. Sarah really should share the first-author spot with Ernazar.  Her contributions to the paper are just flat out amazing. Sarah is in her second year here at Caltech (she came to us from Cardiff), but is already an accomplished junior scientist with several publications in gravitational wave astronomy.  Sarah and I are both members of the LIGO Scientific Collaboration.

What is this paper about? Some Background.

Here is the thing: We don't know how massive stars spin deep inside. Optical Astronomers can measure their surface velocities and it's possible to use asteroseismology to probe the structure and spin of less massive (smaller) stars, but the internal structure and rotation of massive stars that explode in core-collapse supernovae remains a mystery. This is a pity, since rotation may have a big effect on the way the supernova explosion develops once the core has collapsed to a hot newborn proto-neutron star.

When the next core-collapse supernova happens in the Milky Way (this happens every 50-100 years), we may have a way to probe how the cores of massive stars rotate:  Gravitational waves. Gravitational waves are, quite literally, ripples in spacetime. LIGO is hunting for them, and they are emitted by fast large scale asymmetric mass-energy motions that change on short timescales. This is precisely what happens in a rotating massive star: The centrifugal force, which is strongest at the equator deforms the core into an oblate shape, making an oblate proto-neutron star:

This picture shows the proto-neutron star 12 milliseconds after its birth. It's quite aspherical! The contour lines show isodensity contours and the labels are in units of 10x g/ccm. The colors show specific entropy, and yellowish/red regions have been shocked when the neutron star was made. The core of the proto-neutron star is very dense, wasn't shocked and has very low entropy.

At the end of core collapse, the inner core of a massive star moves with about a third of the speed of light. When the core reaches densities comparable to those in an atomic nucleus, the nuclear force suddenly increases the pressure support in the core. This decelerates the inner core within only about a millisecond! It's useful to think of the inner core as a basket ball that falls to the ground (collapse) and elastically bounces back. "Core Bounce" is what supernova theorists call the moment when the inner core reaches nuclear density, gets decelerated, and rebounds into the outer stellar core.

If the core is not spherically symmetric (rotation!), core bounce gives off a burst of gravitational waves:

The blue line is the time of core bounce and the gravitational wave signal is rescaled by distance D.

What did we do?

It has been known for a while (see, e.g., my review from 2009) that the shape of the gravitational wave signal depends quite sensitively on the spin of the inner core. Various authors in the past hinted at the possibility of reading off the core's rotation from the gravitational wave signal, but nobody had demonstrated that it's really possible. In this paper, we quantitatively demonstrated that it is possible to extract information about the core's spin from a gravitational wave signal observed by LIGO.

Ernazar and Alex performed more than a hundred simulations of rotating core collapse with an axisymmetric general-relativistic hydrodynamics code. From these simulations, they generated a gravitational waveform catalog. Sarah turned this catalog into what gravitational-wave astronomers call a numerical template bank for matched filtering. Matched filtering is a method for extracting a signal from noisy data by cross-correlating the data stream (that is, LIGO output in our case) with a known signal model (template). This method is normally used to look for gravitational waves from coalescing binaries of black holes and neutron stars. It was Sarah's idea to try it for rotating massive star collapse and it turns out to work very well!  This is because the signal shape is very regular and determined by relatively few parameters. Sarah injected trial signals that weren't part of the template bank into simulated random LIGO noise and used her matched filtering pipeline to find the best matching template. Sarah then recorded the rotation parameters (e.g., total angular momentum, degree of differential rotation) of the best matching template as the "measured" quantities.

Here is an example result showing how well we can extract the total angular momentum of the inner core:

The x-axis has the true angular momentum of the inner core that corresponds to the injected trial wave signal. The top panel compares this to the measured angular momentum and the bottom panel gives the relative measurement error. Don't worry about the various different types of symbols plotted here (check out the paper if you are interested). The main point here is that we can extract the total angular momentum of the inner core within ~20%! The blue-triangle outlier at very low angular momentum comes from the fact that in slowly spinning cores the wave signal is not dominated by rotation, but by convective motions after bounce. These are stochastic in nature and cannot be predicted precisely.

Detecting the precollapse degree of differential rotation is harder. First of all, it's more difficult to parameterize, since there are many ways in which angular momentum could be distributed inside a massive star's core. We assumed that the core is rotating before collapse with constant angular velocity on cylindrical shells, because the Poincare-Wavre theorem suggests that the electron-degenerate cores of massive stars should rotate in this way (but it is not clear that this is true!). This is a standard assumption in rotating core collapse and we used a rotation law that reduces the angular velocity as a function of cylindrical radius and a single differential-rotation parameter A. We showed that our method can "measure" the correct value of A, provided the core is very rapidly spinning. This is because differential rotation does not affect the wave signal unless the core is fast spinning.

What are the limitations?


Every work has limitations and it's important to mention them:
The biggest potential systematic problems outside of what we control are uncertainties in the equation of state of nuclear matter and the detailed composition of the inner core (the relevant compositional variable is the number of leptons [neutrinos, electrons] per baryon). Variations in these can lead to waveform changes that can spoil the accuracy with which we can extract information about core rotation.  Our simulations also do not include magnetic fields, but others (e.g., Obergaulinger et al. 2006) have shown that the effects of realistic magnetic fields are unlikely to alter the dynamics and gravitational wave signal in the phases that we study.

On the gravitational-wave astronomy side, we considered only a single LIGO detector and assumed optimal source-detector orientation. In reality there will be multiple detectors taking data in coincidence, but a real signal will also be weaker.