Uncovering the quark-gluon plasma

Holographic Evolution with
Dynamical Boundary Gravity

This is an update from the code below that can handle dynamical boundary gravity. Instead of treating the boundary metric as a source the metric follows the Friedman equations. This causes extra complications, as explained in 2109.10411. The code is then able to evolve strongly coupled non-conformal matter in universe with a positive, negative or vanishing cosmological constant. For flat universes this leads at late times to de Sitter, flat or a big crunch geometry.

The code can be found here: dropbox. The code is developed together with Christian Ecker.

What attracts to attractors?

This is an update from the Tehran lectures below with a more refined code. One can look at this folder dropbox, which is the original code used in 1907.08101.

Lectures on numerics in dynamical AdS spacetimes

This page is the result of a week of lectures at the IPM Institute for Research in Fundamental Sciences, in Tehran, Iran (June 2015). I am thankful for the warm hospitality and active participation of the participants during this short course.

The main aim of this page is to present a simple but complete example of an evolving AdS spacetime. This requires several steps, whereby the main ideas are outlined in the short presentation included. These ideas are then worked out in the notebook BIthermalization.nb, which is capable of evolving homogeneous or boost-invariant AdS spacetimes with anisotropy (the only difference in this case is the metric on the boundary).

A few key ingredients are:

  • The notebook is also a package, whereby loading the package is equivalent to evaluating all initialisation cells (conform explanation in the presentation). Before loading the package it is possible to put mode = homogeneous, or mode = BI for the homogeneous or boost-invariant version of the package (BI is standard). This is illustrated in runsimulation.nb.

  • The main command is run, which takes profileB (conform presentation), the initial a4, initial time, end time, grid points and session name as input.

  • The notebook aims to compromise on readability and showing some efficient Mathematica functions. These Mathematica lines sometimes take a bit of effort to write, and hence also to understand their meaning. For this I provide some comments, containing the idea of the code, but usually it is essential to play around with the small parts of the code ('Mathematica surgery').

Also, the notebook can definitely be improved to contain more physics, some of the following can be interesting exercises to implement, in increasing order of difficulty:

  • track the location of the apparent horizon (Sdot = 0)

  • not only save a4, b4 and ksi, but also save the metric functions every few time steps and create an interpolating function of the latter

  • use this function to compute the event horizon and other outgoing null geodesics

  • use the gauge freedom ksi(t) to fix the apparent horizon at r = 1, if it exists (if it doesn't initially this exercise is a bit challenging)

  • optimise the code using Mathematica's Compile routine (such as implemented in i.e. shockwaves notebook on this website)

  • compute other observables, such as entaglement entropy (conform 1506.02658)

  • include other field, such as scalar or Maxwell field (in probe approximation is simpler)

Lastly, as always, comments, feedback and questions are much appreciated!

Tensor networks and MERA package

This page contains the package required for the 7th Mathematica Summer School in Theoretical Physics (MSSTP) and includes a sample MERA for the Ising and Heisenberg model with bond dimension 6 (closely following 0707.1454). Make sure to first extract all files, including the data sub-directory. Most useful will be the ncon function (as in 1402.0939), which provides support for tensor contractions.

Also included are the exercises by Guifre Vidal in the same school, for which the ncon function is essential to construct the tree tensor network.

Update: the ncon function can now also be evaluated using nconp, which then also prints the cost of the contraction, plus a plot of the corresponding network:

Holographic heavy ion collisions at finite coupling

This page accompanies 1610.08976, which studied gravitational shock wave collisions in Einstein-Gauss-Bonnet theory. These collisions model heavy ion collisions at finitely strong coupling, by including inverse coupling constant corrections to the infinite coupling result.

The code can be downloaded here: dropbox. In principle everything is self-contained. Runwide and runnarrow run the wide and narrow simulations (take care to evaluate initialization cells, simulations take about 3 and 14 days respectively). Both require two simulations with different regulators to thereafter allow extrapolation to zero regulator.

The resulting simulation files are then loaded into plotsforpaperfinal. There you can either evaluate all cells (without initialization cells), or first evaluate initialization cells and then make a particular plot. The simulation files are rather large, but please contact me for a link to those files, to avoid rerunning the whole simulations.

Lastly, all equations are included, but the full derivation of the EOM and holographic renormalization has not been updated at the Mathematica level yet.

PS: unfortunately the code does not work with Mathematica 11, due to a bug in that version. (the bug is in Series[x^4, {x, 0, 2}] not expanding to the power asked for; it's not that difficult to work around, but for now I decided to wait till a new version may fix this bug)

Radial Flow project

Radial flow project (accompanying arxiv:1211.2218)

On this page you can find the numerical code, the resulting data and a simple movie displaying the transverse velocities for our models of a nucleus and a fluctuation.

First of all a (big) warning: this notebook and results are not (at all!) meant to be a pedagogical introduction how to do computations as outlined in the paper. Instead, they serve a scientific purpose (data and methods should be publicly available) and more importantly it gives an opportunity to have access and view the results in a very direct way.

Notebook (radial_flow.nb)

The notebook needs these files in the same directory. Alternatively, you can compute the formulas in the files, by putting compute = True; at the top of section Full equations, but this will take about 15 minutes.

Initially you can run everything till section "Do runs". After that you can even do the runs, although they will both take about 30 hours (on a new desktop). Alternatively you can downloadresultnucleus.dat and resultfluc.dat.

Having obtained the data files you can evaluate the Results section, which will reproduce the plots displayed in the article. If you have a look at the section Plots you can find more plot functions (although only the plots appearing in the article are checked on errors). Also, naturally, it should be rather easy to access all relevant data just looking at the formulas plotted.

Comments, questions or criticisms are much appreciated! (e-mail)

PS: at some point I may write a more pedagogical description for these kind of numerics, if it turns out there is demand for this. A helpful start might be an easy introduction to spectral methods, combined with the presentation at Iberian strings on how to solve Einstein equations this way.


Here we show movies of the nucleus and fluctuation model. They show the energy density and the radial velocity in the transverse plane:



PhD Thesis

Download thesis here (book format) or at arXiv (article format)

Cover figure

The cover figure is supposed to represent the story behind this thesis. The upper part is a cartoon of a collision of two nuclei (typically lead or gold). Part of those hit head-on (opague) whereas the rest just flies on. Between the parts of the nuclei that collided a quark-gluon plasma will form, illustrated by (approximate) isothermal contours. One interesting aspect is the somewhat elliptical shape of this region, which induces a (measurable) difference in pressure gradient horizontally and vertically, illustrated by the arrows. This thesis tries to describe this collision using a so-called holographic technique. This technique suggests that such a very quantum mechanical collision may behave similarly to a collision of classical gravitational shocks in a space with one extra dimension, which is illustrated by the direction going down. At the bottom of this extra dimension a black hole horizon will form, with the same temperature as the quark-gluon plasma. This thesis studies the formation of this horizon, and hence can hopefully give some extra understanding of the formation of the quark-gluon plasma.

Computations and movies

PS: unfortunately the code does not work with Mathematica 11, due to a bug in that version. (the bug is in Series[x^4, {x, 0, 2}] not expanding to the power asked for; it's not that difficult to work around, but for now I decided to wait till a new version may fix this bug)

This site contains Mathematica notebooks and resulting simulations for many of the computations presented in my PhD thesis. All computational files are contained in this archive (28MB), while two sample shock wave simulations can be downloaded here (60MB). All files were created with Mathematica 9.0.1; unfortunately they don't seem to be directly compatible with Mathematica 10. It includes the following notebooks:


This is the main (commented) file, with up-to-date efficient computations, deriving all necessary formulae and numerical code. Mathematica groups all initialisation cells into a package shockwaves.m, which can be used both for running, loading and analysing shock wave simulations.


This is a sample file how to run a shock wave simulation, in this case the wide shocks of section 3.2. The basic input parameters are the profile function h(z), the regulator energy densities and the grid points, after which run[] produces the stress-energy tensor. There are various options, which can be found in the "Standard values" section of shockwaves.nb and in Options[run], but usually little has to be changed.


This is a sample file loads the stress-energy tensors produced in the previous sample notebook and produces most plots presented in chapter 3. It also heavily relies on shockwaves.nb, and for a complete overview of all possibilities we advice to scroll through subsection "Make plots" of section "Load and display".


This sample file runs the standard shocks with six different set of gridpoints and plots resulting constraint violations, as presented in subsection 3.2.2.


This small package is used by shockwaves.nb and summarises numerical functions, such as generating the grid, and making differentiation and filter matrices. It has beta-support for adaptive gridding, but this is not functional in the shockwaves package yet.


A pedagogical presentation on spectral methods is included, following the two sample linear and non-linear equations of section 2.2.


This is an old notebook and therefore has a much less efficient numerical interpretation; it is not commented either. We include it for completeness; when evaluating the complete notebook one should be able to solve the homogeneous system described in section 2.2, producing part of the plots presented there.

For gravitational computations in chapter 4 we refer to the radial flow project, which unfortunately is also not fully commented.

Movies in thesis

The thesis contains various movies, which can be viewed by clicking on the links below:

Narrow and wide shock collision in AdS (figure 3.6):

Off-central collision using shock waves: https://www.dropbox.com/s/q0yeu1b7itsvtvd/offcenter.wmv?dl=0