Viscous Gubser BSQ

This repository provides the semi-analytic solution for Gubser Flow with BSQ conserved charges. Currently, there are two equations of state implemented, which can be found in the paper.

The code is organized into fairly stand-alone units (this is because many implementation details are very equation of state specific). Within the plots folder, you will find scripts which generate the figures in the paper; these scripts are not well documented but can serve as examples on how to run the code (details on how to run plotting codes is described below). Also within the plots folder is another folder named for_ccake which contains scripts for generating initial conditions for CCAKE (this is generate_initial_condition.py).

Generate initial conditions

The initial conditions are generated by the generate_initial_conditions.py script. It and the other script are controlled by the run.cfg file. An example of what its contents look like is

tau_0 1.0     # fm/c
tau_f 1.6     # fm/c
tau_step 0.10 # fm/c
x_min -5.0    # fm
x_max 5.0     # fm
x_step 0.025  # fm
y_min -5.0    # fm
y_max 5.0     # fm
y_step 0.025  # fm
temp_0 0.250  # GeV
muB_0 0.050    # GeV
muS_0 0.000    # GeV
muQ_0 0.000    # GeV
ceos_temp_0 1.000  # GeV
ceos_muB_0 1.00    # GeV
ceos_muS_0 1.00    # GeV
ceos_muQ_0 1.00    # GeV
pi_0 0.0      # ---
tolerance 1e-12
output_dir ./data/analytic_solutions

Hence, you can control all the variables for the output, including the directory where the output is being dump to. To run the generation scripts, use

python3 generate_intial_conditions.py
# or
python3 generate_sem-exact_solutions.py

Running plotting code

The files compare_mus_qgp.py, compare_mus_ccake.py and freezeout_surface_qgp.py all can all be run by just typing

python3 <file_name>

compare_mus_qgp.py generates Figures 1 and 2 in the paper, while freezeout_surface_qgp.py generates Figure 3. Figure 5 is generated by the script plot_gubser_figure_4.py, which can be found in the plots/for_ccake directory. Its usage is

python3 plot_gubser_figure_4.py <path_to_analytic_soln> <path_to_numeric_soln>

Lastly, Figure 6 is generated by freezeout_surface_ccake.py, and requires the freezeout surface outputted by CCAKE. Its usage is

python3 freezeout_surface_ccake.py <path_to_CCAKE_freezeout_file>

Running your own code

To write your own script to run the code, follow the following example

from scipy.integrate import odeint
from scipy.interpolate import interp1d
from numpy import concatenate
from numpy import linspace
from numpy import array
from system_massless_qgp import MasslessQGP
from variable_conversions import milne_T
from variable_conversions import milne_mu

# Create an instance of the system you want to solve: MasslessQGP or ConformalPlasma
system = MasslessQGP()

# de Sitter time steps
rhos_backward = linspace(-10, 0, 1000)[::-1]
rhos_forward = linspace(0, 10, 1000)

# Initial conditions: (temperature, chemcical potential, shear pressure)
y0s = array([1.0, 1.0, 0.0])

# Obtain the numerical solutions in de Sitter space (that's how the equations are implemented)
# Here, I assume the initial condition is given at rho = 0
soln_1 = odeint(system.eom.for_scipy, y0s, rhos_forward)    # Forward in de Sitter time 
soln_2 = odeint(system.eom.for_scipy, y0s, rhos_backword)   # Backward in de Sitter time

# Glue forward and backward evolutions together
t_hat = concatenate((soln_1[:, 0][::-1], soln_2[:, 0]))       
mu_hat = concatenate((soln_1[:, 1][::-1], soln_2[:, 1]))
pi_bar_hat = concatenate((soln_1[:, 2][::-1], soln_2[:, 2]))
rhos = concatenate((rhos_1[::-1], rhos_2))

# Create interpolating functions from the solutions
t_interp = interp1d(rhos, t_hat)
mu_interp = interp1d(rhos, mu_hat)
pi_interp = interp1d(rhos, pi_bar_hat)

# Convert temperature and chemcial potential from de Sitter values to Milne values
t_evol = milne_T(tau, xs, 1, t_interp)
mu_evol = milne_mu(tau, xs, 1, mu_interp)

# Milne coordinates, for example, for plotting
tau = 1.0    # fm/c
xs = linspace(-6, 6, 200)

# Convert thermodynamic quantities from de Sitter values to Milne values
e_evol = system.milne_energy(tau, xs, 0.0, 1.0, t_interp, mu_interp)
n_evol = system.milne_number(tau, xs, 0.0, 1.0, t_interp, mu_interp)
s_evol = system.milne_entropy(tau, xs, 0.0, 1.0, t_interp, mu_interp)
pimunu_evol = system.milne_pi(
    tau, xs, 0.0, 1, t_interp, mu_interp, pi_interp, nonzero_xy=True,)

Adding an Equation of State

To add your own equation of state, simply make a copy of system_massless_qgp.py or system_conformal_plasma.py and fill out all the functions with your respective analytic expressions.