# Lattice field computations via recursive numerical integration

Tobias Hartung111Computation-based Science and Technology Research Center, The
Cyprus Institute, 20 Konstantinou Kavafi Street, 2121 Nicosia, Cyprus and
Department of Mathematics, King’s College London, Strand, London WC2R 2LS,
United Kingdom (email: tobias.hartung@desy.de), Karl Jansen222NIC, DESY
Zeuthen, Platanenallee 6, 15738 Zeuthen, Germany (email: karl.jansen@desy.de),
Frances Y. Kuo333School of Mathematics and Statistics, UNSW Sydney, Sydney NSW
2052, Australia (email: f.kuo@unsw.edu.au),
Hernan Leövey444Structured Energy Trading, AXPO Trading & Sales, Parkstrasse
23, 5400 Baden, Germany (email: hernaneugenio.leoevey@axpo.com), Dirk
Nuyens555Department of Computer Science, KU Leuven, Celestijnenlaan 200A, 3001
Leuven, Belgium (email: dirk.nuyens@cs.kuleuven.be), Ian H. Sloan666School of
Mathematics and Statistics, UNSW Sydney, Sydney NSW 2052, Australia (email:
<EMAIL_ADDRESS>
(December 2021
Submitted to the Proceedings of the 38th International Symposium on Lattice
Field Theory, 26th-30th July 2021, Zoom/Gather at Massachusetts Institute of
Technology)

###### Abstract

We investigate the application of efficient recursive numerical integration
strategies to models in lattice gauge theory from quantum field theory. Given
the coupling structure of the physics problems and the group structure within
lattice cubature rules for numerical integration, we show how to approach
these problems efficiently by means of Fast Fourier Transform techniques. In
particular, we consider applications to the quantum mechanical rotor and
compact $\mathit{U}(1)$ lattice gauge theory, where the physical dimensions
are two and three. This proceedings article reviews our results presented in
J. Comput. Phys 443 (2021) 110527 [17].

##  1 Introduction

_Lattice Field Theory_ (LFT) serves as a non-perturbative tool to regulate
Feynman’s path integral, by removing infinities in the infra-red and ultra-
volet. In LFT continuum models in physics are formulated on a finite Euclidean
space-time lattice with lattice spacing $h$. In this way, an ultra-violet cut-
off, $1/h$, is provided through a non-vanishing value of the lattice spacing.
At the same time, the finite volume equips LFTs with an infra-red cut-off
providing thus to a well-defined path integral, see e.g., [14] for an
introduction to lattice field theories. (See [16, 22] for an alternative
approach using the $\zeta$-regularization.)

Another important aspect of Euclidean LFT is that the path integral and
physical observables can be computed numerically within this framework. One
very prominent example of a LFT is _Lattice Quantum Chromodynamics_ (LQCD)
which is the lattice version of _Quantum Chromodynamics_ (QCD), the theory of
the interaction between quarks and gluons. LQCD has been extremely successful
in computing e.g., the low lying hadron spectrum, non-perturbative matrix
elements and form factors, fundamental parameters of QCD and non-zero
temperature physics, see e.g., the _Flavour Lattice Averaging Group_ (FLAG)
review [5].

The main tool for simulating LQCD are _Markov Chain Monte Carlo_ (MCMC)
methods such as Hybrid Monte Carlo, see [25] for an overview. The employed
methods and algorithms have been immensely improved over the last years,
leading to large factors of speedup in the numerical calculations. In
addition, the structure of LQCD and of many other LFTs is of local nature,
coupling essentially only nearest neighbours on the lattice. This allows for
massively parallel simulations on state of the art supercomputers with
(hundred) thousands of processors, which led to computations on lattices of
size close to physical relevance, meaning that the first two quark
generations, i.e., the up, down, strange and charm quarks, assume their
physical values as determined from experiments.

Despite this enormous success, the MCMC approach is unable to address problems
which lead to a complex integrand in the path integral, a situation which is
referred to as sign problem [30]. These problems are of fundamental nature in
the context of high energy physics. In particular, they are related to the
questions: Why is there more matter than anti-matter in the universe? Why is
there more experimentally detected CP violation than the standard model
predicts? What are the physics of topological systems?

In addition, MCMC methods suffer from the problem of auto-correlations which
is an inherent property of MCMC algorithms. This problem becomes severe when
the continuum limit in LQCD is taken, i.e., when the lattice spacing is shrunk
towards zero, see [2] and the references therein. In particular, it has been
noticed [28] that the topological charge in LQCD starts to freeze in Monte
Carlo simulations, making it thus very hard to take the continuum limit.

The above sketched shortcomings of MCMC methods has motivated researchers
working in LFTs to constantly look for alternative approaches. The very large
numerical demand of lattice simulations has triggered new approaches to
accelerate LFT calculations even more and to achieve a very high accuracy in
evaluating the path integral. The sign problem has led to new developments in
using numerical techniques such as tensor networks [6] or tensor
renormalization group methods [1].

In this proceedings article we report on a particular attempt to address the
above mentioned problems by using efficient recursive integration techniques,
combined with a Fourier analysis, following our original work in [17]. The
idea is illustrated and applied to an abelian $\mathit{U}(1)$ compact lattice
gauge theory in two dimensions. By applying _lattice cubature rules_ (see
e.g., [12, 29]), the underlying group structure of this cubature allows for
efficient computation when plugged into integrands exhibiting special
structure as in the physics problems considered here. Although we will
demonstrate the method here in two dimensions only, we see the prospect to
generalize our method to higher dimensions and also to non-abelian groups in
our future work.

To be more concrete, we will initially consider $L$-dimensional integrals (the
“dimensionality” here refers to the number of integration variables rather
than the space/time dimensions) of the form

$\displaystyle\int_{D^{L}}\prod_{i=0}^{L-1}f_{i}\big{(}x_{i},x_{i+1}\big{)}\,\mathrm{d}\boldsymbol{x},\qquad\mbox{with
$\boldsymbol{x}=(x_{0},\ldots,x_{L-1})$ and $x_{L}\equiv x_{0}$,}$ (1)

where each variable $x_{i}$ belongs to a bounded domain $D\subset{\mathbb{R}}$
and each function $f_{i}$ depends only on two consecutive variables $x_{i}$
and $x_{i+1}$ – which is called a first order coupling – and we implicitly
assume periodic boundary conditions. 1 sketches first order as well as third
order coupling.

(a) Coupling of first order in 1D $(r=1)$$x_{0}$$x_{i}$$x_{i+1}$$x_{L}\equiv
x_{0}$(b) Coupling of third order in 1D
$(r=3)$$x_{0}$$x_{i}$$x_{i+1}$$x_{i+2}$$x_{i+3}$$x_{L}\equiv x_{0}$ Figure 1:
Different couplings in LQCD (Figure taken from [17, Fig. 1])

As already said above, we propose to solve lattice systems with such a
coupling structure by approximating the involved integrals employing recursive
numerical integration methods, where the lattice cubature rules used belong to
a special family among _Quasi-Monte Carlo_ rules (see e.g., [12, 13, 21, 24,
26, 27, 29]). The work we describe here is part of an ongoing effort by some
of us to apply alternative to MCMC mathematical methods to tackle lattice
models; previous works employ the Quasi-Monte Carlo approach [23],
polynomially exact integration rules [4, 18], and a first applications of the
recursive numerical integration technique to address quantum mechanical models
[3, 17]. Quasi-Monte Carlo methods have also been considered in [9, 11] for
multiloop calculations in perturbation theory.

In the next section we will provide the physical models considered. In the
subsequent sections we will outline the efficient methods and the
corresponding computational complexities. Executable Julia codes and numerical
results can be found in the original work [17].

##  2 Description of physical models

In this section we will introduce the quantum mechanical rotor [3, 7, 8] and a
$2$-dimensional compact $\mathit{U}(1)$ lattice gauge theory, see e.g., [14]
for an introduction into lattice field theories.

###  2.1 Quantum rotor

The quantum rotor is a quantum mechanical model which describes a particle
with mass $m_{0}$ moving on a circle with radius $r_{0}$, see e.g., [3, 7, 8].
We can supply the particle with a moment of inertia $I=m_{0}r_{0}^{2}$. The
quantum rotor can be formulated on a time lattice with lattice spacing $h=T/L$
between neighbouring lattice points, where $T$ is the final time and $L$ is
the number of lattice points. The variables of the system are taken to be
angles $\phi_{i}\in D=[-\pi,\pi)$ defined at each time lattice point. In the
continuum, the action of the system is given by

$\displaystyle
S(\phi)=\int_{0}^{T}\frac{I}{2}\Big{(}\frac{\mathrm{d}\phi}{\mathrm{d}t}\Big{)}^{2}\mathrm{d}t.$
(2)

The discretized version of the action of the quantum rotor, using
$\frac{1}{2}(\frac{\mathrm{d}\phi}{\mathrm{d}t})^{2}\approx\frac{1}{2}(\frac{\phi_{i+1}-\phi_{i}}{h})^{2}\approx\frac{1-\cos(\phi_{i+1}-\phi_{i})}{h^{2}}$,
is then given by

$\displaystyle
S[\boldsymbol{\phi}]\,=\,\frac{I}{h^{2}}\sum_{i=0}^{L-1}\big{(}1-\cos(\phi_{i+1}-\phi_{i})\big{)}$

and we impose periodic boundary conditions. The choice of the action will
allow us to employ a combination of techniques as detailed below to make use
of the Fast Fourier transform (FFT), reducing significantly the computational
cost of the problem. As an observable we will take

$\displaystyle
O[\boldsymbol{\phi}]\,=\,\cos(\phi_{k+1}-\phi_{k})\quad\mbox{for any $k$},$

which is evaluated through the path integral as

$\displaystyle\langle
O[\boldsymbol{\phi}]\rangle\,=\,\frac{\int_{D^{L}}\cos\big{(}\phi_{k+1}-\phi_{k})\,\exp(\beta\sum_{i=0}^{L-1}\cos(\phi_{i+1}-\phi_{i})\big{)}\,\mathrm{d}\boldsymbol{\phi}}{\int_{D^{L}}\exp\big{(}\beta\sum_{i=0}^{L-1}\cos(\phi_{i+1}-\phi_{i})\big{)}\,\mathrm{d}\boldsymbol{\phi}},$
(3)

where $\beta=\frac{I}{h^{2}}=\frac{IL^{2}}{T^{2}}$. In order to evaluate the
above integrals we will in the following convert their domain into the unit
cube $[0,1)^{L}$ since the integration domain $[0,1)^{L}$ is the classical
domain where Quasi-Monte Carlo rules are defined.

###  2.2 Quantum compact abelian gauge theory

Gauge theories are at the heart of the standard model of high energy physics.
They are constructed such that they exhibit local gauge invariance, i.e.,
physical observables and the action of gauge models are invariant under local
changes of the field variables defining a gauge theory. Using gauge invariant
theories in constructing the standard model has been extremely successful. In
particular, with the discovery of the Higgs boson, all particles predicted by
the standard model have been identified experimentally. Thus, we now have a
complete microscopic description of the interaction of the fundamental
particles building all matter.

There are a number of non-perturbative aspects of gauge theories, such as
confinement, topological effects and the existence of glueballs. The main
tool, as formulated by Wilson [31] to understand these non-perturbative
phenomena is Lattice Gauge Theory (LGT). On the lattice, the fundamental gauge
fields are group valued and mainly taken from the abelian group
$\mathit{U}(1)$, or the non-abelian ones, $\mathit{U}(N)$ or $\mathit{SU}(N)$
with $N\geq 2$.

Ideally, a study of gauge theories in three space and one time dimension
($3+1$ dimensions) would be preferable. However, this is presently not
affordable within our approach and we therefore resort to a lower
($1+1$)-dimensional model, the ($1+1$)-dimensional compact $\mathit{U}(1)$
lattice gauge theory. The path integral of this model reads

$\displaystyle\int_{D^{2L^{2}}}\exp\bigg{(}\beta\sum_{i=0}^{L-1}\sum_{j=0}^{L-1}\cos\Big{(}\phi^{a}_{i,j}+\phi^{b}_{i+1,j}-\phi^{a}_{i,j+1}-\phi^{b}_{i,j}\Big{)}\bigg{)}\,\mathrm{d}\boldsymbol{\phi}.$
(4)

Here $D=[-\pi,\pi]$, and we have parametric periodicity where all indices
should be taken modulo $L$. A typical observable to be considered is the
plaquette expectation value, see 2(a),

$\displaystyle\cos\Big{(}\phi^{a}_{0,0}+\phi^{b}_{1,0}-\phi^{a}_{0,1}-\phi^{b}_{0,0}\Big{)},$

which we also study here.

(a) 2D plaquette
$P_{i,j}^{a,b}$$\phi_{i,j}^{a}$$-\phi_{i,j+1}^{a}$$-\phi_{i,j}^{b}$$\phi_{i+1,j}^{b}$$(i,j)$$(i+1,j+1)$(b)
3D plaquettes  $P_{i,j,k}^{a,b}$,  $P_{i,j,k}^{a,c}$, and
$P_{i,j,k}^{b,c}$,$(i,j,k)$$(i+1,j+1,k+1)$ Figure 2: How to construct a
plaquette in two and three dimensions
(Figure taken from [17, Fig. 2]) In the two dimensional case, we construct the
plaquette $P_{i,j}^{a,b}$ at lattice point $(i,j)$ by carrying the field along
the path $(i,j)\to(i+1,j)\to(i+1,j+1)\to(i,j+1)\to(i,j)$. The contribution of
the angles $\phi$ is of positive sign if moving right or up, and the
contribution is of negative sign if moving left or down. We index left/right
movements with a superscript $a$ and up/down movements with superscript $b$.
For the three dimensional case, we add an extra forward/backward movement,
which is denoted by the superscript $c$. This forms three independent smallest
loops that can be taken: the $(a,b)$-plane (black solid line), the
$(a,c)$-plane (red dotted line), and the $(b,c)$-plane (blue dashed line). One
can see that each of these planes contribute equally with a “2-dimensional”
plaquette. All visualized links are straight lines, they are only bent for
better visualization.

The model above uses the plaquette – the smallest Wilson loop – to define the
action. Alternative forms of the action, employing larger Wilson loops can be
considered and are actually used in lattice simulations as they can improve
the extrapolation behavior to the continuum limit, i.e., the limit where the
lattice spacing is sent to zero. The only requirement is that these Wilson
loops are gauge invariant and that the corresponding lattice action converges
to the continuum action when the lattice spacing is sent to zero. Our method
generalizes straightforwardly to such extended Wilson loops and in our
original work [17] there is a detailed discussion of such a situation.

The construction of a $\mathit{U}(1)$ compact lattice gauge theory can be
generalized to three dimensions as illustrated in 2(b). The path integral then
reads

$\displaystyle\int_{D^{3L^{3}}}\exp\bigg{(}\beta\sum_{i=0}^{L-1}\sum_{j=0}^{L-1}\sum_{k=0}^{L-1}\bigg{[}\cos\Big{(}\phi^{a}_{i,j,k}-\phi^{a}_{i,j+1,k}-\phi^{b}_{i,j,k}+\phi^{b}_{i+1,j,k}\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\cos\Big{(}\phi^{c}_{i,j,k}-\phi^{c}_{i+1,j,k}-\phi^{a}_{i,j,k}+\phi^{a}_{i,j,k+1}\Big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+\cos\Big{(}\phi^{b}_{i,j,k}-\phi^{b}_{i,j,k+1}-\phi^{c}_{i,j,k}+\phi^{c}_{i,j+1,k}\Big{)}\bigg{]}\bigg{)}\,\mathrm{d}\boldsymbol{\phi},$

while for the observable we again take the plaquette expectation value

$O[\boldsymbol{\phi}]\,=\,\cos\Big{(}\phi^{a}_{0,0,0}-\phi^{a}_{0,1,0}-\phi^{b}_{0,0,0}+\phi^{b}_{1,0,0}\Big{)}.$

In a more general way, the plaquette variable at a point $(i,j)$ can be
written as

$\displaystyle
P_{i,j}^{a,b}=\Phi_{i,j}^{a}\Phi_{i+1,j}^{b}(\Phi_{i,j+1}^{a})^{-1}(\Phi_{i,j}^{b})^{-1}\equiv
e^{\mathrm{i}\left(\phi_{i,j}^{a}+\phi_{i+1,j}^{b}-\phi_{i,j+1}^{a}-\phi_{i,j}^{b}\right)}.$
(5)

This generalized form allows us to formulate the path integral and the gauge
invariant observables by taking $\Phi_{i,j,\ldots}^{\alpha}$ from
$\mathit{U}(N)$ or $\mathit{SU}(N)$ and taking the trace as well as real parts
of (5). For the physically relevant cases of the standard model, namely the
weak interaction, or QCD (quantum chromodynamics or strong nuclear force), the
relevant groups are $\mathit{SU}(2)$ and $\mathit{SU}(3)$, respectively.

##  3 Recursive numerical integration and complexity results

We present in this section a strategy to approximate the integral (1) using
iterated integration and then replacing each of the iterated integrals by a
numerical quadrature/cubature. We refer to this method as recursive numerical
integration. The method first uses the fact that the integral (1) can be
rewritten as

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D}\cdots\int_{D}f_{0}\big{(}x_{0},x_{1}\big{)}f_{1}\big{(}x_{1},x_{2}\big{)}f_{2}\big{(}x_{2},x_{3}\big{)}\cdots
f_{L-1}\big{(}x_{L-1},x_{0}\big{)}\,\mathrm{d}x_{0}\cdots\mathrm{d}x_{L-1}$
$\displaystyle\,=\,\int_{D}\bigg{[}\int_{D}\cdots\left(\int_{D}\left(\int_{D}f_{0}\big{(}x_{0},x_{1}\big{)}f_{1}\big{(}x_{1},x_{2}\big{)}\,\mathrm{d}x_{1}\right)f_{2}\big{(}x_{2},x_{3}\big{)}\,\mathrm{d}x_{2}\right)\cdots
f_{L-1}\big{(}x_{L-1},x_{0}\big{)}\,\mathrm{d}x_{L-1}\bigg{]}\,\mathrm{d}x_{0},$
(6)

with $\boldsymbol{x}=(x_{0},\ldots,x_{L-1})$ and $x_{L}\equiv x_{0}$. Due to
the periodicity assumption of the indices of the variables, i.e, $x_{L}\equiv
x_{0}$, we can start the iterated integration by choosing any variable of
preference. The integrand presented here exhibits what is called a first order
coupling, which means that each variable $x_{i}$ is coupled via the physical
model to its nearest neighboring variables $x_{i-1}$ and $x_{i+1}$ in a way
that only these three variables are present in the iterated integration of the
variable $x_{i}$.

Recursive integration has been considered in e.g., [10, 15, 19, 20], but not
for integrands with parametric periodicity $x_{L}\equiv x_{0}$.

If we have an observable function that (i) depends only on one variable, or
(ii) depends on two consecutive variables as in the factors $f_{i}$, or (iii)
depends on many variables but takes the form of a product of first order
couplings as for the integrand in (1), we will also arrive at an integral that
can be expressed iteratively in the same manner as in (3). For all physical
models considered in this proceedings article, we have observable functions
that affect at most two consecutive variables (or can be decomposed into a
direct sum of such functions) and therefore we can condense the effect of the
observable function into a new integrand factor which, without loss of
generality, can be included with the index $0$ and be written as
$f_{0}^{\star}\big{(}x_{0},x_{1}\big{)}:=O[x_{0},x_{1}]f_{0}\big{(}x_{0},x_{1}\big{)}$.
In the coming sections, when the observable function is included in the
integration, we will skip the ⋆ notation and just mention that the integrand
$f_{0}$ may be different from the other integrands $f_{i}$ involved.

When the iterated integrals are replaced by the same approximating quadrature
rule with points $t_{0},\ldots,t_{n-1}\in D$ and weights
$w_{0},\ldots,w_{n-1}\in{\mathbb{R}}$, the resulting approximation scheme
takes the form

$\displaystyle{\mathcal{Q}}$
$\displaystyle\,=\,\sum_{k_{0}=0}^{n-1}w_{k_{0}}\bigg{[}\sum_{k_{L-1}=0}^{n-1}w_{k_{L-1}}\cdots\bigg{(}\sum_{k_{2}=0}^{n-1}w_{k_{2}}\bigg{(}\sum_{k_{1}=0}^{n-1}w_{k_{1}}f_{0}\big{(}t_{k_{0}},t_{k_{1}}\big{)}f_{1}\big{(}t_{k_{1}},t_{k_{2}}\big{)}\bigg{)}f_{2}\big{(}t_{k_{2}},t_{k_{3}}\big{)}\bigg{)}$
$\displaystyle\hskip 284.52756pt\cdots
f_{L-1}\big{(}t_{k_{L-1}},t_{k_{0}}\big{)}\bigg{]}.$ (7)

The recursive approximation procedure is then carried out as follows. Define
$M_{i}$ to be the $n\times n$ matrix with entries

$\displaystyle(M_{i})_{p,q}\,=\,f_{i}(t_{p},t_{q})\qquad\mbox{for}\quad
p,q=0,\ldots,{n-1},$ (8)

and let $W$ denote the $n\times n$ diagonal matrix with the weights
$w_{0},\ldots,w_{n-1}$ on the diagonal. Now we can express the innermost sum
in (3) as

$\displaystyle\sum_{k_{1}=0}^{n-1}w_{k_{1}}f_{0}\big{(}t_{k_{0}},t_{k_{1}}\big{)}f_{1}\big{(}t_{k_{1}},t_{k_{2}}\big{)}$
$\displaystyle\,=\,\sum_{k_{1}=0}^{n-1}(M_{0})_{k_{0},k_{1}}w_{k_{1}}^{1/2}\,w_{k_{1}}^{1/2}(M_{1})_{k_{1},k_{2}}$
$\displaystyle\,=\,\sum_{k_{1}=0}^{n-1}(M_{0}W^{1/2})_{k_{0},k_{1}}\,(W^{1/2}M_{1})_{k_{1},k_{2}}\,=\,(M_{0}WM_{1})_{k_{0},k_{2}}.$

Further we have

$\displaystyle\sum_{k_{2}=0}^{n-1}w_{k_{2}}\bigg{(}\sum_{k_{1}=0}^{n-1}w_{k_{1}}f_{0}\big{(}t_{k_{0}},t_{k_{1}}\big{)}f_{1}\big{(}t_{k_{1}},t_{k_{2}}\big{)}\bigg{)}f_{2}\big{(}t_{k_{2}},t_{k_{3}}\big{)}$
$\displaystyle\,=\,\sum_{k_{2}=0}^{n-1}w_{k_{2}}(M_{0}WM_{1})_{k_{0},k_{2}}(M_{2})_{k_{2},k_{3}}$
$\displaystyle\,=\,(M_{0}WM_{1}WM_{2})_{k_{0},k_{3}}.$

This leads to

$\displaystyle{\mathcal{Q}}$
$\displaystyle\,=\,\sum_{k_{0}=0}^{n-1}w_{k_{0}}(M_{0}WM_{1}WM_{2}W\cdots
M_{L-1})_{k_{0},k_{0}}$
$\displaystyle\,=\,\sum_{k_{0}=0}^{n-1}(W^{1/2}M_{0}WM_{1}WM_{2}W\cdots
M_{L-1}W^{1/2})_{k_{0},k_{0}}$ $\displaystyle\,=\,{\rm
trace}(B),\qquad\mbox{with}\qquad B\,=\,W^{1/2}M_{0}WM_{1}WM_{2}W\cdots
M_{L-1}W^{1/2}.$ (9)

The equation (3) defines what we call the recursive numerical integration
approach. It shows a clear advantage over the naive approach that could be
used for approximation in (3) since the naive approach can be interpreted as a
product rule over the $L$ dimensional domain. If the quadrature rules has $n$
points, this implies that the naive approach would have a cost equal to
$n^{L}$, i.e., exponential in the dimension $L$. On the contrary, applying the
recursive numerical integration approach we have an evaluation procedure for
the approximation of the integral that reduces the complexity significantly
and takes advantage of the particular structure of the matrices $M_{i}$. The
matrix structure is of course a consequence of the structure of the underlying
integrand factors $f_{i}$ and the chosen quadrature rule.

1 reviews the cases considered in the work [17] for different possibilities of
matrices $M_{i}$ one can find in LQCD applications.

There are particular conditions on the integrand and quadrature rule that lead
to great improvements in terms of computational cost. In particular, if

  1. 1.

each function $f_{i}$ depends only on the difference of the two arguments,
i.e., $f_{i}(u,v)=\kappa_{i}(v-u)$ for some function
$\kappa_{i}:D\to{\mathbb{R}}$, and

  2. 2.

each function $\kappa_{i}$ is _periodic_ , and

  3. 3.

we have equally spaced points with equal weights $1/n$ (i.e., we have the
rectangle rule),

then the matrix $M_{i}$ is _circulant_ and therefore FFT can be used to obtain
the eigenvalues, see Scenarios (A5)–(A7) in 1.

If we replace the domain $D$ in (1) by an $s$-dimensional domain $D^{s}$,

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D^{s}}\cdots\int_{D^{s}}\prod_{i=0}^{L-1}f_{i}\big{(}\boldsymbol{x}_{i},\boldsymbol{x}_{i+1}\big{)}\,\mathrm{d}\boldsymbol{x}_{0}\cdots\mathrm{d}\boldsymbol{x}_{L-1},$
(10)

where $\boldsymbol{x}_{i}=(x_{i,0},\ldots,x_{i,s-1})\in D^{s}$, with

$x_{i,j}\equiv x_{i\bmod L,\;j\bmod s}\qquad\mbox{for all}\quad
i,j\in{\mathbb{N}},$

then the one-dimensional quadrature rule in (3) becomes an $s$-dimensional
cubature rule with points $\boldsymbol{t}_{0},\ldots,\boldsymbol{t}_{n-1}\in
D^{s}$ and weights $w_{0},\ldots,w_{n-1}\in{\mathbb{R}}$, and the matrices
$M_{i}$ in (8) become

$\displaystyle(M_{i})_{p,q}\,=\,f_{i}(\boldsymbol{t}_{p},\boldsymbol{t}_{q})\qquad\mbox{for}\quad
p,q=0,\ldots,{n-1}.$

Scenarios (A0)–(A4) from 1 apply again in this case. The following conditions
are sufficient to ensure circulant matrices for the more favorable Scenarios
(A5)–(A7):

  1. 1.

each function $f_{i}$ depends only on the difference of the two arguments,
i.e.,
$f_{i}(\boldsymbol{u},\boldsymbol{v})=\kappa_{i}(\boldsymbol{v}-\boldsymbol{u})$
for some function $\kappa_{i}:D^{s}\to{\mathbb{R}}$, and

  2. 2.

each function $\kappa_{i}$ is _periodic with respect to each of the $s$
components_, and

  3. 3.

we have a _lattice cubature rule_ with points

$\boldsymbol{t}_{k}=\frac{k\boldsymbol{z}\bmod n}{n}\qquad\mbox{for}\quad
k=0,\ldots,n-1,$

and equal weights $1/n$.

A lattice cubature rule has an additive group structure. This means that the
difference of two lattice points is another lattice point. Since the lattice
cubature rule has equal weights $1/n$, we obtain circulant matrices $M_{i}$.
This particular fact is the main motivation for favoring lattice cubature
rules above all other cubature rules. The cost in all scenarios is
_independent of $s$_. Nevertheless, the error is ${\mathcal{O}}(n^{-\alpha})$,
where $\alpha$ is determined by the cubature rule and the implied constant may
depend on $s$.

Table 1: Recursive numerical integration for first order couplings (Table
taken from [17, Table 1])

  * •

$M_{i}$ is the $n\times n$ matrix of $f_{i}$ at quadrature points.

  * •

$W$ is an $n\times n$ diagonal matrix with quadrature weights on the diagonal.

  * •

eig returns a diagonal matrix of eigenvalues.

  * •

fft takes the first column of a circulant matrix and returns a diagonal matrix
of eigenvalues.

  * •

The quadrature error for all cases is ${\mathcal{O}}(n^{-\alpha})$, with
$\alpha$ given by the performance of the quadrature rule.

  * •

The scenarios and strategies extend to an $L$-fold product of $s$-dimensional
integrals with the quadrature rule replaced by an $s$-dimensional cubature
rule. The cost remains independent of $s$. The error is again
${\mathcal{O}}(n^{-\alpha})$, with $\alpha$ given by the performance of the
cubature rule, and with an implied constant dependent on $s$.

Scenario | Strategy | Cost
---|---|---
(A0) naive implementation | $\begin{array}[]{l}{\mathcal{Q}}=\mbox{direct product calculation}\end{array}$ | $n^{L}$
(A1) recursive integration | $\begin{array}[]{l}B=W^{1/2}M_{0}WM_{1}W\cdots M_{L-1}W^{1/2}\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}B_{k,k}\end{array}$ | $L\,n^{3}$
(A2) $M_{i}=M$ | $\begin{array}[]{l}A=W^{1/2}MW^{1/2}\\\ B=A^{L}\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}B_{k,k}\end{array}$ | $\log(L)\,n^{3}$
(A3) $M_{i}=M$ diagonalizable | $\begin{array}[]{l}A=W^{1/2}MW^{1/2}\\\ \Lambda={\tt eig}(A)\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}\Lambda_{k,k}^{L}\end{array}$ | $n^{3}$
(A4) $M_{i}=M$ except $M_{0}$ | $\begin{array}[]{l}A=W^{1/2}MW^{1/2}\\\ B=W^{1/2}M_{0}W^{1/2}A^{L-1}\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}B_{k,k}\end{array}$ | $\log(L)\,n^{3}$
(A5) $M_{i}$ circulant | $\begin{array}[]{l}\Lambda_{i}={\tt fft}(M_{i}/n)\mbox{ for each $i$}\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}\prod_{i=0}^{L-1}(\Lambda_{i})_{k,k}\end{array}$ | $L\,n\log(n)$
(A6) $M_{i}=M$ circulant | $\begin{array}[]{l}\Lambda={\tt fft}(M/n)\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}\Lambda_{k,k}^{L}\end{array}$ | $n\log(n)$
| (A7) $M_{i}=M$ except $M_{0}$
---
all circulant
$\begin{array}[]{l}\Lambda_{0}={\tt fft}(M_{0}/n)\\\ \Lambda={\tt fft}(M/n)\\\ {\mathcal{Q}}=\sum_{k=0}^{n-1}(\Lambda_{0})_{k,k}\,\Lambda_{k,k}^{L-1}\end{array}$ | $n\log(n)$

These strategies have been extended in [17] to higher order couplings of order
$r$, of the form

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D^{L}}\prod_{i=0}^{L-1}f_{i}\big{(}x_{i},x_{i+1},\ldots,x_{i+r}\big{)}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\,=\,\int_{D}\cdots\int_{D}f_{0}\big{(}x_{0},x_{1},\ldots,x_{r}\big{)}f_{1}\big{(}x_{1},x_{2},\ldots,x_{r+1}\big{)}\cdots
f_{r}\big{(}x_{r},x_{r+1},\ldots,x_{2r}\big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdots
f_{L-1}\big{(}x_{L-1},x_{0},x_{1},\ldots,x_{r-1}\big{)}\,\mathrm{d}x_{0}\cdots\mathrm{d}x_{L-1}.$

The trick is to group successive $r$ functions into a new factor to form a
product of $L/r$ factors (assuming for simplicity here that $L$ is a multiple
of $r$), and then apply the recursive strategies to these new factors. This
yields analogous results to Scenarios (A1)–(A7), which are denoted by
Scenarios (B1)–(B7) in [17, Table 2]. In all scenarios the error is again
${\mathcal{O}}(n^{-\alpha})$, with $\alpha$ determined by the
quadrature/cubature rule, and the implied constant now depends on $r$.

##  4 Applications

In this section we consider applications of recursive numerical integration to
the quantum rotor and the 2D compact $\mathit{U}(1)$ LGT problem. The
interested reader can find numerical experiments for these two models in our
recent work [17].

###  4.1 The Quantum rotor

For the quantum rotor, the integrals for the numerator and denominator of the
normalized path integral ratio (3) are of the form

$\int_{D^{L}}\prod_{i=0}^{L-1}f_{i}(x_{i+1}-x_{i})\,\mathrm{d}\boldsymbol{x},$

where we have, after a change of variables, $D=[0,1]$ and

$f_{i}(x)=f(x)=\exp(\beta\cos(2\pi x))\qquad\mbox{for all }i=0,\ldots,L-1.$

The only exception is that in the numerator integral we replace $f_{0}$ by

$f_{0}(x)=\cos(2\pi x)\,\exp(\beta\cos(2\pi x)).$

With a slight abuse of notation comparing with (3), here we have the (special)
case $f_{i}(u,v)=\kappa_{i}(v-u)\equiv f_{i}(v-u)$, which means that each
integrand factor $f_{i}(u,v)$ can be considered in fact as a function of one
single variable by taking the difference of the two arguments $u,v$. The
resulting functions $\kappa_{i}\equiv f_{i}$ are periodic and therefore we
know that for the rectangle quadrature rule we have Scenario (A7) in 1, and
thus the numerical integration cost becomes ${\mathcal{O}}(n\log(n))$.

The interesting case of higher order couplings for the quantum rotor arise
when we consider higher order finite difference approximations of the physical
discretization of the action in continuum (2). For example, the central
difference formula $(-x_{i+2}+8x_{i+1}-8x_{i-1}+x_{i-2})/(12h)$ of order
$h^{4}$ leads to (now with $\beta=IL^{2}/(144T^{2})$)

$\int_{D^{L}}\prod_{i=0}^{L-1}f\big{(}-x_{i+2}+8x_{i+1}-8x_{i-1}+x_{i-2}\big{)}\,\mathrm{d}\boldsymbol{x},$

which has order $r=4$. This case can be solved using $4$-dimensional lattice
cubature rules combined with recursive numerical integration, following [17,
Scenario (B4) in Table 2].

###  4.2 The 2D compact $\mathit{U}(1)$ lattice gauge theory model

Now we consider the model (4) which can be expressed in the generic form

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D^{L^{2}}}\int_{D^{L^{2}}}\prod_{i=0}^{L-1}\prod_{j=0}^{L-1}f_{i,j}\left(x^{a}_{i,j}-x^{a}_{i,j+1}-x^{b}_{i,j}+x^{b}_{i+1,j}\right)\,\mathrm{d}\boldsymbol{x}^{a}\,\mathrm{d}\boldsymbol{x}^{b}.$

After a change of variables, we have $D=[0,1]$ and

$f_{i,j}(x)=f(x)=\exp(\beta\cos(2\pi x))\qquad\mbox{for all
}i,j=0,\ldots,L-1.$

To include the observable, in the numerator integral we will replace $f_{0,0}$
by

$f_{0,0}(x)=\cos(2\pi x)\,\exp(\beta\cos(2\pi x)).$

All these functions are clearly periodic. Below we will outline three stages
to simplify the problem.

#### STAGE ONE

By rearranging the variables in the $a$-direction and the $b$-direction, we
can write

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D^{L^{2}}}\prod_{i=0}^{L-1}\bigg{(}\underbrace{\int_{D^{L}}\prod_{j=0}^{L-1}f_{i,j}\left(x^{a}_{i,j}-x^{a}_{i,j+1}-x^{b}_{i,j}+x^{b}_{i+1,j}\right)\,\mathrm{d}\boldsymbol{x}^{a}_{i}}_{\mbox{$=:\,g_{i}\left(\boldsymbol{x}^{b}_{i+1}-\boldsymbol{x}^{b}_{i}\right)$}}\bigg{)}\,\mathrm{d}\boldsymbol{x}^{b},$

where we used the fact that each factor over the index $i$ depends only on
$\boldsymbol{x}^{a}_{i}$, and wrote the integral over $\boldsymbol{x}^{a}\in
D^{L^{2}}$ as a product of integrals over
$\boldsymbol{x}^{a}_{i}=(x^{a}_{i,0},\ldots,x^{a}_{i,L-1})\in D^{L}$. This
leads to

$\displaystyle{\mathcal{I}}$
$\displaystyle\,=\,\int_{D^{L}}\cdots\int_{D^{L}}\prod_{i=0}^{L-1}g_{i}\left(\boldsymbol{y}_{i+1}-\boldsymbol{y}_{i}\right)\,\mathrm{d}\boldsymbol{y}_{0}\cdots\mathrm{d}\boldsymbol{y}_{L-1},$
(11)

where

$\displaystyle g_{i}(\boldsymbol{y})$
$\displaystyle\,:=\,\int_{D^{L}}\prod_{j=0}^{L-1}f_{i,j}\left(x_{j}-x_{j+1}+y_{j}\right)\,\mathrm{d}\boldsymbol{x}.$
(12)

Thus we have obtained a nested integration problem where for the outer
integral (11) we have an integrand with first order couplings of the form (10)
(with $s$ replaced by $L$), and for the inner integral (12) we have first
order couplings given by (3) at each input $i$ and $\boldsymbol{y}$. Then, the
results in 1 apply for both the inner and outer integrals. If we use an
$n$-point rectangle rule for the inner integral and an $N$-point lattice
cubature rule for the outer integral, and if the functions $f_{i,j}$ are
periodic (thus so are the functions $g_{i}$), then we will fall under
Scenarios (A5)–(A7) for the inner and outer integrals. Furthermore, if all
functions $f_{i,j}$ are the same, then the final cost is of order

$N\,\log(N)+N\,n\,\log(n)\,.$

The cost is independent of $L$, while the error is of order
$N^{-\alpha}+n^{-\alpha}$, with $\alpha$ given by the smoothness of the
functions and the underlying lattice rule, with an implied error constant that
may depend exponentially on $L$.

#### STAGE TWO

The next two lemmas show that the problem can be simplified further for this
particular model with periodic functions $f_{i,j}$. The proofs can be found in
[17].

######  Lemma 1 ([17])

Assume the functions $f_{i,j}$ are periodic. Then the inner integral (12)
simplifies to

$\displaystyle
g_{i}(\boldsymbol{y})\,=\,g_{i}\Big{(}\textstyle\sum_{j=0}^{L-1}y_{j},0,\ldots,0\Big{)}\,=\,g_{i}\Big{(}\textstyle\sum_{j=0}^{L-1}y_{j},\boldsymbol{0}\Big{)},$
(13)

that is, $g_{i}(\boldsymbol{y})$ depends only on the sum of the components of
$\boldsymbol{y}$.

######  Lemma 2 ([17])

Assume the functions $f_{i,j}$ are periodic. Then the outer integral (11)
simplifies to

$\displaystyle{\mathcal{I}}\,=\,\int_{D^{L}}\prod_{i=0}^{L-1}g_{i}\big{(}y_{i+1}-y_{i},\boldsymbol{0}\big{)}\,\mathrm{d}\boldsymbol{y}.$
(14)

The last transformation in Lemma 2 leaves the outer integral (14) in the form
of (3) and so there is no longer a need to use a lattice cubature rule. By
taking now an $n$-point rectangle rule for both the inner and outer integrals,
Scenario (A7) applies in both cases and the cost becomes

$n^{2}\,\log(n)\,.$

The cost is again independent of $L$. The approximation error is then
${\mathcal{O}}(n^{-\alpha})$, where $\alpha$ depends on the smoothness of the
functions, but importantly, the implied constant no longer depends on $L$.

#### STAGE THREE

An alternative approach based on Fourier series was used in [17] to simplify
the expression even further, as shown in the following theorem. The proof can
be found in [17].

######  Theorem 1 ([17])

Suppose that the functions $f_{i,j}$ are periodic and have absolutely
convergent Fourier series. Define $\mu_{i+jL}:=f_{i,j}$ for
$i,j=0,\ldots,L-1$. Then the integral (11), with inner integral (12),
simplifies to

$\displaystyle{\mathcal{I}}\,=\,\int_{D^{L^{2}}}\prod_{k=0}^{L^{2}-1}\mu_{k}\big{(}x_{k+1}-x_{k}\big{)}\,\mathrm{d}\boldsymbol{x},$
(15)

where now the parametric periodicity is to be taken modulo $L^{2}$, i.e.,
$x_{k}\equiv x_{k\bmod L^{2}}$.

Theorem 1 shows that, instead of nested integrals, we now have a single
integral (15) with dimensionality $L^{2}$ of the form (3), with $L$ replaced
by $L^{2}$. We are again in Scenario (A7) with an $n$-point rectangle rule and
the cost is only of order

$n\,\log(n)\,,$

and the error is ${\mathcal{O}}(n^{-\alpha})$. Both the cost and the error
bound are independent of $L$.

##  5 Summary

This proceedings article reviews our recent work published in [17]. We
presented very high-dimensional integrals exhibiting a special strucure, that
arise from models such as the quantum rotor and the 2D compact $U(1)$ LGT. The
numerical methods presented here and in [17] extend the ones presented in [3,
18]. In particular, here and in [17] we investigated in detail the advantages
of applying FFT and gave sufficient conditions in order to meet the
requirements of FFT application. The latter conditions may apply to other LQCD
models than the ones presented in this work. Executable Julia codes were
provided in [17] together with some numerical results.

The use of lattice cubature rules can be beneficial and even essential when we
have an $L$-fold product of $s$-dimensional integrals, in order to meet the
sufficient conditions to apply FFT. These improvements are not limited to
systems with first order couplings. As shown in [17] they can be extended to
systems with higher order couplings. In the latter case, the problems can
still be truly high dimensional.

The application examples of the quantum rotor and 2D compact $U(1)$ LGT are
encouraging us to look more deeply and to further investigate potential
advantages on truly difficult problems in 3D and 4D compact $\mathit{SU}(N)$
LGT. We are currently investigating the latter problems using recursive
numerical integration techniques.

## Acknowledgments

We gratefully acknowledge financial support from the Australian Research
Council (ARC) under grant DP210100831 and the Research Foundation Flanders
(FWO) under grant G091920N.

## References

  * [1] S. Akiyama, Y. Kuramashi, and Y. Yoshimura, Quantum field theories with tensor renormalization group, arXiv:2111.04240 [hep-lat].
  * [2] D. Albandea, P. Hernández, A. Ramos, and F. Romero-López, Topological sampling through windings, Eur. Phys. J. C 81 (2021), 873.
  * [3] A. Ammon, A. Genz, T. Hartung, K. Jansen, H. Leövey, and J. Volmer, On the efficient numerical solution of lattice systems with low-order couplings, Comput. Phys. Commun. 198 (2016), 71–81.
  * [4] A. Ammon, T. Hartung, K. Jansen, H. Leövey, and J. Volmer, Overcoming the sign problem in one-dimensional QCD by new integration rules with polynomial exactness, Phys. Rev. D 94 (2016), 114508.
  * [5] S. Aoki et al. [Flavour Lattice Averaging Group], FLAG Review 2019: Flavour Lattice Averaging Group (FLAG), Eur. Phys. J. C 80 (2020), 113.
  * [6] M. C. Bañuls, K. Cichy, Review on novel methods for lattice gauge theories, Rept. Prog. Phys. 83 (2020), 024401.
  * [7] W. Bietenholz, R. Brower, S. Chandrasekharan, and U. J. Wiese, Perfect lattice topology: The Quantum rotor as a test case, Phys. Lett. B 407 (1997), 283.
  * [8] W. Bietenholz, U. Gerber, M. Pepe, and U.-J. Wiese, Topological Lattice Actions, JHEP 1012 (2010), 020.
  * [9] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, and J. Schlenk, A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec, Comput. Phys. Commun. 240 (2019) 120–137.
  * [10] P. Craig, A new reconstruction of multivariate normal orthant probabilities, J. R. Statist. Soc. B 70 (2008), 227–243.
  * [11] E. de Doncker, A. Almulihi, and F. Yuasa, High-speed evaluation of loop integrals using lattice rules, J. Phys. Conf. Ser. 1085 (2018), no. 5, 052005\.
  * [12] J. Dick, F. Y. Kuo, and I. H. Sloan, High-dimensional integration: the Quasi-Monte Carlo way, Acta Numer. 22 (2013), 133–288.
  * [13] J. Dick and F. Pillichshammer, Digital Nets and Sequences, Cambridge University Press, Cambridge, 2010.
  * [14] C. Gattringer and C. B. Lang, Quantum chromodynamics on the lattice, Lect. Notes Phys. 788 (2010), 1.
  * [15] A. Genz and D. K. Kahaner, The numerical evaluation of certain multivariate normal integrals, J. Comput. Appl. Math. 16 (1986), 255–258.
  * [16] T. Hartung and K. Jansen, Zeta-regularized vacuum expectation values, J. Math. Phys. 60 (2019), 093504.
  * [17] T. Hartung, K. Jansen, F. Y. Kuo, H. Leövey, D. Nuyens, and I. H. Sloan, Lattice meets lattice: Application of lattice cubature to models in lattice gauge theory, J. Comput. Phys. 443 (2021), 110527.
  * [18] T. Hartung, K. Jansen, H. Leövey, and J. Volmer, Avoiding the sign-problem in lattice field theory, in Monte Carlo and Quasi-Monte Carlo Methods 2018 (B. Tuffin and P. L’Ecuyer, eds), Springer Proceedings in Mathematics & Statistics, vol 324, 2020, pp. 231–249.
  * [19] A. J. Hayter, Recursive integration methodologies with statistical applications, J. Statist. Plann. Inference 136 (2006), 2284–2296.
  * [20] A. J. Hayter, Recursive integration methodologies with applications to the evaluation of multivariate normal probabilities, J. Stat. Theory Pract. 5 (2011), 563–589.
  * [21] F. J. Hickernell, Lattice rules: How well do they measure up?, in: Random and Quasi-Random Point Sets (P. Hellekalek and G. Larcher, eds.), Springer, Berlin, 1998, pp. 109–166.
  * [22] K. Jansen and T. Hartung, Zeta-regularized vacuum expectation values from quantum computing simulations, PoS LATTICE 2019 (2020), 363.
  * [23] K. Jansen, H. Leövey, A. Ammon, A. Griewank, and M. Müller-Preussker, Quasi-Monte Carlo methods for lattice systems: a first look, Comput. Phys. Commun. 185 (2014), 948–959.
  * [24] C. Lemieux, Monte Carlo and Quasi-Monte Carlo Sampling, Springer, New York, 2009.
  * [25] M. Lüscher, Computational Strategies in Lattice QCD, arXiv:1002.4232 [hep-lat].
  * [26] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, 1992.
  * [27] D. Nuyens, The construction of good lattice rules and polynomial lattice rules, in: Uniform Distribution and Quasi-Monte Carlo Methods (P. Kritzer, H. Niederreiter, F. Pillichshammer, A. Winterhof, eds.), Radon Series on Computational and Applied Mathematics Vol. 15, De Gruyter, 2014, pp. 223–256.
  * [28] S. Schaefer, R. Sommer, and F. Virotta [ALPHA], Critical slowing down and error analysis in lattice QCD simulations, Nucl. Phys. B 845 (2011), 93–119.
  * [29] I. H. Sloan and S. Joe, Lattice Methods for Multiple Integration, Oxford University Press, Oxford, 1994.
  * [30] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94 (2005), 170201.
  * [31] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10 (1974), 2445–2459.


# Nonlinear waves in a hot, viscous and nonextensive quark gluon plasma

Golam Sarwar<EMAIL_ADDRESS>School of Physical Sciences, National
Institute of Science Education and Research, HBNI, Jatni-752050, Odisha, India
Department of Physics,University of Calcutta, 92, A.P.C. Road, Kolkata-700009,
India  Md Hasanujjaman<EMAIL_ADDRESS>Department of Physics, Darjeeling
Government College, Darjeeling- 734101, India  Trambak Bhattacharyya
<EMAIL_ADDRESS>Bogoliubov Laboratory of Theoretical Physics,
Joint Institute for Nuclear Research, Dubna 141980, Moscow Region, Russian
Federation  Mahfuzur Rahaman<EMAIL_ADDRESS>Variable Energy
Cyclotron Centre, 1/AF Bidhan Nagar, Kolkata- 700064, India  Homi Bhabha
National Institute, Training School Complex, Mumbai - 400085, India  Abhijit
Bhattacharyya<EMAIL_ADDRESS>Department of Physics,University of
Calcutta, 92, A.P.C. Road, Kolkata-700009, India  Jan-e Alam<EMAIL_ADDRESS>Variable Energy Cyclotron Centre, 1/AF Bidhan Nagar, Kolkata- 700064, India
Homi Bhabha National Institute, Training School Complex, Mumbai - 400085,
India

###### Abstract

The effects of the nonextensive statistics on the nonlinear propagation of
perturbations have been studied within the scope of relativistic second order
dissipative hydrodynamics and the non-extensive equation of state. We have
shown that the equations which describe the propagation of nonlinear waves
under such situation are KdV-like (Korteweg–De Vries). Apart from their
preserved solitonic behaviour the dissipative nature of these waves are also
observed. The waves with larger amplitude and width dissipate less and and
propagate faster. The waves deplete more with smaller values of Tsallis
parameter, q and lower temperature of the medium. For vanishingly small
transport coefficients the nonlinear waves show breaking nature. These
findings suggest that the nature of the propagation of the nonlinear waves may
serve as a good probe to discern the extensive and non-extensive thermodynamic
nature of a fluid such as the quark gluon plasma produced in relativistic
nuclear collisions.

###### pacs:

12.38.Mh, 12.39.-x, 11.30.Rd, 11.30.Er

##  I Introduction

The propagation of the perturbations in physical quantities of fluids have
been widely used to unveil the thermodynamic state of the fluid. Depending on
the magnitude of the perturbation it has been categorized as linear Vladimir ;
Lick  or nonlinear respectively if the perturbation is small or comparable to
its average value. In the context of the quark gluon plasma (QGP) formed in
relativistic heavy ion collisions (RHICs), the response of the QGP to the
perturbations may help in characterizing the medium. The QCD fluid formed in
these violent collisions in different collider experiments at Large Hadron
Collider (LHC) and Relativistic Heavy Ion Collider (RHIC) transform to
hadronic matter via a cross over transition or an intermediary coexisting
phase depending on the value of temperature ($T$) and baryonic chemical
potential ($\mu_{B}$). Finally, these hadrons are detected experimentally.
Apart from these hadrons penetrating probes such as, photons and lepton pairs
emanated from each space time point of the medium are also detected to shed
light on the early hot and dense phase of QGP. The locally equilibrated system
of QGP responds to its high internal pressure by rapid expansion against the
vacuum.

The pressure gradient in QGP can be anisotropic and the degree of anisotropy
depends on the value of the impact parameter of the collisions. An anisotropic
interacting system of quarks and gluons created after the collisions will
develop a pressure gradient of different magnitude along different directions.
This anisotropic pressure gradient will result in different expansion rate
which in turn will induce anisotropy in momentum distribution of the hadrons
originating from the system on hadronization. Therefore, the momentum
anisotropy produced due to the initial spatial anisotropy could be interpreted
as a linear for small initial spatial anisotropy. However, nonlinear effects
will be crucial if the anisotropy in initial geometry is large due to
collisions at large impact parameter. Moreover, quarks and gluons produced in
the early stage of collision with high momenta do not equilibrate but
propagate as jets through the thermal medium formed due the rescattering of
the low momentum quarks and gluons. These jets while propagating deposit
energy in the medium Shuryaknlw ; Betz . If the magnitude of the energy
density resulting from the jet-medium interaction is comparable to the energy
density of the medium then the effects of jet propagation will be treated as
nonlinear. On the other hand smaller value of the deposited energy density
allows the study of the jet-medium interaction within the domain of linear
response. The propagation of linear perturbations in QGP has been studied
extensively by several authors Shuryak1 ; Shuryak2 ; Sarwar:2015mma ;
Rafiei:2016zxk ; hasan1 ; hasan2 ; kunihiro ; fluidity . However, the
nonlinear aspects of the perturbations have been studied only by a few authors
Raha1 ; Raha2 ; Raha3 ; Raha4 ; Fogaca1 ; Fogaca2 ; Sarwar1 . In Refs. Raha1 ;
Raha2  the background medium is considered as an ideal fluid and found a
solitonic behavior of the perturbations. Fogaca et al in Refs Fogaca1 ;
Fogaca2  went further and introduced the shear viscous effects in the medium
and considered the perturbations up to second order. They too found the
solitonic behaviour of the perturbations. In these previous studies, the local
equilibrium is treated within the scope of Gibbs-Boltzmann (GB) extensive
distribution. However, particle spectra from small systems formed in $p-p$
collision are found to be well described by non-extensive equilibrium
distribution suggested by Tsallis Tsallis ; Tsallisbook ; TsCMS ; TsALICE ;
Biro09 ; Cleymans13 ; mcdprd ; TsallisTaylor ; tsallisraa ; Grigoryan17 ;
bcmmp ; abcp . In fact the description of data is better Wong ; Azami ; TB ;
Cleymans ; Zheng  with Tsallis distribution than the GB distribution.

The appearance of the entropic index, $q\neq 1$ (also called nonextensive
parameter) is attributed to the local fluctuations in the system Wilk00 . It
is important to recall at this point that in the limit of $q\rightarrow 1$ the
Tsallis distribution approached the GB distribution. In principle, $q$ can be
determined from the dynamics of the system. For some cases like non-linear
Fokker-Planck equation tsallispre , Boltzmann lattice model lattpre , cold
atoms in optical lattices coldprl  etc the $q$ is known analytically in terms
of microscopic or mesoscopic quantities. A relation between the entropic index
$q$ and the basic parameters of the QCD has been established in in DMM1 ; DMM2
and show that this value of $q$ can be used to reproduce the transverse
momentum spectra of hadrons produced in relativistic proton+proton collisions
DMM2020 . The value of the parameter $q$ of a many body system has been
expressed in terms of the number ($N$) of particles as $q=1+2N^{-1}/3$ and
shown that the two particle correlation becomes weaker as
$N\rightarrow\infty$, i.e. when the system approach the domain of molecular
chaos LD  which exposes the importance of non-extensitivity of a system with
finite number of particles e.g.. This has relevance for system formed in high
energy collisions of protons and nuclei. It has been shown in Wilk00  that in
a system with fluctuating temperature zones, $\beta$ is the inverse of average
temperature and $q$ is related to the relative variance in temperature. The
interplay between the parameter $q$ and the QCD dynamics has been exposed in
Ref.s RY1 ; RY2 . Although in relativistic nuclear collisions the $q$ is
treated as a parameter to fit the transverse momentum spectra of hadrons which
can in principle vary from 0 to $\infty$. However, for the phenomenological
studies the values of $q$ are generally greater than one and from the
thermodynamic considerations it has an upper bound of 4/3 BhattaCleMog . The
descriptions of different experimental observables that carries the effect of
fluctuations require $q\neq 1$ as argued in Refs. Navarra1 ; Navarra2 ; Wilk ;
Wilk00 ; Wilk01 ; Wilk09 ; Biro05 . The value of the $q$ varies with collision
centrality, energy, and system size Wilk . It is reported that the system size
alone is not sufficient to explain the values of nonextensive parameter
Sarwar2 .

Moreover, there are other observables advocate for choosing the non-extensive
local thermal equilibrium. The system size dependence of quantum fluctuations
have been studied theoretically in Refs. Arpan1 ; Arpan2 ; Arpan3  and these
fluctuations may be represented by the $q$ parameter. Whatever be the
underlying reason, the appearance of the $q-$parameter in particle spectra for
small systems draws serious consideration to the fact that a similar situation
may arise in some local regions(of similar size) in bigger systems. Because
the spatial inhomogeneity created in a larger system (produced in nuclear
collisions) may be of the size of the entire system produced in a $p-p$
collision. This appeals that the consideration of local equilibrium of such
systems reserves more for non-extensivity. At least, to strike out this
possibility, there should be a way to distinguish whether non-extensive local
equilibrium is the underlying nature of a fluid medium or an extensive local
equilibrium is the correct situation for the fluid medium. This requires
investigations of the effect of non-extensivity on the response of the fluid.

To understand the distinct effect of the non-extensivity from extensivity in
the local equilibrium of the fluid, one considers two aspects to respect the
non extensivity: a) the equation of state to be of a non-extensive form, and
b) the fluid dynamic equations should also be of a non-extensive type, i.e.,
to say that the fluid velocity must also be defined by taking into account the
non-extensivity, which will clearly be different from the extensive one
Osada:2008sw .

Recently the study of the propagation of non-linear waves has considered only
the first modification from the extensive one, where the non-linear wave
equation for first-order perturbation is derived for the ideal fluid
background of fixed thermodynamic quantities TrambakNLW . The authors found
the breaking wave solution with non-dissipative propagation of the nonlinear
waves. However, it has been shown in Ref. Osada:2008sw  that the ideal non-
extensive hydrodynamics ($q$\- hydrodynamics) is equivalent to extensive
dissipative hydrodynamics ($d$\- hydrodynamics). The corresponding dissipation
should result from the non-extensive nature of fluid velocity field
Osada:2008sw . Therefore, even in ideal non-extensive background, one expects
dissipative propagations of perturbations. Moreover, the local equilibrium
leads to inherent gradients of fluid dynamic fields which leads to transport
phenomenon. Therefore, it is important to investigate the dissipative
propagation of the nonlinear waves in the non-extensive background to
understand the role of non-extensivity on the propagation of perturbations.

In this work we proceed to investigate the propagation of nonlinear waves in a
QGP with non-extensive equation of state (EoS) and viscosities appearing from
extensive fluid. The aim is to find how the conclusion changes with the
background (either extensive or non-extensive), and how $q$ affects the
propagation of the nonlinear waves. The required equations for propagation(of
nonlinear waves) are derived from Israel- Stewart (IS) hydrodynamics and the
order of the perturbations are taken up to second order.

The paper is organized as follows: In the next three sections, i.e  in
Sec.III, Sec.IV, and Sec.V, we derive the required equations for nonlinear
waves from IS theory with all the relevant transport coefficients. Then in
Sec.V we discuss the EoS for the non-extensive fluid. Finally we discuss the
results in Sec.VI and summarize the findings in Sec.VII. In this work we have
used natural units ($\hbar=c=1,k_{B}=1$) and the Minkowski metric to be
$g^{\mu\nu}=(1,-1,-1,-1)$.

##  II Hydrodynamic equations

The non-extensive version of ideal and first order dissipative hydrodynamics
have been investigated in Ref. Osada:2008sw  and Ref Biro:2011bq
respectively. It is found in Ref. Biro:2011bq  that in the $q-$ hydrodynamics,
the main tensorial structure of energy momentum tensor remains same as the
$d$\- hydrodynamics whereas changes appears only in thermodynamic quantities
and EoS. Here we use this observations and derive the nonlinear wave equation
from same tensorial structure of second order IS theory, with thermodynamic
quantities and equation of state from non-extensive thermodynamics. The fluid
velocity field considered here is of non-extensive nature.

Relativistic viscous hydrodynamics is a useful tool to describe the space-time
evolution of the fluid created in relativistic heavy ion collisions(RHICs).
Depending on the magnitude of the expansion gradient, the order of the theory
is being anticipated. Navier-Stokes (NS) theory deals with first-order
gradients of hydrodynamic fields, usually designated as the first-order
viscous hydrodynamics Eckart ; LL . The solution of the NS equation is acausal
and unstable. However, there are recently developed theory to show that the
first-order hydrodynamics can be stable as well as causal Bemfica1 ; Bemfica2
; Bemfica3 ; Kovtun ; Das ; Arpan .

Whereas, second-order hydrodynamics takes into account the second order
gradient in hydrodynamic fields and the relaxation effects were introduced to
attain causality and stability. Several prescriptions of the second-order
hydrodynamics can be found in Refs. Israel ; Hiscock ; Dennicol ; Romatschke ;
Van ; Muronga ; Tsumura . In this work, the second-order Israel-Stewart-Muller
hydrodynamics is used to obtain the required evolution equations for nonlinear
waves. Out of different choices of frames of reference for defining fluid
velocity, we use the the Eckart frame of reference. Eckart frame represents a
local rest frame where the net charge dissipation is vanishing but the energy
dissipation is non-vanishing Eckart . In Eckart frame, the energy-momentum
tensor(EMT) and particle current can be written as:

$\displaystyle T^{\mu\nu}$ $\displaystyle=$ $\displaystyle\epsilon
u^{\mu}u^{\nu}-p\Delta^{\mu\nu}+\Delta T^{\mu\nu}$ (1) $\displaystyle J^{\mu}$
$\displaystyle=$ $\displaystyle nu^{\mu}$ (2)

where $\epsilon$ is the local energy density field, $u^{\mu}$ is the fluid
four velocity, $p$ is the local pressure, and $\Delta T^{\mu\nu}$ is the
correction to the EMT due to dissipation. The $u^{\mu}$ satisfies the
condition $u^{\mu}u_{\mu}=1$ implying $u^{\mu}\partial_{\nu}u_{\mu}=0$. The
projection operator normal to $u^{\mu}$ is defined as
$\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ such that
$\Delta^{\mu\nu}u_{\nu}=0$. This gives $\Delta^{\mu\nu}\Delta_{\mu\nu}=3$. And
the symmetric traceless projection normal to $u^{\mu}$ is defined as
$\Delta^{\mu\nu}_{\alpha\beta}=\frac{1}{2}(\Delta^{\nu}_{\alpha}\Delta^{\mu}_{\beta}+\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}-\frac{2}{3}\Delta^{\mu\nu}\Delta_{\alpha\beta})$.

The dissipative part in $T^{\mu\nu}$ can be decomposed in terms of scalar,
vector and tensor as:

$\displaystyle\Delta
T^{\mu\nu}=-\Pi\Delta^{\mu\nu}+q^{\mu}u^{\nu}+q^{\nu}u^{\mu}+\pi^{\mu\nu}$ (3)

The complete energy momentum tensor in Eckart frame reads as:

$\displaystyle T^{\mu\nu}=\epsilon
u^{\mu}u^{\nu}-(p+\Pi)\Delta^{\mu\nu}+q^{\mu}u^{\nu}+q^{\nu}u^{\mu}+\pi^{\mu\nu}$
(4)

In the fluid rest frame the vector and tensor form of dissipations are
considered to be non-existent that is $u_{\mu}q^{\mu}=0$,
$u_{\mu}\pi^{\mu\nu}=0$. The scalar dissipation ($\Pi$) corresponds to the
non-equilibrium pressure perturbation which can not be made to relate to
energy density using the EoS. It is related to volume expansion of the fluid,
triggers small non-equilibrium perturbation, which tends bring a new
equilibrium state through this dissipative correction. The vector
dissipations, $q^{\mu}$, $\pi^{\mu\nu}$ corresponds to dissipative energy flow
and the shear stress tensor respectively. The relations between energy
density, pressure and the dissipative fluxes and EMT are given by the
following relations:

$\displaystyle u_{\mu}T^{\mu\nu}u_{\nu}$ $\displaystyle=$
$\displaystyle\epsilon,\,\,\,\,q_{\alpha}=u_{\mu}T^{\mu\nu}\Delta_{\nu\alpha}=u_{\mu}\Delta
T^{\mu\nu}\Delta_{\nu\alpha},$ $\displaystyle p+\Pi$ $\displaystyle=$
$\displaystyle-\frac{1}{3}\Delta_{\mu\nu}T^{\mu\nu},\,\,\,\,\Pi=-\frac{1}{3}\Delta_{\mu\nu}\Delta
T^{\mu\nu},\,\,\,\,u_{\mu}\Delta T^{\mu\nu}=q^{\nu}$ (5)

The conservation laws for EMT and the net charge (here the net baryon number)
density are given by

$\displaystyle\partial_{\mu}T^{\mu\nu}$ $\displaystyle=$ $\displaystyle 0$ (6)
$\displaystyle\partial_{\mu}J^{\mu}$ $\displaystyle=$ $\displaystyle 0$ (7)

The dissipative fluxes Israel  are given by

$\displaystyle\Pi$ $\displaystyle=$
$\displaystyle-\frac{1}{3}\zeta\Big{[}\partial_{\mu}u^{\mu}+\beta_{0}D\Pi-\tilde{\alpha_{0}}\partial_{\mu}q^{\mu}\Big{]}$
$\displaystyle\pi^{\lambda\mu}$ $\displaystyle=$
$\displaystyle-2\eta\Delta^{\lambda\mu\alpha\beta}\Big{[}\partial_{\alpha}u_{\beta}+\beta_{2}D\pi_{\alpha\beta}-\tilde{\alpha_{1}}\partial_{\alpha}q_{\beta}\Big{]}$
$\displaystyle q^{\lambda}$ $\displaystyle=$ $\displaystyle-\kappa
T\Delta^{\lambda\mu}\Big{[}\frac{1}{T}\partial_{\mu}T+Du_{\mu}+\tilde{\beta_{1}}D{q_{\mu}}-\tilde{\alpha_{0}}\partial_{\mu}\Pi-\tilde{\alpha_{1}}\partial_{\nu}\pi^{\nu}_{\mu}\Big{]}$
(8)

Here, $\eta$, $\zeta$, $\kappa$ are the coefficient of shear viscosity, bulk
viscosity and thermal conductivity respectively. $D\equiv
u^{\mu}\partial_{\mu}$, is the co-moving derivative and in the local rest
frame (LRF) it is simply $D\Pi=\dot{\Pi}$, represents the time derivative.
Here $\beta_{0},\tilde{\beta_{1}},\beta_{2}$ are relaxation coefficients,
$\tilde{\alpha_{0}}$ and $\tilde{\alpha_{1}}$ are coupling coefficients.
$\beta_{0},\tilde{\beta_{1}},\beta_{2}$ are related to the bulk pressure
($\tau_{\Pi}$), the heat flux ($\tau_{q}$) and the shear tensor
($\tau_{\pi}$), defined as Muronga

$\tau_{\Pi}=\zeta\beta_{0},\,\,\,\,\tau_{q}=k_{B}T\beta_{1},\,\,\,\,\tau_{\pi}=2\eta\beta_{2}$
(9)

where $\tilde{\alpha_{0}}$ and $\tilde{\alpha_{1}}$ are related to the
relaxation lengths which couple the heat flux and bulk pressure ($l_{\Pi
q},l_{q\Pi}$). The heat flux and shear tensor $(l_{q\pi},l_{\pi q})$ are
defined as:

$l_{\Pi
q}=\zeta\alpha_{0},\,\,\,\,l_{q\Pi}=k_{B}T\alpha_{0},\,\,\,\,l_{q\pi}=k_{B}T\alpha_{1},\,\,\,\,l_{\pi
q}=2\eta\alpha_{1}$ (10)

The relaxation and coupling coefficients are given by Israel

$\displaystyle\alpha_{0}$ $\displaystyle=$
$\displaystyle(D_{41}D_{20}-D_{31}D_{30})\Lambda\phi\Omega J_{21}J_{31}$
$\displaystyle\alpha_{1}$ $\displaystyle=$
$\displaystyle(J_{41}J_{42}-J_{31}J_{52})\Lambda\phi J_{21}J_{31}$
$\displaystyle\beta_{0}$ $\displaystyle=$
$\displaystyle\frac{3\beta}{\phi^{2}\Omega^{2}}[5J_{52}-\frac{3}{D_{20}}\\{J_{31}(J_{31}J_{30}-J_{41}J_{20})+J_{41}(J_{41}J_{10}-J_{31}J_{20})\\}]$
$\displaystyle\beta_{1}$ $\displaystyle=$
$\displaystyle\frac{D_{41}}{\Lambda^{2}nmJ_{21}J_{31}}$
$\displaystyle\beta_{2}$ $\displaystyle=$ $\displaystyle\frac{\beta
J_{52}}{2\phi^{2}}$ (11)

where,

$\displaystyle D_{rs}$ $\displaystyle=$ $\displaystyle
J_{r+1,s}J_{r-1,s}-(J_{rs})^{2},\,\,\,\,\phi=T(\epsilon+P),\,\,\,\,\psi=\frac{\epsilon+P}{n_{i}m_{i}},$
$\displaystyle\Lambda$ $\displaystyle=$ $\displaystyle
1+5(\frac{\phi}{nm})-\psi^{2},\,\,\,\,\Omega=3\Big{(}\frac{\partial
ln\phi}{\partial lnn}\Big{)}-5$ (12)

where $r$, $s$ are integers, $m_{i}$ and $n_{i}$ are respectively the mass and
density of $i$-th type particle here. $\phi,\psi$ and $\Lambda$ are calculated
by using their relations with $\epsilon$, $P$ and $n_{i}$. We take the current
quark mass of up and down flavours as 10 MeV. $J_{rs}$ is defined as,

$\displaystyle J_{rs}=\frac{A_{0}}{(2s+1)!!}\int_{0}^{\infty}N\Delta
Sinh^{2(s+1)}\mathscr{R}Cosh^{r-2s}\mathscr{R}d\mathscr{R}$ (13)

and

$\displaystyle N=\frac{1}{exp_{q}(\beta Cosh\mathscr{R}-\alpha)-\epsilon}$
(14)

where $q$-deformed exponential is defined as

$\displaystyle exp_{q}(x)=\big{[}1+(1-q)x\big{]}^{1/(1-q)}$ (15)

and

$\displaystyle\Delta=1+\epsilon N,\,\,\,\,A_{0}=4\pi m_{i}^{3}$ (16)

##  III 1-D flow equations

In this section we derive the equations to study the propagation of nonlinear
perturbation in a fluid in $(1+1)$ dimension. The evolution equations of
energy density and fluid velocity are obtained from the Eqs.(6) and (7) by
taking projection along the directions parallel and perpendicular to $u^{\mu}$
as:

$\displaystyle D\epsilon=-(\epsilon+P)\theta+(\partial_{\mu}u_{\nu})\Delta
T^{\nu\mu}-(\partial_{\mu}q^{\mu})$ (17)
$\displaystyle(\epsilon+p)Du^{\alpha}=\Delta^{\alpha\mu}\partial_{\mu}-\Delta^{\alpha}_{\nu}\partial_{\mu}\Delta
T^{\nu\mu}$ (18)

where $u^{\mu}=(\gamma,\gamma v,0,0),\,\,\,\,\gamma=1/\sqrt{1-v^{2}}$ and

$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\partial_{\mu}u^{\mu}=v\gamma^{3}\frac{\partial v}{\partial
t}+\gamma^{3}\frac{\partial v}{\partial x}$ (19)
$\displaystyle\partial_{\mu}q^{\mu}$ $\displaystyle=$ $\displaystyle
q^{x}\frac{\partial v}{\partial t}+v\frac{\partial q^{x}}{\partial
t}+\frac{\partial q^{x}}{\partial x}$ (20)

Here,

$\displaystyle
u_{\mu}q^{\nu}=0,\,\,\,u_{\mu}\pi^{\mu\nu}=0,\,\,\,g_{\mu\nu}\pi^{\mu\nu}=0,\,\,\,\pi^{xy}=\pi^{zx}=\pi^{yz}=0,\,\,\,\pi^{yy}=\pi^{zz}=\pi^{T}$
(21)

which leads to

$\displaystyle\pi^{00}=v^{2}\pi^{xx},\,\,\,\pi^{xt}=\pi^{tx}=v^{2}\pi^{xx},\,\,\,\pi^{T}=\pi^{tx}=\frac{(v^{2}-1)}{2}\pi^{xx}$
(22)

and

$\displaystyle q^{t}=vq^{x},\,\,\,q^{y}=q^{z}=q^{T}$ (23)

Now from the energy-momentum conservation equations, we get

$\displaystyle\gamma v\big{(}\frac{\partial\epsilon}{\partial
t}+\frac{\partial\epsilon}{\partial
x}\big{)}=\\{(h+\Pi)v\gamma^{3}+v\gamma\pi^{xx}+q^{x}\\}\frac{\partial
v}{\partial t}+\\{(h+\Pi)\gamma^{3}+\gamma\pi^{xx}\\}\frac{\partial
v}{\partial x}+v\frac{\partial q^{x}}{\partial t}+\frac{\partial
q^{x}}{\partial x}$ (24)
$\displaystyle\\{(h+\Pi)v\gamma^{4}+v\gamma^{2}(1-2v^{2})v\pi^{xx}+\gamma^{3}(1+v^{2})vq^{x}\\}\frac{\partial
v}{\partial t}+\\{(h+\Pi)v\gamma^{4}-\gamma^{2}v^{2}\pi^{xx}$
$\displaystyle+\gamma v(1+2\gamma^{2})q^{x}\\}\frac{\partial v}{\partial
x}=v^{2}\big{(}-\gamma^{2}\frac{\partial p_{\zeta}}{\partial
t}-\frac{\gamma^{2}}{v}\frac{\partial p_{\zeta}}{\partial
x}-\frac{\partial\pi^{xx}}{\partial
t}-\frac{1}{v}\frac{\partial\pi^{xx}}{\partial
x}-\frac{\gamma}{v}\frac{\partial q^{x}}{\partial t}-\gamma\frac{\partial
q^{x}}{\partial x}$ (25)

where $h=\epsilon+p$ is the enthalpy density. We define $p+\Pi=p_{\zeta}$.
From the equation governing the number conservation we get

$\displaystyle\frac{\partial n}{\partial t}+v\frac{\partial n}{\partial
x}=-n\gamma^{2}(v\frac{\partial v}{\partial t}+\frac{\partial v}{\partial x})$
(26)

The other three equation originating from dissipative fluxes are,

$\displaystyle\Pi+\zeta\beta_{0}(\gamma\frac{\partial\Pi}{\partial
t}+v\frac{\partial\Pi}{\partial
x})=-\zeta[\\{v\gamma^{3}-\alpha_{0}q^{x}\\}\frac{\partial v}{\partial
t}+\\{\gamma^{3}\\}\frac{\partial v}{\partial x}-\alpha_{0}v\frac{\partial
q^{x}}{\partial t}-\alpha_{0}\frac{\partial q^{x}}{\partial x}]$ (27)
$\displaystyle q^{x}+\kappa T\beta_{1}(\gamma\frac{\partial q^{x}}{\partial
t}+v\frac{\partial q^{x}}{\partial x})=\kappa
T[\\{-\gamma^{4}+(\beta_{1}\gamma^{3}v)q^{x}+\alpha_{1}(2-\gamma^{2})\pi^{xx}\\}\frac{\partial
v}{\partial t}+\\{-v\gamma^{4}+(\beta_{1}\gamma^{3}v^{2})q^{x}$
$\displaystyle-\alpha_{1}\gamma^{2}v\pi^{xx}\\}\frac{\partial v}{\partial
x}+\alpha_{1}v\frac{\partial\pi^{xx}}{\partial
t}+\alpha_{1}\frac{\partial\pi^{xx}}{\partial
x}-\alpha_{0}v\gamma^{2}\frac{\partial\Pi}{\partial
t}-\alpha_{0}\gamma^{2}\frac{\partial\Pi}{\partial
x}+\frac{1}{T}v\gamma^{2}\frac{\partial T}{\partial
t}+\frac{1}{T}\gamma^{2}\frac{\partial T}{\partial x}]$ (28)
$\displaystyle\pi^{xx}+\frac{4}{3}\eta\beta_{2}(\gamma\frac{\partial\pi^{xx}}{\partial
t}+v\frac{\partial\pi^{xx}}{\partial
x})=\frac{4}{3}\eta[\\{\gamma^{5}+(\alpha_{1}\gamma^{4}v)q^{x}+(2\beta_{2}\gamma^{3})\pi^{xx}\\}\frac{\partial
v}{\partial t}+\\{\gamma^{5}+(\alpha_{1}\gamma^{4}v)q^{x}$
$\displaystyle+(2\beta_{2}\gamma^{3}v^{2})\pi^{xx}\\}\frac{\partial
v}{\partial x}-\alpha_{1}\gamma^{2}v\frac{\partial q^{x}}{\partial
t}-\alpha_{1}\gamma^{2}\frac{\partial q^{x}}{\partial x}]$ (29)

##  IV Derivation of nonlinear equations

The above hydrodynamic equations from IS causal theory is used to derive the
nonlinear equations. We have adapted Reductive Perturbative Method (RPM) RPM1
; RPM2 ; RPM3 ; RPM4 . To proceed further we need to define ’stretched co-
ordinates’ as

$\displaystyle X=\frac{\sigma^{1/2}}{L}(x-c_{s}t),\text{and}\hskip
11.38092ptY=\frac{\sigma^{3/2}}{L}(c_{s}t)$ (30)

where $L$ is the characteristic length and $c_{s}$ is the speed of sound.
Therefore,

$\displaystyle\frac{\partial}{\partial
x}=\frac{\sigma^{1/2}}{L}\frac{\partial}{\partial X},\text{and}\hskip
11.38092pt\frac{\partial}{\partial
t}=-c_{s}\frac{\sigma^{1/2}}{L}\frac{\partial}{\partial
X}+c_{s}\frac{\sigma^{3/2}}{L}\frac{\partial}{\partial Y}$ (31)

where, $\sigma$ is the expansion parameter. The coordinate $X$ is measured
from the frame of propagating sound waves, whereas the $Y$ represents fast
moving coordinate. The RPM technique is devised to preserve the structural
form of the parent equation. In the equations for perturbations for different
order of $\sigma$ i.e breaking wave, Berger’s, KdV equation etc. Along with
this, we perform the simultaneous series expansion of hydrodynamic quantities,
in terms of the expansion parameter $\sigma$ to get,

$\displaystyle\hat{\epsilon}$ $\displaystyle=$
$\displaystyle\frac{\epsilon}{\epsilon_{0}}=1+\sigma\epsilon_{1}+\sigma^{2}\epsilon_{2}+\sigma^{3}\epsilon_{3}+...,\,\,\,\hat{p}=\frac{p}{p_{0}}=1+\sigma
p_{1}+\sigma^{2}p_{2}+\sigma^{3}p_{3}+...,$ $\displaystyle\hat{v}$
$\displaystyle=$ $\displaystyle\frac{v}{c_{s}}=\sigma
v_{1}+\sigma^{2}v_{2}+\sigma^{3}v_{3}+...,\,\,\,\hat{\Pi}=\frac{\Pi}{p_{0}}=\sigma\Pi_{1}+\sigma^{2}\Pi_{2}+\sigma^{3}\Pi_{3}+...,$
$\displaystyle\hat{q^{x}}$ $\displaystyle=$
$\displaystyle\frac{q^{x}}{\epsilon_{0}}=\sigma
q^{x}_{1}+\sigma^{2}q^{x}_{2}+\sigma^{3}q^{x}_{3}+...,\,\,\,\hat{\pi^{xx}}=\frac{\pi^{xx}}{p_{0}}=\sigma\pi^{xx}_{1}+\sigma^{2}\pi^{xx}_{2}+\sigma^{3}\pi^{xx}_{3}+...$
(32)

We rewrite the equations by collecting different order of $\sigma$ and keep up
to third order. Finally we revert from $(X,Y)\to(t,x)$ to get the following
equations for the perturbation of the energy density ($\hat{\epsilon}$) as:

$\displaystyle\frac{\partial\hat{\epsilon_{1}}}{\partial
t}+\left[1+(1-c_{s}^{2})\frac{\epsilon_{0}}{\epsilon_{0}+p_{0}}\hat{\epsilon_{1}}\right]c_{s}\frac{\partial\hat{\epsilon_{1}}}{\partial
x}-\left[\frac{1}{2(\epsilon_{0}+p_{0})}(\zeta+\frac{4}{3}\eta)\right]\frac{\partial^{2}\hat{\epsilon_{1}}}{\partial
x^{2}}=0,$ (33)

and

$\displaystyle\frac{\partial\hat{\epsilon_{2}}}{\partial
t}+\mathscr{S}_{1}\frac{\partial\hat{\epsilon_{2}}}{\partial
x}+\mathscr{S}_{2}\frac{\partial\hat{\epsilon_{1}}}{\partial
x}+\mathscr{S}_{3}\frac{\partial^{2}\hat{\epsilon_{1}}}{\partial
x^{2}}+\mathscr{S}_{4}\frac{\partial^{3}\hat{\epsilon_{1}}}{\partial
x^{3}}+\mathscr{S}_{5}\frac{\partial^{2}\hat{\epsilon_{2}}}{\partial x^{2}}=0$
(34)

where $\hat{\epsilon_{i}}=\sigma^{i}\epsilon_{i}$ for $i=1,2,..$. The
coefficients $\mathscr{S}_{i}$’s for $i=1$ to $5$ are given by,

$\displaystyle\mathscr{S}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{0}+p_{0}}\Big{[}c_{s}\\{\epsilon_{0}\left(1-\hat{\epsilon}_{1}\left(c_{s}^{2}-1\right)\right)+p_{0}\\}\Big{]};$
$\displaystyle\mathscr{S}_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon_{0}+p_{0}}\Big{[}\epsilon_{0}c_{s}\\{c_{s}^{2}-1\\}\\{\epsilon_{0}\left(\left(2c_{s}^{2}+1\right)\hat{\epsilon}_{1}{}^{2}-\hat{\epsilon}_{2}\right)-p_{0}\hat{\epsilon}_{2}\\}\Big{]};$
$\displaystyle\mathscr{S}_{3}$ $\displaystyle=$
$\displaystyle\frac{1}{12(\epsilon_{0}+p_{0})^{2}}
\Big{[}\epsilon_{0}\hat{\epsilon}_{1}\\{3c_{s}^{2}(7\zeta+8\eta)+3\zeta+4\eta\\}\Big{]};$
$\displaystyle\mathscr{S}_{4}$ $\displaystyle=$
$\displaystyle-\frac{1}{72c_{s}c_{V}(\epsilon_{0}+p_{0})^{2}}\Big{[}4c_{V}c_{s}^{2}\\{3\kappa
T\epsilon_{0}(3\alpha_{0}\zeta+4\alpha_{1}\eta)+3\kappa
T(3\zeta+4\eta)+3\kappa Tp_{0}(3\alpha_{0}\zeta+4\alpha_{1}\eta)$
$\displaystyle+(\epsilon_{0}+p_{0})(9\beta_{0}\zeta^{2}+16\beta_{2}\eta^{2})\\}-c_{V}(3\zeta+4\eta)^{2}+12\kappa\\{\epsilon_{0}+p_{0}\\}\\{\epsilon_{0}(3\alpha_{0}\zeta+4\alpha_{1}\eta)+3\zeta+4\eta$
$\displaystyle+p_{0}(3\alpha_{0}\zeta+4\alpha_{1}\eta)\\}\Big{]};$
$\displaystyle\mathscr{S}_{5}$ $\displaystyle=$
$\displaystyle-\frac{3\zeta+4\eta}{6\left(p_{0}+\epsilon_{0}\right)}$ (35)

##  V Tsallis MIT bag equation of state

We consider that quarks and gluons in the viscous plasma medium follow the
quantum Tsallis fermionic (‘$f$’) and bosonic (‘$b$’) single particle
distributions given by tsallisM ,

$\displaystyle
n_{fq}=\frac{1}{\left[1+(q-1)\frac{E_{p}-\mu}{T}\right]^{\frac{q}{q-1}}+1}$
$\displaystyle
n_{bq}=\frac{1}{\left[1+(q-1)\frac{E_{p}-\mu}{T}\right]^{\frac{q}{q-1}}-1},$
(36)

where $E_{p}=\sqrt{p^{2}+m^{2}}$ is the single particle energy of a particle
of mass $m$, $\mu$ is the chemical potential (assumed zero in this paper), $q$
is the Tsallis parameter, and $T$ is the Tsallis temperature. For a gas of
quarks and gluons, the thermodynamic variables like the energy density
$\epsilon$ and pressure $P$ in the MIT bag model Weisskopf  are given by,

$\displaystyle\epsilon_{\mathrm{bag}}$ $\displaystyle=$
$\displaystyle\mathscr{B}+\epsilon_{b}+2\epsilon_{f},$ (37) $\displaystyle
P_{\mathrm{bag}}$ $\displaystyle=$ $\displaystyle-\mathscr{B}+P_{b}+2P_{f},$
(38)

where $\mathscr{B}$ is the bag constant. The thermodynamic variables
$\epsilon$ and $P$ are given by,

$\displaystyle\epsilon_{i}=g\int\frac{d^{3}p}{(2\pi)^{3}}E_{p}n_{i},~{}P_{i}=g\int\frac{d^{3}p}{(2\pi)^{3}}\frac{p^{2}}{3E_{p}}n_{i},$
(39)

where $g$ is the degeneracy factor and $i=b,f$.

The factor 2 in front of the fermionic parts in Eqs. (37), (38) takes into
account both particles and anti-particles. The net baryonic number density
$\rho_{\mathrm{B}}$ is assumed to be zero. The fermionic thermodynamic
variables have been evaluated up to $\mathscr{O}(m^{2}T^{2})$, where $m$ is
the mass of the fermions. It is assumed that the quark-gluon plasma consists
of the up quark, down quark and gluons. It has been verified that for a wide
range of the $q$ and $T$ parameter values relevant for phenomenological
studies of high-energy collisions, the $\mathscr{O}(m^{2}T^{2})$ approximated
results overlap with the exact results when mass is small ($\sim$10 MeV). We
write down the expressions for the energy density and pressure in the non-
extensive bag model TsallisMIT1 ; TsallisMIT2 ,

$\displaystyle\epsilon_{\mathrm{bag}}(x,t)$ $\displaystyle=$
$\displaystyle\mathscr{B}+\epsilon_{b,1}T^{4}+2(\epsilon_{f,1}T^{4}+\epsilon_{f,2}m^{2}T^{2}),$
(40) $\displaystyle=$
$\displaystyle\mathscr{B}+\epsilon_{\mathrm{bag},1}T^{4}+\epsilon_{\mathrm{bag},2}m^{2}T^{2}$
$\displaystyle P_{\mathrm{bag}}(x,t)$ $\displaystyle=$
$\displaystyle-\mathscr{B}+P_{b,1}T^{4}+2(P_{f,1}T^{4}+P_{f,2}m^{2}T^{2})$
(41) $\displaystyle=$
$\displaystyle-\mathscr{B}+P_{\mathrm{bag},1}T^{4}+P_{\mathrm{bag},2}m^{2}T^{2}.$

When written out long-hand, $\epsilon_{i,\ell}$ and $P_{i,\ell}$ ($i=f,b$;
$\ell=1,2$) are given by TrambakNLW ,

$\displaystyle\epsilon_{b,1}$ $\displaystyle=$
$\displaystyle\frac{g}{2\pi^{2}(q-1)^{3}q}\left[3\psi^{(0)}\left(\frac{3}{q}-2\right)+\psi^{(0)}\left(\frac{1}{q}\right)\right.\left.-3\psi^{(0)}\left(\frac{2}{q}-1\right)-\psi^{(0)}\left(\frac{4}{q}-3\right)\right],$
(42) $\displaystyle\epsilon_{f,1}$ $\displaystyle=$
$\displaystyle\frac{g}{2\pi^{2}(q-1)^{3}q}\left[3\Phi\left(-1,1,\frac{2}{q}-1\right)-3\Phi\left(-1,1,\frac{3}{q}-2\right)+\Phi\left(-1,1,\frac{4}{q}-3\right)\right.$
(43) $\displaystyle\left.-\Phi\left(-1,1,\frac{1}{q}\right)\right],$
$\displaystyle\epsilon_{f,2}$ $\displaystyle=$
$\displaystyle-\frac{g}{4\pi^{2}(q-1)q}\left[\Phi\left(-1,1,\frac{2}{q}-1\right)\right.\left.-\Phi\left(-1,1,\frac{1}{q}\right)\right],$
(44) $\displaystyle P_{b,1}$ $\displaystyle=$
$\displaystyle\frac{\epsilon_{b,1}}{3},~{}P_{f,1}=\frac{\epsilon_{f,1}}{3},~{}P_{f,2}=\epsilon_{f,2}.$

In terms of the above variables, $\epsilon_{\rm{bag},\ell}$ and
$P_{\rm{bag},\ell}$ ($\ell=1,2$) are given by,

$\displaystyle\epsilon_{\mathrm{bag},1}$ $\displaystyle=$
$\displaystyle\epsilon_{b,1}+2\epsilon_{f,1},~{}P_{\mathrm{bag},1}=P_{b,1}+2P_{f,1},$
$\displaystyle~{}\epsilon_{\mathrm{bag},2}$ $\displaystyle=$ $\displaystyle
2\epsilon_{f,2},~{}P_{\mathrm{bag},2}=2P_{f,2}.$ (46)

##  VI Results and discussion

We present here the results on the fate of the nonlinear perturbations in the
QGP fluid at zero baryonic chemical potential which leads the vanishing
thermal conductivity of the medium. We solve the equations describing the
evolution of $\hat{\epsilon_{i}}$ for the initial profile of the perturbations
given by:

$\displaystyle\hat{\epsilon}_{1}$ $\displaystyle=$ $\displaystyle
A_{1}\big{[}\mathrm{sech}\frac{(x-10)}{B_{1}}\big{]}^{2}$ (47)
$\displaystyle\hat{\epsilon}_{2}$ $\displaystyle=$ $\displaystyle
A_{2}\big{[}\mathrm{sech}\frac{(x-10)}{B_{2}}\big{]}^{2}$ (48)

Figure 1: (Color online) a) $\hat{\epsilon}_{1}$ as a function of $x$ at
different time for $T=200MeV$ and $q=1.08$, $A_{1}=0.65,B_{1}=1$ fm. b)
$\hat{\epsilon}_{2}$ as a function of $x$ at different time for $T=200MeV$ and
$q=1.08$, $A_{2}=0.35,B_{2}=1$ fm. c)
$\hat{\epsilon}=1+\hat{\epsilon}_{1}+\hat{\epsilon}_{2}$ as a function of $x$
at different time for $T=200MeV$ and $q=1.08$, $A_{1}=0.65,B_{1}=1$ fm,
$A_{2}=0.35,B_{2}=1$ fm.

Firstly, we investigate the effects of non-extensivity only through the EoS
keeping value of the shear viscosity at the KSS bound KSS , $\eta/s=1/4\pi$.
The bulk viscosity is taken as $\zeta/s=1/4\pi$ and $\kappa=0$. Fig.1(a) and
Fig.1(b) show the time evolution of the spatial shape of the first-order
$(\hat{\epsilon_{1}})$ and the second-order $(\hat{\epsilon_{2}})$
perturbations respectively. The initial profile is given by the Eq.(48). In
Fig.1(a), the initial magnitude $(A)$ and width $(B)$ of the first-order
perturbation is taken as $A=0.65$,$B=1(fm)$, whereas in Fig.1(b), the initial
magnitude of the second-order perturbation is taken as $A=0.35$, and width is
kept same as the first-order i.e $B=1(fm)$. The value of $T$ and $q$ of the
background are taken as $T=200MeV$ and $q=1.08$. Comparison of results
depicted in Fig.1(a) and Fig.1(b) indicates that the dissipation of the
second-order perturbation is less as compared to first order. This can be
understood from the Eq.(34), where $\hat{\epsilon_{2}}$ gets contributions
from $\hat{\epsilon_{1}}$ which contains terms with third-order derivatives in
space, known as dispersive term, responsible for preserving the shape of
nonlinear wave propagations. Fig.1(c) shows the evolution of the energy
density up to second order scaled by the background energy density. Though the
actual QGP system produced in heavy-ion collisions has a lifetime of the order
of $10-15fm$, we have considered the time evolution up to $80fm$ to better
understand late time behaviour. It is evident from the results that with time,
amplitude decreases and width increases indicating towards the possibility of
complete dissipation, unlike the case studied earlier in TrambakNLW . However,
up to $15fm$, the dissipation of the nonlinear wave is not very drastic, and
the change in amplitude is faster in the beginning than later times and the
propagation speed is smaller at latter times. The slowing down of dissipation
can be understood from a comparison of the results of Fig.1(c), Fig.2,
Fig.3(a), and Fig.3(b). In Fig.1(c) and Fig.2, the amplitudes of the initial
perturbations are different with same value of width, whereas in Fig.1(c), the
magnitude is taken as $A=2$ and in Fig.2 it is considered to be $A=4$. It is
observed that with larger amplitude, the speed of propagation is larger, which
is clearly shown in Fig.2($t=80fm$ provides for better visibility). This is
due to amplitude dependent propagation of nonlinear waves, a well known
feature of the nonlinear propagation. So at latter times, with the decrease of
amplitude due to dissipation, the propagation gets slower in all the cases.

Figure 2: (Color online) $\hat{\epsilon}$ as a function of $x$ at different
time for $T=200MeV$ and $q=1.08$ (same as Fig.1(c) but with higher magnitude),
$A_{1}=1.65,B_{1}=1,A_{2}=1.35,B_{2}=1$. We see that with larger amplitude of
the perturbation the propagation speed is larger.

Figure 3: (Color online) a) $\hat{\epsilon}$ as a function of $x$ at different
time for $T=200MeV$ and $q=1.08$ (same as Fig.1(c) but with larger width),
$A_{1}=0.65,B_{1}=2,A_{2}=0.35,B_{2}=2$, We observe that more the width less
the dissipation of the wave. Also we find the propagation speed of the waves
is larger. b) Same as Fig.1(c), and
$A_{1}=0.65,B_{1}=0.5,A_{2}=0.35,B_{2}=0.5$, dissipation is more and the speed
of the propagation is less.

To compare the results displayed in Fig.3(a) and Fig.3(b) the widths are taken
as $B=2fm$ and $B=0.5fm$ respectively with the same value of the amplitude,
$A=2$. It is observed that the increase in width of initial perturbation
causes less dissipation of the perturbation with time for the same background
fluid conditions. The propagation speed is found to be greater for
perturbation with larger width. This is because, the perturbation with larger
width dissipates less, and thus due to amplitude dependent speed of
propagation, it propagates more in space. Now comparing Fig.2 with Fig.3 we
find similar propagation behaviour. For the increase in either amplitude or
the width, the damping is lesser, but we see a bit more shift in the peak of
the perturbation (at $80fm$) in results shown in Fig.2 than the results
displayed in Fig.3(a).

Now we move to discuss the role of the background fluid inscribed through the
$q$ and $T$ on the perturbations in Fig.4 and Fig.5. The results in Fig.4 show
the fate of the same initial perturbation as of Fig.1(c), after $t=10$ fm for
three different background temperatures with a fixed value of $q$. It shows
that with the larger background temperature the dissipation is less but the
stiffness of the front is more for the perturbation in the direction of
propagation in the dissipative background. However, for ideal background, the
magnitude as well as tilting remains unchanged with variation of background
temperature. This is in contrast to previous work on ideal hydrodynamic
background at fixed background energy density with varying temperature of the
propagating profile but in the present work we consider the background energy
density as a function of the temperature.

Figure 4: (Color online) $\hat{\epsilon}$ as a function of $x$ at different
time for $t=10$ fm and $q=1.08$, $A_{1}=0.65,B_{1}=1,A_{2}=0.35,B_{2}=1$. Here
we see that with larger values of temperature of the medium less the
dissipation but tilting of the wave is more along x axis. Figure 5: (Color
online) $\hat{\epsilon}$ as a function of $x$ at different time for $t=10$ fm
and $T=200MeV$, $A_{1}=0.65,B_{1}=1,A_{2}=0.35,B_{2}=1$. Similar as Fig.4,
where increasing $q$ values give less dissipation and tilting along x axis is
more for smaller values of $q$. Figure 6: (Color online) $\hat{\epsilon}$ as a
function of $x$ at different time for $t=10$ fm and $T=200MeV$,
$A_{1}=0.65,B_{1}=1,A_{2}=0.35,B_{2}=1$. Here we considered $\eta=\zeta\to 0$.
We observe the waves tend to loose their localization with time. Figure 7:
(Color online) A comparison on propagation of nonlinear waves is shown for
extensive vs non-extensive background .

Fig.5 shows the fate of the same initial perturbation as of Fig.1(c) after
$10$ fm time for three different values of $q$ with fixed temperature at
$T=200MeV$. It is observed that higher the $q$ values lower is the dissipation
and more is the stiffness of the front due to accumulation. This is in the
presence of the effect of dissipation or within the scope of $q$-hydrodynamics
which imitate the dissipative hydrodynamics. This is in contrast to the ideal
background case where tilting (equivalent to increase in stiffness of the
front) is higher for lower $q$ values, as reported in the earlier work
TrambakNLW  and the same result can be found by setting the transport
coefficient to zero here with the analytic solution of the resulting
equations. Moreover, Fig.6 shows that, with almost vanishing transport
coefficient, the waves loose their localization(breaking waves) and the shape
preserving nature no longer survives. It is also to be noted from results
shown in Fig.5 that for $q$ dependence only through EoS, the propagation speed
of the peak increases with $q$.

So far we have discussed the effects of non-extensivity on the propagation of
nonlinear waves only through the EoS. Since the EoS and transport coefficients
control the bulk evolution of a system, therefore, for a thorough treatment of
the propagation of nonlinear waves in a non-extensive
background($q$-background), one should consider the effects of non-
extensivity, not only through the EoS($q$-EoS) but also on the transport
coefficients($q$-viscosities) which are estimated for $q$\- background (see
Biro:2011bq ; Kadam:2015xsa  for details). The results with the
$q$-viscosities along with the $q$-EoS are presented in Fig.7 where a
comparison with the results with only the $q$-EoS is made. It is found that
the dissipation of nonlinear waves is more when non-extensivity is introduced
via the $q$-EoS and the $q$-viscosities together. Interestingly, here, the
propagation speed of the peak of the perturbations decrease with the
increasing value of $q$. Along with reducing the dissipation, the higher $q$
also reduces the speed of propagation. This means that introduction of non-
extensivity favours the non-dissipative (or anti-dissipative) slowing down of
nonlinear perturbations. It has been shown that for an extensive background,
the slowing down is caused by the dissipative effects Sarwar1 . This may be
useful as a distinguishing aspect of the extensive and the non-extensive
background. Interestingly, we also find that the dissipation is more for
smaller $q$-values (compare with Fig. 5). With the q-viscosities along with
q-EoS, the relative dissipation with $q$ is more than that of the q-EoS only.
This is in contrast to the expectation from the notion that the non-extensive
distribution is treated as a non-equilibrium deviation from the extensive one.
In such a scenario the dissipative correction is interpreted as the result of
arising due to $q\neq 1$. These findings are opposite to that expectation.
Rather, it indicates that $q$-equilibrium should not always be interpreted as
a dissipative non-equilibrium form of extensive statistics.

For a $q$-background, the higher $q$ is accompanied by less dissipation in the
system. This can be understood from the fact that the enhanced correlation in
a system can be represented by higher $q$-values LD ; Abe1990 . The
correlations enhance the fluctuations, which can also be represented by a
value of $q\neq 1$ Wilk00 . The correlations opposes the entropy generation in
a system and hence dissipation. Therefore, a higher value of $q$ parameter
representing a higher degree of correlation is expected to result in less
dissipation. In terms of the wave equations, this fact is manifested as the
variation of viscosity to entropy ratio($\eta/s$) with $q$ Biro:2011bq . The
$\eta/s$ decreases with increasing $q$ due to greater rate of change of non-
extensive entropy with $q$ than of the bulk and shear viscosity, thereby
reducing dissipation.

However, it is important to note that higher $q$ causes less dissipation,
through the second and higher order space derivatives in the wave equations,
where the effects of dissipation is present through the coefficients of the
derivatives. Without the presence of the dissipative terms, the $q$
alone(through $q$-EoS only) can not produce such dissipative nature, i.e, $q$
can only have a role in dissipation in the presence of explicit dissipative
terms. That is why without considering the $q$-viscosities(or any dissipative
effects arising from $d$-hydrodynamics), the authors in TrambakNLW  do not
find any dissipation with $q$-EoS alone. Also they have considered only the
first order of energy perturbation and found the equation to be of the
breaking-wave type, causing a tilting (breaking) of nonlinear waves without
any dissipation. Whereas, in this work we have obtained KdV like equation by
considering perturbations up to second-order with inclusion of all the
relevant transport coefficients and found to produce prominent dissipative
effects without much of the tilting.

In such scenario, one may be tempted to think that, perhaps, less dissipation
is the unique signature of presence of finite $q$(increasing $q$-values).
However, it may be noted that the less dissipation may be accounted through
the smaller values of the transport coefficients. So there may be an ambiguity
behind the observed smaller dissipation in a system, whether it is non-
extensive system with larger $q$ or it is extensive with lower viscosity.
However this ambiguity can be removed if other distinct effects are present
for these two separate situations. The role of higher $q$ as the cause of less
dissipation can be mimicked by reducing the values of the transport
coefficients, but it may possess other distinct effects which will help in
removing the ambiguity. In the case of nonlinear wave propagation, the nature
of the nonlinear wave will be the relevant distinct feature. With less
dissipation in extensive background (i.e without $q$-background) the nonlinear
wave will not preserve the solitonic nature simultaneously due to breaking of
the wave(see Fig.6), whereas the $q$ makes dissipation less with preserving
the solitonic nature(see Fig.7). Therefore, the causes for less dissipation
with Tsallis background can be distinguished with the propagation of a
nonlinear wave. This may pave the way for investigation whether the viscosity
of the QGP is really low or it has $q$-thermalized local equilibrium with not
so low transport coefficients. Since $p-p$ collision shows $q$-distributed
particle spectra, therefore, one would argue for consideration of
$q$-thermalized local equilibrium. In that case it is very important to remove
the ambiguity about viscosities of QGP measured by fitting the results with
models. This needs a close attention to look after. Another important point is
to note that, when the ideal $q$-hydrodynamics is interpreted as equivalent to
$d$-hydrodynamics, one should not take $q$ as a cause of dissipation in
general sense. Because, $q$ does not necessarily represents dissipation,
especially when the background is in $q$ equilibrium.

##  VII Summary and conclusion

In summary, we have derived the equations for the propagation of nonlinear
waves involving dissipative transport coefficients like shear viscosity, bulk
viscosity, and thermal conductivity. The effect of the $q-$statistics is
incorporated by taking into account the $q-$EoS and $q$-viscosities. The
equation for the second order perturbation is similar to the KdV equation
resulting in the solitonic nature of nonlinear waves. We observe dissipation
of the nonlinear waves at both the first and second order. But the dissipative
effect is more prominent for the first order perturbation. For larger
perturbations, we find speed of the propagation to be larger, whereas for
broader perturbations, the dissipation is comparatively less and propagation
is faster. As far as the background temperature of the medium is concerned,
dissipation is less for background with higher temperature. The effect of $q$
is evident from the results shown in Fig.5, where we see more dissipation of
nonlinear waves for smaller values of $q$. This is in line with the
interpretation of $q$ as representative of correlations as well as
fluctuations in a system. As larger $q$ values may imitate the similar
dissipative effects as that of smaller values of the transport coefficients,
one may not treat them in the same footing as they have other features like
the maintenance of the solitonic nature of the nonlinear waves. Moreover, it
is found that the non-extensivity favours the non-dissipative slowing down of
nonlinear perturbations. This is in contrast to the extensive background where
the slowing down comes with the dissipation Sarwar1 . This also may be a
useful as distinguishing aspect between the extensive and the non-extensive
background. So the nature of propagation of a nonlinear wave can be
potentially used to distinguish a non-extensive medium from an extensive one.
This could be especially helpful in situations where momentum distribution of
particles can not distinguish the nature of underlying statistics. This
indistinguishability may occur when there will be an uncertainty on whether
the equilibration is complete or not and when there are other sources of power
law distributions such as radiation from high momentum of particles e.g jets
through the QGP medium, instead of the presence of correlations that causes
non-extensivity of the medium.

##  VIII Acknowledgement

MH, MR would like to thank Department of Higher Education, Govt. of West
Bengal, India. TB acknowledges partial support from the joint project between
the JINR and IFIN-HH.

## References

  * (1) V. Kirtskhalia,“Correct Definition of Sound Speed and Its Consequences in the Tasks of Hydrodynamics”, Journal of Fluids, vol. 2016, Article ID 4519201, 9 pages, 2016.
  * (2) W. Lick, “Nonlinear Wave Propagation in Fluids,” Annual Review of Fluid Mechanics, Vol. 2, Annual Reviews, Inc., Palo Alto, California, 1970.
  * (3) J. C. Solana, E. Shuryak, D. Teaney, J. Phys.: Conf. Ser. 27, 003 (2005).
  * (4) B. Betz, [arXiv:0910.4114 [nucl-th]].
  * (5) E. Shuryak, Phys. Rev. C 80, 069902 (2009).
  * (6) P. Staig, E. Shuryak, Phys. Rev. C 84, 034908 (2011), Phys. Rev. C 84, 044912 (2011).
  * (7) G. Sarwar and J. Alam, Int. J. Mod. Phys. A 33 (2018) no.08, 1850040 doi:10.1142/S0217751X18500409 [arXiv:1503.06019 [nucl-th]].
  * (8) A. Rafiei, K. Javidan, Phys. Rev. C 94, 034904 (2016).
  * (9) M. Hasanujjaman, M. Rahman, A. Bhattacharyya, J. Alam, Phys. Rev. C 102, 034910 (2020).
  * (10) M. Hasanujjaman, G. Sarwar, M. Rahman, A. Bhattacharyya, J. Alam, arXiv:2008.03931v2.
  * (11) Y. Minami, T. Kunihiro, Prog. Th. Phys. 122, 881 (2009).
  * (12) M. Rahaman and J. Alam, Phys. Rev. C 97 (2018) 054906; arXiv: 1712.09175 [nucl-th].
  * (13) S. Raha, K. Wehrberger, R.M. Weiner, Nucl. Phys. A 433, 427 (1984).
  * (14) G.N. Fowler, S. Raha, N. Stelte, R.M. Weiner, Phys. Lett. B 115, 286 (1982).
  * (15) S. Raha and R. M. Weiner, Phys. Rev. Lett. 50, 407 (1983).
  * (16) E. F. Hefter, S. Raha and R. M. Weiner, Phys. Rev. C 32, 2201 (1985).
  * (17) D. A. Fogaca, L. G. Ferreira Filho, F. S. Navarra, Phys. Rev. C 81, 055211 (2010).
  * (18) D.A. Fogaca, H. Marrochio, F.S. Navarra, J. Noronha, Nucl. Phys. A 934, 18 (2015).
  * (19) G. Sarwar, M. Hasanujjaman, M. Rahaman, A. Bhattacharyya and J. Alam, Phys. Lett. B 820 (2021), 136583: [arXiv:2012.12668 [nucl-th]].
  * (20) C Tsallis, J. Stat. Phys. 52, 479 (1988).
  * (21) C Tsallis, _Introduction to non extensive statistical mechanics_ (Springer Science+Business Media, 2009).
  * (22) V. Khachatryan et al. (CMS collaboration), Phys. Rev. Lett. 105, 022002 (2010).
  * (23) S. Acharya et al. (ALICE collaboration), Phys. Rev. C 97, 024615 (2018).
  * (24) T. S. Biró, G. Purcsel, K. Ürmössy, Eur. Phys. J. A 40 325 (2009).
  * (25) J. Cleymans et al., Phys. Lett. B 723 351 (2013).
  * (26) L. Marques, J. Cleymans, and A. Deppman, Phys. Rev. D 91, 054025 (2015).
  * (27) T. Bhattacharyya et al., Eur. Phys. J. A 52 30 (2016).
  * (28) S. Tripathy et al., Eur. Phys. J. A 52, no. 9 289 (2016).
  * (29) S. Grigoryan, Phys. Rev. D 95 056021 (2017).
  * (30) T. Bhattacharyya et al., J. Phys. G 45, no. 5, 055001 (2018).
  * (31) M. D. Azmi, T. Bhattacharyya, J. Cleymans and M. Paradza, J. Phys. G 47, 045001 (2020)
  * (32) C.Y. Wong et al., Phys. Rev. D 91, 114027 (2015).
  * (33) M.D. Azmi, J. Cleymans, Eur. Phys. J. C 75, 430 (2015).
  * (34) T. Bhattacharyya et al., J. Phys. G: Nucl. Part. Phys. 45, 055001 (2018).
  * (35) J. Cleymans, D. Worku, Eur. Phys. J. A 48, 160 (2012), J. Phys. G 39, 025006 (2012).
  * (36) H. Zheng and Lilin Zhu, Advances in High Energy Phys. 2016, 9632126 (2016).
  * (37) G. Wilk, and Z. Włodarczyk, _Interpretation of the Nonextensivity Parameter $q$ in Some Applications of Tsallis Statistics and Lèvy Distributions, Phys. Rev. Lett._ 84 (2000) 2770.
  * (38) C. Tsallis and D.J. Bukman, _Anomalous diffusion in the presence of external forces: Exact time-dependent solutions and their thermostatistical basis, Phys. Rev. E_ 54 (1996) R2197(R).
  * (39) B.M. Boghosian, P.J. Love, P.V. Coveney, I.V. Karlin, S. Succi, and J. Yepez, _Galilean-invariant lattice-Boltzmann models with H theorem, Phys. Rev. E_ 68 (2003) 025103(R).
  * (40) P. Douglas, S. Bergamini, and F. Renzoni, _Tunable Tsallis Distributions in Dissipative Optical Lattices, Phys. Rev. Lett._ 96 (2006) 110601.
  * (41) A. Deppman, E. Megias and D. P. Menezes, _Fractal structure of Yang-Mills fields_ , Phys. Scripta 95 (2020) 094006; arXiv: 2002.12667 [hep-th].
  * (42) A. Deppman, E. Megías and D. P. Menezes, _Fractals, nonextensive statistics, and QCD_ Phys.Rev.D 101 (2020) 034019.
  * (43) A. Deppman, E. Megias and D. P. Menezes,_Fractal Structures of Yang-Mills Fields and Non Extensive Statistics: Applications to High Energy Physics_ ,
  * (44) J. A. S. Lima and A. Deppman, _Tsallis meets Boltzmann: q-index for a finite ideal gas and its thermodynamic limit_ , Phys. Rev. E 101 (2020) 040102(R).
  * (45) J. Rozynek and G. Wilk, _An example of the interplay of nonextensivity and dynamics in the description of QCD matter”_ , Eur.Phys.J.A 52 (2016) 294.
  * (46) J. Rożynek and G. Wilk, _Nonextensive quasiparticle description of QCD matter_ , Symmetry 11 (2019) 401.
  * (47) T. Bhattacharyya, J. Cleymans, and S. Mogliacci, _Analytic results for the Tsallis thermodynamic variables, Phys. Rev. D_ 94 (2016) 094026.
  * (48) F. S. Navarra, O. V. Utyuzh, G. Wilk, and Z. Włodarczyk, Phys. Rev. D 67, 114002 (2003); Physica A 344, 568 (2004); Nukleonika 49, (Supplement 2), S19 (2004).
  * (49) F. S. Navarra, O. V. Utyuzh, G. Wilk, and Z. Włodarczyk, Physica A 340, 467 (2004).
  * (50) M. Biyajima, M. Kaneyama, T. Mizoguchi, and G. Wilk, Eur. Phys. J. C 40, 243 (2005); M. Biyajima, T. Mizoguchi, N. Nakajima, N. Suzuki, and G. Wilk, ibid. 48, 597 (2006).
  * (51) G. Wilk and Z. Wlodarczyk, Chaos Solitons Fractals 13, 581-594 (2002): [arXiv:hep-ph/0004250 [hep-ph]]..
  * (52) G. Wilk, and Z. Włodarczyk, Phys. Rev. C 79 054903 (2009).
  * (53) T.S. Biró, and A. Jakovác, Phys. Rev. Lett. 94 132302 (2005).
  * (54) S. Deb, G. Sarwar, R. Sahoo and J. Alam, Eur. Phys. J. A 57 (2021): [arXiv:1909.02837 [hep-ph]].
  * (55) A. Das, W. Florkowski, R. Ryblewski and R. Singh: [arXiv:2012.05662 [hep-ph]].
  * (56) A. Das, W. Florkowski, R. Ryblewski and R. Singh, Phys. Rev. D 103 (2021): [arXiv:2103.01013 [nucl-th]].
  * (57) A. Das, W. Florkowski, R. Ryblewski and R. Singh, [arXiv:2105.02125 [nucl-th]].
  * (58) T. Osada and G. Wilk, Phys. Rev. C 77, 044903 (2008) [erratum: Phys. Rev. C 78, 069903 (2008)]: [arXiv:0710.1905 [nucl-th]].
  * (59) T. Bhattacharyya and A. Mukherjee, Eur. Phys. J. C 80, 656 (2020): [arXiv:2003.10692 [nucl-th]].
  * (60) T. S. Biro and E. Molnar, Phys. Rev. C 85, 024905 (2012): [arXiv:1109.2482 [nucl-th]].
  * (61) C. Eckart, Phys. Rev. 58, 919 (1940).
  * (62) L. D. Landau and E. M. Lifshitz, Fluid Mechanics (AddisonWesley, Boston, 1959).
  * (63) F. S. Bemfica, M. M. Disconzi and J. Noronha, Phys. Rev. D 98, 104064 (2018).
  * (64) F. S. Bemfica, M. M. Disconzi and J. Noronha, Phys. Rev. D 100, 104020 (2019).
  * (65) F. S. Bemfica, M. M. Disconzi and J. Noronha: [arXiv:2009.11388 [gr-qc]].
  * (66) P. Kovtun, JHEP 10, 034 (2019).
  * (67) A. Das, W. Florkowski, J. Noronha and R. Ryblewski, Phys. Lett. B 806, 135525 (2020).
  * (68) A. Das, W. Florkowski, R. Ryblewski, Phys. Rev. D 102, 031501 (2020), Phys. Rev. D 103, 014011 (2021).
  * (69) W. Israel, J. M. Stewart, Ann. Phys. (NY) 118, 341 (1979).
  * (70) W. A. Hiscock, L. Lindblom, Annals of Physics, 151, 466-496 (1983).
  * (71) T. Koide, G. S. Denicol, P. Mota and T. Kodama, Phys. Rev. C 75, 034909 (2007): [arXiv:hep-ph/0609117 [hep-ph]].
  * (72) P. Romatschke, Int. J. Mod. Phys. E 19, 1-53 (2010): [arXiv:0902.3663 [hep-ph]].
  * (73) P. Van and T. S. Biro, Eur. Phys. J. ST 155, 201-212 (2008): [arXiv:0704.2039 [nucl-th]].
  * (74) K.Tsumura and T.Kunihiro, Phys. Lett. B 690, 255-260 (2010): [arXiv:0906.0079 [hep-ph]].
  * (75) A. Muronga, Phys. Rev. C 69, 03490 (2004); Phys. Rev. Lett. 88, 062302 (2002); Phys. Rev. Lett. 89, 159901 (2002) (erratum).
  * (76) H. Washimi, T. Taniuti, Phys. Rev. Lett. 17, 996 (1966).
  * (77) R.C. Davidson, “Methods in Nonlinear Plasma Theory”, Academic Press, New York an London, (1972).
  * (78) H. Leblond, J. Phys. B: At. Mol. Opt. Phys. 41, 043001 (2008).
  * (79) R. A. Kraenkel, M.A. Manna, J.C. Montero, J.G. Pereira, arXiv:patt-sol/9509003v1
  * (80) M. Rahaman, T. Bhattacharyya and J. Alam, Int. J. Mod. A36 (2021) 2150154; arXiv: 1906.02893 [hep-ph].
  * (81) A. Chodos, R. L. Jaffe, K. Johnson, C. B. Thorn, and V. F. Weisskopf, Phys. Rev. D 9, 3471 (1974).
  * (82) A. Lavagno, D. Pigato, and P. Quarati, J. Phys. G: Nuclear and Particle Physics 37, 11 (2010).
  * (83) P. H. G. Cardoso, T. N. da Silva, A. Deppman and D. P. Menezes, Eur. Phys. J. A 53, 191 (2017).
  * (84) P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 94, 111601 (2005): [arXiv:hep-th/0405231 [hep-th]].
  * (85) G. P. Kadam and H. Mishra, Phys. Rev. C 92, no.3, 035203 (2015): [arXiv:1506.04613 [hep-ph]].
  * (86) S. Abe, Physica A, 269,403(1990).