# Fluidity of the system produced in relativistic pp and heavy-ion collisions:
Hadron resonance gas model approach

Ronald Scaria1 Dushmanta Sahu1 Captain R. Singh1 Raghunath Sahoo1,2
111Corresponding author<EMAIL_ADDRESS>Jan-e Alam3 1Department of
Physics, Indian Institute of Technology Indore, Simrol, Indore 453552, India
2CERN, CH 1211, Geneva 23, Switzerland  3Variable Energy Cyclotron Centre,
1/AF, Bidhan Nagar, Kolkata, India

###### Abstract

We have estimated the dimensionless parameters such as Reynolds number ($Re$),
Knudsen number ($Kn$) and Mach number ($Ma$) for a multi-hadron system by
using the excluded volume hadron resonance gas (EVHRG) model along with
Hagedorn mass spectrum to include higher resonances in the system. The size
dependence of these parameters indicate that the system formed in
proton+proton collisions may achieve thermal equilibrium making it unsuitable
as a benchmark to analyze the properties of the system produced in heavy ion
collisions at similar energies. While the magnitude of $Kn$ can be used to
study the degree of thermalization, the variations of $Re$ and $Ma$ with
temperature ($T$) and baryonic chemical potential ($\mu_{B}$) assist to
understand the change in the nature of the flow in the system - from laminar
to turbulent and from subsonic to supersonic, respectively.

##  I Introduction

Relativistic heavy-ion collisions aim to create extreme conditions of
temperature ($T$) and baryonic chemical potential ($\mu_{B}$) so as to melt
the hadrons to produce a deconfined state of quarks and gluons, called Quark-
Gluon Plasma (QGP). Although initially it was envisaged that such a state of
matter would behave as an ideal gas, the analysis of the results of the
nuclear collisions at the Relativistic Heavy-Ion Collider (RHIC) and Large
Hadron Collider (LHC) energies reveal that the created matter behaves like
strongly interacting liquid with lowest shear viscosity ($\eta$) to entropy
density ($s$) ratio. The study of experimental data within the framework of
relativistic hydrodynamics shows that the $\eta/s$ of QGP is close to
$1/4\pi$, which is the lowest AdS/CFT bound obtained in Ref. Kovtun:2004de ,
indicating the creation of a “perfect fluid” in relativistic heavy-ion
collisions at the RHIC Arsene2005 ; Back2005 ; Adams2005 ; Adcox2005 . So far,
the results from proton+proton (pp) collisions were used as a baseline for
establishing QCD medium formation in heavy-ion collisions at RHIC and LHC
energies. However, recent observations of strangeness enhancement,
collectivity, long-range particle correlations etc. in pp and p+Pb collisions
at the LHC energies Adam2017 ; Khachatryan2017  opened a new domain of
studies. These observations advocate for the possibility of thermalization and
collectivity in the small system produced in pp collisions. Final state
multiplicity dependent event characterization and collision centrality
selection even in such collisions emerged out with intense search for QGP-like
properties in high-multiplicity pp collisions. Although the formation of QGP-
droplets is yet to be established, a detailed characterization of the system
formed in pp collisions at the LHC through various theoretical models and
experimental searches using new observables are underway. Combining pp, p+Pb
and heavy-ion collisions at the LHC energies, final state event multiplicity
has become one of the unique observables to describe the system size and
collision energy dependence of various observables. In this direction, one
tries to find a threshold in the final state multiplicity Velicanu:2011zz ;
Campanini:2011bj ; Thakur:2017kpv ; Sahu:2020nbu  after which QGP-like
properties and the corresponding theoretical models to understand a QCD medium
formation, including those to study thermalization and hydrodynamic models
could be employed.

The nature of the hydrodynamic flow (laminar or turbulent) can be classified
based on the magnitude of the Reynolds number ($Re$). The Mach number ($Ma$)
helps in determining whether the flow is supersonic or subsonic and the
Knudsen number ($Kn$) can be an indicator of thermalization. Fluid with high
value of $Re$ may show turbulent behavior. As
$Re\sim\frac{inertialforce}{viscousforce}\sim 1/\eta$, systems with low
viscosity are likely to manifest turbulent nature. At RHIC energy for a system
size of 5 fm and temperature around $T=300$ MeV, $Re$ varies between 10$-$100
for a fluid velocity in the range 0.1$-$1 (in the unit of $c=1$) Deng:2016gyh
. $Kn$ is the ratio between mean free path of the system and the system size.
Therefore, $Kn$ value tending to zero implies a system with a high degree of
thermalization. A large value of $Kn$ suggests a system that is far from
equilibrium and in such a case, the application of hydrodynamics is not
allowed. For Au$+$Au collision at $\sqrt{s_{\rm NN}}$ = 200 GeV, the Knudsen
number was estimated from a hydrodynamic model, in which the elliptic flow as
a function of centrality suggests the value of $Kn$ $\simeq$ 0.3 for central
collisions and 0.5 for semi-central collisions Bhalerao:2005mm ;
Drescher:2007cd . Similarly, Mach number is defined as the ratio of the flow
velocity to the speed of sound in the system. For $Ma$ $<<$ 1, the fluid
density is almost uniform and consequently suggests incompressible flow. Thus,
Mach number can help in characterizing the compressibility of the fluid. These
three parameters together define the applicability of hydrodynamics and give
information about properties like compressibility and viscosity of the fluid
Ollitrault2008 . These observables provide a unique way to observe the
deviation of a fluid (system) from that of the perfect fluid under the limits:
$Kn$ $<<$ 1, $Ma$ $<<$ 1 and $Re$ $>>$ 1\.

To understand the thermodynamics of the system produced in the ultra-
relativistic collisions, inputs from lattice quantum chromodynamics (LQCD) has
been extremely useful. However, it has limitations, such that it can work
accurately for zero or low baryochemical potential ($\mu_{B}/T\rightarrow 0$).
In this context models such as Hadron Resonance Gas (HRG) can be useful to
analyze data without having restrictions on the value of baryochemical
potential. The HRG model has been used successfully to study the
thermodynamics of the hadronic phase that appears in the course of evolution
of the system produced in ultra-relativistic collisions. In this article, we
have used the HRG model with excluded volume effect along with a Hagedorn mass
spectrum to study the thermal behavior of the hadronic phase formed in
relativistic collisions Andronic2012 ; Vovchenko2015 ; Rischke1991 . This
article is organized as follows. Detailed descriptions of the formalism and
methodology to calculate the thermodynamic properties are given in Sec. II. In
Sec. III, we present the results along with a detailed discussion. Finally, we
summarize our study in Sec. IV.

##  II Formalism

###  II.1 Hadron Resonance Gas with Excluded Volume

The hadron resonance gas model describes the hadronic phase of the matter by
assuming that the system is composed of free hadrons of different species. It
treats various hadrons and their resonances as point particles. However, the
existence of repulsive interactions among the hadrons has already been known
from the nucleon-nucleon scattering experiments. The excluded volume hadron
resonance gas (EVHRG) model takes such repulsive interactions into account. It
reduces the volume available for the hadrons to move by the volume that they
occupy. The EVHRG model estimations for various thermodynamic observables
match the LQCD estimations up to $T\sim$ 140 MeV. It is also essential to take
into account the effect of the Hagedorn mass spectrum, which includes the yet
unmeasured higher masses of the hadrons, as they can have a significant
contribution to the Equation of State (EoS) near the critical temperature
($T_{c}$).

In a Grand Canonical Ensemble (GCE) of ideal HRG (IHRG) formalism, the
partition function for the $i^{th}$ particle can be written as Andronic2012 ;

$lnZ^{id}_{i}=\pm\frac{Vg_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\
ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\}$ (1)

where $g_{i}$ and $E_{i}=\sqrt{p^{2}+m_{i}^{2}}$ are the degeneracy and energy
of the $i^{th}$ hadron, respectively. The corresponding chemical potential is
given by,

$\mu_{i}=B_{i}\mu_{B}+S_{i}\mu_{S}+Q_{i}\mu_{Q},$ (2)

where $B_{i}$, $S_{i}$ and $Q_{i}$ are the baryon number, strangeness and
electric charge of the $i^{th}$ hadron. The $\mu_{B}$, $\mu_{S}$, $\mu_{Q}$
are the baryon chemical potential, strangeness chemical potential and chemical
potential related to the electric charge, respectively. From the partition
function one can obtain the pressure $P_{i}$, energy density
$\varepsilon_{i}$, number density $n_{i}$ and entropy density $s_{i}$ as
follows:

$P^{id}_{i}(T,\mu_{i})=\pm\frac{Tg_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\
ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\}$ (3)
$\varepsilon^{id}_{i}(T,\mu_{i})=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\frac{E_{i}\
p^{2}dp}{\exp[(E_{i}-\mu_{i})/T]\pm 1}$ (4)
$n^{id}_{i}(T,\mu_{i})=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\frac{p^{2}dp}{\exp[(E_{i}-\mu_{i})/T]\pm
1}$ (5) $\displaystyle s^{id}_{i}(T,\mu_{i})=$
$\displaystyle\pm\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\Big{[}\ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\}$
$\displaystyle\pm\frac{(E_{i}-\mu_{i})/T}{\exp[(E_{i}-\mu_{i})/T]\pm
1}\Big{]}$ (6)

The effect of higher states and resonances is included into IHRG by making use
of the following parametrization for the Hagedorn mass spectrum Hagedorn1968 ;

$\rho(m)=C\frac{\theta(m-M)}{(m^{2}+m_{0}^{2})^{a}}\exp\Big{(}\frac{m}{T_{H}}\Big{)},$
(7)

where $M$ is the mass cut-off for selected hadrons and $T_{H}$ is the Hagedorn
temperature chosen as 180 MeV for our analysis. For other parameters, we have
chosen $C=0.05$ GeV3/2, $m_{0}=0.5$ GeV and $a=5/4$ Vovchenko2015 . The
inclusion of Hagedorn states leads us to adopt Boltzmann approximation. Thus,
the previously defined thermodynamic properties are modified as;

$\displaystyle P^{H}(T,\mu)=$ $\displaystyle\frac{T}{2\pi^{2}}\int
dm\int_{0}^{\infty}p^{2}dp\ \exp\Big{(}-\frac{E-\mu}{T}\Big{)}$
$\displaystyle\times\Big{[}\sum_{i}g_{i}\delta(m-m_{i})+\rho(m)\Big{]}$ (8)
$\displaystyle\varepsilon^{H}(T,\mu)=$ $\displaystyle\frac{1}{2\pi^{2}}\int
dm\int_{0}^{\infty}p^{2}E\ dp\ \exp\Big{(}-\frac{E-\mu}{T}\Big{)}$
$\displaystyle\times\Big{[}\sum_{i}g_{i}\delta(m-m_{i})+\rho(m)\Big{]}$ (9)
$\displaystyle n^{H}(T,\mu)=$ $\displaystyle\frac{1}{2\pi^{2}}\int
dm\int_{0}^{\infty}p^{2}dp\ \exp\Big{(}-\frac{E-\mu}{T}\Big{)}$
$\displaystyle\times\Big{[}\sum_{i}g_{i}\delta(m-m_{i})+\rho(m)\Big{]}$ (10)
$\displaystyle s^{H}(T,\mu)=$ $\displaystyle\frac{1}{2\pi^{2}}\int
dm\int_{0}^{\infty}p^{2}dp\ \exp\Big{(}-\frac{E-\mu}{T}\Big{)}$
$\displaystyle\times\Big{(}1+\frac{E-\mu}{T}\Big{)}\Big{[}\sum_{i}g_{i}\delta(m-m_{i})+\rho(m)\Big{]}$
(11)

Repulsive interactions at short distances are introduced into IHRG by using
thermodynamically consistent excluded volume approach Rischke1991 . This, in
essence is a van der Waals correction. The system volume $V$ is now replaced
by the available volume $V_{av}$,

$V\rightarrow V_{av}=V-\sum_{i}v_{i}N_{i},$ (12)

where $v_{i}=16\pi r_{i}^{3}/3$ is the excluded volume parameter to take care
of two overlapping hadrons (also called as eigen volume) Andronic2012 ;
Vovchenko2015  and $N_{i}$ is the number of particles. From here on, we
restrict ourselves to the case of equal volume parameter $v_{i}$ for all
hadrons, $v=16\pi r_{h}^{3}/3$. Use of Eq. 12 in grand canonical ensemble
changes pressure through an iterative process Andronic2012 ; Vovchenko2015 ,

$P^{ev}(T,\mu)=\kappa P^{H}(T,\mu),$ (13)

where, $\kappa$ is the excluded volume suppression factor given by,

$\kappa=exp\bigg{(}-\frac{vP^{ev}}{T}\bigg{)}.$ (14)

Accordingly, the thermodynamic entities like $\varepsilon$, $n$ and $s$ are
modified as follows;

$\varepsilon^{ev}(T,\mu)=\frac{\kappa\varepsilon^{H}(T,\mu)}{1+\kappa
vn^{H}(T,\mu)}$ (15) $n^{ev}(T,\mu)=\frac{\kappa n^{H}(T,\mu)}{1+\kappa
vn^{H}(T,\mu)}$ (16) $s^{ev}(T,\mu)=\frac{\kappa s^{H}(T,\mu)}{1+\kappa
vn^{H}(T,\mu)}.$ (17)

###  II.2 $Re$, $Kn$ and $Ma$ Dynamics in Hadron Gas

The transport properties of hadron gas can be investigated by using the
Boltzmann Transport Equation (BTE) which is given as,

$\frac{\partial f_{p}}{\partial t}+v^{i}_{p}\frac{\partial f_{p}}{\partial
x^{i}}+F^{i}_{p}\frac{\partial f_{p}}{\partial p^{i}}=I(f_{p}),$ (18)

where $v^{i}_{p}$ is the velocity of the $i$th hadron and $F^{i}_{p}$ is the
external force acting on it. $I(f_{p})$ is defined as the collision integral
which gives the rate of change of distribution function due to collisions. The
thermal equilibrium in a system is achieved and maintained by the collisions
among the constituents. Under the relaxation time approximation the collision
integral can be written as:

$I(f_{p})=-\frac{(f_{p}-f_{p}^{0})}{\tau(E_{p})},$ (19)

here $\tau(E_{p})$ is the particle dependent relaxation time and $f_{p}^{0}$
is the equilibrium distribution function defined as;

$f_{p}^{0}=\exp(-\frac{E_{i}-\mu}{T}).$ (20)

For a hadron gas at finite chemical potential, the shear viscosity can be
obtained by using the following expression Kadam2015 ,

$\eta=\sum_{i}\frac{1}{15T}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{4}}{E_{i}^{2}}(\tau_{i}f^{0}_{i}+\bar{\tau}_{i}\bar{f}^{0}_{i})$
(21)

where $\tau_{i}$ is the relaxation time of $i^{th}$ hadron and $f^{0}_{i}$ is
its corresponding distribution function. The particle dependent average
relaxation time is given by,

$\tilde{\tau}^{-1}_{i}=\sum_{j}n_{j}\langle\sigma_{ij}v_{ij}\rangle,$ (22)

where $n_{j}$ is the number density of $j^{th}$ hadronic species. Furthermore,
considering hadrons as hard spheres of equal radius $r_{h}$ (having constant
cross section $\sigma=4\pi r_{h}^{2}$), the thermal average of total cross
section times relative velocity i.e. $\langle\sigma v\rangle$ can be
calculated as Kadam2015 ; Cannoni2014 ; Gondolo1991 ,

$\langle\sigma_{ij}v_{ij}\rangle=\sigma\frac{\int
d^{3}p_{i}d^{3}p_{j}v_{ij}f^{0}_{i}f^{0}_{j}}{\int
d^{3}p_{i}d^{3}p_{j}f^{0}_{i}f^{0}_{j}}.$ (23)

Now, rewriting the momentum space volume elements in terms of scattering angle
$\theta$ Kadam2015 ; Tiwari2018 ;

$d^{3}p_{i}d^{3}p_{j}=8\pi^{2}p_{i}p_{j}dE_{i}dE_{j}d\cos\theta,$ (24)

the numerator in Eq. 23 becomes

$\displaystyle d^{3}p_{i}d^{3}p_{j}v_{ij}f^{0}_{i}f^{0}_{j}=\sigma\int
8\pi^{2}p_{i}p_{j}dE_{i}dE_{j}d\cos\theta f^{0}_{i}f^{0}_{j}$
$\displaystyle\times\frac{\sqrt{(E_{i}E_{j}-p_{i}p_{j}cos\
\theta)^{2}-(m_{i}m_{j})^{2}}}{E_{i}E_{j}-p_{i}p_{j}\cos\theta}$ (25)

and the denominator reads as:

$d^{3}p_{i}d^{3}p_{j}f^{0}_{i}f^{0}_{j}=\int
8\pi^{2}p_{i}p_{j}dE_{i}dE_{j}d\cos\theta f^{0}_{i}f^{0}_{j}.$ (26)

Now $\langle\sigma_{ij}v_{ij}\rangle$ takes the following form

$\langle\sigma_{ij}v_{ij}\rangle=\frac{\sigma\int
p_{i}p_{j}dE_{i}dE_{j}d\cos\theta
f^{0}_{i}f^{0}_{j}\times\frac{\sqrt{(E_{i}E_{j}-p_{i}p_{j}\cos\theta)^{2}-(m_{i}m_{j})^{2}}}{E_{i}E_{j}-p_{i}p_{j}\cos\theta}}{\int
p_{i}p_{j}dE_{i}dE_{j}d\cos\theta f^{0}_{i}f^{0}_{j}}$ (27)

$E_{i}$ and $E_{j}$ are integrated in the limit $m_{i}$ to $\infty$ and
$m_{j}$ to $\infty$, respectively. The limit of integration of $\cos\theta$
runs from -1 to +1. Finally, the relaxation time is calculated using Eqs. 22
and 27.

The use of Boltzmann states cannot rigorously describe massive Hagedorn
states. This makes estimation of their contribution to shear viscosity
difficult. These states being highly unstable, decay rapidly and it is
reasonable to assume that their presence affects the mean free paths of other
hadrons in the medium. They also significantly affect the thermodynamics of
hadronic matter close to the QCD critical temperature. Hence their
contribution cannot be ignored. The contribution from these states to the
shear viscosity can be estimated by using the approximation that the mean free
path of such resonances is inversely related to their decay widths
Kadamahep2019 ; Noronha2012 ; Noronha2009 . The decay width of Hagedorn states
is obtained by a linear fit to decay widths of all resonances in particle data
book. Assuming that these particles have hard-core radius $r_{h}$, the
excluded volume approximation of shear viscosity contribution of Hagedorn
states is given by,

$\displaystyle\eta_{H}=\frac{5T^{3/2}}{128\pi^{5/2}n^{H}(T)}\int m^{5/2}dm\
\rho(m)K_{\frac{5}{2}}\left(\frac{m}{T}\right).$ (28)

Here $n^{H}(T)$ is the number density of Hagedorn states in the medium and
$K_{n}$ is the modified Bessel function of the second kind.

The Reynolds number, which varies inversely to shear viscosity can be written
as,

$Re=\frac{D\langle v\rangle T}{\eta/s^{ev}},$ (29)

where $\langle
v\rangle=m^{-1}\sum_{i=1}^{m}{\int\frac{d^{3}p}{(2\pi)^{3}}v_{i}f^{0}_{i}}/{\int\frac{d^{3}p}{(2\pi)^{3}}f^{0}_{i}}$
is the fluid velocity and $D$ is the system size. The Knudsen number is
written as,

$Kn=\lambda/D,$ (30)

where $\lambda=\frac{1}{\sqrt{2}n^{ev}\sigma}$ is the mean free path of the
system.

Finally, the Mach number is written as,

$Ma=\langle v\rangle/c_{s},$ (31)

where

$c_{s}^{2}=\frac{\frac{\partial P}{\partial T}+\frac{\partial
P}{\partial\mu_{B}}\frac{d\mu_{B}}{dT}}{\frac{\partial\varepsilon}{\partial
T}+\frac{\partial\varepsilon}{\partial\mu_{B}}\frac{d\mu_{B}}{dT}},$ (32)

is the speed of sound in the system which reduces to
$c_{s}^{2}=(\frac{\partial P}{\partial\varepsilon})$ at vanishing baryon
chemical potential. Here,

$\frac{d\mu_{B}}{dT}=\frac{s\frac{\partial n}{\partial T}-n\frac{\partial
s}{\partial T}}{n\frac{\partial s}{\partial\mu_{B}}-s\frac{\partial
n}{\partial\mu_{B}}}.$ (33)

With the help of the formulae given above, we now move forward to estimate the
dimensionless parameters $Re$, $Kn$ and $Ma$.

##  III Results and discussion

We include all the hadrons and resonances up to a mass of 2.25 GeV listed in
the particle data book Zyla2020  along with the Hagedorn spectrum. The
thermodynamic properties are dependent on the eigen volume of hadrons which
have to be chosen appropriately. It is dependent only on the unknown parameter
$r_{h}$. The hadron-hadron interactions are not repulsive for an $r_{h}$ value
comparable to the charge radii. Hence, it was argued that $r_{h}$ should be
taken according to the hard-core volume known from nucleon-nucleon scattering.
The repulsive interactions are mediated by $\omega$ mesons and the range of an
interaction is inversely proportional to the mass of the mediator. Based upon
the above argument, we have chosen a hard-core radius, $r_{h}=0.3$ fm, to
incorporate the repulsive interaction. Also there is not much deviation in
charge radii between baryons and mesons, therefore, uniform hard-core radius
($r_{h}=0.3$ fm) is considered for all hadrons Munzinger1999 ; Kadam2015 ;
Tiwari2018 .

The effect of system size on volume and particle numbers, considered infinite
in the thermodynamic limit, can be implemented by providing a lower momentum
cut-off to the integral over momentum space Redlich2016 ; Bhattacharyya2015 ;
Sarkar2017 . The finite-size effect is introduced by using lower limit of
momentum $p_{\text{cutoff}}$(MeV) = 197$\pi/D$(fm) Sarkar2017 . Here $D=2R$ is
characteristic system size, R being system radius. We chose few representative
radii, $R=$ 3, 5, 7 and 10 fm to cover system properties from high
multiplicity pp to ultra-relativistic heavy-ion collisions Sarkar2017 ;
Chatterjee2015 . This choice of system size goes inline with the recent
observations of finite hadronic phase lifetime observed in pp and heavy-ion
collisions Sahu:2019tch ; ALICE:2019xyr .

Figure 1: Comparison of $\eta/s$ obtained for $R=10.0$ fm with the estimations
of other models at $\mu_{B}=0.0$ GeV Gorenstein2008 ; Denicol2013 ; Kadam2019
; Kadamahep2019 .

To validate our study, we have compared our results with others in Fig. 1. It
shows the dependence of $\eta/s$ on system temperature for a fixed system size
at zero chemical potential ($\mu_{B}=0.0$ GeV). The AdS/CFT lower limit to
$\eta/s$ (also known as KSS limit), which is $1/4\pi$ Kovtun:2004de , is also
shown in the figure. It is observed that the $\eta/s$ in the given temperature
range of the hadronic phase, shows a near perfect fluid behavior towards
higher temperature. Our study is in good agreement with the values of $\eta/s$
obtained by Gorenstein et al., where the authors considered EVHRG within
relativistic molecular kinetic theory Gorenstein2008 . We have also compared
our $\eta/s$ values with the results obtained from Chapman-Enskog theory,
relativistic mean-field HRG and EVHRG formalism with Hagedorn spectrum
Denicol2013 ; Kadam2019 ; Kadamahep2019 .

Figure 2: The square of the speed of sound ($c_{s}^{2}$) as a function of
temperature ($T$). This is compared with lattice data Borsanyi2013 .

Further, the dependence of the speed of sound on temperature is compared with
lattice QCD data at $\mu_{B}=0.0$ GeV in Fig. 2. It is observed that
$c_{s}^{2}$ obtained in EVHRG formalism is in agreement with the LQCD data
except at the high-temperature regime Borsanyi2013 , where the HRG model is
not a good approximation for the description of the system. The addition of
the Hagedorn spectrum accentuates this deviation because of the rapid increase
in energy density as compared to pressure on approaching the Hagedorn
temperature.

Figure 3: The $c_{s}^{2}$ as a function of temperature and $\mu_{B}$ for
$R=10.0$ fm.

In Fig. 3 the variation of $c_{s}^{2}$ is explored with temperature ($T$) and
chemical potential ($\mu_{B}$). We have observed that increasing temperature
at very low chemical potential ($\mu_{B}=0-0.2$ GeV) gives the similar results
as discussed in Fig. 2. However, at $\mu_{B}>0.2$ GeV with increasing $T$
speed of sound (or $c_{s}^{2}$) decreases initially up to $T\sim 0.140$ GeV
and then start increasing with temperature. As $\mu_{B}$ increases the
$c_{s}^{2}$ dip point on the $T$-axis shifts towards the lower temperature
end. This smooth change in the behaviour of $c_{s}^{2}$ with high $\mu_{B}$
mimics the phase transition at $T<T_{c}\approx 0.160$ GeV. The sudden drop in
the $c_{s}^{2}$ at high $T$ is caused by the same effect of the Hagedorn
states as seen in Fig. 2.

Figure 4: Reynolds number as a function of temperature for different values of
$R$ at $\mu_{B}$ = 0 GeV (left panel) and different values of $\mu_{B}$ at
$R=10.0$ fm, which is the case of heavy-ion collisions (right panel).

Reynolds number, $Re$, is an important parameter in hydrodynamics that
contains information regarding instabilities in the fluid. A very high
Reynolds number is usually associated with the turbulence in the fluid. The
inverse of Reynolds number, $Re^{-1}$, is generally seen as a proxy for
$\eta/s$ in ultra-relativistic heavy-ion collisions as both of them are
directly related. As a result of this relation, $Re>>1$ implies inviscid flow.
We present the dependence of Reynolds number on system size in Fig. 4. We note
that it increases with temperature for all values of $R$, indicating that the
system becomes less viscous at higher temperatures and approaches a ’perfect
fluid’ limit. Fig. 4 also explores the variation of Reynolds number with
chemical potential (right panel). We have chosen five different values of
chemical potentials; $\mu_{B}=$ 0.0, 0.025, 0.200, 0.436 and 0.630 GeV, which
corresponds to LHC energies, RHIC at $\sqrt{s_{\rm NN}}=200$ GeV and
$\sqrt{s_{\rm NN}}=19.6$ GeV, RHIC/FAIR at $\sqrt{s_{\rm NN}}=7.7$ GeV and
NICA at $\sqrt{s_{\rm NN}}=3$ GeV respectively Tawfik2016 ; Munzinger2001 ;
Cleymans2006 ; Khuntia2019 . The system radius is fixed at 10.0 fm, which
roughly corresponds to ultra-relativistic heavy-ion collisions. At low
temperatures, $Re$ increases with $\mu_{B}$ indicating that viscosity of the
system decreases with increasing $\mu_{B}$. As it appears, at higher
temperature the Reynolds number, $Re$, becomes almost independent of the
baryochemical potential of the system.

Figure 5: Knudsen number, Kn, as a function of temperature for different
values of system size, $R$ at $\mu_{B}=0.0$ GeV (left panel) and for different
values of $\mu_{B}$ at $R=10.0$ fm, which is the case of larger systems (right
panel).

Large values for the Knudsen number, $Kn$, imply that the system is far from
thermodynamic equilibrium, while small values, tending to zero, indicate a
high degree of thermalization within the system. Fig. 5 shows the dependence
of Knudsen number with system radius. We see a lower value of $Kn$ even at low
temperatures for large $R$. This indicates that larger systems have a higher
tendency to achieve thermodynamic equilibrium. $Kn$ is small at high
temperatures and does not vary much with system radius. This behavior is
expected because at lower temperatures the mean free path ($\lambda$) is
larger, indicating less interaction among the hadrons. The variation of $Kn$
also indicates that the system formed at high temperatures is near to the
equilibrium and therefore, we can apply hydrodynamics to study the evolution
of the system. The variation of Knudsen number as a function of temperature
for different chemical potentials is also shown in the right panel of Fig. 5.
We see that $Kn$ decreases with increase in $\mu_{B}$.

Figure 6: Mach number, Ma, as a function of temperature for different values
of system size, $R$ at $\mu_{B}=0.0$ GeV (left panel) and for different values
of $\mu_{B}$ at $R=10.0$ fm, for larger systems (right panel).

Supersonic flows significantly impact two-particle correlations because of the
formation of Mach cones. The interest in Mach cones formed in ultra-
relativistic heavy-ion collisions was sparked from the two-particle
correlations observed at RHIC Aggarwal2010 ; Adare2010 . However, subsequent
measurements at the LHC energies do not support these observations Aad2014 ;
Chatrchyan2014 . It is argued that the structure observed might be an artifact
of incorrect background considerations Nattrass2016 . Although there is no
direct experimental evidence for their existence, theoretical interest in Mach
cones still persists as its disappearance could hint at a possible QCD
critical point Sarwar2021 . Mach number, $Ma$, being the ratio of
characteristic flow velocity to the velocity of sound, could tell about the
possible formation of Mach cones in the fluid. If $Ma<<1$, the density is
almost uniform throughout the system and the fluid is incompressible. Fig. 6
depicts the variation of Mach number with system radius. We observe that $Ma$
is independent of system radii at a particular temperature. As temperature
increases, the system flow changes from subsonic to supersonic, indicating
possible Mach cone formation.

It was observed that as $\eta/s$ increases, the typical Mach cone structure
smears out and vanishes Bouras2014 ; Bouras2011 . Right panel of Fig. 6
suggests the dependence of Mach number on $\mu_{B}$. We see that it increases
with increasing $\mu_{B}$ and remains above the supersonic limit at all
temperatures for the highest value of $\mu_{B}$ considered here. This might be
an indication of the early formation of a rapidly expanding system with shock
waves at higher $\mu_{B}$ values. From Fig. 6, it can also be inferred that
since Mach number is of the order of unity for all temperatures and at all
values of $\mu_{B}$, the system produced in ultra-relativistic collisions is
compressible Ollitrault2008 .

##  IV Summary

In summary, we have estimated various relevant parameters like Reynolds
number, Knudsen number and Mach number which assist in characterizing many-
particle system. The dependence of these parameters on the system size,
temperature, chemical potential for hadrons produced in ultra-relativistic
nuclear collisions have been investigated. The obtained values ($Kn<<1$,
$Ma\sim 1$ and $Re>>1$) indicate the occurrence of compressible inviscid flows
at high temperatures close to the phase transition region ($T\sim 150-170$
MeV). The degree of thermalization of hadron gas estimated is comparable over
different system sizes, indicating the applicability of hydrodynamics in
interpreting the results from high multiplicity pp to heavy-ion collisions.
This puts into question the applicability of pp collisions as a baseline for
heavy-ion collisions at the high beam energies presently accessible at
different particle accelerators.

## Acknowledgment

This research work has been carried out under financial supports from DAE-
BRNS, the Government of India, Project No. 58/14/29/2019-BRNS of Raghunath
Sahoo. For the research fellowship, Ronald Scaria acknowledges CSIR, Govt. of
India. CRS and RS acknowledge the financial supports under the above BRNS
project. Further R.S. acknowledges the financial supports under the CERN
Scientific Associateship, CERN, Geneva, Switzerland. The authors acknowledge
the Tier-3 computing facility in the experimental high-energy physics
laboratory of IIT Indore supported by the ALICE project.

## References

  * (1) P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005).
  * (2) I. Arsene et al. (BRAHMS Collaboration), Nucl. Phys. A 757, 1 (2005).
  * (3) B. B. Back et al. (PHOBOS Collaboration), Nucl. Phys. A 757, 28 (2005).
  * (4) J. Adams et al. (STAR Collaboration), Nucl. Phys. A 757, 102 (2005).
  * (5) K. Adcox et al. (PHENIX Collaboration), Nucl. Phys. A 757, 184 (2005).
  * (6) J. Adam et al. (ALICE Collaboration), Nat. Phys. 13, 535 (2017).
  * (7) V. Khachatryan et al. (CMS Collaboration), Phys. Lett. B 765, 193 (2017).
  * (8) D. Velicanu [CMS Collaboration], J. Phys. G 38, 124051 (2011).
  * (9) R. Campanini and G. Ferri, Phys. Lett. B 703, 237 (2011).
  * (10) D. Thakur, S. De, R. Sahoo and S. Dansana, Phys. Rev. D 97, 094002 (2018).
  * (11) D. Sahu, S. Tripathy, R. Sahoo and S. K. Tiwari, [arXiv:2001.01252 [hep-ph]].
  * (12) W. T. Deng and X. G. Huang, Phys. Rev. C 93, 064907 (2016).
  * (13) R. S. Bhalerao, J. P. Blaizot, N. Borghini and J. Y. Ollitrault, Phys. Lett. B 627, 49 (2005).
  * (14) H. J. Drescher, A. Dumitru, C. Gombeaud and J. Y. Ollitrault, Phys. Rev. C 76, 024905 (2007).
  * (15) J-Y Ollitrault, Eur. J. Phys. 29, 275 (2008).
  * (16) A. Andronic, P. Braun-Munzinger, J. Stachel, and M. Winn, Phys. Lett. B 718, 80 (2012).
  * (17) V. Vovchenko, D. V. Anchishkin, and M. I. Gorenstein, Phys. Rev. C 91, 024905 (2015).
  * (18) D. H. Rischke, M. I. Gorenstein, H. Stöcker, W. Greiner, Z. Phys. C 51, 485 (1991).
  * (19) R. Hagedorn and J. Ranft, Nuovo Cim. Suppl. 6, 169 (1968).
  * (20) G. P. Kadam and H. Mishra, Phys. Rev. C 92, 035203 (2015).
  * (21) Mirco Cannoni, Phys. Rev. D 89, 103533 (2014).
  * (22) P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145 (1991).
  * (23) S. K. Tiwari, S. Tripathy, R. Sahoo and N. Kakati, Eur. Phys. J. C 78, 938 (2018).
  * (24) G. Kadam and S. Pawar, Adv. High Energy Phys. 2019, 6795041 (2019).
  * (25) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. C 86, 024913 (2012).
  * (26) J. Noronha-Hostler, J. Noronha, and C. Greiner, Phys. Rev. Lett. 103, 172302 (2009).
  * (27) P.A. Zyla et al. (Particle Data Group), PTEP 2020, no.8, 083C01 (2020).
  * (28) P. Braun-Munzinger, I. Heppe, J. Stachel, Phys. Lett. B 465, 15 (1999).
  * (29) K. Redlich and K. Zalewski, arXiv:1611.3746 (2016).
  * (30) A. Bhattacharyya, R. Ray, S. Samanta and S. Sur, Phys. Rev. C 91, 041901(R) (2015).
  * (31) N. Sarkar and P. Ghosh, Phys. Rev. C 96, 044901 (2017).
  * (32) S. Chatterjee, S. Das, L. Kumar, D. Mishra, B. Mohanty, R. Sahoo and N. Sharma, Adv. High Energy Phys. 2015, 349013 (2015).
  * (33) D. Sahu, S. Tripathy, G. S. Pradhan and R. Sahoo, Phys. Rev. C 101, 014902 (2020).
  * (34) S. Acharya et al. [ALICE], Phys. Lett. B 802, 135225 (2020).
  * (35) M. Gorenstein, M. Hauer, O. Moroz, Phys.Rev.C 77, 024911 (2008).
  * (36) G.S. Denicol, C. Gale, S. Jeon, and J. Noronha, Phys. Rev. C 88, 064901 (2013).
  * (37) G. Kadam and H. Mishra, Phys. Rev. D 100, 074015 (2019).
  * (38) S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, K. K. Szabó, Phys. Lett. B 730, 99 (2014).
  * (39) N. A. Tawfik, L. I. Abou-Salem, A. G. Shalaby, M. Hanafy, A. Sorin, O. Rogachevsky and W. Scheinast, Eur. Phys. J. A 52, 324 (2016).
  * (40) P. Braun-Munzinger, D. Magestro, K. Redlich and J. Stachel, Phys. Lett. B 518, 41 (2001).
  * (41) J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
  * (42) A. Khuntia, S. K. Tiwari, P. Sharma, R. Sahoo and T. K. Nayak, Phys. Rev. C 100, 014910 (2019).
  * (43) M. M. Aggarwal et al. (STAR collaboration), Phys. Rev. C 82, 024912 (2010).
  * (44) A. Adare et al. (PHENIX collaboration), Phys. Rev. Lett. 104, 252301 (2010).
  * (45) G. Aad et al. (ATLAS collaboration), Phys. Lett. B. 739, 320 (2014).
  * (46) S. Chatrchyan et al. (CMS collaboration), Phys. Rev. C. 90, 024908 (2014).
  * (47) C. Nattrass, N. Sharma, J. Mazer, M. Stuart and A. Bejnood, Phys. Rev. C. 94, 011901 (2016).
  * (48) G. Sarwar, M. Hasanujjaman, M. Rahaman, A. Bhattacharyya and Jan-e Alam, Phys. Lett. B 820, 136583 (2021).
  * (49) I. Bouras, A. El, O. Fochler, H. Niemi, Z. Xu, C. Greiner, Phys. Lett. B 728, 156 (2014).
  * (50) I Bouras et al, J. Phys.: Conf. Ser. 270, 012012 (2011).


# A Phase Lag Between Disc and Corona in GRMHD Simulations of Precessing
Tilted Accretion Discs

M. Liska1, C. Hesp1, A. Tchekhovskoy2, A. Ingram3, M. van der Klis1 & S.B.
Markoff1
1Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science
Park 904, 1098 XH Amsterdam, The Netherlands
2Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA),
Physics & Astronomy, Northwestern University, Evanston, IL 60202, USA
3Department of Physics, Astrophysics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford, OX1 3RH, UK E-mail<EMAIL_ADDRESS>
(Accepted. Received; in original form)

###### Abstract

In the course of its evolution, a black hole (BH) accretes gas from a wide
range of directions. Given a random accretion event, the typical angular
momentum of an accretion disc would be tilted by $\sim$60∘ relative to the BH
spin. Misalignment causes the disc to precess at a rate that increases with BH
spin and tilt angle. We present the first general-relativistic
magnetohydrodynamic (GRMHD) simulations spanning a full precession period of
highly tilted (60∘), moderately thin ($h/r=0.1$) accretion discs around a
rapidly spinning ($a\simeq 0.9$) BH. While the disc and jets precess in phase,
we find that the corona, sandwiched between the two, lags behind by $\gtrsim
10^{\circ}$. For spectral models of BH accretion, the implication is that hard
non-thermal (corona) emission lags behind the softer (disc) emission, thus
potentially explaining some properties of the hard energy lags seen in Type-C
low frequency quasi-periodic oscillations in X-Ray binaries. While strong jets
are unaffected by this disc-corona lag, weak jets stall when encountering the
lagging corona at distances $r\sim 100$ black hole radii. This interaction may
quench large-scale jet formation.

###### keywords:

accretion, accretion discs – black hole physics – MHD – galaxies: jets –
methods: numerical

††pagerange: A Phase Lag Between Disc and Corona in GRMHD Simulations of
Precessing Tilted Accretion Discs–References††pubyear: 2017

##  1 Introduction

Evidence is growing that many accreting black holes (BHs) have accretion discs
that are misaligned relative to the BH spin. Because supermassive BH spin
magnitude and direction are set by the long history of randomly oriented gas
accretion events and galaxy mergers, misaligned accretion can naturally arise
(Volonteri et al., 2005; King et al., 2005). Indeed, misalignment between
inner and outer disc is found in AGN maser observations (see e.g. Herrnstein
et al. 2005; Caproni et al. 2006, 2007; Greene et al. 2013). For X-ray
binaries (XRBs) the dynamical evolution of the system, which may include
asymmetric supernova kicks and 3-body interactions in stellar clusters, can
also lead to substantial misalignment. Indeed, the jets in XRB GRO J1655–40
are offset by $15^{\circ}$ from the binary plane (Hjellming & Rupen 1995;
Greene et al. 2001; Maccarone 2002).

In view of these observations, it is important to understand how accretion
discs respond to the misalignment relative to the BH spin. An important effect
here is general relativistic (GR) frame dragging, which is associated with the
BH spin and known to induce nodal Lense & Thirring (1918) precession of orbits
inclined to the BH equatorial plane. As the precession frequency decreases
with distance from the BH roughly as $\Omega_{LT}\propto a/r^{3}$, frame
dragging is predicted to warp tilted accretion discs. When the disc is
geometrically thick, i.e. the dimensionless scale height is larger than the
viscosity parameter, $h/r>\alpha$, as would be expected in all low luminosity
AGN and hard state XRBs (Narayan & Yi, 1994), viscous diffusion of the warp
plays a minimal role (Papaloizou & Pringle, 1983) and the warp is communicated
radially through bending waves traveling at approximately half the speed of
sound (Papaloizou & Lin, 1995). When the tilt (measured in radians) is small
compared to $h/r$, a bending wave dominated disc accretes misaligned and forms
smooth radial oscillations in tilt closer to the BH (Ivanov & Illarionov,
1997; Lubow et al., 2002; Fragile et al., 2007; Liska et al., 2018). Since a
larger non-linear tilt, exceeding the disc scale height, is expected in many
systems (King et al., 2005), it is crucial to understand how a disc responds
to such a large tilt.

Due to the absence of a fully non-linear analytic theory in the bending wave
regime, general relativistic magnetohydrodynamic (GRMHD) simulations are an
excellent tool to gain insight into such problems since they include magnetic
fields, which generate the turbulent viscosity (Balbus & Hawley, 1991, 1998)
driving accretion, and they naturally reproduce LT precession. Solid-body
precession of tilted accretion discs has been suggested (Stella & Vietri,
1998; Fragile et al., 2007) as the origin of low frequency quasi-periodic
oscillations (QPOs) observed in XRB lightcurves (e.g., van der Klis 1989).
However, solid-body precession is inconsistent with the observed energy
dependence of QPO frequency (e.g. van den Eijnden et al. 2016; van den Eijnden
et al. 2017), suggesting that a more complex geometry, involving differential
precision between various disc and/or jet components, may be present.

Whereas the relevant timescale in the problem is comparable to the precession
period, GRMHD simulations of such a long duration have usually been
prohibitively expensive. Thick discs, $h/r\sim 0.3$, were shown to viscously
spread and their precession to stall within $40^{\circ}$ (Liska et al., 2018).
Thinner discs have only been simulated at a relatively small tilt of around
$10^{\circ}$ in both the viscous (Liska et al., 2018) and bending wave
(Morales Teixeira et al., 2014) regimes, impeding the study of precession.

Here, we present the first GRMHD simulations of highly tilted accretion discs
in the bending wave, $h/r\gtrsim\alpha$, regime around rapidly spinning BHs
that span the full precession cycle. We describe our numerical setup n section
2, present our results in section 3, and conclude in section 4.

##  2 Numerical models

We use our recently developed massively parallel 3D GRMHD code H-AMR
(pronounced “hammer”, Liska et al. 2018) accelerated by Graphical Processing
Units (GPUs). H-AMR is based on a 2D serial open-source code, HARM2D (Gammie
et al., 2003; Noble et al., 2006), which we rewrote and improved with a
staggered grid for constrained transport of magnetic fields (Gardiner & Stone,
2005), adaptive mesh refinement (AMR, not used adaptively here), an advanced
treatment of the polar region, and a locally adaptive time step. We use a
spherical-polar grid of resolution $1372\times 480\times 738$, which is
uniform in $\log r$ and $\phi$ and approximately uniform in $\theta$. It
resolves the disc scale height, $h/r\sim 0.1$, by $12{-}15$ cells in the $r$-,
$\theta$\- and $\phi$-dimensions. To prevent the cells from becoming elongated
near the polar singularity, we reduce the $\phi$-resolution to $368$ cells at
$15^{\circ}<|\theta|<30^{\circ}$, $184$ cells at
$7.5^{\circ}<|\theta|<15^{\circ}$ and $92$ cells at
$0^{\circ}<|\theta|<7.5^{\circ}$ around both poles. We use outflow boundary
conditions in the $r$-direction, with the inner radial boundary inside the BH
event horizon and the outer boundary at $10^{5}r_{g}$, where $r_{g}=GM/c^{2}$
is the gravitational radius, such that both of the boundaries are causally
disconnected from the flow. We apply transmissive boundary conditions in the
$\theta$-direction that ensure free passage of the plasma through the polar
singularity (Liska et al., 2018).

We initialize the simulations with a torus in hydrostatic equilibrium
(Fishbone & Moncrief, 1976) around a BH with spin $a=0.9375$ in a Kerr-Schild
foliation. We place the torus inner edge at $r_{\rm in}=12.5r_{g}$ and density
maximum, scaled to (initial) $\rho_{\rm max}=1$, at $r_{\rm max}=25r_{g}$. We
use an ideal gas equation of state, $p_{g}=(\Gamma-1)u_{g}$, where $p_{g}$ and
$u_{g}$ are gas thermal pressure and energy densities, and non-relativistic
$\Gamma=5/3$. We tilt the initial torus relative to the BH spin axis (and
grid) by an angle $\mathcal{T}_{\rm init}=60\degree$ (see Liska et al. 2018
for details). We cool the disc to its target thickness, $h/r=0.1$, on the
Keplerian timescale using a source term (Noble et al., 2009; Liska et al.,
2018).

Figure 1: The mass accretion rate ($\dot{M}$) in (a,b) and magnetic fluxes in
the disc ($\Phi_{max}$) and on the BH ($\Phi_{BH}$) in (c,d) as function of
time. While in (e,f) the radiative ($\eta_{rad}$) and wind ($\eta_{wind}$)
efficiency are similar for both models, the jet efficiency ($\eta_{jet}$) is
an order of magnitude higher in S-T60. This can be explained by the presence
of large scale magnetic flux in S-T60.

We consider two models, W-T60 and S-T60, with weak and strong jets,
respectively, whose discs are both tilted by $60^{\circ}$. W-T60 features a
small field loop concentrated near the BH with vector potential
$A_{\phi}\propto(\rho-0.05)$, while S-T60 has a much larger poloidal field
loop, described by $A_{\phi}\propto(\rho-0.05)^{2}r^{3}$. We initially
normalize the field strength by setting $\max p_{\rm g}/\max p_{\rm b}=5$ for
W-T60 and $\max p_{\rm g}/\max p_{\rm b}=100$ for S-T60. The size of the loop
turns out to be more important than the normalization, and S-T60 ends up with
an order of magnitude larger magnetic flux than W-T60.

Figure 2: A colour map of the fluid-frame rest mass density $\log(\rho)$ for
models S-T60 (left panel) and W-T60 (right panel) at $t=4.5\times 10^{4}$
$r_{g}/c$. Magenta lines show the corona-jet boundary, and black lines the
disc-corona boundary. The strong jet in model S-T60 (left panel ) readily
accelerates to relativistic speeds (Fig. 3). However, since the corona in
model W-T60 precesses slower than the disc, it becomes misaligned relative to
the disc and obstructs jet propagation. This causes the jet to stall beyond
$\sim 10^{2}r_{g}$ (Fig. 3).

##  3 Results

Differential rotation shears up magnetic field lines, activating the
magnetorotational instability (MRI; Balbus & Hawley 1991, 1998) and leading to
turbulence, which causes accretion onto the BH. We study the system at late
time, $t\gtrsim 10^{4}r_{g}/c$, after the disc has cooled down from to its
target thickness. Figure 1(a,b) shows that the mass accretion rate, $\dot{M}$,
is similar for both models. Figure 1(c,d) shows the positive (there is no
negative) poloidal magnetic flux on the BH, $\Phi_{BH}=0.5\int_{r=r_{\rm
H}}\left|B^{r}\right|dA_{\theta\phi}$ and in the disc,
$\Phi_{\rm,ax}=\max_{r}\Phi_{p}(r)$ with
$\Phi_{p}(r)=\max_{\theta}\int_{0}^{\theta}B^{r}dA_{\theta\phi}$. Compared to
W-T60, model S-T60 contains an order of magnitude more flux, most of which
threads the event horizon. For both models 3D animations are available (see SI
or this YouTube playlist).

To quantify the conversion efficiency of gravitational potential energy as
matter falls into the black hole we break down in Fig. 1(e,f) the energy
outflow efficiency $\eta_{tot}$, defined as the net energy outflow rate
divided by the BH rest-mass energy accretion rate, $\dot{M}c^{2}$, into
radiative, $\eta_{rad}$, mechanical outflow, $\eta_{wind}$, and jet,
$\eta_{jet}$, constituents that we define below. The radiative efficiency,
$\eta_{rad}$, approximated by integrating the total radiative cooling rate
from the inner photon orbit at $r=1.43r_{g}$ to $100r_{g}$ hovers around the
canonical Novikov & Thorne (1973) value of $18\%$. The jet efficiency,
$\eta_{jet}$, corresponding to the mechanical energy carried by the jets,
reaches peak values of around $0.1\%$ in W-T60 and $10\%$ in S-T60. The wind
efficiency, $\eta_{wind}$, or the mechanical outflow efficiency outside of the
jet, reaches around $20\%$ for both models. The sum $\eta_{tot}\approx 40\%$.

Figure 2 shows a colour map of the density $\rho$ in a vertical slice for
models W-T60 and S-T60. We define the disc-corona boundary at
$\rho=0.025\times\rho_{max}$ and the corona-jet boundary at $p_{b}/\rho
c^{2}=5$. While the disc shapes appear to be similar in both models, the
corona (yellow-green-cyan regions) in model W-T60 shows an asymmetric
elongation in the upper-left and lower-right directions, as seen in Fig. 2(a).
This elongation forms an obstacle for the weak jet, causing it to bend and
eventually stall within a few $\times 100r_{g}$. On the other hand, the strong
jet in S-T60 propagates unhindered through the corona, pushing it aside. The
Lorentz factors, $\gamma$, and mass/energy fluxes, $F_{M}$/$F_{E}$, for both
jets are shown in Fig. 3: while the jet in model S-T60 (solid lines)
accelerates up to $\gamma\sim 7$ and conserves its mass and energy flux, the
jet in model W-T60 (dashed lines) loses energy and mass due to the interaction
with the corona and fails to accelerate to relativistic speeds. Below we show
that this misalignment between corona and jet arises through different rates
of precession resulting in a phase lag.

Figure 4(a,b), shows the evolution of tilt and precession angles (as
calculated in Liska et al. 2018) for the disc, corona and jets, allowing us to
quantitatively study the phase lags. The precession angle increases quasi-
monotonically at approximately the same rate in both models, reaching
$200^{\circ}-400^{\circ}$ by the simulation’s end. This differs from our
previous work (Liska et al., 2018), where the disc stopped precessing after
$40^{\circ}$ due to viscous expansion of the thicker disc. The phase lag
between disc and corona, $\Delta\mathcal{P}$, in Fig. 4(e,f) initially peaks
at $\Delta\mathcal{P}\sim 30^{\circ}-50^{\circ}$ before stabilizing between
$\Delta\mathcal{P}\sim 10^{\circ}-20^{\circ}$. The minimum in both models
coincides with a peak in jet power (Fig. 1(e,f)), suggesting jet (which is
aligned with the inner disc) and corona try to push each other into alignment.
In both models the disc-corona-jet system aligns as a whole with the BH on the
viscous timescale, from the initial $\mathcal{T}=60^{\circ}$ to
$\mathcal{T}\sim 30^{\circ}$ at $t=10^{5}r_{g}/c$. This could be due to
misaligned angular momentum mixing and cancellation between annuli at
different phases of precession (see Fig. 5(c,d) and Liska et al. 2018; see
also Sorathia et al. 2013). This is different from Bardeen & Petterson (1975)
alignment, which occurs on time-scales much shorter than the viscous timescale
and is expected to affect the inner regions of the disc.

Figure 5(a–d) shows the tilt and precession angles of the disc, corona, and
jet as function of radius from the BH. A single radial tilt oscillation
centered at $10r_{g}$ develops as predicted by analytic theory for the bending
wave regime valid at low tilts (Ivanov & Illarionov, 1997; Lubow et al.,
2002). Furthermore, the disc, corona and jet are aligned with each other up to
the disc’s outer edge located at $\sim 100r_{g}$. Beyond $\sim 200r_{g}$, the
corona’s tilt and precession angle oscillate wildly. However, since the
density of the corona is very low there, this region does not contribute
significantly to the corona’s radially averaged tilt and precession angle in
Fig. 4(a–d).

Figure 5(e),(f) shows that both models remain gas pressure dominated
($\beta>1$) down to the event horizon. The Maxwell
($\alpha_{M}=b^{r}b^{\phi}/(p_{g}+p_{b})$, $b^{\mu}$ is the magnetic 4-vector)
and Reynolds ($\alpha_{R}=\rho u^{r}u^{\phi}/(p_{g}+p_{b})$, $u^{\mu}$ is the
4-velocity) stresses in Fig. 5(g,h), calculated in a coordinate system aligned
with the local angular momentum, are similar. However, their sum does not
match the effective viscosity ($\alpha_{eff}=v_{r}v_{k}/c_{s}^{2}$, where
$v_{k}$ is the Keplerian 3-velocity and $c_{s}$ is the sound speed),
indicating that large scale torques drive outflows and remove the disc’s
angular momentum. Due to the rapid mass infall caused by high $\alpha_{eff}$
ions and electrons are expected to decouple into a two-temperature plasma
within the inner disc (Esin et al., 1997), which may make the disc hotter and
the spectrum harder than predicted for a standard Shakura & Sunyaev (1973)
accretion disc (see also Liska et al. 2018).

Figure 3: Radial profiles at $t=4.0\times 10^{4}r_{g}/c$ of the Lorentz
factor, $\gamma$, in (a), and, in arbitrary units, the mass/energy flux,
$F_{M}/F_{E}$, through the jet in (b). The low Lorentz factor ($\gamma\sim 1$)
and drop in mass and energy fluxes at $r\gtrsim(50{-}100)r_{g}$ in model W-T60
(dashed line) indicate that a weak jet stalls and disrupts, because it runs
into the lagging corona. In contrast, acceleration to a relativistic Lorentz
factor ($\gamma\sim 7$) and near-constancy of mass and energy fluxes in model
S-T60 (solid line) indicate that a stronger jet survives the interaction with
the corona and escapes to large radii.

##  4 Discussion and Conclusion

In this work we present the first GRMHD simulations of highly tilted ($60$
degrees), moderately thin ($h/r\approx 0.1$) accretion discs around rapidly
spinning BHs ($a=0.9375$). The 3 times smaller disc thickness compared to our
previous work (Liska et al., 2018) prevents the disc from viscously spreading
on the Lense-Thirring timescale, allowing the disc to complete 0.6-1.4
precession cycles, depending on the model.

The precession-induced quasi-periodic variability in XRB lightcurves could
potentially be exploited to measure BH spin (Stella & Vietri, 1998; Ingram et
al., 2009). However, the precession period is expected to be determined by
both the BH spin and the (initial) disc size and density profile, making such
measurements impossible if the disc’s morphology is poorly constrained. It
would be interesting if a precessing disc formed self-consistently, with a
relatively well-defined size and density profile. One way this could happen is
by tearing up of a large misaligned accretion disc, as seen in hydrodynamic
smoothed particle hydrodynamics (SPH) simulations: the disc tears up at a
radius where the Lense-Thirring torques, which try to tear the disc apart,
overwhelm the viscous torques, which hold the disc together (Nixon et al.,
2012; Nealon et al., 2015). Even though both our GRMHD and those SPH
simulations are in the bending wave dominated regime, neither of our
magnetized discs tore up. This discrepancy might be due to the larger
thickness of our discs (Liska et al. 2019B in prep.).

Figure 4: The angular momentum averaged tilt angle, $\mathcal{T}$, in (a,b)
and precession angle, $\mathcal{P}$, in (c,d) as functions of time, for the
disc, corona, and jet. Model S-T60 and W-T60 both align with the BH over time,
while their precession angle grows semi-monotonically. The phase lag,
$\Delta\mathcal{P}$, between disc and corona in (e,f) first peaks, but seems
to stabilize later.

We found that the magnetic pressure dominated disc wind, which we refer to as
the corona, lags the disc in phase by $10{-}40^{\circ}$ (Fig. 4c). This can
help to interpret the phase lags of Type-C QPOs observed in XRBs (Wijnands et
al., 1999; Reig et al., 2000; Qu et al., 2010; Pahari et al., 2013; van den
Eijnden et al., 2016; van den Eijnden et al., 2017). More specifically, van
den Eijnden et al. (2016) recently used a time-resolved cross-spectral
analysis of GRS 1915+105 and found that the phase lags systematically change
on short timescale in a way indicating that for observations with
$\nu_{QPO}>2$ Hz for $5$-$10$ cycles the QPO in the hard band is slightly
faster than that in the soft band, resulting in a gradually increasing soft
lag, and for observations with $\nu_{QPO}<2$ Hz the QPO in the hard band is
slightly slower than that in the soft band, resulting in a gradually
increasing hard lag. The increase in soft lag for $\nu_{QPO}>2$ Hz was
suggested to result from differential precession, since inner regions are
expected to have a harder spectrum and to precess faster than outer regions
(Stella & Vietri, 1998). However, the increase in hard lag for $\nu_{QPO}<2$
Hz over a QPO period remained a mystery. The evolving coronal lag in our
simulations (i.e. the increase in $\Delta\mathcal{P}$) may provide a natural
solution to this puzzle: the corona, which emits a harder spectrum than the
disc (via inverse Compton scattering), precesses more slowly than the disc,
resulting in the QPO being faster in the soft band. Translating the coronal
lag into actual X-Ray lightcurves requires a fuller treatment of radiation
including Compton scattering, taking into account the complex geometry of this
problem. However, with the present data, the timescale for the phase lag to
build up and its total size seem to be able to match observations. The model
may also be able explain why for $\nu_{QPO}>2$ Hz the phase lag turns soft. As
$\nu_{QPO}$ increases, one would expect the precessing inner disc to become
smaller (Stella & Vietri, 1998; Ingram et al., 2009). Since a smaller disc is
hotter (Shakura & Sunyaev, 1973), its spectrum will become harder and, if the
hardness of the corona does not increase by the same amount, the phase lag
will gradually move from hard to soft.

While the strong jet in S-T60 propagates unhindered, the weak jet in W-T60
disrupts and stalls within a few $\sim 10^{2}r_{g}$, obstructed by the lagging
corona (Figs. 2, 3). However, in reality a precessing jet of any arbitrary
strength may become truncated if the ambient medium is massive enough. Such
jets may transfer significant amounts of energy and angular momentum to this
ambient medium and, by doing so, contribute to weakly collimated matter-
dominated sub-relativistic outflows. In addition, the jet-corona interaction
might lead to reconnection and acceleration of (non-) thermal particles,
possibly producing flares across the electromagnetic spectrum. However, due to
the, presumably, small size of such truncated jets, it may be hard to directly
detect them.

We have shown in previous work (Liska et al., 2018) that BH-launched jets
(Blandford & Znajek, 1977) around gas pressure supported discs are aligned
with the disc, not the BH as was shown for magnetic pressure supported discs
(McKinney et al., 2013). However, our current work demonstrates that on larger
scales the orientation of the jet will be determined by an intricate force
balance between the jet and outer surrounding medium. If the jet contains a
large amount of energy relative to this surrounding medium, the jet will
remain aligned with the inner (precessing) flow. On the other hand, if the jet
is weak, it will be torqued into alignment with the surrounding medium and, as
in model W-T60, may stall. This result suggests caution in the interpretation
of the observed $15^{\circ}$ misalignment between the orbital plane of XRB GRO
J1655-40 and its jet (Hjellming & Rupen, 1995; Orosz & Bailyn, 1997; Greene et
al., 2001; Maccarone, 2002) as a tilt between disc and BH spin.

Figure 5: The tilt angle $\mathcal{T}$ in (a,b) and precession $\mathcal{P}$
angle in (c,d) as function of radius. Strong flux model S-T60 (left panel) is
averaged between $3.5\times 10^{4}<t<4.0\times 10^{4}$ $r_{g}/c$ and weak flux
model W-T60 (right panel) is averaged between $5.5\times 10^{4}<t<6.0\times
10^{4}$ $r_{g}/c$. The disc does not align with the BH (Bardeen & Petterson,
1975), and instead develops a single radial tilt oscillation peaking at
$r\approx 10$ $r_{g}$, which is consistent with analytic theory for a
geometrically thick, $h/r>\alpha$, disc (Ivanov & Illarionov, 1997; Lubow et
al., 2002). The density, $\rho$, and gas-to-magnetic pressure ratio, $\beta$,
in (e,f) are similar for both models. In (g,h) the Maxwell, $\alpha_{M}$, and
Reynolds, $\alpha_{R}$, stresses in the disc increase steeply within the inner
disc ($r<10$ $r_{g}$), while their sum poorly matches the effective viscosity
$\alpha_{eff}$. This may be explained by the presence of outflows carrying
away significant portions of the orbital angular momentum (e.g. Liska et al.
(2018)).

##  5 Acknowledgments

We thank Marta Volonteri for useful feedback. This research was made possible
by NSF PRAC award no. 1615281 and OAC-1811605 at the Blue Waters sustained-
petascale computing project. ML and MK were supported by the NWO Spinoza
Prize, AI by the Royal Society URF, AT thanks Northwestern University for
support and SM by the NWO VICI grant (no. 639.043.513).

##  6 Supporting Information

Additional Supporting Information may be found in the online version of this
article: movie files. See our YouTube playlist for 3D visualizations of all
models.

## References

  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Caproni et al. (2007) Caproni A., Abraham Z., Livio M., Mosquera Cuesta H. J., 2007, MNRAS, 379, 135
  * Caproni et al. (2006) Caproni A., Abraham Z., Mosquera Cuesta H. J., 2006, ApJ, 638, 120
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fishbone & Moncrief (1976) Fishbone L. G., Moncrief V., 1976, ApJ, 207, 962
  * Fragile et al. (2007) Fragile P. C., Blaes O. M., Anninos P., Salmonson J. D., 2007, ApJ, 668, 417
  * Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
  * Gardiner & Stone (2005) Gardiner T. A., Stone J. M., 2005, Journal of Computational Physics, 205, 509
  * Greene et al. (2001) Greene J., Bailyn C. D., Orosz J. A., 2001, ApJ, 554, 1290
  * Greene et al. (2013) Greene J. E., Seth A., den Brok M., Braatz J. A., Henkel C., Sun A.-L., Peng C. Y., Kuo C.-Y., Impellizzeri C. M. V., Lo K. Y., 2013, ApJ, 771, 121
  * Herrnstein et al. (2005) Herrnstein J. R., Moran J. M., Greenhill L. J., Trotter A. S., 2005, ApJ, 629, 719
  * Hjellming & Rupen (1995) Hjellming R. M., Rupen M. P., 1995, Nature, 375, 464
  * Ingram et al. (2009) Ingram A., Done C., Fragile P. C., 2009, MNRAS, 397, L101
  * Ivanov & Illarionov (1997) Ivanov P. B., Illarionov A. F., 1997, MNRAS, 285, 394
  * King et al. (2005) King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MNRAS, 363, 49
  * Lense & Thirring (1918) Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liska et al. (2018) Liska M., Tchekhovskoy A., Ingram A., van der Klis M., 2018, ArXiv e-prints, p. arXiv:1810.00883
  * Lubow et al. (2002) Lubow S. H., Ogilvie G. I., Pringle J. E., 2002, MNRAS, 337, 706
  * Maccarone (2002) Maccarone T. J., 2002, MNRAS, 336, 1371
  * McKinney et al. (2013) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2013, Science, 339, 49
  * Morales Teixeira et al. (2014) Morales Teixeira D., Fragile P. C., Zhuravlev V. V., Ivanov P. B., 2014, ApJ, 796, 103
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, ApJ, 428, L13
  * Nealon et al. (2015) Nealon R., Price D. J., Nixon C. J., 2015, MNRAS, 448, 1526
  * Nixon et al. (2012) Nixon C., King A., Price D., Frank J., 2012, ApJ, 757, L24
  * Noble et al. (2006) Noble S. C., Gammie C. F., McKinney J. C., Del Zanna L., 2006, ApJ, 641, 626
  * Noble et al. (2009) Noble S. C., Krolik J. H., Hawley J. F., 2009, ApJ, 692, 411
  * Novikov & Thorne (1973) Novikov I. D., Thorne K. S., 1973, in Dewitt C., Dewitt B. S., eds, Black Holes (Les Astres Occlus) Astrophysics of black holes.. pp 343–450
  * Orosz & Bailyn (1997) Orosz J. A., Bailyn C. D., 1997, ApJ, 477, 876
  * Pahari et al. (2013) Pahari M., Neilsen J., Yadav J. S., Misra R., Uttley P., 2013, The Astrophysical Journal, 778, 136
  * Papaloizou & Lin (1995) Papaloizou J. C. B., Lin D. N. C., 1995, ARA&A, 33, 505
  * Papaloizou & Pringle (1983) Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181
  * Qu et al. (2010) Qu J. L., Lu F. J., Lu Y., Song L. M., Zhang S., Ding G. Q., Wang J. M., 2010, The Astrophysical Journal, 710, 836
  * Reig et al. (2000) Reig P., Belloni T., van der Klis M., Méndez M., Kylafis N. D., Ford E. C., 2000, The Astrophysical Journal, 541, 883
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Sorathia et al. (2013) Sorathia K. A., Krolik J. H., Hawley J. F., 2013, ApJ, 777, 21
  * Stella & Vietri (1998) Stella L., Vietri M., 1998, ApJ, 492, L59
  * van den Eijnden et al. (2016) van den Eijnden J., Ingram A., Uttley P., 2016, MNRAS, 458, 3655
  * van den Eijnden et al. (2017) van den Eijnden J., Ingram A., Uttley P., Motta S. E., Belloni T. M., Gardenier D. W., 2017, MNRAS, 464, 2643
  * van der Klis (1989) van der Klis M., 1989, Annual Review of Astronomy and Astrophysics, 27, 517
  * Volonteri et al. (2005) Volonteri M., Madau P., Quataert E., Rees M. J., 2005, ApJ, 620, 69
  * Wijnands et al. (1999) Wijnands R., Homan J., van der Klis M., 1999, ApJ, 526, L33