aainstitutetext: School of Mathematical Sciences, Queen Mary University of
London, Mile End Road, London E1 4NS, United Kingdombbinstitutetext:
Department of Physics and Helsinki Institute of Physics, PL 64, FI-00014
University of Helsinki, Finland

# Evolutions in first-order viscous hydrodynamics

Hans Bantilan b Yago Bea a Pau Figueras<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Motivated by the physics of the quark-gluon plasma created in heavy-ion
collision experiments, we use holography to study the regime of applicability
of various theories of relativistic viscous hydrodynamics. Using the
microscopic description provided by holography of a system that relaxes to
equilibrium, we obtain initial data with which we perform real-time evolutions
in 2+1 dimensional conformal fluids using the first-order viscous relativistic
hydrodynamics theory of Bemfica, Disconzi, Noronha and Kovtun (BDNK), BRSSS
and ideal hydrodynamics. By initializing the hydrodynamics codes at different
times, we can check the constitutive relations and assess the predictive power
and accuracy of each of these theories.

##  1 Introduction

Relativistic hydrodynamics is a very general theory that provides the
effective description of the out-of-equilibrium real-time dynamics of many
physical systems of interest. In fact, hydrodynamics is better understood as a
low energy effective field theory (EFT) constructed as a derivative expansion
near a local thermodynamic equilibrium. As such, the fundamental variables are
the basic quantities describing the equilibrium solution, such as energy
density and velocities, and their gradients. The zeroth order piece in the
gradient expansion corresponds to ideal hydrodynamics.

The equations of ideal hydrodynamics are hyperbolic, having a well-posed
initial value problem for general initial conditions, and are amenable to
numerical simulations. However, there are important physical situations where
dissipative effects may play a relevant role, such as in the Quark-Gluon
Plasma (QGP) Romatschke:2017ejr ; Busza:2018rrf  or in certain astrophysical
systems Alford:2017rxf ; Chabanov:2021dee ; Shibata:2017jyf ;
Fujibayashi:2017puw ; Fujibayashi:2020jfr . However, it is well-known that if
first-order viscous terms are included in the relativistic hydrodynamic
expansion (in the Landau frame), the resulting equations of motion are
parabolic and hence incompatible with the causality postulates of relativity
Hiscock:1985zz .

A well-known approach to address this issue of lack of hyperbolicity and
causality is exemplified by the work of Müller, Israel and Stewart (MIS)
Muller:1967zza ; Israel:1979wp ; Israel:1976tn . The strategy consists in
including second order terms and introduce new variables with the respective
equations of motion, such that the overall system of equations is hyperbolic
and causal. The new theory should have the same IR, i.e., hydrodynamic regime,
but different UV properties.111It is not clear that such modification of the
theory correctly describes the physics in situations where there is a strong
UV cascade, such as in a 3+1 turbulent flow. The MIS formulation is not
unique; there are several formulations which include extra variables and
equations to recover hyperbolicity: BRSSS Baier:2007ix , DNMR Denicol:2012cn ,
divergence type theories Geroch:1990bw ; Lehner:2017yes , etc.. As discussed
in Geroch:1995bx ; Geroch:2001xs , in principle all these theories should be
equivalent in the IR.

In recent years it has been realized that it is possible to write the
equations of first-order viscous hydrodynamics in a way such that they are
hyperbolic, without the need of adding extra variables or equations. When
including first order terms in the hydrodynamic expansion, the basic
hydrodynamic variables present an ambiguity when defined out of equilibrium,
and fixing this ambiguity is known as choosing a ‘frame’. Traditionally, it
was common to use only natural frames such as the Landau frame. The insight of
BDNK Bemfica:2017wps ; Bemfica:2019knx ; Kovtun:2019hdm  was to notice that
the choice of frame affects the hyperbolicity of the equations of motion and a
certain frame can be chosen such that the latter are hyperbolic.222See
Armas:2022wvb  for a recent generalization to magnetohydrodynamics.

MIS-based hydrodynamic codes have been extensively used for many years to
study dissipative effects. Given that there is now a new theory of viscous
relativistic hydrodynamics, it is natural to explore it. The first ones to do
so were the authors of Pandya:2021ief , who studied the evolution of smooth
and non-smooth initial data using BDNK in effectively 1+1 dimensions. The
purpose of the present paper is to carry out the first studies of BDNK
focusing on applications to the QGP. More concretely, in this article we
report on studies of dynamical evolutions of the BDNK equations in a 2+1
dimensional conformal uncharged fluid. In order to acquire a better
understanding of BDNK as a theory of viscous hydrodynamics, we quantitatively
compare it to ideal and BRSSS Baier:2007ix  hydrodynamics, as well as to the
microscopic theory provided by holography.

The QGP created in heavy-ion collision experiments is initially far from
equilibrium and it subsequently relaxes to a hydrodynamic regime. Hydrodynamic
numerical codes are used to provide effective descriptions of these
evolutions, where the underlying microscopic theory is Quantum Chromodynamics
(QCD). Motivated by this picture, in this paper we use numerical general
relativity and holography to obtain a first principles description of the
microscopic real-time dynamics of the stress tensor in a strongly coupled
field theory of a system which relaxes to equilibrium. The microscopic
solution is valid for all times, covering both the far-from-equilibrium and
near equilibrium hydrodynamic regimes. Whilst previous studies of the
applicability of hydrodynamics in holographic systems have focused on checking
the constitutive relations, in this paper we go beyond the state of the art
and use the microscopic data to initialize the hydrodynamic codes at different
times 333Previous studies in a one-dimensional setting were performed in
Heller:2014wfa .. In this way, by comparing with the UV-complete solution
provided by holography, we verify both the constitutive relations and the
predictivity of each hydrodynamic theory for a certain class of initial
conditions. We use units $c=G_{4}=1$.

When our article was nearing completion, Pandya:2022pif  appeared on arXiv,
extending the numerical studies of the BDNK equations.

##  2 Hydrodynamics: the equations

We use three sets of hydrodynamic evolution equations: ideal hydrodynamics,
BRSSS and BDNK. We write these equations for a $2+1$ dimensional conformal
fluid in Minkowski spacetime. Conformal symmetry fixes the equation of state:

$p=\frac{\epsilon}{2}\,,$ (1)

where $\epsilon$ is energy density and $p$ is pressure. The temperature $T$ is
related to the energy density as $\epsilon=\frac{2}{3}A\,T^{3}$, where $A$ is
an arbitrary constant that in our simulations we fix to $A=4\pi^{2}/9$.444For
a holographic fluid in $d$ spacetime dimensions, this constant $A$ is related
to the bulk’s Newton’s constant of gravitation $G_{d+1}$ as
$A=\frac{(4\pi)^{d}}{16\pi\,G_{d+1}d^{d-1}}$. Considering the 3-velocity of
the fluid $u^{\mu}$, normalized such that $u^{2}=-1$, we define
$\dot{\epsilon}\equiv u^{\mu}\nabla_{\mu}\epsilon$,
$\nabla_{\perp}^{\mu}\equiv\Delta^{\mu\nu}\nabla_{\nu}$, $\nabla\cdot
u\equiv\nabla_{\rho}u^{\rho}$ , $\Delta^{\mu\nu}\equiv
g^{\mu\nu}+u^{\mu}\,u^{\nu}$,
$\Omega^{\mu\nu}\equiv\left(\nabla^{\mu}_{\perp}u^{\nu}-\nabla^{\nu}_{\perp}u^{\mu}\right)/2$,
$\sigma^{\mu\nu}\equiv 2\nabla^{\langle\mu}u^{\nu\rangle}$, where

$\displaystyle
A^{\langle\mu\nu\rangle}\equiv\frac{1}{2}\Delta^{\mu\alpha}\Delta^{\nu\beta}(A_{\alpha\beta}+A_{\beta\alpha})-\frac{1}{2}\Delta^{\mu\nu}\Delta^{\alpha\beta}A_{\alpha\beta}\,,$
(2)

is symmetric, traceless and transverse to $u_{\mu}$. The constitutive
relations of relativistic hydrodynamics up to second order in the derivative
expansion in the Landau frame can be written as Baier:2007ix ,

$T^{\mu\nu}=\epsilon\,u^{\mu}\,u^{\nu}+p\,\Delta^{\mu\nu}+\Pi^{\mu\nu}\,,$
(3a) with
$\displaystyle\Pi^{\mu\nu}=-\eta\,\sigma^{\mu\nu}+\eta\,\tau_{\pi}\left(\dot{\sigma}^{\langle\mu\nu\rangle}+\frac{1}{2}\,\sigma^{\mu\nu}\,\nabla\cdot
u\right)$
$\displaystyle+{\lambda_{1}}\,{\sigma^{\langle\mu}}_{\rho}\sigma^{\nu\rangle\rho}+{\lambda_{2}}\,{\sigma^{\langle\mu}}_{\rho}\Omega^{\nu\rangle\rho}+\lambda_{3}\,{\Omega^{\langle\mu}}_{\rho}\Omega^{\nu\rangle\rho}\,,$
(3b)

where $\eta$ is the shear viscosity and $\tau_{\pi}$, $\lambda_{1}$,
$\lambda_{2}$, $\lambda_{3}$ are second order transport coefficients. The
constitutive relations of ideal hydrodynamics are given by (3a) with
$\Pi^{\mu\nu}=0$.

The conservation of the stress-energy tensor

$\nabla_{\mu}T^{\mu\nu}=0\,,$ (4)

provides the evolution equations for the dynamical variables $\epsilon$ and
the independent components of the velocity vector $u^{\mu}$. The equations of
ideal hydrodynamics are obtained by plugging (3a) with $\Pi^{\mu\nu}=0$ into
(4).

The equations of BRSSS Baier:2007ix  are obtained from (3) and (4) by
promoting $\Pi^{\mu\nu}$ to a new variable; then, using the first order
relation $\Pi^{\mu\nu}=-\eta\sigma^{\mu\nu}$, equation (3b) becomes an
independent evolution equation for $\Pi^{\mu\nu}$:

$\displaystyle\Pi^{\mu\nu}=-\eta\sigma^{\mu\nu}-\tau_{\pi}\left(\dot{\Pi}^{\langle\mu\nu\rangle}+\frac{3}{2}\Pi^{\mu\nu}\nabla\cdot
u\right)$
$\displaystyle+\frac{\lambda_{1}}{\eta^{2}}{\Pi^{\langle\mu}}_{\rho}\Pi^{\nu\rangle\rho}-\frac{\lambda_{2}}{\eta}{\Pi^{\langle\mu}}_{\rho}\Omega^{\nu\rangle\rho}+\lambda_{3}{\Omega^{\langle\mu}}_{\rho}\Omega^{\nu\rangle\rho}\,.$
(5)

We now present the BDNK equations. Given a timelike unit vector $u^{\mu}$, a
symmetric tensor can be decomposed as

$T^{\mu\nu}=\mathcal{E}\,u^{\mu}\,u^{\nu}+\mathcal{P}\,\Delta^{\mu\nu}+(\mathcal{Q}^{\mu}\,u^{\nu}+\mathcal{Q}^{\nu}\,u^{\mu})+\mathcal{T}^{\mu\nu}\,,$
(6)

where

$\displaystyle\mathcal{E}$ $\displaystyle\equiv
u_{\mu}\,u_{\nu}\,T^{\mu\nu}\,,$ (7a) $\displaystyle\mathcal{P}$
$\displaystyle\equiv\frac{1}{2}\Delta_{\mu\nu}T^{\mu\nu}\,,$ (7b)
$\displaystyle\mathcal{Q}_{\mu}$
$\displaystyle\equiv-\Delta_{\mu\alpha}u_{\beta}T^{\alpha\beta}\,,$ (7c)
$\displaystyle\mathcal{T}_{\mu\nu}$
$\displaystyle\equiv\frac{1}{2}\left(\Delta_{\mu\alpha}\Delta_{\nu\beta}+\Delta_{\nu\alpha}\Delta_{\mu\beta}-\Delta_{\mu\nu}\Delta_{\alpha\beta}\right)T^{\alpha\beta}\,.$
(7d)

For a conformal fluid, the expansion in derivatives of each component to first
order is

$\displaystyle\mathcal{E}$
$\displaystyle=\epsilon+2\pi_{2}\left(\frac{2}{3}\frac{\dot{\epsilon}}{\epsilon}+\,\nabla\cdot
u\right)\,,$ (8a) $\displaystyle\mathcal{P}$
$\displaystyle=p+\pi_{2}\,\left(\frac{2}{3}\frac{\dot{\epsilon}}{\epsilon}+\,\nabla\cdot
u\right)\,,$ (8b) $\displaystyle\mathcal{Q}^{\mu}$
$\displaystyle=\theta_{1}\left(\dot{u}^{\mu}+\frac{1}{3}\frac{\nabla_{\perp}^{\mu}\epsilon}{\epsilon}\right)\,,$
(8c) $\displaystyle\mathcal{T}^{\mu\nu}$
$\displaystyle=-\eta\,\sigma^{\mu\nu}\,,$ (8d)

where $\pi_{2}$, $\theta_{1}$ and $\eta$ are transport coefficients. Conformal
symmetry fixes the dimensions of $\pi_{2}$ and $\theta_{1}$, and we write them
as a constant times the shear viscosity $\eta$

$\pi_{2}=a_{2}\,\eta\,,~{}~{}~{}~{}\theta_{1}=a_{1}\,\eta\,.$ (9)

The stress tensor, to first order in the derivative expansion, is then

$\displaystyle T^{\mu\nu}=$
$\displaystyle\left[\epsilon+2\,a_{2}\eta\left(\frac{2}{3}\frac{\dot{\epsilon}}{\epsilon}+\nabla\cdot
u\right)\right]\left(u^{\mu}\,u^{\nu}+\frac{\Delta^{\mu\nu}}{2}\right)$ (10)
$\displaystyle+a_{1}\eta\left[\left(\dot{u}^{\mu}+\frac{1}{3}\frac{\nabla_{\perp}^{\mu}\epsilon}{\epsilon}\right)u^{\nu}+(\mu\leftrightarrow\nu)\right]$
$\displaystyle-\eta\,\sigma^{\mu\nu}+O(\partial^{2})\,.$

The BDNK equations are obtained from plugging (10) into the conservation
equation (4).

The constants $a_{1}$ and $a_{2}$ specify the frame. The Landau frame
corresponds to $a_{1}=a_{2}=0$, by which we recover (3) up to first order
terms. The BDNK equations are hyperbolic iff Bemfica:2019hok

$a_{2}>1\,,~{}~{}~{}~{}a_{1}>\frac{4\,a_{2}}{a_{2}-1}\,.$ (11)

The Landau frame lies outside this region, and there is a gap between the
causal frames and the Landau frame.

The frame change to Landau frame is straightforwadly obtained from (10)

$\displaystyle\epsilon_{\text{Landau}}$
$\displaystyle=\epsilon+2\,a_{2}\,\eta\left(\frac{2}{3}\frac{\dot{\epsilon}}{\epsilon}+\nabla\cdot
u\right)\,,$ (12a) $\displaystyle u^{\mu}_{\text{Landau}}$
$\displaystyle=u^{\mu}+\frac{2\,a_{1}\,\eta}{3\epsilon}\left(\dot{u}^{\mu}+\frac{1}{3}\frac{\nabla_{\perp}^{\mu}\epsilon}{\epsilon}\right)\,.$
(12b)

In our simulations we use the transport coefficients of the microscopic
holographic theory that will be presented in the next Section Herzog:2002fn ;
Natsuume:2007ty ; VanRaamsdonk:2008fp :

$\displaystyle\eta$ $\displaystyle=$ $\displaystyle\frac{s}{4\pi}~{},$ (13a)
$\displaystyle\tau_{\pi}$ $\displaystyle=$ $\displaystyle\frac{1}{24\pi
T}\left(18-9\ln 3+\sqrt{3}\pi\right)~{},$ (13b) $\displaystyle\lambda_{1}$
$\displaystyle=$ $\displaystyle\frac{3s}{32\pi^{2}T}~{},$ (13c)
$\displaystyle\lambda_{2}$ $\displaystyle=$ $\displaystyle\frac{2\pi
T}{27}\left(\sqrt{3}\pi-9\ln 3\right)~{},$ (13d) $\displaystyle\lambda_{3}$
$\displaystyle=$ $\displaystyle 0~{}.$ (13e)

For details on the implementation of the hydrodynamic equations in the
numerical codes and convergence tests see Appendix A.

##  3 Holography: microscopic evolution

Our model is Einstein’s gravity with negative cosmological constant coupled to
a massless scalar field in four spacetime dimensions.555The scalar field is
used to form a black brane by prompt collapse of the initial data. This prompt
colapse happens immediately after initialization of the code and this
timescale is well separated from the timescales relevant for the physics
studied in this paper. Thus, after the immediate collapse the scalar field is
completely negligible and the evolution will be well described by pure
gravity. This pure gravity theory can be obtained as a consistent truncation
in a set of top-down theories which includes the ABJM theory Aharony:2008ug .
By holography, it describes the decoupled sector of the stress tensor on the
conformal field theory side. The real-time quantum dynamics on the field
theory maps to the bulk classical dynamics of Einstein’s gravity in anti-de
Sitter (AdS) space coupled to a massless scalar field. We use numerical
relativity techniques to solve for the classical dynamics of the bulk fields
employing the same code as in Bantilan:2020pay ; Bantilan:2020xas . For more
details of our implementation see these two references and Appendix B.

We consider the gravitational collapse of massless scalar field in the 3+1
dimensional Poincaré patch of AdS using time symmetric initial data. The
scalar field has an initial deformed Gaussian profile along the boundary
directions and is localized in the AdS radial direction. The initial data is
“strong” in the sense that there is a trapped surface on the initial Cauchy
surface. Therefore, one can think of our initial data as corresponding to a
highly deformed black brane.

From the dual CFT point of view, the initial state corresponds to a large
localized perturbation on top of a homogeneous plasma at temperature
$\overline{T}$. The initial energy density profile on the field theory is
shown in Fig. 1 (top). The initial data is rotationally symmetric and has
vanishing initial velocities (since the initial data is time-symmetric), even
though our code does not assume rotational symmetry along the boundary
directions. We have also considered initial data with small deformations that
break the rotational symmetry along the boundary directions. However, the non-
rotationally symmetric modes decay exponentially and by the time the system
approaches the hydrodynamic regime, those modes have already decayed. We
measure all quantities in units of $\overline{T}$ or
$\overline{\mathcal{E}}\equiv\frac{2}{3}A\overline{T}^{3}$ and we use
Cartesian coordinates $\\{t,x,y\\}$ to label the boundary directions.

The time evolution of the system is shown in Fig. 1 in three representative
snapshots at times $t\overline{T}=0,0.08,0.16$. The initial peak explodes and
the system expands and disperses away, and at late times it relaxes to a
homogeneous solution.

Figure 1: Snapshots of the energy density in the lab frame, $T_{tt}$, of the
holographic solution at times $t\overline{T}=0,0.08,0.16$, from top to bottom.

##  4 Hydrodynamics: evolutions

In the following we study the applicability of hydrodynamics in this system.
We do this in two ways: first we check the constitutive relations pointwise in
spacetime and, second, we perform time evolution of the hydrodynamic equations
using the holographic solution as initial data.

We start by evaluating the constitutive relations of hydrodynamics (3)
pointwise in spacetime using the holographic data to calculate the various
terms in this expression. In the Landau frame, gradient corrections in the
hydrodynamic expansion are transverse to the velocity and, together with
conformal symmetry, this implies that they vanish at the center of a
spherically symmetric system. Away from the center however the gradients are
non-trivial; we choose $\\{x,y\\}\overline{T}\simeq\\{0.12,0\\}$ as a
representative point where to evaluate the constitutive relations.

In Fig. 2 we compare the size of the 1st and 2nd order corrections to the
ideal term in the derivative expansion of the stress tensor by plotting
$T_{xx}^{1st}/T_{xx}^{id}$ and $T_{xx}^{2nd}/T_{xx}^{id}$ as functions of
time. This figure indicates that, according to the constitutive relations, the
system is initially far from equilibrium, as second order gradients are as
large as $40\%$ compared to the ideal terms. The size of the derivative
corrections quickly decays and the system has hydrodynamized at times around
$t\overline{T}\simeq 0.5$, if our criterium is that the ratios
$T_{xx}^{1st}/T_{xx}^{id}$ and $T_{xx}^{2nd}/T_{xx}^{id}$ are smaller than
$2\%$.

Figure 2: Ratios $T_{xx}^{n}/T_{xx}^{ideal}$, $n=1,2$, of the constitutive
relations (3) evaluated at $\\{x,y\\}\overline{T}\simeq\\{0.12,0\\}$ as a
function of time for the holographic solution. The vertical dashed lines show
the times at which we initialize the hydrodynamics codes $t\overline{T}\simeq
0.019,0.16,0.58$.

We now consider the evolutions of the hydrodynamic equations of the various
theories using initial data from the holographic solution. We start by
diagonalizing the holographic stress tensor and obtaining the local energy
density and velocities (these are the energy and velocities in the Landau
frame), which will be used to initialize the hydro codes. In the particular
case of BRSSS, we have to also specify the dissipative part of the stress
tensor $\Pi^{\mu\nu}$ from the microscopic holographic data. But notice that
the use of time symmetric initial data implies that at $t=0$ the dissipative
part of the stress tensor vanishes identically and so in the subsequent
evolution it is fairly small and hence one could initialize the code using
very small initial data for the viscous tensor $\Pi^{\mu\nu}$. However,
initially the second order terms are large and therefore, alternatively, one
could initialize de code using the dissipative part of the stress tensor
computed with the first and second order consitutive relations (3), which is
not small. These two possible ways of initializing the code match in the
hydrodynamic regime (up to 3rd order terms), but under these far-from-
equilibrium conditions they differ. We choose the second option, which
captures the initial presence of large gradients.

For the evolutions in the BDNK theory, we consider the following causal
frames:

$\displaystyle\text{Frame}\ 1:$ $\displaystyle\quad\quad\quad a_{1}=10,\quad
a_{2}=10\,,$ (14a) $\displaystyle\text{Frame}\ 2:$
$\displaystyle\quad\quad\quad a_{1}=6,\quad\,\ a_{2}=4\,,$ (14b)

which satisfy the hyperbolicity conditions (11). We obtain the energy and
velocities in the causal frame from the Landau frame by inverting the
expressions (12) and neglecting second order terms. The inverted expressions
are the following:

$\displaystyle\epsilon$
$\displaystyle=\epsilon_{\text{Landau}}-2\,a_{2}\,\eta\left(\frac{2}{3}\frac{\dot{\epsilon}_{\text{Landau}}}{\epsilon_{\text{Landau}}}+\nabla\cdot
u_{\text{Landau}}\right)\,,$ (15a) $\displaystyle u^{\mu}$
$\displaystyle=u^{\mu}_{\text{Landau}}-\frac{2\,a_{1}\,\eta}{3\epsilon}\left(\dot{u}_{\text{Landau}}^{\mu}+\frac{1}{3}\frac{\nabla_{\perp}^{\mu}\epsilon_{\text{Landau}}}{\epsilon_{\text{Landau}}}\right)\,,$
(15b)

Where we have used the fact that in the first order terms we can replace
$\epsilon=\epsilon_{\text{Landau}}$, $u^{\mu}=u^{\mu}_{\text{Landau}}$ as the
difference is second order. Recall that the evolution equations of BDNK are of
second order and hence one also has to specify the time derivatives of the
evolved variables in the initial data. We compute these by using the
holographic data to obtain the constitutive relations in the corresponding
causal frame as functions of time using (15), and then calculating the time
derivatives in the causal frame. Other procedures might be equivalent in the
hydrodynamic regime up to second order terms, but in the far from equilibrium
regime they will generically differ.

We start the discussion by initializing the hydrodynamic codes using the
holographic data at $t\overline{T}\simeq 0.019$, when the system is still far
from equilibrium according to the constitutive relations. Fig. 3 (top) shows
the evolution of the energy density in the lab frame at a representative off-
center point of the domain $\\{x,y\\}\overline{T}=\\{0.17,0\\}$ obtained using
ideal (dashed green), BRSSS (dashed blue), BDNK in frame 2 (purple curves) and
the holographic solution (solid black curve). For the BDNK evolutions, we plot
two different quatities: the result of obtaining $T_{tt}$ from the evolved
variables using (10) (solid line), and the result of evaluating only the 0th
order terms (ideal tems) in (10) (dotted line). The reason for doing this is
to show explicitly the size of the first order terms compared to the ideal
terms. As this figure shows, the three hydrodynamic theories that we consider
initialized at $t\overline{T}\simeq 0.019$ exhibit a similar level of (large)
disagreement with the microscopic theory throughout the evolution. The fact
that none of the theories of hydrodynamics provides an accurate description of
the evolution of the fluid confirms that at this initial time, the system is
still very far from the hydrodynamic regime. This is expected since the
gradients are large initially and each theory has a different UV completion.
Therefore, there is no reason why these theories should agree with each other
away from the hydrodynamic regime.

Figure 3: Energy density in the lab frame, $T_{tt}$, at an off center location
$\\{x,y\\}\overline{T}=\\{0.17,0\\}$ as a function of time. Black continuous
line corresponds to the microscopic (holographic) solution. We include the
results of the hydrodynamic evolutions initialised with holographic data at
$t\overline{T}\simeq 0.019,0.16,0.58$, from top to bottom. For the BDNK
evolutions, we also include the result of using only the 0th order terms (the
ideal part) in (10), in dotted lines. The other components of the stress
energy tensor behave in a quantitatively similar manner.

We continue the discussion by considering hydrodynamic evolutions with
holographic initial data at a later time $t\overline{T}\simeq 0.16$, when the
size of the gradients is smaller but the system has not hydrodynamized yet
according to the constitutive relations, see Fig. 3 (middle panel). Indeed, at
this initial time the constitutive relations indicate that the first and
second order terms are comparable, see the second vertical line in Fig. 2. We
find that all hydrodynamic evolutions still exhibit large deviations from the
microscopic theory which confirms that the system has not reached the
hydrodynamic regime yet. Again, the deviations of BRSSS and BDNK from the
microscopic theory are comparable.

Finally, we consider hydrodynamic evolutions with holographic initial data at
a later time $t\overline{T}\simeq 0.58$, when the size of the gradients is
small, see Fig. 3 (bottom panel). At $t\overline{T}\simeq 0.58$, by checking
the constitutive relations (3) we find that the first order gradients have a
maximum value in the whole domain of the order of $~{}1\%$, i.e.,
$T_{xx}^{1st}/T_{xx}^{id}\sim 1\%$, and similarly for the second order
gradients, i.e., $T_{xx}^{2nd}/T_{xx}^{id}\sim 1\%$. These values are small,
and we expect the hydrodynamic evolutions to follow the microscopic solution
when initialized at this time. Our results shown in the bottom panel of Fig. 3
confirm this expectation. This is true within the numerical errors of the
holographic solution, which are inherited by the hydrodynamic evolutions
through the initial data. We have checked that the energy density in the lab
frame for the hydrodynamic evolutions differs by less than $\sim 2\%$ for
BRSSS, BDNK and ideal hydrodynamics, compared to the microscopic evolution in
whole the domain and at all times beyond $t\overline{T}\simeq 0.58$. These
results also confirm that the various theories of hydrodynamics provide
equivalent descriptions of the system, further supporting the applicability of
hydrodynamics.

Fig. 3 (bottom) suggests that ideal hydrodynamics provides a better
description of the system at late times than viscous hydrodynamics. The reason
might be the following. Our particular choice of initial data implies that the
holographic system at $t=0$ is exactly described by the constitutive relations
of ideal hydrodynamics (3a) with $\Pi^{\mu\nu}=0$ (even if second order terms
are large and the system is not within the regime of hydrodynamics) 666The
argument of why the system at t=0 is exactly described by the constitutive
relations of ideal hydrodynamics is as follows. Rotational symmetry and
vanishing of the initial velocities imply that the stress tensor is diagonal
(in cartesian coordinates) at t=0, and that $T_{xx}=T_{yy}$. Conformal
symmetry implies that the trace is vanishing, $-T_{tt}+T_{xx}+T_{yy}=0$,
obtaining $T_{xx}=T_{tt}/2$, which is the equation of state (1). Thus, eq.
(3a) with $\Pi^{\mu\nu}=0$ exactly describes the microscopic holographic
tensor at t=0. . This aspect may partially survive through the far from
equilibrium region explaining why ideal hydrodynamics provides a better
description of the microscopic system than viscous hydrodynamics around the
hydrodynamization time. We think that this might be a particularity of the
specific system under consideration, and we wonder if this could be a generic
feature of systems sharing similar initial conditions.

See Appendix C for more details and further examples.

##  5 Discussion

We have used holography to obtain, from first principles, the real-time
quantum dynamics of a large-$N$ strongly coupled conformal field theory
initially far from equilibrium which relaxes to a hydrodynamic regime. This
constitutes a microscopic solution that we used to test the applicability of
hydrodynamics. This solution is a very simple toy model that captures one
important aspect of the physics of the QGP created in heavy-ion collision
experiments: the initial far-from-equilibrium conditions and subsequent
relaxation to a hydrodynamic regime.

We have considered the applicability of three theories of hydrodynamics,
namely ideal hydrodynamics, BRSSS and the newly formulated BDNK theory. For
each of these theories, we evolve the equations of motion numerically using
the microscopic solution as initial data at different times to explore the
non-linear and far-from-equilibrium regimes. We assess the applicability of
hydrodynamics by comparing both the constitutive relations, as previously done
in the literature, and the evolution of the fluid as predicted by each theory.
These comparisons are shown in Fig. 3, and we outline the main conclusions in
the following.

The predictions of the different hydrodynamic theories differ from the
microscopic solution and from each other in the far from equilibrium region,
$t\overline{T}\simeq 0.019,0.16$ in our case. This is expected since each
theory has a different a UV completion and, in this regime, the latter is
relevant. On the other hand, at late times, $t\overline{T}\simeq 0.58$ in our
case, the gradients are small and hydrodynamics provides a good description of
the system in the sense that the three theories that we have considered agree
well with the microscopic evolution and with each other.

In this article we have studied BDNK as a theory of relativistic viscous
hydrodynamics having in mind applications to the QGP. BDNK is particularly
interesting because its only evolution variables are the thermodynamic
quantities such as energy density and fluid velocities. To initialize the BDNK
equations, we have to consider data for the fundamental thermodynamic
variables and their first time derivatives in the causal frame obtained by
inverting (12) and neglecting higher order terms. This is only justified in
the hydrodynamic regime. We performed evolutions in different causal frames;
only when the system is in the hydrodynamic regime we find that the physics is
not affected by the choice of causal frame, up to second order terms. Away
from the hydrodynamic regime, evolutions carried out in different frames in
general differ.

In the QGP created in heavy-ion collisions the hydrodynamic evolutions provide
a good description of the system even if gradients are not very small, which
sometimes is known as ’unreasonable success of hydrodynamics’. In our case, we
do not find such an unreasonable success, but instead we found that
hydrodynamics applies precisely when it should. The theories that we
considered and the initial data are very different from the QGP case, so there
is a priori no reason why the success should happen also in our scenario.

One may wonder which causal formulation of viscous hydrodynamics, BRSSS or
BDNK, provides a better description of the system. However, this may not be
the right question to ask given that all theories of hydrodynamics should be
equivalent in the regime of validity of hydrodynamics. Therefore, the question
that one would really like to answer is which of these theories (if any!)
provides a better (i.e., more in accordance to the microscopic theory)
description of the system slightly outside the regime of hydrodynamics. It is
likely that the answer depends on the details of the system that one wants to
model and the initial conditions. In our case, where we have full control of
the microscopic theory, we did not find that either BRSSS or BDNK provides a
better description of the evolution outside the hydrodynamic regime. More
studies are needed to address this question.

Hydrodynamics can be defined by considering a specific well-posed theory
subject to suitable initial conditions. We may wonder if the constitutive
relations evaluated pointwise in spacetime on the solution obtained by solving
the initial value problem in hydrodynamics and those obtained from the
microscopic solution may provide a different answer regarding the regime of
applicability of hydrodynamics. This is an important question since in
practical applications in general one does not have a theoretical control of
the microscopic theory, and our solutions provide a concrete example where
this question can be addressed. Above we have defined the hydrodynamization
time using the second approach, obtaining $t\overline{T}\simeq 0.5$.
Alternatively, we could define the hydrodynamization time using the first
approach, defined as the time at which, if we initialize the hydrodynamic code
with microscopic data, the stress tensors along the rest of the evolution
differ less than a $2\%$ with the holographic one. For this case we obtain
hydrodynamization times compatible with $t\overline{T}\simeq 0.5$. Thus, both
approaches provide a compatible answer. Having two different theories of
relativistic viscous hydrodynamics, i.e., BDNK and BRSSS, allows us to use a
third criterium, which might be useful in practical applications where the
microscopic theory cannot be solved. If initialising both theories at a given
time using suitable initial conditions, the evolutions differ by less than
some prescribed small amount at all times, then the system is in the
hydrodynamic regime. The caveat is that one has to be careful to ensure that
the initial data for both theories is the “same”.

With this work we initiate a program to study evolutions of the BDNK equations
and to test the applicability of causal viscous hydrodynamic theories by
comparing with microscopic holographic evolutions. Possible extensions include
non-conformal and charged theories and initial data that models heavy-ion
collisions. We hope that our work, together with Pandya:2021ief , provide the
first steps towards the implementation of the BDNK equations to describe
relevant physical systems like the QGP or neutron star mergers.

### Acknowledgements

We thank Anton Faedo and Jorge Noronha for useful discussions. Y.B. and P.F.
are financially supported by the European Research Council grant ERC-2014-StG
639022-NewNGR. Y.B. is also supported by the Academy of Finland grant no.
333609. P.F. is also supported by a Royal Society University Research
Fellowship (Grant No. UF140319 and URF\R\201026). We thank the MareNostrum
supercomputer at the BSC (activity Id FI-2020-1-0007) for significant
computational resources.

##  Appendix A Implementation of hydrodynamic codes and tests

In our hydrodynamic codes we work directly with the primitive variables,
namely energy density and velocities of the fluid, and the viscosity tensor in
the case of BRSSS. This implies that our code cannot deal with shocks. That
said, we point out that for the class of initial conditions that we have
considered in this article, we did not observe the formation of shocks or
steep features in the fluid flows that would require the use high resolution
shock capturing techniques.

For the BRSSS case, we follow Green:2013zba  and we impose the constraints
$u^{\mu}u_{\mu}=-1$, $u^{\mu}\Pi_{\mu\nu}=0$ and
$\Pi^{\mu}_{\phantom{\mu}\mu}=0$ algebraically to reduce the number of
dynamical variables to 5, which we take to be
$\mathcal{U}=\\{\rho,u_{x},u_{y},\Pi_{xx},\Pi_{xy}\\}$. Then we use
Mathematica to solve the hydrodynamics equations of motion in terms of the
time derivatives of $\mathcal{U}$ and write the them as

$\partial_{t}\mathcal{U}=\mathcal{F}(\mathcal{U},\partial_{i}\mathcal{U})\,.$
(16)

The ideal case can be recovered from the BRSSS case by setting to zero all the
transport coefficients and $\Pi_{xx}=\Pi_{xy}=0$, as well as their
derivatives. In our code, we discretize the spatial derivatives in (16) using
6th order finite difference stencils and we integrate them forward in time
using a standard RK4 time-integrator. For simplicity, we imposed periodic
boundary conditions; therefore, in our hydrodynamic simulations we have to
choose a large enough domain to avoid boundary effects during the duration of
the simulations. In practice this is not a problem since the hydrodynamic
simulations are very cheap.

The equations of motion for BDNK are second order in time and space. In this
case, we implement the constraint $u^{\mu}u_{\mu}=-1$ algebraically to reduce
the number of independent variables to 3; we choose
$\mathcal{U}=\\{\rho,u_{x},u_{y}\\}$. Again, we use Mathematica to solve for
the second time derivatives of the dynamical variables, and re-write the
equations of motion as

$\partial_{t}^{2}\mathcal{U}=\mathcal{G}(\mathcal{U},\partial_{t}\mathcal{U},\partial_{i}\mathcal{U},\partial^{2}_{t,i}\mathcal{U},\partial^{2}_{i,j}\mathcal{U})\,.$
(17)

One could write (17) as a first order in time system of equations in the
obvious way, e.g., $\mathcal{P}\equiv\partial_{t}\mathcal{U}$, and use a
standard integrator such as RK4. However, we found that the resulting system
was numerically unstable; it is not clear to us what is the origin of this
instability. We did not attempt to write (17) as a fully first order system by
further defining $\mathcal{V}_{i}\equiv\partial_{i}\mathcal{U}$.

As in the BRSSS and Ideal cases, in BDNK we discretise the spatial derivatives
using 6th order finite differences and we impose periodic boundary conditions.
To proceed, we implemented two different time integrators.

First, we implemented an implicit second order in time scheme by discretizing
the time derivatives of a given variable $f$ at the time $t^{n}=n\,\Delta t$
as:

$\displaystyle\partial_{t}\mathcal{U}$
$\displaystyle\to\frac{\mathcal{U}^{n+1}_{i}-\mathcal{U}^{n-1}_{i}}{2\,\Delta
t}\,,$ (18) $\displaystyle\partial^{2}_{t}\mathcal{U}$
$\displaystyle\to\frac{\mathcal{U}^{n+1}_{i}-2\,\mathcal{U}^{n}_{i}+\mathcal{U}^{n-1}_{i}}{\Delta
t^{2}}\,,$

where $i$ denotes a collection of indices that labels a given spatial grid
point. In this way, for known $\mathcal{U}^{n}_{i}$ and
$\mathcal{U}^{n-1}_{i}$, the discrete equations of motion at time $t^{n}$
become a non-linear algebraic system for the values of the variables at the
next time level, $\mathcal{U}^{n+1}_{i}$, which can be solved using standard
techniques; in our case, for the BDNK equations we used a Newton-Raphson
algorithm. This method is robust and it works well in practice, and it has
been successfully used in simulations of black hole binary mergers
Pretorius:2004jg  and in our holographic simulations Bantilan:2012vu ;
Bantilan:2020pay ; Bantilan:2020xas . However, for high spatial resolutions,
it becomes considerably slower and memory-demanding than an explicit time
integrator such as RK4. Therefore, for practical applications that one could
easily run on a laptop, such as the ones presented here, it would be desirable
to have a stable explicit time integrator for BDNK. We now turn to this.

We have successfully implemented an explicit time integrator for the BDNK
equations of motion, written as second order in time, i.e., (17). The
analogues of the Runge-Kutta methods for second order equations are known as
Runge-Kutta-Nyström Generalized (RKNG) methods. We will review them here since
they may not be as well-known as the standard first order methods. Here we
follow Fine:1987ul  since we have implemented their RKNG34 scheme. This method
is competitive in terms of both accuracy and efficiency with respect to RK
methods of similar orders, see Fine:1987ul  for detailed comparisons.

Consider a system of $N$ second order ordinary differential equations:

$y^{\prime\prime}(x)=f(x,y(x),y^{\prime}(x))\,,$ (19)

on an interval $I=[x_{0},x_{F}]$ with initial conditions

$y(x_{0})=y_{0}\,,\quad y^{\prime}(x_{0})=y^{\prime}_{0}\,.$ (20)

The generalization of this method to PDEs is straightforward, just as in the
standard RK4 case. An explicit RKNG method of $s$ stages allows to compute the
approximations $y_{n+1}$ and $y_{n+1}^{\prime}$ of the solution $y(x_{n+1})$
and its derivative $y^{\prime}(x_{n+1})$ at $x_{n+1}\in I$ from their values
in the previous steps, $y_{n}$ and $y_{n}^{\prime}$, as follows:

$\displaystyle y_{n+1}$
$\displaystyle=y_{n}+h\,y_{n}^{\prime}+h^{2}\sum_{i=1}^{s}b_{i}\,f_{i}\,,$
(21) $\displaystyle y^{\prime}_{n+1}$
$\displaystyle=y^{\prime}_{n}+h\sum_{i=1}^{s}b^{\prime}_{i}\,f_{i}\,,$ (22)

where

$\displaystyle f_{1}$ $\displaystyle=f(x_{n},y_{n},y^{\prime}_{n})\,,$ (23)
$\displaystyle f_{i}$
$\displaystyle=f\bigg{(}x_{n}+c_{i}\,h,\,y_{n}+c_{i}\,h\,y^{\prime}_{n}+h^{2}\sum_{j=1}^{i-1}a_{ij}\,f_{j},$
$\displaystyle\hskip
28.45274pty_{n}^{\prime}+h\sum_{j=1}^{i-1}a^{\prime}_{ij}f_{j}\bigg{)}\,,\quad
i=2,\ldots,s\,,$ (24)

and $h=x_{n+1}-x_{n}$. Here $c_{i}$, $b_{i}$ and $b_{i}^{\prime}$ are
$s$-dimensional constant vectors, and $a_{ij}$ and $a^{\prime}_{ij}$ are lower
triangular $s\times s$ constant matrices; together, they constitute the
parameters that define the method. An RKGN method is said to have order $p$ if
the local convergence order for both the approximation $y_{n+1}$ that of the
derivative $y^{\prime}_{n+1}$ is $p$,

$\displaystyle y_{n+1}-y(x_{n+1})$
$\displaystyle=T_{p+1}\,h^{p+1}+O(h^{p+2})\,,$ (25) $\displaystyle
y^{\prime}_{n+1}-y^{\prime}(x_{n+1})$
$\displaystyle=T^{\prime}_{p+1}\,h^{p+1}+O(h^{p+2})\,,$ (26)

where $T_{p+1}$ and $T^{\prime}_{q+1}$ are the principal error functions of
the solution and the derivative respectively and $h$ is the grid spacing. See
Fine:1987ul  for more details. As in other multistep methods, such as standard
Runge-Kutta methods, the global truncation error, i.e., the accumulated error,
is of order $O(h^{p})$.

Ref. Fine:1987ul  gives a detailed derivation of a RKNG45 and a RKNG34
methods; here we only give the values of the parameters in (21)–(24) of the
5-stages RKNG34 scheme, which is a fourth order method, since this is what we
implemented in our BDNK code as this method is competitive with the standard
RK4 method. For the RKNG34 method of Fine:1987ul , one has:

$\displaystyle b_{i}$
$\displaystyle=\left(\tfrac{19}{180},0,\tfrac{63}{200},\tfrac{16}{225},\tfrac{1}{120}\right)\,,$
(27) $\displaystyle b^{\prime}_{i}$
$\displaystyle=\left(\tfrac{1}{9},0,\tfrac{9}{20},\tfrac{16}{45},\tfrac{1}{12}\right)\,,$
(28) $\displaystyle c_{i}$
$\displaystyle=\left(0,\tfrac{2}{9},\tfrac{1}{3},\tfrac{3}{4},1\right)\,,$
(29)

with

$\displaystyle a_{21}$ $\displaystyle=\tfrac{2}{81}\,,$ $\displaystyle a_{31}$
$\displaystyle=\tfrac{1}{36}\,,\,\,a_{32}=\tfrac{1}{36}\,,$ (30)
$\displaystyle a_{41}$
$\displaystyle=\tfrac{9}{128}\,,\,\,a_{42}=0\,,\,\,a_{43}=\tfrac{27}{128}\,,$
$\displaystyle a_{51}$
$\displaystyle=\tfrac{11}{60}\,,\,\,a_{52}=-\tfrac{3}{20}\,,\,\,a_{53}=\tfrac{9}{25}\,,\,\,a_{54}=\tfrac{8}{75}\,,$

and

$\displaystyle a^{\prime}_{21}$ $\displaystyle=\tfrac{2}{9}\,,$ $\displaystyle
a^{\prime}_{31}$
$\displaystyle=\tfrac{1}{12}\,,\,\,a^{\prime}_{32}=\tfrac{1}{4}\,,$ (31)
$\displaystyle a^{\prime}_{41}$
$\displaystyle=\tfrac{69}{128}\,,\,\,a^{\prime}_{42}=-\tfrac{243}{128}\,,\,\,a^{\prime}_{43}=\tfrac{135}{64}\,,$
$\displaystyle a^{\prime}_{51}$
$\displaystyle=-\tfrac{17}{12}\,,\,\,a^{\prime}_{52}=-\tfrac{27}{4}\,,\,\,a^{\prime}_{53}=-\tfrac{27}{5}\,,\,\,a^{\prime}_{54}=\tfrac{16}{15}\,,$

and the remaining components of $a_{ij}$ and $a^{\prime}_{ij}$ all vanish.

We study the convergence of our BDNK, BRSSS and ideal hydrodynamics codes by
considering the evolution of a Gaussian initial profile for the energy density
of the form,

$\rho=\rho_{0}\bigg{(}1+\tfrac{A}{R\sqrt{2\pi}}\,e^{-\frac{(r-r_{0})^{2}}{2R^{2}}}\bigg{)}\,,$
(32)

centred at $r=r_{0}$, with width $R$ and amplitude $A$, on top of a background
energy density $\rho_{0}$. The initial values for all the remaining variables
are chosen to be zero. For the runs shown below, we chose a computational
domain of size $L=10$, $\rho_{0}=0.5$, $A=0.5$, $R=1$ and we centred the
Gaussian at the centre of our grid. For the BDNK simulations, we choose a
frame with $a_{1}=a_{2}=10$. We monitor the convergence rate $Q_{N}(t)$,

$Q_{N}(t)=\frac{1}{\ln
2}\,\ln\left(\frac{||u_{4h}(t)-u_{2h}(t)||_{1}}{||u_{2h}(t)-u_{h}(t)||_{1}}\right)\,,$
(33)

where $||\cdot||_{1}$ denotes the numerical $L^{1}$-norm taken over our
computational domain, and $h$ is the grid spacing. To compute these norms, we
sum over all the evolution variables.

Figure 4: Convergence tests for our ideal (top), BRSSS (middle) and BDNK
(bottom) codes. After an initial transient, the convergence rate approaches
(and even surpasses) 4, as one would expect given our differencing and time
integration schemes.

The results of our convergence tests are displayed in Fig. 4 for the three
theories that we consider in this article. As this figure shows, the
convergence rate is close to 4 in all cases. Even though in our code we use
6th order spatial differences, this figure shows that the dominant error comes
from the time integration. This justifies the use of 6th order Kreiss-Oliger
dissipation. Furthermore, the convergence rate that we find is consistent with
the order of the time integrators that we used in our code, namely RKNG34 for
BDNK and RK4 for BRSSS and ideal hydrodynamics.

##  Appendix B Numerical relativity details

In this appendix we summarize the basic aspects of our numerical relativity
simulations. We use the same code as in Bantilan:2020pay ; Bantilan:2020xas ,
so more details can be found in these references.

We solve the Einstein equations in AdS coupled to a massless scalar field in
3+1 dimensions using generalized harmonic coordinates,

$\displaystyle-\frac{1}{2}g^{\gamma\delta}g_{\alpha\beta,\gamma\delta}-{g^{\gamma\delta}}_{,(\alpha}g_{\beta)\gamma,\delta}-H_{(\alpha,\beta)}+H_{\gamma}{\Gamma^{\gamma}}_{\alpha\beta}$
$\displaystyle-{\Gamma^{\gamma}}_{\delta\alpha}{\Gamma^{\delta}}_{\gamma\beta}-\kappa\left(2n_{(\alpha}C_{\beta)}-(1+P)g_{\alpha\beta}n^{\gamma}C_{\gamma}\right)$
$\displaystyle=\frac{2}{d-2}\Lambda
g_{\alpha\beta}+8\pi\left(T_{\alpha\beta}-\frac{1}{d-2}{T^{\gamma}}_{\gamma}g_{\alpha\beta}\right),$

where $n_{\alpha}=-\partial_{\alpha}t$ is the timelike, future-directed unit
1-form normal to slices of constant $t$, and $\kappa$ and $P$ are the damping
parameters. Here $C^{\alpha}\equiv H^{\alpha}-\Box x^{\alpha}$ are the
constraints and $H_{\alpha}$ are the source functions, which we can freely
prescribe. We choose both the parameters and the source functions as in
Bantilan:2020pay ; Bantilan:2020xas .

For convenience, we couple gravity in AdS to a massless scalar field, since
this allows us to easily construct far-from-equilibrium initial data by
prescribing a suitable profile for the scalar field. Therefore, in addition to
the metric, we also evolve a massless scalar field, with equation of motion

$\Box\phi=0\,,$ (34)

and stress tensor

$T_{\alpha\beta}=\partial_{\alpha}\phi\partial_{\beta}\phi-
g_{\alpha\beta}\frac{1}{2}g^{\gamma\delta}\partial_{\gamma}\phi\partial_{\delta}\phi\,.$
(35)

To carry out the evolution, we write the general spacetime metric
$g_{\mu\nu}=\hat{g}_{\mu\nu}+\bar{g}_{\mu\nu}$, where $\hat{g}_{\mu\nu}$ is
the metric of the Poincaré patch of AdS4 (with the AdS length $L=1$) written
as

$\hat{g}=\frac{\rho^{4}}{(1-\rho^{2})^{2}}\left(-dt^{2}+\tfrac{4}{\rho^{6}}\,d\rho^{2}+dx^{2}+dy^{2}\right),$
(36)

where $z=(1-\rho^{2})/\rho^{2}$ is the usual AdS radial coordinate, and $x$,
$y$ are the spatial boundary directions, which we take to be Cartesian
directions with infinite range. In practice we compactify these directions as
$x=\tan(\frac{\pi}{2}\bar{x}_{1})$, $y=\tan(\frac{\pi}{2}\bar{x}_{2})$.
$\bar{g}_{\mu\nu}$ is a general deviation away from AdS4, not necessarily
small, and it satisfies Dirichlet boundary conditions at the AdS boundary.
Similarly, we write the scalar field as

$\phi=(1-\rho^{2})^{2}\bar{\phi}\,,$ (37)

and the evolved variable $\bar{\phi}$ vanishes at the boundary of AdS4.

Following Bantilan:2020pay ; Bantilan:2020xas , we write the source functions
as

$H_{\mu}=\hat{H}_{\mu}+(1-\rho^{2})\bar{H}_{\mu}\,,$ (38)

where $\hat{H}_{\mu}$ are the source functions for AdS4 in the coordinates of
equation (36) and $\bar{H}_{\mu}$ are the actual evolved source functions; the
power of $(1-\rho^{2})$ in (38) has been chosen so that
$\bar{H}_{\mu}\big{|}_{\rho=1}=0$ at the boundary of AdS. Following the
prescription of Bantilan:2020xas , we fix our gauge such that near the
boundary of AdS, we have

$\displaystyle\bar{H}^{(1)}_{t}$
$\displaystyle=\frac{3}{2}\bar{g}^{(1)}_{\text{$t$$\rho$}}\,,$
$\displaystyle\bar{H}^{(1)}_{\rho}$
$\displaystyle=\frac{3}{2}\bar{g}^{(1)}_{\rho\rho}\,,$
$\displaystyle\bar{H}^{(1)}_{x}$ $\displaystyle=\frac{3}{2}\bar{g}^{(1)}_{\rho
x}\,,$ $\displaystyle\bar{H}^{(1)}_{y}$
$\displaystyle=\frac{3}{2}\bar{g}^{(1)}_{\rho y}\,,$

where

$\displaystyle\bar{H}_{\mu}$
$\displaystyle=(1-\rho^{2})\bar{H}^{(1)}_{\mu}+O((1-\rho^{2})^{2})\,,$ (39)
$\displaystyle\bar{g}_{\mu\nu}$
$\displaystyle=(1-\rho^{2})\bar{g}^{(1)}_{\mu\nu}+O((1-\rho^{2})^{2})\,,$

near the boundary. The source functions are set to zero in the interior of the
spacetime. We choose the same smooth transition functions as Bantilan:2020xas
to fix the source functions everywhere in the spacetime.

We construct time symmetric initial data by prescribing a Gaussian profile for
the initial scalar field,

$\phi=(1-\rho^{2})^{3}\,A\,e^{-\tilde{r}(\rho,\bar{x}_{1},\bar{x}_{2})^{2}/\Delta^{2}},$
(40)

where $A$ and $\Delta$ control the amplitude and the width of the Gaussian
respectively, and

$\tilde{r}(\rho,x,y)^{2}=\tfrac{(\rho-\rho_{0})^{2}}{a_{\rho}^{2}}+\tfrac{\bar{x}_{1}}{\alpha_{1}^{2}}+\tfrac{\bar{x}_{2}}{\alpha_{2}^{2}}\,.$
(41)

Here $\rho_{0}$ allows us to localize the initial Gaussian in the AdS radial
direction, while the constant parameters
$\\{a_{\rho},\alpha_{1},\alpha_{2}\\}$ control the shape along the boundary
directions. Having specified an initial profile for the scalar field, the
Hamiltonian constraint is solved as in Bantilan:2020xas . We choose “strong"
data so that there is a trapped surface, with planar topology, in the initial
data slice.

We extract the stress energy tensor of the boundary CFT as in Bantilan:2020xas
. We find,

$\displaystyle\langle T_{tt}\rangle$
$\displaystyle=\frac{1}{32\,\pi}\left(6\,\bar{g}^{(1)}_{xx}+6\,g^{(2)}_{yy}+\bar{g}^{(1)}_{\rho\rho}\right)\,,$
(42) $\displaystyle\langle T_{tx}\rangle$
$\displaystyle=\frac{3}{16\,\pi}\,\bar{g}^{(1)}_{tx}\,,$ (43)
$\displaystyle\langle T_{ty}\rangle$
$\displaystyle=\frac{3}{16\,\pi}\,\bar{g}^{(1)}_{ty}\,,$ (44)
$\displaystyle\langle T_{xx}\rangle$
$\displaystyle=-\frac{1}{32\,\pi}\left(-6\,\bar{g}^{(1)}_{tt}+6\,g^{(2)}_{yy}+\bar{g}^{(1)}_{\rho\rho}\right)\,,$
(45) $\displaystyle\langle T_{xy}\rangle$
$\displaystyle=\frac{3}{16\,\pi}\,\bar{g}^{(1)}_{xy}\,,$ (46)
$\displaystyle\langle T_{yy}\rangle$
$\displaystyle=-\frac{1}{32\,\pi}\left(-6\,\bar{g}^{(1)}_{tt}+6\,g^{(2)}_{xx}+\bar{g}^{(1)}_{\rho\rho}\right)\,.$
(47)

Fig. 5 shows a convergence test for the trace of the boundary stress tensor of
the evolution presented in the main text, which shows first order convergence
with the number of points along the holographic direction $N_{\rho}$.

Figure 5: Convergence test for the trace. (Top) Trace at the center of the
domain $\\{x,y\\}\overline{T}\simeq\\{0,0\\}$ as a function of time, for three
different resolutions along the holographic direction $N_{\rho}=257,513,1025$.
The data for the $N_{\rho}=257$ case is only available up to
$t\overline{T}\simeq 0.17$. (Bottom) Trace as a function of $1/N_{\rho}$ at a
fixed time $t\overline{T}\simeq 0.14$, at the center of the domain. The dashed
line is a straight line fit that shows first order convergence.

A convergence test of the conservation equation of the stress tensor for our
code was presented in Bantilan:2012vu , and we do not perform this convergence
test in our case.

Figure 6: Temporal component of the stress tensor conservation equation,
$\partial_{\mu}T^{\mu t}$, at the center of the domain as a function of time.

Fig. 6 shows the temporal component of the conservation equation for the
stress tensor, $\partial_{\mu}T^{\mu t}$, for our highest resolution case
$N_{\rho}\times N_{x}\times N_{y}=1025\times 129\times 129$. We present the
result at the center of the domain because the numerical error is maximal
there. For the other components, $\partial_{\mu}T^{\mu x}$,
$\partial_{\mu}T^{\mu y}$, the error is smaller.

We use the numerical error of the conservation of the stress tensor, Fig. 6,
to provide an estimate of the numerical error in the first order terms of the
hydrodynamic expansion. Schematically, the first order terms are given by the
shear viscosity times derivatives of the velocity, so we estimate the error as

$\frac{T^{1st}_{xx}}{T^{id}_{xx}}\sim\frac{\eta\,\partial_{x}v_{x}}{\overline{\mathcal{E}}/2}\sim\frac{\eta\,\partial_{\mu}T^{\mu
t}/(2\overline{\mathcal{E}})}{\overline{\mathcal{E}}/2}\,.$ (48)

We estimate $T^{id}_{xx}\sim\overline{\mathcal{E}}/2$ and
$\partial_{x}v_{x}\sim\partial_{\mu}T^{\mu t}/(2\overline{\mathcal{E}})$. In
the last expression, by dividing by $\overline{\mathcal{E}}$ we obtain the
correct units, and the factor of two is because at the center the
$\partial_{x}T^{xt}$ and $\partial_{y}T^{yt}$ terms contribute equally and
from the actual data we check that the $\partial_{t}T^{tt}$ term is much
smaller. For times around the hydrodynamization time $t\overline{T}\simeq
0.5$, the right hand side of (48) gives an error of the order

$\frac{T^{1st}_{xx}}{T^{id}_{xx}}\sim 2\%\,.$ (49)

A numerical error of the order of $\sim 2\%$ is of the same order of the first
and second order terms around times $t\overline{T}\simeq 0.5$, see Fig. 2. So,
all the statements in this paper about the applicability of hydrodynamics for
quantities beyond this time must be considered only up to numerical errors.

The change from the Landau frame to the causal frame is specially affected by
the error in the conservation equation of the stress tensor. This is because
the first order terms in the change of frame expression are precisely the same
terms as in the ideal conservation equation of the stress tensor. For this
reason, we subtract by hand the error of the conservation of the stress tensor
in the change of frame to obtain the initial data for the BDNK evolutions.

##  Appendix C Details on the dynamical evolutions

In this appendix we provide more details about the dynamical evolutions
presented in the main text.

Fig. 2 shows the ratios $T_{xx}^{n}/T_{xx}^{ideal}$, $n=1,2$, of the
constitutive relations (3) evaluated at
$\\{x,y\\}\overline{T}\simeq\\{0.12,0\\}$ as a function of time for the
holographic solution. This shows explicitly the validity of the hydrodynamic
expansion. Alternatively, the constitutive relations can also be compared with
the microscopic stress tensor values, to assess if the constitutive relations
describe the microscopic evolution. This is shown in Fig. 7. On general
grounds, we would expect that if the system is in the hydrodynamic regime,
then the gradient corrections should improve the ideal description. However,
in our system, at the hydrodynamization time $t\overline{T}\simeq 0.5$ first
and second order terms do not improve the ideal description of the microscopic
solution. As stated in the main text, this might be explained by our choice
initial data. We expect that at later times viscous descriptions would improve
the ideal result. However, at later times the numerical error (discussed in
appendix B) does not allow us to make conclusive statements in this regard.
Also, possibly for the same reason, ideal terms provide a good description of
the system even before the hydrodynamization time $t\overline{T}\simeq 0.5$.
The microscopic $T_{xx}$ presented in Fig. 7 presented a numerical oscillation
of a cero mode (homogeneous). This numerical oscillation converges to zero
similarly to Fig. 5, and we removed it by hand in Fig. 7.

Figure 7: $T_{xx}$ component of the holographic stress tensor at
$\\{x,y\\}\overline{T}\simeq\\{0.12,0\\}$ as a function of time, solid black
line. We include the values of $T_{xx}$ obtained by evaluating the
constitutive relations (3) up to ideal terms (dashed green), first order terms
(dotted blue) and second order terms (dashed red).

Figure 8: Energy density in the lab frame, $T_{tt}$, at the center of the
domain $\\{x,y\\}\overline{T}=\\{0,0\\}$ as a function of time. Black
continuous line corresponds to the microscopic holographic solution. We
include the results of the hydrodynamic evolutions initialised with
holographic data at $t\overline{T}\simeq 0.019,0.16,0.58$, from top to bottom.
For the BDNK evolutions, we include the result of using only the 0th order
terms (the ideal part) in (10), in dotted lines. We also include the
evolutions of BDNK in the frame 1.

Fig. 8 shows the same evolutions as in Fig. 3 but at an off center location
$\\{x,y\\}\overline{T}=\\{0.17,0\\}$. The conclusions are similar to those
obtained from Fig. 3. If the gradients are large (top and middle), the
hydrodynamic evolutions provide a description that differs from the
microscopic one, and are also different among them, as all theories have
different UV behaviors. At late times (bottom), when gradients are small, the
hydrodynamic theories provide a better description of the system, matching the
microscopic theory within a maximum deviation in the domain of a $\sim 2\%$ in
the energy density.

Fig. 8 also includes the evolutions of the BDNK equations in frame 1, and we
discuss now the differences of evolving the BDNK equations in the different
causal frames (14). First, recall that $a_{1}$ and $a_{2}$ are multiplying
first order terms in (10). This means that the larger these values, the larger
are those first order terms. So, for frame 1 these terms are larger than for
frame 2. This is the reason why in the main text we focus on the frame 2, for
which the values for $a_{1}$ and $a_{2}$ are closer to the smallest ones
allowed by causality. The size of the constants $a_{1}$ and $a_{2}$ will be
relevant when changing from the Landau frame to the causal frame, when
evolving the equations, and also when going back to the Landau frame. We
observe that the size of these first order terms (the ones proportional to
$a_{1}$ and $a_{2}$) in the microscopic solution is typically small at all
times, even far from equilibrium. On the other hand, in the BDNK simulations
in the far from equilibrium regime, these terms become large, even of order
one, see Fig. 8 (top and middle). Moreover, the BDNK evolutions in the
different frames in the far from equilibrium region may provide considerably
different results, see Fig. 8 (middle). Finally, in Fig. 8 (bottom) we see
that the difference between frame 1 and 2 is compatible with second order
terms, in accordance with the fact that the system has hydrodynamized.

## References

  * (1) P. Romatschke and U. Romatschke, _Relativistic Fluid Dynamics In and Out of Equilibrium_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (5, 2019), 10.1017/9781108651998, [1712.05815].
  * (2) W. Busza, K. Rajagopal and W. van der Schee, _Heavy Ion Collisions: The Big Picture, and the Big Questions_ , _Ann. Rev. Nucl. Part. Sci._ 68 (2018) 339 [1802.04801].
  * (3) M.G. Alford, L. Bovard, M. Hanauske, L. Rezzolla and K. Schwenzer, _Viscous Dissipation and Heat Conduction in Binary Neutron-Star Mergers_ , _Phys. Rev. Lett._ 120 (2018) 041101 [1707.09475].
  * (4) M. Chabanov, L. Rezzolla and D.H. Rischke, _General-relativistic hydrodynamics of non-perfect fluids: 3+1 conservative formulation and application to viscous black-hole accretion_ , 2102.10419.
  * (5) M. Shibata, K. Kiuchi and Y.-i. Sekiguchi, _General relativistic viscous hydrodynamics of differentially rotating neutron stars_ , _Phys. Rev. D_ 95 (2017) 083005 [1703.10303].
  * (6) S. Fujibayashi, K. Kiuchi, N. Nishimura, Y. Sekiguchi and M. Shibata, _Mass Ejection from the Remnant of a Binary Neutron Star Merger: Viscous-Radiation Hydrodynamics Study_ , _Astrophys. J._ 860 (2018) 64 [1711.02093].
  * (7) S. Fujibayashi, M. Shibata, S. Wanajo, K. Kiuchi, K. Kyutoku and Y. Sekiguchi, _Viscous evolution of a massive disk surrounding stellar-mass black holes in full general relativity_ , _Phys. Rev. D_ 102 (2020) 123014 [2009.03895].
  * (8) W.A. Hiscock and L. Lindblom, _Generic instabilities in first-order dissipative relativistic fluid theories_ , _Phys. Rev. D_ 31 (1985) 725.
  * (9) I. Muller, _Zum Paradoxon der Warmeleitungstheorie_ , _Z. Phys._ 198 (1967) 329.
  * (10) W. Israel and J.M. Stewart, _Transient relativistic thermodynamics and kinetic theory_ , _Annals Phys._ 118 (1979) 341.
  * (11) W. Israel, _Nonstationary irreversible thermodynamics: A Causal relativistic theory_ , _Annals Phys._ 100 (1976) 310.
  * (12) R. Baier, P. Romatschke, D.T. Son, A.O. Starinets and M.A. Stephanov, _Relativistic viscous hydrodynamics, conformal invariance, and holography_ , _JHEP_ 04 (2008) 100 [0712.2451].
  * (13) G.S. Denicol, H. Niemi, E. Molnar and D.H. Rischke, _Derivation of transient relativistic fluid dynamics from the Boltzmann equation_ , _Phys. Rev. D_ 85 (2012) 114047 [1202.4551].
  * (14) R.P. Geroch and L. Lindblom, _Dissipative relativistic fluid theories of divergence type_ , _Phys. Rev. D_ 41 (1990) 1855.
  * (15) L. Lehner, O.A. Reula and M.E. Rubio, _Hyperbolic theory of relativistic conformal dissipative fluids_ , _Phys. Rev. D_ 97 (2018) 024013 [1710.08033].
  * (16) R.P. Geroch, _Relativistic theories of dissipative fluids_ , _J. Math. Phys._ 36 (1995) 4226.
  * (17) R.P. Geroch, _On hyperbolic ’theories’ of relativistic dissipative fluids_ , gr-qc/0103112.
  * (18) F.S. Bemfica, M.M. Disconzi and J. Noronha, _Causality and existence of solutions of relativistic viscous fluid dynamics with gravity_ , _Phys. Rev. D_ 98 (2018) 104064 [1708.06255].
  * (19) F.S. Bemfica, M.M. Disconzi and J. Noronha, _Nonlinear Causality of General First-Order Relativistic Viscous Hydrodynamics_ , _Phys. Rev. D_ 100 (2019) 104020 [1907.12695].
  * (20) P. Kovtun, _First-order relativistic hydrodynamics is stable_ , _JHEP_ 10 (2019) 034 [1907.08191].
  * (21) J. Armas and F. Camilloni, _A stable and causal model of magnetohydrodynamics_ , 2201.06847.
  * (22) A. Pandya and F. Pretorius, _Numerical exploration of first-order relativistic hydrodynamics_ , _Phys. Rev. D_ 104 (2021) 023015 [2104.00804].
  * (23) M.P. Heller, R.A. Janik, M. Spaliński and P. Witaszczyk, _Coupling hydrodynamics to nonequilibrium degrees of freedom in strongly interacting quark-gluon plasma_ , _Phys. Rev. Lett._ 113 (2014) 261601 [1409.5087].
  * (24) A. Pandya, E.R. Most and F. Pretorius, _Conservative finite volume scheme for first-order viscous relativistic hydrodynamics_ , 2201.12317.
  * (25) F.S. Bemfica, M.M. Disconzi, C. Rodriguez and Y. Shao, _Local well-posedness in Sobolev spaces for first-order conformal causal relativistic viscous hydrodynamics_ , 1911.02504.
  * (26) C.P. Herzog, _The Hydrodynamics of M theory_ , _JHEP_ 12 (2002) 026 [hep-th/0210126].
  * (27) M. Natsuume and T. Okamura, _Causal hydrodynamics of gauge theory plasmas from AdS/CFT duality_ , _Phys. Rev. D_ 77 (2008) 066014 [0712.2916].
  * (28) M. Van Raamsdonk, _Black Hole Dynamics From Atmospheric Science_ , _JHEP_ 05 (2008) 106 [0802.3224].
  * (29) O. Aharony, O. Bergman, D.L. Jafferis and J. Maldacena, _N=6 superconformal Chern-Simons-matter theories, M2-branes and their gravity duals_ , _JHEP_ 10 (2008) 091 [0806.1218].
  * (30) H. Bantilan, P. Figueras and D. Mateos, _Real-time Dynamics of Plasma Balls from Holography_ , _Phys. Rev. Lett._ 124 (2020) 191601 [2001.05476].
  * (31) H. Bantilan, P. Figueras and L. Rossi, _Cauchy Evolution of Asymptotically Global AdS Spacetimes with No Symmetries_ , _Phys. Rev. D_ 103 (2021) 086006 [2011.12970].
  * (32) S.R. Green, F. Carrasco and L. Lehner, _Holographic Path to the Turbulent Side of Gravity_ , _Phys. Rev. X_ 4 (2014) 011001 [1309.7940].
  * (33) F. Pretorius, _Numerical Relativity Using a Generalized Harmonic Decomposition_ , _Class. Quant. Grav._ 22 (2005) 425.
  * (34) H. Bantilan, F. Pretorius and S.S. Gubser, _Simulation of Asymptotically AdS5 Spacetimes with a Generalized Harmonic Evolution Scheme_ , _Phys. Rev. D_ 85 (2012) 084038 [1201.2132].
  * (35) J.M. Fine, _Low order practical runge-kutta-nyström methods_ , _Computing_ 38 (1987) 281.


# The effect of three matters on KSS bound

Mehdi Sadeghi
Department of Physics, School of Sciences,
Ayatollah Boroujerdi University, Boroujerd, Iran Email:
<EMAIL_ADDRESS>
###### Abstract

In this paper we introduce the black brane solutions in AdS space in
4-dimensional (4D) Einstein-Gauss-Bonnet-Yang-Mills theory in the presence of
string cloud and quintessence. Shear viscosity to entropy density ratio is
computed via fluid-gravity duality, as a transport coefficient for this model.

PACS numbers: 11.10.Jj, 11.10.Wx, 11.15.Pg, 11.25.Tq

Keywords: AdS/CFT duality, Fluid-Gravity duality, Shear viscosity, Green-Kubo
formula

##  1 Introduction

The AdS/CFT duality [1]-[4] originates from string theory and provides a new
perspective on quantum gravity. In its most general form, it is known as a
gauge/gravity duality. This duality relates weakly coupled gravitational
theories in $(D+1)$-dimensional AdS space-time with strongly coupled conformal
field theories (CFTs) defined on the $D$-dimensional boundary of AdS. In
particular, it gives support to study the transport coefficients from
hydrodynamics to the quark-gluon plasma formed at relativistic heavy-ion
collisions [5]-[8]. Shear viscosity $\eta$ [9]-[12] as one of these
coefficients is the most well-known coefficient calculated by this duality,
especially fluid-gravity duality[13]-[16].
The Kovtun-Son-Starinets (KSS) bound states that the ratio $\eta/s$ has a
lower bound, $\frac{\eta}{s}\geq\frac{1}{4\,\pi}$, for all relativistic
quantum field theories [11],[12] and can be interpreted as the Heisenberg
uncertainty principle [14],[17] where $\eta$ and s are the shear viscosity and
entropy density, respectively. However, this conjecture violates for some
theories like the Gauss-Bonnet gravity [18], the Gauss-Bonnet gravity in
presence of $U(1)$ gauge field [19],[20], Horndeski theory [21], massive
gravity for $c_{i}<0$ [22], scalar-tensor gravity[23]-[25] and anisotropic
black brane [26].
Observational cosmology shows the accelerated expansion of the Universe
[27],[28]. Therefore, the Universe at present is becoming filled with a
strange substance. This outcome is confirmed by the measurement of the Cosmic
Microwave Background (CMB) using PLANCK space Satellite [29]. It is believed
that this expansion is due to a substantial negative pressure arising from the
dark energy that is up to 70$\%$ of the total energy of the Universe. There
are two sources of this negative pressure: one is cosmological constant and
the other is the so called quintessence [30],[31] which acts as a repulsive
force against the gravity. Quintessence is described by a scalar field and a
parameter $\omega$, defined as the ratio of the pressure to the energy density
of the dark energy. Black hole solution surrounded by the quintessence
introduced in [41].
String theory in which the fundamental ingredient of our universe is assumed
to be one-dimensional strings, instead of particles, tries to merge the
gravitational theory and quantum mechanics. An extension of this idea is to
consider a cloud of strings [33]-[39] and study its possible measurable
effects on long range gravitational fields of some sources, shuch as black
holes. Letelier [33] has tried to extend the idea of a string to a black hole,
possibly to reveal some properties of the blac khole. Ghosh and et al.
obtained a generalization for third-order Lovelock gravity [34], and
Herscovich and et al. [35] obtained the solution for Einstein-Gauss-Bonnet
theory in the Letelier spacetime.
General relativity (GR) provides the standard description of gravity and
Gauss-Bonnet gravity is an extension gravity theory. Recently a new proposal
has been offered wherein GB term is made in $4D$ dimentions[40]. In that ,the
GB coupling is scaled as $\alpha\to\frac{\alpha}{D-4}$ and thereby canceling
out $(D-4)$ factor in the equation, and then taking the limit $(D\to 4)$. This
results into an effective equation in $4D$ which is in fact the Einstein-
Gauss-Bonnet (EGB) equation written in $D=4$. Then it could be solved in
spacetime with some specific symmetries for different situations, black holes
and cosmology.
In this paper, we consider $4D$ Einstein-Gauss-Bonnet-Yang-Mills theory
(EGBYM) in the presence of string cloud and quintessence. Then, by calculation
of $\frac{\eta}{s}$ via fluid-gravity duality we study the effects of
quintessence, string cloud and Yang-Mills charge on the field theory dual
side.

##  2 $4D$ AdS Einstein-Gauss-Bonnet-Yang-Mills in presence of string cloud
and quintessence black brane

The action of AdS Einstein-Gauss-Bonnet-Yang-Mills in the presence of string
cloud and quintessence is:

$I=\frac{1}{16\pi}\int{d^{D}x\sqrt{-g}\Bigg{[}R-2\Lambda+\frac{\alpha^{\prime}}{D-4}\mathcal{G}+F^{(a)}_{\mu\alpha}F^{(a)\mu\alpha}\Bigg{]}}+S_{\text{cs}}+S_{\text{quint}},$
(1)

where $R$ is the scalar curvature, $\Lambda=\frac{-(D-1)(D-2)}{2l^{2}}$ ,
$\alpha^{\prime}$ is a (positive) Gauss-Bonnet coupling constant with
dimension $(\text{length})^{2}$ ,
$\mathcal{G}=R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$
,
$F^{(a)}_{\mu\nu}=\partial_{\mu}A^{(a)}_{\nu}-\partial_{\nu}A^{(a)}_{\mu}-i[A^{(a)}_{\mu},A^{(a)}_{\nu}]$
is the Cartan subalgebra of $SU(2)$ Yang-Mills field strength tensor in which
the gauge coupling constant is 1, $A_{\nu}$’s are the Cartan subalgebra of the
$SU(2)$ gauge group Yang-Mills potentials.
By considering the following ansatz for the metric,

$ds^{2}=-f(r)N^{2}(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega^{2}_{2,\kappa},$
(2)

where $\kappa=-1,0,1$ and where $\Omega_{D-2}$ is the volume of unit
$(D-2)$-sphere.
$f(r)$ is defined by a new variable $\phi(r)$,

$f(r)=\kappa-r^{2}\phi(r).$ (3)

The action of cloud of string is as follows,

$S_{\text{cs}}=\int_{\Sigma}P\sqrt{-\gamma}d\tau d\sigma,$ (4)

where $P$ is non-negative constant related to string tension,
$(\tau,\sigma)=(\lambda^{0},\lambda^{1})$ is a coordiantes of worldsheet of
string with $\lambda^{0}$ and $\lambda^{1}$ are timelike and spacelike
parameters respectively, $\gamma_{ab}$ is the induced metric on the world
sheet is as follows,

$\gamma_{ab}=g_{\mu\nu}(x)\frac{\partial
x^{\mu}}{\partial\lambda^{a}}\frac{\partial x^{\nu}}{\partial\lambda^{b}},$
(5)

and $\gamma=\text{det}\gamma_{ab}$ is the determinant of the induced metric.
The bivector of string worldsheet $\Sigma$ is given by,

${\Sigma}^{\mu\nu}={\epsilon}^{ab}\frac{\partial
x^{\mu}}{\partial\lambda^{a}}\frac{\partial x^{\nu}}{\partial\lambda^{b}},$
(6)

where ${\epsilon}^{ab}$ is the Levi-Civita tensor in two dimensions,
$\epsilon^{01}=-\epsilon^{10}=1$.
Therefore, the action of cloud of string (4) can be written as,

$S_{\text{cs}}=P\int_{\Sigma}\sqrt{-\frac{1}{2}{\Sigma}^{\mu\nu}{\Sigma}_{\mu\nu}}d\tau
d\sigma,$ (7)

The non-zero components of the bivector $\Sigma$ is
$\Sigma^{tr}=-\Sigma^{rt}=-\frac{a}{Pr^{D-2}}$. Where $a$ is an integration
constant which is related to the cloud of string [33]-[39].
The action of quintessence is introduce by Kiselev [41] and written in terms
of scalar field as follows,

$S_{\text{quint}}=\frac{1}{16\pi
G}\int{d^{D}x\sqrt{-g}\Big{[}-\frac{1}{2}(\vec{\nabla}\phi)^{2}-V(\phi)\Big{]}}.$
(8)

The energy-momentum tensor of the quintessence dark energy in $D$ dimensions
can be described by [42] ,

$\displaystyle T_{t}\,^{t}$
$\displaystyle=T_{r}\,^{r}=\rho=\frac{\omega\alpha(D-1)(D-2)}{4r^{(D-1)(\omega+1)}},$
(9) $\displaystyle T_{x_{i}}\,^{x_{i}}$
$\displaystyle=\frac{\rho}{D-2}((D-1)\omega+1),\,\,i=(1,...,D-2),$ (10)

where $\omega$ ($-1<\omega<-\frac{D-3}{D-1}$) and $\alpha$ are the
quintessential state parameter and the positive normalization constant related
to the density of quintessence, respectively.
By variation of the action with respect to $A^{(a)}_{\mu}$ we have,

$\displaystyle\nabla_{\mu}F^{(a)\mu\nu}=0,$ (11)

For solving Eq.(11), we consider the ansatz for the gauge field as
follows[43],

${\bf{A}}^{(a)}=\frac{i}{2}h(r)dt\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix},$
(12)

so,

$h(r)=C_{2}+Q\int^{r}\frac{1}{u^{D-2}}du,$ (13)
$h^{\prime}(r)=\frac{Q}{r^{D-2}},$ (14)

By plugging the ansatz 2 into the action Eq.(1) yields,

$\displaystyle I=\frac{\Omega_{D-2}(D-2)}{16\pi}$
$\displaystyle\int{dtdrN(r)\Bigg{[}r^{D-1}\phi\bigg{(}1+\alpha^{\prime}(D-3)\phi\bigg{)}+\frac{r^{D-1}}{l^{2}}+\frac{2Q^{2}r^{3-D}}{(D-3)(D-2)}-\frac{2ar}{(D-2)}-\frac{\alpha}{r^{\omega(D-1)}}\Bigg{]}^{\prime}},$
(15)

in which ′ denotes derivative with respect to $r$. Equation of motion is given
by variation of $\delta N(r)$ [44]. The function $\phi$ is given by solving
for the roots of a quadratic polynomial,

$\phi+\alpha^{\prime}(D-3)\phi^{2}=\frac{16\pi
M}{(D-2)r^{D-1}\Omega_{D-2}}-\frac{1}{l^{2}}-\frac{2Q^{2}r^{4-2D}}{(D-3)(D-2)}-\frac{2a}{(D-2)r^{D-2}}+\frac{\alpha}{r^{(\omega+1)(D-1)}},$
(16)

where $M$ is a constant.
By solving Eq.(16) and plugging $\phi$ into Eq.(3) we will have,

$\begin{split}&f(r)=\kappa-\frac{r^{2}}{2\alpha^{\prime}(D-3)}\Bigg{[}-1\pm\\\
&\sqrt{1-4(D-3)\alpha^{\prime}\bigg{(}\frac{1}{l^{2}}+\frac{2Q^{2}r^{4-2D}}{(D-3)(D-2)}-\frac{16\pi
Mr^{1-D}}{(D-2)\Omega_{D-2}}+\frac{2ar^{2-D}}{D-2}+\frac{\alpha}{r^{(D-1)(\omega+1)}}\bigg{)}}\Bigg{]}.\end{split}$
(17)

The line elements of planar black branes can be written as,

$ds^{2}=-H(r)N(r)^{2}dt^{2}+\frac{dr^{2}}{H(r)}+\frac{r^{2}}{l^{2}}\sum_{i=1}^{D-2}dx_{i}^{2}.$
(18)

That’s easy to show $N(r)$ is constant by variation of $\phi(r)$ from Eq.(2).

$\begin{split}&H(r)=\frac{r^{2}}{2\alpha^{\prime}(D-3)}\Bigg{[}1-\\\
&\sqrt{1-4(D-3)\alpha^{\prime}\bigg{(}\frac{1}{l^{2}}+\frac{2Q^{2}r^{4-2D}}{(D-3)(D-2)}-\frac{16\pi
Mr^{1-D}}{(D-2)\Omega_{D-2}}+\frac{2ar^{2-D}}{D-2}+\frac{\alpha}{r^{(D-1)(\omega+1)}}\bigg{)}}\Bigg{]},\end{split}$
(19)

where $M$ and $Q$ are integration constants proportional to the mass and
charge of the black hole respectively given by the following formulas:

$\begin{split}&M=\frac{(D-2)\Omega_{D-2}}{16\pi}m,\\\
&Q^{2}=\frac{(D-2)(D-3)}{2}q^{2}.\\\ \end{split}$ (20)

By defining new parameters as below,

$\begin{split}&\frac{(D-3)\alpha^{\prime}}{l^{2}}=\lambda_{gb},\\\
&\frac{2a}{D-2}=A,\end{split}$ (21)

and substituting them into the Eq.(19) we will have,

$H(r)=\frac{r^{2}}{2\lambda_{gb}l^{2}}\Bigg{[}1-\sqrt{1-4\lambda_{gb}\bigg{(}1+\frac{q^{2}l^{2}}{r^{2D-4}}-\frac{ml^{2}}{r^{D-1}}+\frac{Al^{2}}{r^{D-2}}+\frac{\alpha
l^{2}}{r^{(D-1)(\omega+1)}}\bigg{)}}\Bigg{]}.$ (22)

Black hole has an event horizon, so by applying the condition $H(r_{+})=0$ we
can find $m$ as follows,

$m=\frac{r_{+}^{D-1}}{l^{2}}+\frac{q^{2}}{r_{+}^{D-3}}+Ar_{+}+\frac{\alpha
l^{2}}{r_{+}^{(D-1)\omega}}.$ (23)

By plugging $m$ in Eq.(22)

$H(r)=\frac{r^{2}}{2\lambda_{gb}l^{2}}\Bigg{[}1-\sqrt{1-4\lambda_{gb}\Bigg{(}1-(\frac{r_{+}}{r})^{D-1}+l^{2}q^{2}\bigg{(}\frac{1}{r^{2D-4}}-(\frac{r_{+}}{r})^{D-1}\frac{1}{r_{+}^{2D-4}}\bigg{)}+\mathcal{P}\bigg{)}\Bigg{)}}\Bigg{]},$
(24)

where,

$\mathcal{P}=Al^{2}\bigg{(}\frac{1}{r^{D-2}}-(\frac{r_{+}}{r})^{D-1}\frac{1}{r_{+}^{D-2}}\bigg{)}+\frac{\alpha
l^{2}}{r^{(D-1)}}\bigg{(}\frac{1}{r^{(D-1)\omega}}-\frac{1}{r_{+}^{(D-1)\omega}}\bigg{)}.$
(25)

In AdS/CFT correspondence, the speed of light in the boundary CFT is simply
$c=1$. So that, for the black brane solution in the asymptotic region we have
$\lim_{r\to 0}{N(r)^{2}H(r)}=1$ to recover a causal boundary. By applying this
criterion we will have,

$N^{2}=\frac{1+\sqrt{1-4\lambda_{gb}}}{2}.$ (26)

The temperature and the Hawking-Bekenstein entropy density follow as,

$\displaystyle
T=\frac{1}{4\pi\sqrt{g_{rr}g_{tt}}}\partial_{r}g_{tt}|_{r=r_{+}}=\frac{Nr_{+}}{4\pi
l^{2}}\Big{[}(D-1)+\frac{Al^{2}}{r_{+}^{D-2}}+\frac{(D-3)l^{2}q^{2}}{r_{+}^{2D-4}}-\frac{(D-1)\alpha\omega}{r_{+}^{(D-1)(\omega+1)}}\Big{]},$
(27) $\displaystyle s=\frac{4\pi}{V}\int
d^{3}x\sqrt{-g}=4\pi\left(\frac{r_{+}}{l}\right)^{3}.$ (28)

##  3 $\frac{\eta}{s}$ of this solution

By introducing new variables $z=\frac{r}{r_{+}}$ ,
$\omega=\frac{l^{2}}{r_{+}}\tilde{\omega}$ ,
$k_{3}=\frac{l^{2}}{r_{+}^{2}}\tilde{k}_{3}$, ,
$\tilde{H}(z)=\frac{l^{2}}{r_{+}^{2}}H$, The line element in terms of new
variables is as follows,

$ds^{2}=-N^{2}\tilde{H}dt^{2}+\frac{dz^{2}}{r_{+}^{2}\tilde{H}}+\frac{r_{+}^{2}z^{2}}{l^{2}}\sum_{i=1}^{D-2}dx_{i}^{2},$
(29)

where

$\tilde{H}(z)=\frac{z^{2}}{2\lambda_{gb}}\Bigg{[}1-\sqrt{1-4\lambda_{gb}\Bigg{(}1-\frac{1}{z^{D-1}}+\frac{l^{2}q^{2}}{r_{+}^{2D-4}}\bigg{(}\frac{1}{z^{2D-4}}-\frac{1}{z^{D-1}}\bigg{)}+\mathcal{P}\Bigg{)}}\Bigg{]},$
(30)
$\mathcal{P}=\frac{Al^{2}}{r_{+}^{D-2}}\bigg{(}\frac{1}{z^{D-2}}-\frac{1}{z^{D-1}}\bigg{)}+\frac{\alpha
l^{2}}{r_{+}^{(D-1)(\omega+1)}}\bigg{(}\frac{1}{z^{(D-1)\omega}}-\frac{1}{z^{(D-1)(\omega+1)}}\bigg{)}.$
(31)

To calculate the shear viscosity, we perturb the background metric by
$h^{y}_{x}$. We call the perturbed part of metric $\psi=h^{y}_{x}$ [45, 46].

$ds^{2}=-N^{2}\tilde{H}dt^{2}+\frac{dz^{2}}{u^{2}\tilde{H}}+\frac{z^{2}}{u^{2}l^{2}}\Bigg{(}\sum_{i=1}^{D-2}dx_{i}^{2}+2\psi(t,\vec{x},z)dx_{1}dx_{2}\Bigg{)},$
(32)

where $u=\frac{1}{r_{+}^{2}}$.
Using Fourier decomposition,

$\psi(t,x_{i},z)=\int{\frac{d^{D-1}k}{(2\pi)^{D-1}}e^{-i\tilde{\omega}t+\tilde{k_{i}}x_{i}}}\psi(k,z).$
(33)

We substitute the perturbed metric Eq.(32) in Eq.(1) and expand the action up
to second order of $\psi$. Finally, the equation of motion for $\psi(z)$ can
be obtained by variation of the perturbed action with respect to $\psi$ as
follows,

$(K\psi^{\prime})^{\prime}+\omega^{2}K_{1}\psi-k_{i}^{2}K_{2}\psi=0,$ (34)

where

$\displaystyle K=\frac{1}{16\pi}\sqrt{-g}g^{zz}g^{x_{i}x_{i}}K_{3}(z),$
$\displaystyle K_{1}=-\frac{1}{16\pi}\sqrt{-g}g^{tt}K_{3}(z),$ $\displaystyle
K_{2}=\sqrt{-g}g^{x_{i}x_{i}}K_{4}(z),$ $\displaystyle
K_{3}(z)=1-\frac{2\lambda_{gb}}{D-3}[z^{-1}\tilde{H}^{\prime}+z^{-2}(D-5)\tilde{H}],$
$\displaystyle
K_{4}(z)=1-\frac{2\lambda_{gb}}{(D-3)(D-4)}[\tilde{H}^{\prime\prime}+(D-5)(D-6)z^{-2}\tilde{H}+2(D-5)z^{-1}\tilde{H}^{\prime}],$

in which ′ denotes derivative with respect to $z$. The factors $(D-5)$ and
$(D-6)$ in the expression of $K_{4}(z)$ is for $D>5$ of the Gauss-Bonnet
theory.
The Green function is,

$G_{x_{i}x_{j},x_{i}x_{j}}=\frac{K\psi^{\prime}}{\psi}\,.$ (36)

The shear viscosity is calculated by Green-Kubo formula as the following,

$\eta_{x_{i}x_{j},x_{i}x_{j}}=\frac{-G_{x_{i}x_{j},x_{i}x_{j}}}{i\omega},$
(37)

Now we write Eq.(34) in terms of the shear viscosity,

$\partial_{z}\eta_{x_{i}x_{j},x_{i}x_{j}}=(\frac{\eta^{2}_{x_{i}x_{j},x_{i}x_{j}}}{K}-K_{1})+\frac{i}{\omega}K_{2}k_{i}^{2},$
(38)

and then we can compute the shear viscosity by requiring horizon regularity as
follows,

$\eta_{x_{i}x_{j},x_{i}x_{j}}=\sqrt{KK_{1}}\Bigg{|}_{z=1}=\frac{1}{16\pi}\frac{r_{+}^{D-2}}{l^{D-2}}\left(1-\frac{2\lambda}{D-3}\tilde{H}^{\prime}(1)\right),$
(39)

The ratio of the shear viscosity to the entropy density for 4D charged black
hole solutions in Gauss-Bonnet gravity in the presence of string cloud and
quintessence is then

$\frac{\eta}{s}=\frac{1}{4\pi}[1-\frac{2\lambda}{D-3}\tilde{H}^{\prime}(1)],$
(40)
$\frac{\eta}{s}=\frac{1}{4\pi}\Bigg{[}1-\frac{2\lambda}{D-3}\Bigg{(}(D-1)-(D-3)l^{2}q^{2}r_{+}^{4-2D}+Al^{2}r_{+}^{2-D}+\alpha
l^{2}(D-1)r_{+}^{(1-D)(1+\omega)}\Bigg{)}\Bigg{]}.$ (41)

For $\lambda\to 0$ this value is $\frac{\eta}{s}=\frac{1}{4\pi}$. It means
that the KSS bound is saturated for Einstein-Hilbert gravity.
For $q=A=\alpha=0$ we will have,

$\frac{\eta}{s}=\frac{1}{4\pi}\big{[}1-2\frac{(D-1)}{(D-3)}\lambda\big{]},$
(42)

which is consistent with the literature[47].
In the limit $D\to 4$, we obtain

$\frac{\eta}{s}=\frac{1}{4\pi}\Bigg{[}1-2\lambda\Bigg{(}3-l^{2}q^{2}r_{+}^{-4}+Al^{2}r_{+}^{-2}+3\alpha
l^{2}r_{+}^{-3(1+\omega)}\Bigg{)}\Bigg{]}.$ (43)

For $A=\alpha=0$ the result is the same as [48],

$\frac{\eta}{s}=\frac{1}{4\pi}\Bigg{[}1-2\lambda\Big{(}3-l^{2}q^{2}r_{+}^{-4}\Big{)}\Bigg{]}.$
(44)

##  4 Conclusion

In summary, we have investigated the effect of string cloud , quintessence and
Yang-Mills charge in $4D$ Gauss-Bonnet gravity on the field theory dual side
by the fluid-gravity duality. Our result shows that theses three matters are
affected in dual of Gauss-Bonnet gravity but they have no contribution in dual
of Einstein-Hilbert gravity. Therefore, the KSS conjecture [12] is saturated
in the limit $\lambda\to 0$ . This conjecture tells us that the ratio $\eta/s$
has a lower bound, $\frac{\eta}{s}\geq\frac{\hbar}{4\,\pi\,k_{B}}$, for all
relativistic quantum field theories at finite temperature without chemical
potential [14] and can be interpreted as the Heisenberg uncertainty principle
[11]. String cloud [49] and quintessence [22] in Einstein massive gravity does
not contribute to the $\frac{\eta}{s}$. However, this conjecture violates for
higher derivative theories of gravity like the AdS Einstein-Gauss-Bonnet-Yang-
Mills in the presence of string cloud and quintessence as we confirmed in
section 3 . Our outcome shows that the coupling of this model is scaled as
$\lambda\to\lambda\Big{(}(D-1)-(D-3)l^{2}q^{2}r_{+}^{4-2D}+Al^{2}r_{+}^{2-D}+\alpha
l^{2}(D-1)r_{+}^{(1-D)(1+\omega)}\Big{)}$.

Acknowledgment Author would like to thank Xian-Hui Ge , M.H. Abbasi ,and the
referees of MPLA for useful comments and suggestions.

Data Availability  This manuscript has no associated data.

## References

  * [1] J. M. Maldacena, Int. J. Theor. Phys. 38 (1999) 1113 [Adv. Theor. Math. Phys. 2 (1998) 231] [hep-th/9711200].
  * [2] S. S. Gubser, I. R. Klebanov and A. M. Polyakov, Phys. Lett. B 428, 105-114 (1998) doi:10.1016/S0370-2693(98)00377-3 [arXiv:hep-th/9802109 [hep-th]].
  * [3] O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri and Y. Oz, Phys. Rept. 323 (2000), 183-386 doi:10.1016/S0370-1573(99)00083-6 [arXiv:hep-th/9905111 [hep-th]].
  * [4] E. Witten, Adv. Theor. Math. Phys. 2, 253-291 (1998) doi:10.4310/ATMP.1998.v2.n2.a2 [arXiv:hep-th/9802150 [hep-th]].
  * [5] W. J. Jiang, H. S. Liu, H. Lu and C. N. Pope, JHEP 07 (2017), 084 doi:10.1007/JHEP07(2017)084 [arXiv:1703.00922 [hep-th]].
  * [6] M. Baggioli and W. J. Li, JHEP 07 (2017), 055 doi:10.1007/JHEP07(2017)055 [arXiv:1705.01766 [hep-th]].
  * [7] E. Shuryak, Prog. Part. Nucl. Phys. 53, 273-303 (2004) doi:10.1016/j.ppnp.2004.02.025 [arXiv:hep-ph/0312227 [hep-ph]].
  * [8] S. A. Hartnoll, Class. Quant. Grav. 26 (2009), 224002 doi:10.1088/0264-9381/26/22/224002 [arXiv:0903.3246 [hep-th]].
  * [9] D.T. Son, Nuclear Physics B (Proc. Suppl.) 192–193 (2009) 113–118.
  * [10] G. Policastro, D. T. Son and A. O. Starinets,JHEP 0209, 043 (2002) [hep-th/0205052].
  * [11] D. T. Son and A. O. Starinets, Ann. Rev. Nucl. Part. Sci. 57 (2007), 95-118 doi:10.1146/annurev.nucl.57.090506.123120 [arXiv:0704.0240 [hep-th]].
  * [12] P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94 (2005), 111601 doi:10.1103/PhysRevLett.94.111601 [arXiv:hep-th/0405231 [hep-th]].
  * [13] S. Bhattacharyya, V. E. Hubeny, S. Minwalla and M. Rangamani, JHEP 02 (2008), 045 doi:10.1088/1126-6708/2008/02/045 [arXiv:0712.2456 [hep-th]].
  * [14] G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 87 (2001), 081601 doi:10.1103/PhysRevLett.87.081601 [arXiv:hep-th/0104066 [hep-th]].
  * [15] M. Rangamani, Class. Quant. Grav. 26 (2009), 224003 doi:10.1088/0264-9381/26/22/224003 [arXiv:0905.4352 [hep-th]].
  * [16] J. Bhattacharya, S. Bhattacharyya, S. Minwalla and A. Yarom, JHEP 05 (2014), 147 doi:10.1007/JHEP05(2014)147 [arXiv:1105.3733 [hep-th]].
  * [17] P. Kovtun, J. Phys. A 45 (2012), 473001 doi:10.1088/1751-8113/45/47/473001 [arXiv:1205.5040 [hep-th]].
  * [18] M. Brigante, H. Liu, R. C. Myers, S. Shenker and S. Yaida, Phys. Rev. Lett. 100 (2008), 191601 doi:10.1103/PhysRevLett.100.191601 [arXiv:0802.3318 [hep-th]].
  * [19] X. H. Ge, Y. Matsuo, F. W. Shu, S. J. Sin and T. Tsukioka, JHEP 10, 009 (2008) doi:10.1088/1126-6708/2008/10/009 [arXiv:0808.2354 [hep-th]].
  * [20] X. H. Ge and S. J. Sin, JHEP 05, 051 (2009) doi:10.1088/1126-6708/2009/05/051 [arXiv:0903.2527 [hep-th]].
  * [21] X. H. Feng, H. S. Liu, H. Lü and C. N. Pope, JHEP 11 (2015), 176 doi:10.1007/JHEP11(2015)176 [arXiv:1509.07142 [hep-th]].
  * [22] M. Sadeghi, Mod. Phys. Lett. A 36 (2021) no.28, 2150202 doi:10.1142/S0217732321502023 [arXiv:2007.09688 [hep-th]].
  * [23] M. Bravo-Gaete and F. F. Santos, Eur. Phys. J. C 82, no.2, 101 (2022) doi:10.1140/epjc/s10052-022-10064-y [arXiv:2010.10942 [hep-th]].
  * [24] M. Bravo-Gaete and M. M. Stetsko, Phys. Rev. D 105, no.2, 024038 (2022) doi:10.1103/PhysRevD.105.024038 [arXiv:2111.10925 [hep-th]].
  * [25] M. Bravo-Gaete, F. F. Santos and H. Boschi-Filho, [arXiv:2201.07961 [hep-th]].
  * [26] K. A. Mamo, JHEP 10, 070 (2012) doi:10.1007/JHEP10(2012)070 [arXiv:1205.1797 [hep-th]].
  * [27] A. G. Riess et al. [Supernova Search Team], Astron. J. 116 (1998), 1009-1038 doi:10.1086/300499 [arXiv:astro-ph/9805201 [astro-ph]].
  * [28] S. Perlmutter et al. [Supernova Cosmology Project], Astrophys. J. 517, 565-586 (1999) doi:10.1086/307221 [arXiv:astro-ph/9812133 [astro-ph]].
  * [29] P. A. R. Ade et al. [Planck], Astron. Astrophys. 571 (2014), A16 doi:10.1051/0004-6361/201321591 [arXiv:1303.5076 [astro-ph.CO]].
  * [30] L. M. Wang, R. R. Caldwell, J. P. Ostriker and P. J. Steinhardt, Astrophys. J. 530 (2000), 17-35 doi:10.1086/308331 [arXiv:astro-ph/9901388 [astro-ph]].
  * [31] S. Tsujikawa, Class. Quant. Grav. 30 (2013), 214003 doi:10.1088/0264-9381/30/21/214003 [arXiv:1304.1961 [gr-qc]].
  * [32] V. V. Kiselev, Class. Quant. Grav. 20 (2003), 1187-1198 doi:10.1088/0264-9381/20/6/310.
  * [33] P. S. Letelier, Phys. Rev. D 20 (1979), 1294-1302 doi:10.1103/PhysRevD.20.1294
  * [34] S. G. Ghosh, U. Papnoi and S. D. Maharaj, Phys. Rev. D 90, no. 4, 044068 (2014) doi:10.1103/PhysRevD.90.044068 [arXiv:1408.4611 [gr-qc]].
  * [35] E. Herscovich and M. G. Richarte, Phys. Lett. B 689 (2010), 192-200 doi:10.1016/j.physletb.2010.04.065 [arXiv:1004.3754 [hep-th]].
  * [36] T. H. Lee, D. Baboolal and S. G. Ghosh, Eur. Phys. J. C 75, no. 7, 297 (2015) doi:10.1140/epjc/s10052-015-3515-5 [arXiv:1409.2615 [gr-qc]].
  * [37] S. H. Mazharimousavi and M. Halilsoy, Eur. Phys. J. C 76, no. 2, 95 (2016) doi:10.1140/epjc/s10052-016-3954-7 [arXiv:1511.00603 [gr-qc]].
  * [38] J. P. Morais Graça, G. I. Salako and V. B. Bezerra, Int. J. Mod. Phys. D 26, no. 10, 1750113 (2017) doi:10.1142/S0218271817501139 [arXiv:1604.04734 [gr-qc]].
  * [39] H. Ranjbari, M. Sadeghi, M. Ghanaatian and G. Forozani, Eur. Phys. J. C 80, no.1, 17 (2020) doi:10.1140/epjc/s10052-019-7592-8 [arXiv:1911.10803 [hep-th]].
  * [40] D. Glavan and C. Lin, Phys. Rev. Lett. 124, no.8, 081301 (2020) doi:10.1103/PhysRevLett.124.081301 [arXiv:1905.03601 [gr-qc]].
  * [41] V. V. Kiselev, Class. Quant. Grav. 20, 1187-1198 (2003) doi:10.1088/0264-9381/20/6/310 [arXiv:gr-qc/0210040 [gr-qc]].
  * [42] S. Chen, B. Wang and R. Su, Phys. Rev. D 77 (2008), 124011 doi:10.1103/PhysRevD.77.124011 [arXiv:0801.2053 [gr-qc]].
  * [43] B. L. Shepherd and E. Winstanley, Phys. Rev. D 93, no. 6, 064064 (2016) doi:10.1103/PhysRevD.93.064064 [arXiv:1512.03010 [gr-qc]].
  * [44] R. C. Myers and B. Robinson, JHEP 08, 067 (2010) doi:10.1007/JHEP08(2010)067 [arXiv:1003.5357 [gr-qc]].
  * [45] G. Policastro, D. T. Son and A. O. Starinets,JHEP 0209, 043 (2002) [hep-th/0205052].
  * [46] P. Kovtun, D. T. Son and A. O. Starinets, JHEP 0310, 064 (2003) [hep-th/0309213].
  * [47] M. Brigante, H. Liu, R. C. Myers, S. Shenker and S. Yaida, Phys. Rev. D 77, 126006 (2008) doi:10.1103/PhysRevD.77.126006 [arXiv:0712.0805 [hep-th]].
  * [48] X. H. Ge and S. J. Sin, Eur. Phys. J. C 80, no.8, 695 (2020) doi:10.1140/epjc/s10052-020-8288-9 [arXiv:2004.12191 [hep-th]].
  * [49] M. Sadeghi and H. Ranjbari, Class. Quant. Grav. 36, no.20, 205012 (2019) doi:10.1088/1361-6382/ab436e [arXiv:1905.12856 [hep-th]].