Holographic EoS Code — C++ Version

0. Directories and files.

0.1. Main components

The source code for the main components of the holographic equation of state, including headers, can be found under ./src/.

The makefile can be found in the root directory, at ./Makefile.

Compiled objects are divided into auxiliary objects, that are built into ./obj, and executables, which are built under ./bin/.

1. How this code works

1.1. Overview

The present code provides a holographic model of the QCD equation of state as a function of temperature $T$ and baryon chemical potential $\mu_B$. This model is calibrated to reproduce lattice-QCD data at $\mu_B =0$ and presents a good agreement with extrapolations to finite $\mu_B$. Because of its holographic nature, it is expected to be more reliable at moderate to high temperatures.

This model predicts a first-order phase transition at large densities, with a critical endpoint at ($T$, $\mu_B$) $\approx$ (103, 597) MeV for the polynomial hyperbolic model.

1.2. Physics model

The present model of the QCD equation of state relies on a bottom-up construction based on the gauge-gravity duality. In short, the equation of state is calculated for a theory living in the 3+1-dimensional boundary of a Einstein–Maxwell–Dilaton (EMD - that is, gravity + scalar field + gauge field) theory in a 5-dimensional asymptotically anti-de Sitter (AdS) spacetime.

In this description, one identifies the properties of a charged black-hole in the bulk theory with the thermodynamics of a strongly-coupled fluid in a 3+1-dimensional slice of the AdS spacetime. In order to extract meaningful quantities from this model, it is necessary to connect the infrared physics on the black-hole horizon at the AdS radius $r=0$ to the ultraviolet scalings at AdS radius $r\to\infty$, where the AdS spacetime geometry must be recovered.

The construction of the model is presented in Refs. [1] and [2] and is inspired by earlier models such as the one in Ref. [3]. Since the publication of [2], two other groups have also proposed similar models capable of reproducing lattice data, and which will also be made available in the present code: [4] and [5].

1.2.1. Parametrizations

These models are specified by two functions:

the potential $V(\phi)$ for the dilaton (scalar) field $\phi$ in the action for the 5-dimensional theory,
the dilaton-gauge coupling $f(\phi)$ that multiplies the gauge-field kinetic term $f(\phi)F^{\mu\nu}F_{\mu\nu}/4$ in the same action,

where it is important that $V(0)=-12$ and $f(0)=1$.

Besides the parameters defining $V(\phi)$ and $f(\phi)$, one must specify the gravitational constant in five dimensions $\kappa_5^2=8\pi G_5$ and the energy scale $\Lambda$, with default values

\[\kappa_5^2=8\pi G_5 = 8\pi(0.46),\; \Lambda = 1058.83\; \text{MeV} .\]

These parameters and functions, as their derivatives, are implemented in the BhModel class.

In its current version, this code supports the following parametrizations for modeling these functions. Other parametric forms can be included by editing BhModel.h and BhModel.cpp.

Polynomial Hyperbolic Model

Implemented in class PolyHyperbolicModel. Reproduces the model of [2] for $c_3=0$,

\[V(\phi) = -12 \cosh(\gamma \phi) + b_2 \phi^2 + b_4 \phi^4 + b_6 \phi^6 ,\]

$f(\phi) =\frac{\operatorname{sech}\,(c_{1}\phi+c_{2}\phi^{2}+c_{3}\phi^{3})}{1+d_{1}} + \frac{d_{1}}{1+d_{1}}\operatorname{sech}\,(d_{2}\phi) .$

MUSES Parametric Model

Implemented in class MusesParametricModel, with the following functional forms:

\[V(\phi) = -12\cosh\left[\left(\frac{\gamma_1\,\Delta\phi_V^2 + \gamma_2 \,\phi^2}{\Delta \phi_V^2 + \phi^2}\right) \phi\right],\]

\[f(\phi) = 1 - (1-A_1) \left[\frac{1}{2} + \frac{1}{2}\tanh\left(\frac{\phi - \phi_1}{\delta \phi_1}\right)\right] - A_1\left[\frac{1}{2} + \frac{1}{2}\tanh\left(\frac{\phi - \phi_2}{\delta \phi_2}\right)\right].\]

1.3. Initial conditions: phi0 and Phi1

To calculate the equation of state, one must solve a set of 4 second-order ordinary differential equations, describing the evolution of the dilaton $\phi(r)$. the gauge field $A^\mu(r)$ and the metric $g_{\mu\nu}(r)$ in the AdS radius coordinate $r$. An ansatz is used to characterize these degrees of freedom in terms of 4 dynamical fields: $\phi(r)$, $\Phi(r)$, $h(r)$, and $A(r)$.

Initial conditions for these fields and their derivatives must be provided at the black hole horizon which is taken as r=0 for the numerical coordinates. By using a convenient gauge, they are reduced to $\phi_0=\phi(0)$ and $\Phi_1=\Phi'(0)$, aka phi0 and Phi1_raw in the code. For each value of phi0, there is a maximum allowed value of Phi1, Phi1Max, depending on the particular model. This maximum value is given by the method BhModel.get_Phi1_max(double phi0) in the BhModel class. For that reason, we often use the normalized variable Phi1_norm = Phi1_raw/Phi1Max instead of Phi1_raw. Both variables are members of the BlackHole class.

Given a pair phi0, Phi1_norm, the function set_ics(phi0, Phi1N) in the BlackHole class initializes the fields and first deriatives to their values at a small radius $r = r_0 > 0$, from where the equations of motion can be integrated. Initial conditions at $r=r_0$ are determined by substituting a Taylor expansion for the different fields into equations (3) to (7) in Ref. [2]. By starting numerical integration at $r_0>0$, one avoids numerical issues that take place exactly at the balck-hole horizon. Function BlackHole::set_ics(phi0, Phi1N) also calculates the conserved black-hole charges, which can be used as a sanity check for the numerical solutions.

It is only after solving the equations of motion that one finds which point $(T,\mu_B)$ correspond to a given (phi0, Phi1) pair. The mapping between (phi0, Phi1) and $(T,\mu_B)$ is not bijective. To each point in the phase-coexistence region of the $(T,\mu_B)$ phase diagram, there correspond three distinct (phi0, Phi1) pairs, corresponding to the stable, metastable and unstable phases. In the equilibrium equation of state, these phases must be matched according to the Maxwell construction. Lines of constant phi0 and varying Phi1 in the $(T,\mu_B)$ phase diagram are shown in the figure below:

Lines of constant phi0 in the T, mu phase diagram

The files Tmu_utils.h and Tmu_utils.cpp provide root finders to solve for $:raw-latex:phi_0 $ and $\Phi1$ as functions of $T$ and $\mu_B$.

1.4. Equations of motion

To calculate the equation of state, one must solve equations (3) to (6) in Ref. [2]. In practice, we reduce these 4 second-order differential equations to 8 first-order differential equations by treating the derivatives of each field $\phi'(r)$, $\Phi'(r)$, $h'(r)$ and $A'(r)$ as independent variables , which follow their own equations of motion. These equations are solved by the function solve(), in the BlackHole Class.

Numerical integration is performed by using an adaptive stepper. The stepper is initialized using the Runge-Kutta Dormand-Prince 5 method in Lib Boost ODEint, where it starts at a small radius $r = r_0 > 0$ and proceeds up to a radius INITIAL_RADIUS. For each loop, the step is adaptively adjusted, and the function BlackHole::convergence_condition() is called to test for numerical convergence at each step. The integration continues unless the MAX_RADIUS is reached or the integration failed the numerical convergence test for any step. INITIAL_RADIUS and MAX_RADIUS are defined in blackhole_engineering.h.

An extra ordinary differential equation of first order remains, which provides a constraint for the system, imposed when calculating initial conditions for the different fields in BlackHole::SetICs(phi0,Phi1N). This constraint is presented in Eq. (7) in Ref. [2] and is also used as a sanity check for the solutions in function BlackHole::check_constraint(), called by function BlackHole::do_physics().

1.5. Asymptotic expansion and thermodynamics

After numerical integration, function solve() calls function do_physics in the BlackHole Class to calculate the temperature $T$, the chemical potential $\mu_B$, the entropy density $s$ and the baryon density $\rho$. These are calculated from equations (11) to (14) in Ref. [2]. These expressions depend on the field values near the horizon, the prefactor $\phi_A$ in the scaling of $\phi(r)\sim \phi_A e^{-\nu A(r)}$, and the exponent exponent $\nu$, calculated by function nu() in Class BhModel.

1.5.1. Relaxational method

The most challenging part in estimating the different thermodynamical quantities is finding the prefactor $\phi_A$ in the scaling of $\phi(r)$. Because $A(r)$ grows linearly in the scaling region, a direct evaluation of $\phi_A\approx e^{\nu A(r)}\phi(r)$ is numerically unstable: numerical imprecisions in the exponentially small $\phi(r)$ are amplified by the exponentially large $e^{\nu A(r)}$ factor.

Previously in the literature, this problem has been dealt with by fitting the behavior of $\phi(r)$ with a decreasing exponential to extract $\phi_A$, in a carefully chosen window in $r$. However, this procedure often lead to noisy features in the extracted $\phi_A$ [2].

Here, we tackle this issue by implementing one extra differential equation in Class BlackHole::solve():

\[\frac{dC(r)}{dr} = - \Gamma (C(r) - e^{\nu A(r)}\phi(r)),\]

such that $C(r) \to \phi_A$ as $r\to \infty$. Its solution:

\[C(r) = C(r_H) e^{-\Gamma_C(r-r_H)} + \Gamma_C \int^r_{r_H} dr' e^{-\Gamma_C(r-r')}\phi(r) e^{\nu A(r)},\]

shows that the resulting estimate of $\phi_A$ corresponds to a weighted average over an interval $\Delta r \sim 1/\Gamma$.

The relaxation rate $\Gamma$ should be large enough for $C(r)$ to quicly relax to $\phi_A$, but small enough to make the average interval sufficiently large for numerical stability. At this point, this parameter is hardcoded to be REL_RATE = 0.5 in blackhole_engineering.h. Results seem to be stable to changing this parameter by a factor of 2, but this should be incorporated to automated tests as a sanity check.

This relaxational method has been found to yield results similar to the previous fitting procedure, while being faster and numerically more stable. As an independent test of its robustness, we have implemented another extra differential equation in `BlackHole::solve() <#521-Class-BlackHole>`__: $\frac{dD(r)}{dr} = -\Gamma(D(r) - \Phi'(r)\,e^{2A(r)})$, which can be used to extract the factor $\Phi_2^{far}\approx D(r\to \infty)$ in the near-boundary scaling of $\Phi(r) \sim \Phi_0^{far} + \Phi_2^{far}e^{-2 A(r)}$. Because the same $\Phi_2^{far}$ is extracted from equation (17) in Ref. [2], this equation provides a reliable test of both the relaxational method and the choice of REL_RATE. It should also be incorporated to automated tests.

1.5.2. Precision check

After the integration of the black hole is performed, the accuracy is checked by comparing the Ricci scalar and the $\phi_A$ with the desired values.

It is assumed that the EMD fields reach the their vultraviolet behavior corresponding to the $\text{AdS}_5$ geometry, which has a Ricci scalar $R = -20$. Therefore, for a set accuracy $\epsilon_R$, for the calculation to be valid, we require $|R_{calculated} - R_{\text{AdS}_5}| < \epsilon_R$. If the users do not specify the accuracy in the input file, the default value will be used.

As noted in section 1.5.1, ideally $C \approx \phi_A$. Therefore, for a set accuracy $\epsilon_A$, for the calculation to be valid, we require that $|C - \phi e^{\nu A}|/C < \epsilon_A$. If the users do not specify the accuracy in the input file, the default value will be used.

If the accuracy is not reached, the code will increase the radius for integration $r$ to get within the accuracy. If the accuracy is not reached for several iterations, this point will be abandoned and the data related will not enter the final table.

1.5.3. Baryon susceptibility at $\mu_B=0$

As $\mu_B \to 0$, the second order baryon susceptibility is given by

\[\chi_2^B \equiv \left.\frac{\partial \rho_B}{\partial \mu_B}\right|_T = \frac{\rho_B}{\mu_B} .\]

The baryon density $\rho_B$, entropy $s$, and baryon chemical potential $\mu_B$ can be written in terms of the near-boundary EMD fields, as in Ref. [2],

\[\rho_B = - \frac{\Phi_2^{far}}{\kappa^2_5 \phi_A^{3/\nu} \sqrt{h_0^{far}}} \Lambda^3 ,\]

\[s = \frac{2\pi}{\kappa_5^2 \phi_A^{3/\nu}}\Lambda^3 ,\]

\[\mu_B = \frac{\Phi_0^{far}}{\phi_A^{1/\nu} \sqrt{h_0^{far}}} \Lambda ,\]

where $\Lambda$ is the energy scale, $\kappa_5^2 \equiv 8\pi G_5$, and $G_5$ is the five-dimensional gravitational constant.

Also using the Gauss charge and the fact that $\Phi \to 0$ as $r \to r_H$,

\[Q_F = f(\phi) e^{A-B} \Phi ,\]

\[\Phi_0^{far} = \int_{r_H}^\infty dr \Phi ,\]

We can obtain the expression of the dimentionless baryon susceptibility $\hat{\chi}_2^B \equiv \chi_2^B/T^2$ at $\mu_B = 0$,

\[\chi_2^B (\mu_B = 0) = \frac{1}{16\pi^2} \frac{s}{T^3} \left[f(0) \int_{r_H}^\infty dr\; e^{-2A(r)} f(\phi(r))^{-1} \right]^{-1}.\]

The baryon susceptibility at $\mu_B = 0$ is used to match the holographic model with the lattice QCD data to fix the dilaton-gauge coupling $f(\phi)$.

1.6. Integrating the pressure

After the user specified the range of $T$ and $\mu_B$ for the phase diagram, a table of phi0 and Phi1 will be created. The pressure of the point with maximum phi0 and minimum Phi1 will be set to 0, and the pressure of the other points will be generated by numerical integration.

First, fixing Phi1 to be the minimum, we loop the phi0 values to get the pressure values a constant Phi1 line. For each iteration in the phi0 values, we extract the values of $T$, $\mu_B$, $s$, and $n$ for the current phi0 and the next phi0 in line. The change in pressure is then given by the Gibbs–Duhem relation.

\[dp = \left( n(i) \left.\frac{d\mu_B}{dT}\right|_{\Phi_1} + s(i) \right) dT,\]

where i is the index of the phi0 values, and the derivative and difference are calculated such that

\[\left.\frac{d\mu_B}{dT}\right|_{\Phi_1} = \frac{\mu_B(i+1) - \mu_B(i)}{T(i+1) - T(i)} ,\]

\[dT = T(i+1) - T(i) .\]

Then, the pressure is calculated such that

\[p(i+1) = p(i) + dp .\]

To get the pressure for the whole phase diagram, we loop the phi0 values, and for each phi0 value, we loop the Phi1 values. For each iteration, because the rough correspondense of phi0 to $T$ and Phi1 to $\mu_B$, the change in pressure is better calculated by

\[dp = \left( n(i) d\mu_B + s(i) \left.\frac{dT}{d\mu_B}\right|_{\phi_0} \right) d\mu_B\]

where i is the index of the Phi1 values, and the derivative and difference are calculated such that

\[d\mu_B = \mu_B(i+1) - \mu_B(i) .\]

Then, the pressure is calculated such that

\[p(i+1) = p(i) + dp .\]

1.7. Input format

The users need to use a yaml file to speficy the parameters used for the generation of the phase diagram and calculation of the critical point. The current input file has four sections that need to be specified. In practice, the users only need to care about chaninge the parameters in bh_eos/input/holo_eos_input_user.yaml and run through the input validator to check if input file is still valid. The details can be found in section 4 Usage.

The first section is to choose the model and specify the parameters. Here, as an example, the model is chosen to be the polynomial hyperbolic model, and the parameters associated with it is specified as follows in the yaml input file:

model_options:
variant: polynomial_hyperbolic
parameters:
    description: Parameters of the model in question.
    flag_default: no
    Lambda: 1148.121631179041
    kappa2: 11.350245262021625
    gamma: 0.5884489971316446
    b2: 0.301821894469367
    b4: -0.04261907765391285
    b6: 0.0004566058948051932
    c1: -0.06469185351673243
    c2: -0.22517557265138566
    c3: 0.03664412474454205
    d1: 1.7231476382001352
    d2: 472.2764686476346

The second section is the output, which determines what files should be generated. There are a few types of files that can be generated, depending on the extension. The extensions allowed are .dat for csv, .yaml for yaml.:

output_options:
    yaml_output_file: muses_polyhyperbolic_bestfit.yaml
    output_path: muses_polyhyperbolic_bestfit

The third section is for the generation of the phase diagram:

eos_options:
    description: Options for EoS calculation.
    eos_stages:
        description: Which stages of the calculation to perform.
        interpolate_points: true
        maxwell_construction: true

    temperature_options:
        description: Determines which values of the temperature will be calculated for the equation of state.
        T_min: 20
        T_max: 400
        T_step: 2.5
        N_phi0_lines: 200

    chemical_potential_options:
        description: Determines which values of the baryon chemical potential will be calculated for the equation of state.
        mu_B_min: 0
        mu_B_max: 2000
        mu_B_step: 5
        N_Phi1_lines: 1000

The fourth section is for the calculation of the critical points:

critical_point_options:
    description: Options for searching for the critical point position
    yaml_output: critical_point_interp.yaml
    phi0_min: 0.5
    phi0_max: 10
    N_lines: 800
    max_tries: 1000
    initial_Phi1_step: 1e-5
    max_Phi1_step: 0.2
    prec_T: 0.0001
    prec_mu: 0.0001
    acc_T: 0.001
    acc_mu: 0.001

References

[1] R. Critelli, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti and R. Rougemont, ``Critical point in the phase diagram of primordial quark-gluon matter from black hole physics,’’ Phys. Rev. D 96 (2017) no.9, 096026
[2] J. Grefa, J. Noronha, J. Noronha-Hostler, I. Portillo, C. Ratti and R. Rougemont, ``Hot and dense quark-gluon plasma thermodynamics from holographic black holes,’’ Phys. Rev. D 104 (2021) no.3, 034002
[3] O. DeWolfe, S. S. Gubser and C. Rosen, ``A holographic critical point,’’ Phys. Rev. D 83 (2011), 086005
[4] J. Knaute, R. Yaresko and B. Kämpfer, ``Holographic QCD phase diagram with critical point from Einstein–Maxwell–dilaton dynamics,’’ Phys. Lett. B 778 (2018), 419-425
[5] R. G. Cai, S. He, L. Li and Y. X. Wang, ``Probing QCD critical point and induced gravitational wave by black hole physics,’’ arXiv:2201.02004 [hep-th]

2. Dependencies

Lib Boost ODEint (tested with boost 1.6.7).
Lib yaml-cpp
Lib yaml_wrap / lib muses_yaml
[python3] (and some required packages listed in requirements.txt)

So far, the code does not seem to run natively on Windows. Mac users might also face a few issues. These problems should be addressed by using Docker containers.

3. Docker file

The base image with the underlying OS image and compiled Boost library is defined in the /base directory. The following build command uses that image.

docker build -t registry.gitlab.com/nsf-muses/module-holographic-eos/bh_eos/bh-eos:dev .

4. Usage

The user can modify the input in the file bh_eos/input/holo_eos_input_user.yaml. They need to specify either “polynomial_hyperbolic” or “muses_parametric” as the model type. They can also modify the model parameters in the files “holo_musesparametric.yaml” or “holo.polyhyperbolic.yaml”, but it is not recommended. After the user made the changes to input, suppose they are in the home directory, they need to run the following commands to validate the input file: - cd open-API - python3 validate_combine_input.py openapi.yaml’ - cd ..

The executables can be generated by compiling through the make function, with the command in the terminal - make

The compiling can be made faster by parallelization, with the command - make -j# where # denotes the numbers of cores used for compiling.

To clean the compiled objects and recompile, one can use the following command, - make clean - make distclean

The compiled executables will be in the folder “bin”. Suppose one is in the folder of the module “bh_eos”, the phase diagram can be generated by the following command - ./bin/muses_numrelholo.exec ./input/input_full.yaml The command will generate the tables for the phase diagram. By default, all the output files will be saved in the folder “output”.

5. Output files

Note that because we define the zero pressure to be the on the line where $\mu_B=0$. If the user set the minimum $\mu_B!=0$, all the pressure and energy density will be NaN in the table.

5.1 Yaml output file

The user can name this file in the yaml node “yaml_output_file” of the input file. This file is stored in the repo home directory and shows info of the format of the output files.

5.2 EoS output files

The user can specify the folder name where they want to store the EoS output files in the yaml node “output_path” of the input file.

Inside the folder, there are several different output files.

5.2.1 phi0Phi1_lines.csv

This table is directly calculated from the initial conditions given by phi0 Phi1 values.

5.2.2 eos.csv

This table is interpolated based on the output file phi0Phi1_lines.csv and thus has a uniform grid in $T$ and $\mu_B$. It has all the phases: stable, unstable, and meta phases.

5.2.3 stable_eos.csv

This table is obtained by performing a Maxwell construction on the output file eos.csv and thus only has the stable phase.

5.2.4 transition_line.csv

This table records the first-order transition line and relevant thermodynamics given by the EoS. The calculation is not precise enough yet and the first-order transition line does not go all the way to the critical point.

5.2.5 spinodal_lines.csv

This table records the spinodal lines and relevant thermodynamics.