Gravitational waves, neutron stars, and numerical simulations

Ian Hawke

github.com/IanHawke
STAG, University of Southampton

ianhawke.github.io/slides/seds2020

GW170817

Two neutron stars merging seen

in gravitational waves;
in every electromagnetic band;
well localized.

Amazing insight into extreme astrophysics.

Sometimes a breakthrough needs some luck. Three years ago, after thousands of physicists had worked for decades, we got lucky to observe two neutron stars merging 40 Mpc (130 Mly) away. In astrophysics terms this is "close". Even better, some of the resulting signal pointed in just the right direction. But first, let's go through this plot.
Top left we've got the signal in gravitational waves. Usually described as tiny ripples in spacetime, they're generated by very massive objects moving very quickly. When two neutron stars, each more massive than our sun, are orbiting each other at about a third the speed of light, plenty of gravitational waves are produced. The energy that goes into these waves sucks energy from the orbit of the neutron stars, reducing the orbital separation, and hence speeding up the orbital velocity. This is seen as an increase in the frequency of the gravitational waves with time, which goes up until they crash into each other.
The gravitational waves propagate at the speed of light and are generated even before merger, so are the first signal to reach us. However, they're hard to analyse and detect. After the neutron stars smash together they may, as in this case, collapse to a black hole, leaving a hot soup of matter around. When this remaining matter falls in to the black hole it can produce a massive jet of energized particles. This kilonova is one form of gamma-ray burst. For GW170817, this flash of high energy gamma rays was actually the first indication that observers noticed of this violent event. Gravitational wave observers were able to rapidly dig into their data to find it, but the automated algorithms didn't spot it.
As the remaining matter slowly cools, and the jet of energized particles moves towards us, the event has been seen in every single EM band. As EM observatories have much narrower fields of view this means the event has been pinned down to one specific galaxy.
This one observation had a huge impact on the field. Leaving aside the Nobel prize, and the newspaper headlines about "generating all the gold in the universe", this has led to thousands of papers talking about the impact of this one observation on theories of gravity, theories of electromagnetism, theories of matter in extremes of temperature and density, and theories of how these objects formed in the first place. So, we should talk a little more about the objects that we're studying and why.

Andersson

Most of the information that we have about neutron stars comes from observing pulsars. Observations of neutron stars typically put them around 1.4 times the mass of our sun. The most massive neutron stars are around twice the mass of our sun. So we tend to think of them like our sun, but compressed.
However, neutron stars are really compact. The most accurate results come from the NICER satellite which observed burning on the surface of neutron stars. We can also get constraints from pulsars, and from our gravitational wave observation (singular!). These all put the neutron star radius at between 10 and 15 kilometres: I believe that's roughly the distance between your main campus and the Aguada Fort. Crushing the $10^{60}$ particles contained within the sun into such a small volume leads to extremes of density, where matter can behave in ways not seen elsewhere in the universe.
A neutron star may look a little like a "planet". The outer layer has a "crust", just like the Earth. However, this crust is not made of rock. Instead, neutrons have been pushed so close together that they form a lattice, like a metal or a crystal. It's not the weak EM force that pins them together like terrestrial materials, however, but the strong nuclear force. This is probably the toughest "solid" material in the known universe.
Moving further in, the density increases to the point where fundamental particles separate and flow freely. In fact, they may interact quantum mechanically and flow as a superfluid, without any resistance from their neighbours. This is mostly a theoretical prediction, but there are observations of temperature changes in some young neutron stars that is best explained by the superfluid.
Further in again, we reach densities where our best theoretical models have significant uncertainties. It may be that the quarks making up the fundamental particles become deconfined, or other exotic particles are formed.
Finally, it's important to note that the charged fundamental particles like electrons and protons are still around: they do not vanish, even when the forces are dominated by the heavy neutrons. These charged particles will generate strong magnetic fields, as the neutron stars and so compact that they can rotate quickly, and the compression of the magnetic energy raises the peak field strength. A typical neutron star is eight orders of magnitude greater than our sun. This can go up rapidly during a merger, especially, possibly by another eight orders of magnitude!

What do we learn?

Colliders: hot, low density.
NSs: cold (!), high density.

So NSs tell us about

extremes of gravity;
extremes of particle physics.

Results should complement colliders.

Watts et al 2016.

There are two main things we expect to learn from neutron stars.
First, as we've said, neutron stars are really compact. That makes them very extreme objects for gravity. In standard theories of matter and gravity, the only object more compact than a neutron star is a black hole. But black holes are harder to observe: we can now see the gravitational waves as they merge, and we can see hot matter near the horizon, but we don't get any observations in the electromagnetic spectrum through merger itself. With neutron stars we do get EM signals near this extreme event.
As an example, in GW170817, the gamma ray burst first reached us about two seconds after the peak in gravitational waves. Over hundreds of millions of light years this allows us to compare the speed of gravitational waves to the speed of light to incredible precision. It's the tightest constraint to date on how fast information about spacetime propagates through spacetime, and a crucial test of Einstein's theory (which, of course, it passed).
Second, we get an insight into fundamental particles in an entirely novel way. When talking about particle physics experiments we mostly think of the big particle colliders like the LHC, which accelerate particles to near light speed and bash them together. This probes the nature of the fundamental particles, but at low densities (as they can't accelerate much matter) and high temperatures (from the extreme acceleration).
Neutron stars, in contrast, are high densities and low temperatures. This gives us complementary information about the matter, and by piecing together observations we narrow down the acceptable physics.

Summit supercomputer. Photo: Carlos Jones/ORNL

Sun: $\sim 10^{60}$ elementary particles.
Computer resource on Earth: $\sim 10^{25}$ FLOPs.

Averaging

If we can't treat each particle separately, we'll have to average over them in some way. The hydrodynamics approximation averages over them all. Our $10^{60}$ separate particles, each with their own positions, velocities, and properties, are averaged to give a single (vector) field. At each point in spacetime this field will describe the number density of particles, and the direction in which they flow.
This averaging process gives us a small set of differential equations to solve. With some standard approximations, it actually reduces to five PDEs (on top of those describing the evolution of the spacetime). This makes the numerical simulation practical. However, it has also thrown away (or compressed) a lot of physics.
It also means that there is a crucial difference of principle and of scale between vacuum relativity and relativistic hydrodynamics. In the vacuum case we have Einstein's field equations which, in principle, hold on all scales within our simulations and models. In hydrodynamics, however, our model and description breaks down on certain scales. We must keep in mind that as computer power and numerical accuracy increases our results may not get better, as we may hit the limitations of our models.

Characteristics

Nothing moves faster than light;
Cause and effect: only look at the "past light cone".

If a particle follows a path $X(t)$, and $q(X)$ is constant on that path, then $$ \partial_t q(x, t) + \frac{\partial_x X}{\partial_t X} \partial_x q(x, t) = 0. $$

The path $X(t)$ is a characteristic: we trace information back in time to get future solutions.

Having averaged our matter to get a set of differential equations, we need to think about how the solutions of these differential equations behave.
As we're working within Einstein's theory, we know that nothing can propagate faster than light. That means that the only things that can change a quantity "now" must have been able to interact with that quantity. Think of a particle $P_1$ at the current time, and ask if it can interact with a different particle $P_2$. If they can interact then it must be possible to connect, using some curve moving at less than lightspeed, to the spacetime point where the particle $P_1$ is now to some spacetime point where $P_2$ was in the past.
We then assume that interaction carries information: that is, some quantity $q$ that does not change along that path in spacetime $X$ between the two particles. It is the ordinary derivative of the quantity $q$ that does not change. If we sit at a coordinate location and let the path move past us then we have to take partial derivatives to see how the quantity changes in our local coordinates $(x, t)$. This gives us a PDE satisfied by the quantity $q$.
The equations of motion of relativistic matter will have this sort of form: information will propagate at large, but finite, speeds. We want to work backwards from the PDEs to get characteristics along which information does not change. If we knew that, then we might be able to predict future information from our initial conditions.

Conservation

EFEs $G_{ab} = 8 \pi \kappa T_{ab}$ imply $\nabla_a T^{ab} = 0$.

Pick a tetrad, $e_b^{(j)}$ to get $$ \begin{aligned} && \nabla_a \left[ e_b^{(j)} T^{ab} \right] &= \tfrac{1}{\sqrt{-g}} \partial_a \left( \sqrt{-g} e_b^{(j)} T^{ab} \right) = -T^{ab} \nabla_a e_b^{(j)} \\ \implies && \color{red}{\partial_t {\bf q} + \partial_i {\bf f}^{(i)}({\bf q})} &= \color{red}{{\bf s}}. \end{aligned} $$ Balance law form.

Only four equations: need other constituitive equations for, eg, EM, particle number, etc.

We now know the outlines of our problem. We want a model of matter, coupled to GR, that we can write as a set of PDEs. This model will encode the detailed physics we're interested in, at the level we can simulate with our given resources. But before we try and be specific, let's look at what we can say in general.
We know the field equations require total stress energy conservation: the divergence of $T^{ab}$ is zero. We can use this to get four equations of motion. The most useful form for our purposes is to use the standard identity for the 4-divergence of a vector, which can be written as the partial 4-divergence of that vector, weighted by the metric determinant. To get from the stress-energy, which is a two-tensor, to a vector, we contract the free index with one member of a tetrad, which is a set of orthonormal (in our case) vectors that we can link to the coordinates. This then allows us to write stress energy conservation as four PDEs, one for each member of the tetrad.
The key point is that the only derivatives of the matter appear in total divergence form. The term on the right-hand-side which isn't in this form only involves derivatives of the tetrad, which can be written as derivatives of spacetime quantities. We write this abstractly in balance law form: all derivatives of the matter are total derivatives, but these fluxes must balance off against the source terms which come from the geometry.
If we simplify to Minkowski space in standard Cartesian coordinates then the source terms vanish and these balance laws reduce to the conservation laws for energy-momentum. However, we need to be careful of this interpretation (locally): if we write Minkowski spacetime in spherical coordinates then there will be source terms, but that certainly does not mean that energy-momentum is being transferred between the matter and the spacetime.

Shock formation

Advection equation $$ \partial_t q + \partial_x (v q) = 0. $$ Information moves right, speed $v$.

Burgers equation $$ \partial_t q + \tfrac{1}{2} \partial_x q^{2} = 0. $$ Information moves right, speed $q$. Shocks form.

We see that generically the matter will obey conservation laws. We expect the fluxes ${\bf f}$ to depend nonlinearly on the conserved variables ${\bf q}$. This immediately tells us something very important.
The simplest conservation law is the advection equation, where the solution is propagated to the right with constant speed $v$. This is exactly the situation discussed when talking about characteristics: the information is constant along straight lines of slope $v$.
For a nonlinear conservation law we can locally approximate it as an advection-like equation, where the solution propagates with speed $\partial_q f$. This means the propagation speed depends on the data. Depending on the form of the flux function, some gradients will steadily get steeper, and some will get shallower. In the absence of dissipation (and crucially the generic total energy-momentum equations have no dissipative terms) these steepening effects will lead to discontinuities: a generic shock will form.

PDEs and numerics

Finite differences:

Store point values $q_i$;
$$ \partial_x f(q) \to \tfrac{1}{2 \Delta x} \left( f_{i+1} - f_{i-1} \right). $$
Works well for smooth data.
Ok convergence, communication issues.

Finite volumes:

Store cell averages $\hat{q}_i = \tfrac{1}{\Delta x} \int_{x_{i-1/2}}^{x_{i+1/2}} \text{d} x \, q$;
$$ \partial_x f \underset{\int \text{d} x}{\rightarrow} \tfrac{1}{\Delta x} \left( f_{i+1/2} - f_{i-1/2} \right). $$
Works well for discontinuous data.
Slow convergence, communication issues.

Finite elements:

Store modal coefficients $q^{(m)}_i$;
$$ \partial_x f \underset{\int \text{d} x \, \phi(x)}{\rightarrow} M {\bf q} + \dots . $$
Works well for smooth data; tricky for discontinuous.
Good convergence, good communication.

On to numerical approximation. We want to solve continuum partial differential equations on a computer. The solution, $q(x)$, needs - in principle - an infinite amount of information. That's because the continuum solution could, in principle, take totally different unpredictable values at every separate point $x$. This way of thinking is clearly useless when working with a computer with a finite amount of working memory.
Instead we consider ways of approximating the solution $q(x)$. We do this by making assumptions about the functions we're trying to solve for. The more correct assumptions we make, the better our solution will be. However, as soon as we make one inconsistent assumption our solution may be completely wrong. Think of this as a "garbage in, garbage out" type of saying. The solutions for the spacetime that we evolve will be nice and smooth nearly everywhere, so we can make very strong assumptions and get very accurate solutions. However, our solutions for relativistic fluids can have strong shocks, so we need to make much more conservative assumptions, making it harder to get accurate solutions.
The first is the point value approximation, used in finite differencing. In this case the domain is discretized into points, and we assume that we know (or want to compute) the value of the solution at each of these points. The solution is not known between the points. When we want to link the solution at different points - for example, when approximating a derivative - we have to impose the general behaviour of the solution between those points. For example, if we interpolate between points using a straight line then the derivative is the slope of this line.
The second approach is to split the domain into cells, and within each cell to store the average value. This is the finite volume approach. In this case some information is known at every point in spacetime: the average value near that point. However, the exact value is known nowhere. Again, to get the value of the solution at a specific point we have to impose the general behaviour of the solution. When restricted to a given cell this behaviour has to be consistent with the average value within that cell. The finite volume approach implicitly "smears out" the local behaviour of the solution, which is exactly what we want when dealing with discontinuous solutions like shocks.
The final approach has links to the finite volume approach, but is more often referred to as a finite element method. Instead of storing the average value of the solution within a cell, the moments of the solution with respect to some function basis are stored. Typically the zeroth moment would be the cell average of the solution, the first moment linked to its derivative, the second to its curvature and so on. This uses a function basis expansion: think about a Taylor series, or a Fourier series, but truncated after a few terms. This is fundamental to the Discontinuous Galerkin method that will be one of the central methods for the future of our field.

Messy mergers

The merger process

mixes the particles strongly;
generate turbulence and small scale structures;
the high temperatures and mixing lead to strong signals, nucleosynthesis, etc.

Radice et al, 1809.11161

This simulation by Radice illustrates how messy the merger process can be. The figure is illustrating the lepton fraction: that is, how many particles in a local region are electrons. This is very important for computing local properties, as the electrons are so much lighter than the other particles.
As the stars merge we see the matter is reasonably well ordered. However, after merger we see small scale features everywhere. The stars do not so much crash straight into each other as graze past each other. This generates lots of turbulence, leading to small scale features all over the place. This mixes up the particles in ways that is hard to predict. A lot of current working is going into understanding how we can model and analyse these small scale features, and what impact they will have on observed signals.

Structured jets

After the merger

jets form, accelerated by magnetic fields;
but the physics needed is still unclear, and the strength very model dependent;
the magnetic field structure is crucial in creating a kilonova signal.

Ciolfi, 2001.10241; Wright & Hawke, 1906.03150.

Summary

We have discussed

neutrons stars as astrophysical laboratories;
building a fluid model;
numerical evolutions of neutron star models.

There's a lot more to do!