# The $\pi f_{0}(500)$ decay of the $a_{1}(1260)$

R. Molina<EMAIL_ADDRESS>Departamento de Física Teórica and IFIC,
Centro Mixto Universidad de Valencia-CSIC Institutos de Investigación de
Paterna, Aptdo.22085, 46071 Valencia, Spain  M. Döring<EMAIL_ADDRESS>Department of Physics, The George Washington University, Washington, DC 20052,
USA  Theory Center, Thomas Jefferson National Accelerator Facility, Newport
News, Virginia 23606, USA  W. H. Liang<EMAIL_ADDRESS>Department of
Physics, Guangxi Normal University, Guilin 541004, China  Guangxi Key
Laboratory of Nuclear Physics and Technology, Guangxi Normal University,
Guilin 541004, China  E. Oset<EMAIL_ADDRESS>Department of Physics, Guangxi
Normal University, Guilin 541004, China  Departamento de Física Teórica and
IFIC, Centro Mixto Universidad de Valencia-CSIC Institutos de Investigación de
Paterna, Aptdo.22085, 46071 Valencia, Spain

###### Abstract

We evaluate the $a_{1}(1260)\to\pi\sigma(f_{0}(500))$ decay width from the
perspective that the $a_{1}(1260)$ resonance is dynamically generated from the
pseudoscalar-vector interaction and the $\sigma$ arises from the pseudoscalar-
pseudoscalar interaction. A triangle mechanism with $a_{1}(1260)\to\rho\pi$
followed by $\rho\to\pi\pi$ and a fusion of two pions within the loop to
produce the $\sigma$ provides the mechanism for this decay under these
assumptions for the nature of the two resonances. We obtain widths of the
order of $13-22$ MeV. Present experimental results differ substantially from
each other, suggesting that extra efforts should be devoted to the precise
extraction of this important partial decay width, which should provide
valuable information on the nature of the axial vector and scalar meson
resonances and help clarify the role of the $\pi\sigma$ channel in recent
lattice QCD calculations of the $a_{1}$.

††preprint: JLAB-THY-21-3464

##  I Introduction

The nature of the axial vector meson resonances is a subject of continuous
debate. It is well known that quark models, very successful in correlating a
vast amount of data concerning meson spectra, have difficulties reproducing
the properties of axial vector mesons Isgur ; Vijande . On the other hand, the
advent of the chiral unitary approach, combining chiral dynamics with a
unitary approach in coupled channels, showed that the lowest-energy axial
vector mesons can be fairly well reproduced from the interaction of vector
mesons with pseudoscalars Lutz ; Luis ; GengRoca ; RenGeng . The interaction
used between vector-pseudoscalar ($VP$) channels is taken from the chiral
Lagrangian of Ref. Birse . In Ref. Nagahiro , for the particular case of the
$\rho\pi$ interaction, and in Ref. Jorgivan  for the general $VP$ interaction,
it was shown that this Lagrangian can be obtained from the exchange of vector
mesons within the local hidden gauge approach hidden1 ; hidden2 ; hidden4 ;
Nagahiro . This conclusion is in line with the findings in the pseudoscalar-
pseudoscalar ($PP$) interaction Rafael . However,the $PV$ interaction shows a
peculiarity since one can make the transition $PV\to VP$ mediated by a
pseudoscalar exchange of one $P$ from the decay of $V\to PP$. Fortunately,
such contributions are very small compared to those from vector exchange, as
recently shown in Ref. Jorgivan , something corroborated in the study of $D$
decays into three mesons looking at meson final state interactions Nakamura .
Thus, one can rely on the information for the axial vector mesons obtained
within the chiral unitary approach with the chiral Lagrangian of Ref. Birse .
Lattice QCD has also shown interest in describing the axial vector resonances.
Concretely, the $a_{1}(1260)$ was already studied in Ref. Sasa  (see Ref.
RocaFinite  in connection to this work). More detailed work on the
$a_{1}(1260)$ and $b_{1}(1235)$ is shown in Ref. Lang , and the $b_{1}(1235)$
is also investigated along this line considering the coupled channels
$\pi\omega,\pi\phi$ in Ref. Dudek . The first extraction of a three-body
resonance from lattice QCD, using up to three pion operators and a suitable
three-body finite-volume formalism was achieved recently for the example of
the $a_{1}(1260)$ Mai:2021nul . The finite-volume mapping (FVU) relies on
unitarity for three particles Mai:2017vot  leading to the aforementioned
pseudoscalar exchange which develops singularities in the finite volume
Mai:2017bge . The largest source of systematics in the $a_{1}(1260)$ pole
extraction from lattice QCD is the role of the $\pi\sigma$ channel Mai:2021nul
. The prediction of the partial decay width $a_{1}\to\pi\sigma$ in the present
paper is, therefore, essential to assess the role of such sub-dominant
channels in upcoming lattice QCD calculations.

As these examples show, combined work using effective field theories, quark
models, lattice QCD simulations and experiment will ultimately unveil the
nature of the axial vector meson resonances, and work goes along in these
direction. One of the successful strategies pursuing this aim is to look for
particular decay modes of these states which reveal the coupling of the
resonances to some bound meson-meson components. Let us take as an example the
$f_{1}(1285)$, which appears as a $K^{*}\bar{K}+c.c.$ single-channel state,
bound by about $100\textrm{ MeV}$ Lutz ; Luis ; RenGeng . The $f_{1}(1285)$
can naturally decay to $K^{*}\bar{K}$ (due to the $K^{*}$ width), then the
$K^{*}$ decays to $K\pi$ and $K\bar{K}$ fuses into $a_{0}(980)$, which couples
strongly to $K\bar{K}$. This process constitutes a triangle mechanism that
leads to the $f_{1}(1285)\to\pi a_{0}$ decay which is tied to the coupling of
the $f_{1}(1285)$ to the bound $K^{*}\bar{K}$ component. This work is done in
Ref. Aceti  and a good agreement with experiment is found. In a similar way
one can study the $f_{1}(1285)\to\pi K\bar{K}$ decay, and a good agreement
with experiment is found in Ref. AcetiKK .

The triangle mechanisms have been widely used to test molecular components,
with the information encoded in the resonance decay vertex of the first step,
or on the fusion vertex of the last step. In Refs. Zou ; Acetieta  a triangle
mechanism was used to explain the unexpectedly large isospin breaking of the
$\eta(1405)\to\pi^{0}f_{0}(980)$ decay compared to
$\eta(1405)\to\pi^{0}a_{0}(980)$ observed at BES BESeta  (see also following
works of Refs. Achasov ; QiangZhao ). Recent work on triangle mechanisms to
describe different processes can be seen in Refs. He ; HuangGeng ; XieChi ;
Xiemore ; HuangXie ; YuKe ; GengXie .

One particular case of a triangle mechanisms is the one that develops a
triangle singularity (TS). This was early discussed in Refs. Karplus ; Landau
, which showed that a singularity in the triangle diagram emerged when the
three particles in the loop are placed simultaneously on shell, while being
collinear and the decay of the first step and fusion of the second one are
possible at a classical level (Coleman-Norton theorem Norton ). The topic has
experienced a rebirth recently since one finds nowadays many physical examples
of reactions which are interpreted in terms of a TS. A modern reformulation of
the problem, intuitive and practical, is given in Ref. BayarGuo  and a recent
review on the topic is given in Ref. Guorev .

Many examples of TS are discussed in Ref. Guorev . Here we just mention two
examples which have particular relevance. One of them is the TS developed in
the $a_{1}(1260)$ decay into $\pi f_{0}(980)$, related to the topic that we
face in this paper (the $a_{1}(1260)\to\pi\sigma(f_{0}(500))$ decay) which was
originally associated to a new resonance “$a_{1}(1420)$” by the COMPASS
Collaboration COMPASS , but which was soon interpreted in terms of a TS Oka ;
Mikha ; AceDai  and acknowledged as such by the COMPASS Collaboration Alexeev
. The other recent example is the case of the $pp\to\pi^{+}d$ fusion reaction,
which, as shown in Ref. IkenoRaq , has an unexpectedly large strength due to a
TS. This work, together with the one of Ref. RaquelIke , present an
alternative explanation of the peak observed in the $pn\to\pi^{0}\pi^{0}d$
($pn\to\pi^{+}\pi^{-}d$) reaction, so far associated to a dibaryon
“$d^{*}(2380)$” dibaryon .

It is interesting to recall that not all triangle diagrams develop a TS,
actually these are exceptional cases, but it is easy to test if a triangle
diagram develops a TS applying a simple equation (see Eq. (18) of Ref.
BayarGuo ). We shall take advantage of this equation too, because it not only
serves to see if there is a TS, but also to have a feeling of the strength of
the loop by observing how far this equation is from being fulfilled.

In the present work we address the issue of the partial decay width of the
$a_{1}(1260)$ into $\pi\sigma(f_{0}(500))$. This can be seen for instance in
the $a^{+}_{1}(1260)\to\pi^{+}\pi^{0}\pi^{0}$ decay, where the
$\pi^{0}\pi^{0}$ will come from the $\sigma$ meson, hence, this is part of the
three pion decay of the $a_{1}(1260)$. This latter problem has been
extensively studied from different points of view and it is not our purpose to
go through it again. A thorough account and discussions on the different
approaches can be seen in the works of Refs. Misha ; ZhangXie ; DaiRoca , but
in none of those works the $\pi\sigma$ decay mode has been isolated.

Experimentally there are only a few experiments tabulated in the PDG PDG
(although “they are not taken for averages, fits, limts, etc”) and the
branching ratios are

$\Gamma[a_{1}\to\pi f_{0}(500),f_{0}\to\pi\pi]/\Gamma_{\rm tot}=(18.76\pm
4.29\pm 1.48)\times 10^{-2}~{}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Asner}{\@@citephrase{(}}{\@@citephrase{)}}}},$ (1)
$\displaystyle\Gamma[a_{1}\to\pi
f_{0}(500),f_{0}\to\pi\pi]/\Gamma[(\rho\pi)_{S{\rm{-wave}}},\rho\to\pi\pi]=\left\\{\begin{array}[]{l}~{}~{}(6\pm
5)\times 10^{-2}~{}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Salvini}{\@@citephrase{(}}{\@@citephrase{)}}}},\\\ \sim
30\times 10^{-2}~{}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Akhme}{\@@citephrase{(}}{\@@citephrase{)}}}},\\\
~{}~{}(0.3\pm 0.3)\times 10^{-2}~{}\text{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Longa}{\@@citephrase{(}}{\@@citephrase{)}}}},\end{array}\right.$
(5)

and in Ref. Chung  this decay mode is reported as “seen”. As we can see, the
decay ratios are very different, and with large uncertainties. The data of
Ref. Asner  come from the $\tau^{-}\to\nu_{\tau}\pi^{-}\pi^{0}\pi^{0}$ decay
measured at CLEO II. These data should be taken with certain caution since
$m_{\sigma}=860\textrm{ MeV}$ and $\Gamma_{\sigma}=880\textrm{ MeV}$ are used
in the substructure fits, which are very different from present values PDG ;
Colangelo ; sigma . 111From Ref. Colangelo  we find  $m_{\sigma}=470\pm
30\textrm{ MeV},~{}~{}~{}~{}~{}\Gamma_{\sigma}=2\times(295\pm 20)\textrm{
MeV},$ and from Ref. sigma  $m_{\sigma}=449^{+22}_{-16}\textrm{
MeV},~{}~{}~{}~{}~{}\Gamma_{\sigma}=2\times(275\pm 12)\textrm{ MeV}.$ We take
values close to these as $m_{\sigma}=458$ MeV and $\Gamma_{\sigma}=464$ MeV,
coming from the use of Eqs. (35), (36).  Dalitz plots with a more realistic
$\sigma$ resonance were predicted by the three-body framework of Ref.
Kamano:2011ih . That work and Ref. Misha  use two-body amplitudes for the sub-
channels in contrast to earlier work Janssen:1993nj  in this direction.
However, none of these frameworks or related theoretical studies Wagner:2007wy
; Wagner:2008gz  have performed a refit of the CLEO Dalitz distributions Asner
which are among the most precise sources of information on the $a_{1}$ decay
modes. The ALEPH data on the lineshape ALEPH:2005qgp , in contrast, served for
the determination of the pole position of the $a_{1}$ Janssen:1993nj ; Mikha ;
Mai:2021nul  but does not contain information on branching ratios.

In order to evaluate the $a_{1}(1260)\to\pi\sigma$ decay width, we follow the
procedure of Ref. ZhangXie  in the study of
$a_{1}(1260)\to\pi^{+}\pi^{+}\pi^{-}$ decay, where the tree level
$a_{1}(1260)\to\rho\pi$ together with final state interaction of two pions via
a triangle loop mechanism are considered. The formalism follows closely the
one of Refs. Aceti ; AceDai  and is based on the chiral unitary approach for
the axial vector mesons Luis  and for scalar mesons npa . The present work has
to be considered as complementary to the one of Ref. ZhangXie , where the
$\pi^{+}\pi^{-}$ mass distribution is obtained, that is largely dominated by
the $a_{1}\to\rho\pi$ tree level mechanism, and that is shown to be consistent
with the experimental data of ARGUS ARGUS . In the low energy part of the
$\pi^{+}\pi^{-}$ invariant mass spectrum, some extra strength appears which
should be attributed to the $\pi^{+}\sigma$ production but the $\pi\sigma$
decay mode was not isolated. The precise evaluation of this decay mode is the
purpose of the present work.

##  II Formalism

We follow closely the approach of Ref. Aceti  for the $f_{1}(1285)\to\pi
a_{0}(980)$ decay. In that case the $f_{1}(1285)$ decays to $K^{*}\bar{K}$,
the $K^{*}$ to $K\pi$ and then the $K\bar{K}$ fuse to give the $a_{0}(980)$.
In the present case, we have the $a_{1}(1260)$ coupled to $\rho\pi$ and
$\bar{K}^{*}K-K^{*}\bar{K}$, with the largest coupling to $\rho\pi$ Luis . In
view of this, and similar to Ref. ZhangXie , the mechanisms to produce the
$\sigma$ are given by the diagrams of Fig. 1.

Figure 1: Diagrams leading to the $a_{1}(1260)\to\pi\sigma$ decay: (a) through
$a_{1}\to\pi\rho$ in the first step; (b) through $a_{1}\to K^{*}\bar{K}$ in
the first step; (c) through $a_{1}\to\bar{K}^{*}K$ in the first step.

One can see that none of these triangle mechanisms produces a TS BayarGuo ,
yet the diagrams Fig. 1 (b) and Fig. 1 (c) are highly suppressed with respect
to diagram Fig. 1 (a). Indeed, in diagram Fig. 1 (a) one can place the $\rho$,
$\pi$, $\pi$ intermediate states on shell, yet the condition of the Coleman-
Norton theorem Norton  is not satisfied, since, in a collinear decay, the
$\pi$ coming from the $\rho$ decay has a smaller velocity than the one from
$a_{1}$ decay, hence, it cannot catch up with the latter pion to produce the
$\sigma$. The processes in Fig. 1 (b) and Fig. 1 (c) are much weaker since,
first, the $a_{1}\to\bar{K}^{*}K$ coupling is smaller than the
$a_{1}\to\rho\pi$ one Luis , second, because the $\sigma$ couples strongly to
$\pi\pi$ but very weakly to $K\bar{K}$ npa  and, third, because the $K\bar{K}$
state is highly off shell in the loop if we wish to produce the $\sigma$.
Hence, we shall not consider diagrams Fig. 1 (b) and Fig. 1 (c) in our
evaluation.

It is curious to note that, conversely, the diagrams Fig. 1 (b) and Fig. 1 (c)
are those responsible for the $a_{1}(1260)\to\pi f_{0}(980)$ decay, replacing
the $\sigma$ by the $f_{0}(980)$ AceDai . This is because in this case the
diagrams develop a TS, and also the $f_{0}(980)$ coupling to $K\bar{K}$ is
much stronger than the coupling of the $\sigma$ to $K\bar{K}$ nsd . This decay
mode was the one originally assigned to the $a_{1}(1420)$ resonance, now
accepted as a TS Alexeev . Yet, the contribution of this mode to the
$a_{1}(1260)$ width is small, of the order of 1 MeV AceDai  which gives us an
idea of the smallness of the diagrams in Fig. 1 (b) and Fig. 1 (c) in the much
less favored actual situation of $\sigma$ production.

We, thus, proceed with the diagram of Fig. 1(a) to produce $\pi\sigma$ and we
show in detail in Fig. 2 the diagrams to consider.

Figure 2: Detail of diagram (a) of Fig 1 showing explicitly the
$\rho^{+}\pi^{0}$ and $\rho^{0}\pi^{+}$ decays of the $a_{1}$. Particle
momenta are shown in parentheses.

We need the states of the $\sigma$ and $a_{1}^{+}$ in terms of $\pi$ and
$\rho$ multiplets $(-\pi^{+},\pi^{0},\pi^{-})$,
$(-\rho^{+},\rho^{0},\rho^{-})$,

$\displaystyle\ket{\pi\pi,\,I=0}$
$\displaystyle=\frac{-1}{\sqrt{3}}\ket{\pi^{+}\pi^{-}+\pi^{-}\pi^{+}+\pi^{0}\pi^{0}}\
,$ (6) $\displaystyle\ket{\rho\pi,\,I=1,I_{3}=1}$
$\displaystyle=\frac{1}{\sqrt{2}}\ket{\rho^{0}\pi^{+}-\rho^{+}\pi^{0}}\ .$ (7)

We take the coupling $g_{\sigma,\pi\pi}$ of the $\sigma$ to $\pi\pi$ from the
chiral unitary approach npa  with a single channel $\pi\pi$ in $I=0$, but we
also consider the $\pi\pi$, $I=0$ amplitudes in coupled channels. Further
details will be given at the beginning of Sec. III and Appendix A. In the
chiral unitary approach, for convenience, one works with a unitary
normalization where Eq. (6) has a coefficient $1/\sqrt{6}$ instead of
$1/\sqrt{3}$. This $\sqrt{2}$ factor has to be restored at the end for the
external pions so that

$g_{\sigma,\pi\pi}=\sqrt{2}\,(1248-i\,2717)~{}\mathrm{MeV}\ .$ (8)

Similarly, from Ref. Luis , in the $a_{1}$ rest frame, the decay vertex reads

$\displaystyle-it_{a_{1},\rho\pi}$
$\displaystyle=-ig_{a_{1},\rho\pi}\vec{\epsilon}_{a_{1}}\cdot\vec{\epsilon}_{\rho}\
,$ (9)

where $\vec{\epsilon}_{a_{1}},\,\vec{\epsilon}_{\rho}$ are the polarization
vectors of the $a_{1},\,\rho$, with

$g_{a_{1},\rho\pi}=(-3795+i\,2330)~{}\mathrm{MeV}\ .$ (10)

Thus,

$\displaystyle g_{a_{1},\rho^{0}\pi^{+}}=\frac{1}{\sqrt{2}}g_{a_{1},\rho\pi}\
,\quad g_{a_{1},\rho^{+}\pi^{0}}=-\frac{1}{\sqrt{2}}g_{a_{1},\rho\pi}\ .$ (11)

Analogously, we have (see Appendix A)

$\displaystyle g_{\sigma,\pi^{+}\pi^{-}}=\frac{1}{\sqrt{3}}g_{\sigma,\pi\pi}\
,\quad g_{\sigma,\pi^{0}\pi^{0}}=\frac{1}{\sqrt{3}}g_{\sigma,\pi\pi}\ .$ (12)

We also need the $\rho\to\pi\pi$ Lagrangian which is given by

$\displaystyle{\cal L}=-ig\braket{V^{\mu}\left[P,\partial_{\mu}P\right]}\
,\quad g=\frac{M_{V}}{2f}\ ,\quad(M_{V}\approx 800\text{ MeV, }f=93\text{
MeV})\ ,$ (13)

with $V,P$ the SU(3) matrices for the vector mesons and pseudoscalars,
respectively Luis , which gives

$\displaystyle-it_{\rho^{+},\pi^{+}\pi^{0}}$
$\displaystyle=i\sqrt{2}g\left(2\vec{p}_{\pi}+\vec{q}\,\right)\cdot\vec{\epsilon}_{\rho}\
,$ $\displaystyle-it_{\rho^{0},\pi^{+}\pi^{-}}$
$\displaystyle=-i\sqrt{2}g\left(2\vec{p}_{\pi}+\vec{q}\,\right)\cdot\vec{\epsilon}_{\rho}\
,$ (14)

where we keep only the spatial components of the $\rho$ polarization, given
the coupling of Eq. (9) and the fact that the $\epsilon^{0}$ component
provides only a very small contribution even for $\rho$ mesons with the
momenta involved in the $a_{1}\to\rho\pi$ decay (see Appendix A of Ref.
SakaiRamos ).

With these ingredients we obtain the amplitudes corresponding to the diagrams
of Fig. 2 (a1) and (a2):

$\displaystyle-it^{(a1)}$
$\displaystyle=\int\frac{d^{4}q}{(2\pi)^{4}}\,(-i)\,\left(-\frac{1}{\sqrt{2}}\right)\,g_{a_{1},\rho\pi}\,\vec{\epsilon}_{a_{1}}\cdot\vec{\epsilon}_{\rho}\,ig\sqrt{2}\left(2\vec{p}_{\pi}+\vec{q}\,\right)\cdot\vec{\epsilon}_{\rho}\,(-i)\left(-\frac{1}{\sqrt{3}}\right)g_{\sigma,\pi\pi}\,\Pi(q)\
,$ (15) $\displaystyle-it^{(a2)}$
$\displaystyle=\int\frac{d^{4}q}{(2\pi)^{4}}\,(-i)\,\left(\frac{1}{\sqrt{2}}\right)\,g_{a_{1},\rho\pi}\,\vec{\epsilon}_{a_{1}}\cdot\vec{\epsilon}_{\rho}\,ig(-\sqrt{2})\left(2\vec{p}_{\pi}+\vec{q}\,\right)\cdot\vec{\epsilon}_{\rho}\,(-i)\left(-\frac{1}{\sqrt{3}}\right)g_{\sigma,\pi\pi}\,\Pi(q)\
,$ (16)

where $\Pi$ is the product of propagators

$\displaystyle\Pi(q)$
$\displaystyle=\frac{1}{2\omega_{\rho}}\,\frac{i}{P^{0}-q^{0}-\omega_{\rho}+\frac{i\Gamma_{\rho}}{2}}\,\frac{i}{q^{2}-m_{\pi}^{2}+i\epsilon}\,\frac{i}{(P-q-
p_{\pi})^{2}-m_{\pi}^{2}+i\epsilon}\ ,$ (17)

with $\omega_{\rho}^{2}=m_{\rho}^{2}+\vec{q}^{\,\,2}$. In Eq. (17) we have
kept the positive-energy part of the $\rho$ propagator, corresponding to a
heavy particle that will be mostly on shell in the loop, but keep the full
propagators for the two pions in the loop. Separating the pion propagators
into their positive and negative energy parts,

$\displaystyle\frac{1}{q^{2}-m_{\pi}^{2}+i\epsilon}=\frac{1}{2\omega(\vec{q}\,)}\left(\frac{1}{q^{0}-\omega(\vec{q}\,)+i\epsilon}-\frac{1}{q^{0}+\omega(\vec{q}\,)-i\epsilon}\right),$
(18)

where $\omega^{2}(\vec{q}\,)=\vec{q}^{\,\,2}+m_{\pi}^{2}$, we can readily
perform the $q^{0}$ integration of the amplitudes $t^{(a1)}$ and $t^{(a_{2})}$
of Eqs. (15) and (16). We obtain

$\displaystyle t$
$\displaystyle=t^{(a1)}+t^{(a2)}=\frac{-2}{\sqrt{3}}\,g\,g_{a_{1},\rho\pi}\,g_{\sigma,\pi\pi}\,\epsilon_{a_{1},j}\int\frac{d^{3}q}{(2\pi)^{3}}\,F(\vec{q},\vec{p}_{\pi})\,(2\vec{p}_{\pi}+\vec{q}\,)_{j}$
(19)

where

$\displaystyle F(\vec{q},\vec{p}_{\pi})$
$\displaystyle=\frac{1}{2\omega(\vec{q}\,)}\,\frac{1}{2\omega_{\rho}}\,\frac{1}{2\omega(\vec{q}+\vec{p}_{\pi})}\,\bigg{[}\frac{1}{P^{0}-\omega(\vec{q}\,)-\omega_{\rho}(\vec{q}\,)+\frac{i\Gamma_{\rho}}{2}}$
$\displaystyle\times\left(\frac{1}{P^{0}-p_{\pi}^{0}-\omega(\vec{q}\,)-\omega(\vec{q}+\vec{p}_{\pi})+i\epsilon}+\frac{1}{p_{\pi}^{0}-\omega_{\rho}(\vec{q}\,)-\omega(\vec{q}+\vec{p}_{\pi})+\frac{i\Gamma_{\rho}}{2}}\right)$
$\displaystyle+\frac{1}{p_{\pi}^{0}-\omega(\vec{q}+\vec{p}_{\pi})-\omega_{\rho}(\vec{q}\,)+\frac{i\Gamma_{\rho}}{2}}\,\frac{1}{p_{\pi}^{0}-P^{0}-\omega(\vec{q}\,)-\omega(\vec{q}+\vec{p}_{\pi})+i\epsilon}\bigg{]},$
(20)

which, taking into account that

$\displaystyle\int\frac{d^{3}q}{(2\pi)^{3}}\,F(\vec{q},\vec{p}_{\pi})\,q_{i}=p_{\pi,i}\int\frac{d^{3}q}{(2\pi)^{3}}F(\vec{q},\vec{p}_{\pi})\frac{\vec{p}_{\pi}\cdot\vec{q}}{\vec{p}_{\pi}^{\,\,2}}\
,$ (21)

can be cast as

$\displaystyle t$
$\displaystyle=\frac{-2}{\sqrt{3}}\,g\,g_{a_{1},\rho\pi}\,g_{\sigma,\pi\pi}\,\vec{\epsilon}_{a_{1}}\cdot\vec{p}_{\pi}\int\frac{d^{3}q}{(2\pi)^{3}}\,F(\vec{q},\vec{p}_{\pi})\left(2+\frac{\vec{p}_{\pi}\cdot\vec{q}}{\vec{p}_{\pi}^{\,\,2}}\,\right)\
.$ (22)

In Ref. ZhangXie  the integral of Eq. (22) is regularized in terms of a cut
off in $\vec{q}$ of 630 MeV and an additional form factor. We regularize it
with the intrinsic cut offs that stem from the chiral unitary approach.
Indeed, as seen in Ref. DaniJuan , the scattering amplitudes appear with the
factorization

$\displaystyle
t(\vec{q},\vec{q}\,^{\prime})=t\,\theta(q_{\text{max}}-|\vec{q}\,|)\,\theta(q^{\prime}_{\text{max}}-|\vec{q}\,^{\prime}\,|)\
,$ (23)

implying that the vertices are factorized as

$\displaystyle t_{a_{1},\rho\pi}$ $\displaystyle\equiv
g_{a_{1},\rho\pi}\,\theta(q_{\text{max}}-|\vec{q}\,|)\,\vec{\epsilon}_{a_{1}}\cdot\vec{\epsilon}_{\rho}$
$\displaystyle t_{\sigma,\pi\pi}$ $\displaystyle\equiv
g_{\sigma,\pi\pi}\,\theta(q^{\prime}_{\text{max}}-|\vec{q}\,^{\prime}|)\ ,$

and $\vec{q}\,^{\prime}$ is the boosted value of the momentum $\vec{q}$ to the
frame where the $\sigma$ is at rest (see formula in Appendix B). Thus,
finally,

$\displaystyle
t=-\frac{2}{\sqrt{3}}\,g\,g_{a_{1},\rho\pi}\,g_{\sigma\pi\pi}\,\vec{\epsilon}_{a_{1}}\cdot\vec{p}_{\pi}\,\tilde{t}$
(24)

where

$\displaystyle\tilde{t}$
$\displaystyle=\int\frac{d^{3}q}{(2\pi)^{3}}\,F(\vec{q},\vec{p}_{\pi})\left(2+\frac{\vec{p}_{\pi}\cdot\vec{q}}{p_{\pi}^{\,\,2}}\,\right)\,\theta(q_{\text{max}}-|\vec{q}\,|)\,\theta(q^{\prime}_{\text{max}}-|\vec{q}\,^{\prime}|)$
(25)

and, following Ref. Luis , we take $q_{\text{max}}=800$ MeV and
$q^{\prime}_{\text{max}}=750$ MeV (as discussed in Sec. III), by means of
which we obtain the width for $a_{1}\to\pi\sigma$, averaging over the $a_{1}$
initial polarizations,

$\displaystyle\Gamma(M_{I})$
$\displaystyle=\frac{1}{8\pi}\frac{1}{M_{a_{1}}^{2}}\frac{1}{3}\left|\frac{2}{\sqrt{3}}\,g\,g_{a_{1},\rho\pi}\,g_{\sigma,\pi\pi}\right|^{2}\,p_{\pi}^{3}\,|\tilde{t}|^{2}$
(26)

where $M_{I}$ is the invariant mass of the $\sigma$ and
$p_{\pi}\equiv|\vec{p}_{\pi}|=\lambda^{\nicefrac{{1}}{{2}}}\left(M_{a_{1}}^{2},m_{\pi}^{2},M_{I}^{2}\right)/\left(2M_{a_{1}}\right)$
here and in the following. Eq. (26) assumes a $\sigma$ with fixed mass
$M_{I}$, but the $\sigma$ has, in fact, a large mass distribution. We can
formally take this into account (it will done more rigorously below) assuming
a Breit-Wigner mass distribution

$\displaystyle
S_{\sigma}(M_{I})=-\frac{1}{\pi}\,\text{Im}\,\frac{1}{M_{I}^{2}-m_{\sigma}^{2}+im_{\sigma}\Gamma_{\sigma}}$
(27)

and, thus,

$\displaystyle\Gamma=\frac{-1}{18\pi^{2}M_{a_{1}}^{2}}\,g^{2}\,|g_{a_{1},\rho\pi}|^{2}\int
dM_{I}^{2}\,\text{Im}\frac{|g_{\sigma,\pi\pi}|^{2}}{M_{I}^{2}-m_{\sigma}^{2}+im_{\sigma}\Gamma_{\sigma}}|\tilde{t}(M_{I})|^{2}\,p_{\pi}^{3}(M_{I})\
.$ (28)

But now we see that formally
$|g_{\sigma,\pi\pi}|^{2}(M_{I}^{2}-m_{\sigma}^{2}+i\,m_{\sigma}\Gamma_{\sigma})^{-1}$
is the Breit-Wigner representation of the $\pi\pi$ amplitude. Hence, we
replace,

$\displaystyle\frac{|g_{\sigma,\pi\pi}|^{2}}{M_{I}^{2}-m_{\sigma}^{2}+im_{\sigma}\Gamma_{\sigma}}\to
t_{\pi\pi,\pi\pi}^{I=0}(M_{I})\ ,$ (29)

which allows us to use the realistic scattering amplitude
$t_{\pi\pi,\pi\pi}^{I=0}$ from Ref. npa , or from Ref. Liang  where the
coupled channels $\pi^{+}\pi^{-}$, $\pi^{0}\pi^{0}$, $K^{+}K^{-}$,
$K^{0}\bar{K}^{0}$ and $\eta\eta$ are used. Then, as found in Appendix A,

$\displaystyle
t_{\pi\pi,\pi\pi}^{I=0}=3t_{\pi^{+}\pi^{-},\pi^{+}\pi^{-}}^{I=0}\ .$ (30)

Hence, we can write

$\displaystyle\Gamma(M_{a_{1}})=\frac{-1}{18\pi^{2}M_{a_{1}}^{2}}\,g^{2}\,|g_{a_{1},\rho\pi}|^{2}\int
dM_{I}^{2}\,\text{Im}\,[t_{\pi\pi,\pi\pi}^{I=0}(M_{I})]|\tilde{t}(M_{I})|^{2}\,p_{\pi}^{3}(M_{I})\
,$ (31)

where the dependence on $M_{a_{1}}$ is explicitly highlighted. Since the
$a_{1}(1260)$ also has a large width we can do a convolution of the width of
Eq. (31) with the mass distribution of the $a_{1}(1260)$ in analogy to Eq.
(28). We use the most accurate results for the $a_{1}(1260)$ mass and width
from the COMPASS collaboration COMPASSa1 ,

$\displaystyle M_{a_{1}}=(1255\pm 6^{+7}_{-17})~{}\mathrm{MeV}\ ,$
$\displaystyle\Gamma_{a_{1}}=(367\pm 9^{+28}_{-25})~{}\mathrm{MeV}\ ,$ (32)

and obtain

$\displaystyle\Gamma=\frac{1}{N}\int\limits_{(M_{a_{1}}-2\Gamma_{a_{1}})^{2}}^{(M_{a_{1}}+2\Gamma_{a_{1}})^{2}}d\tilde{M}_{a_{1}}^{2}\left(-\frac{1}{\pi}\right)\,\text{Im}\,\frac{1}{\tilde{M}_{a_{1}}^{2}-M_{a_{1}}^{2}+iM_{a_{1}}\Gamma_{a_{1}}}\,\Gamma(\tilde{M}_{a_{1}})$
(33)

with

$\displaystyle
N=\int\limits_{(M_{a_{1}}-2\Gamma_{a_{1}})^{2}}^{(M_{a_{1}}+2\Gamma_{a_{1}})^{2}}d\tilde{M}_{a_{1}}^{2}\left(-\frac{1}{\pi}\right)\,\text{Im}\,\frac{1}{\tilde{M}_{a_{1}}^{2}-M_{a_{1}}^{2}+iM_{a_{1}}\Gamma_{a_{1}}}\
.$ (34)

##  III Results

To calculate the $a_{1}\to\pi\sigma$ width we use the model of Ref. npa  with
a single channel $\ket{\pi\pi,I=0}$ amplitude using unitary normalization
(extra $1/\sqrt{2}$ in Eq. (6)),

$\displaystyle t^{I=0,u}_{\pi\pi,\pi\pi}=\frac{V^{I=0,u}}{1-V^{I=0,u}G}\ ,$
(35)

with

$\displaystyle V^{I=0,u}=-\frac{1}{f^{2}}(M^{2}_{I}-\frac{m^{2}_{\pi}}{2})\ ,$
(36)

and $G$ the $\pi\pi$ loop function regularized with a cut off
$q’_{\mathrm{max}}=750$ MeV, such as to get the $I=0$, $\pi\pi$ phase shifts.
The coupling that we obtain is given in Eq. (8). Thus,

$\displaystyle|g_{\sigma,\pi\pi}^{(u)}|=2.99~{}\mathrm{GeV};~{}|g_{\sigma,\pi\pi}|=\sqrt{2}\,g^{(u)}_{\sigma,\pi\pi}\
,$ (37)

and we have,

$\displaystyle t^{I=0}_{\pi\pi,\pi\pi}=2\,t^{I=0,(u)}_{\pi\pi,\pi\pi}\ .$ (38)

We have also carried out the calculations using the model of Ref. Liang ,
considering explicitly the $\pi^{+}\pi^{-},\pi^{0}\pi^{0}$,
$K^{+}K^{-},K^{0}\bar{K}^{0}$, and $\eta\eta$ channels, and using Eq. (30),
and the results are practically the same except that including the extra
channels demands the use of a different cut off, $q_{\mathrm{max}}=600$ MeV.
The results that we show are those based on Eqs. (35)-(38). We shall estimate
uncertainties at the end using results of different models for $\pi\pi$
scattering.

In Fig. 3 we show first the results obtained with Eq. (26) as a function of
$M_{I}$.

Figure 3: $\Gamma(a_{1}^{+}\to\pi^{+}\sigma)$ from Eq. (26) as a function of
$M_{I}$ for fixed $M_{a_{1}}=1255$ MeV.

As expected for lower $\sigma$ masses $\Gamma$ grows since $p_{\pi}$ is
bigger, and we find a cusp like peak at the two-pion threshold as it should
be. We already see from Fig. 3 that around $M_{I}\simeq 500$ MeV the width is
of the order of $33$ MeV.

Next we use Eq. (31) in terms of $t^{\mathrm{I=0}}_{\pi\pi,\pi\pi}$ of Eqs.
(35)-(38). We show the results in terms of $M_{a_{1}}$ which we plot in Fig.
4.

Figure 4: $\Gamma(a_{1}^{+}\to\pi^{+}\sigma)$ from Eq. (31) as a function of
$M_{a_{1}}$.

We see that $\Gamma$ has been reduced, since around $M_{a_{1}}=1255$ MeV it
has a value about $14$ MeV, around one half what we estimated from Fig. 3. We
should explain this difference, which we do in Fig. 5.

Figure 5: Comparison of the imaginary parts of
$g^{2}_{\sigma}/(M^{2}_{I}-m^{2}_{\sigma}+im_{\sigma}\Gamma_{\sigma})$ and
$t_{\pi\pi,\pi\pi}^{I=0}$.

There we plot the imaginary part of the Breit-Wigner form of the $\sigma$,

$\displaystyle\mathrm{Im}\frac{|g_{\sigma,\pi\pi}|^{2}}{M^{2}_{I}-m^{2}_{\sigma}+im_{\sigma}\Gamma_{\sigma}}\
,$ (39)

versus that of the realistic $t^{I=0}_{\pi\pi,\pi\pi}$ amplitude,
$\mathrm{Im}\,t^{I=0}_{\pi\pi,\pi\pi}$. What we observe is that these two
magnitudes are very different, both in strength and shape. In particular, the
imaginary part of the full $\pi\pi$ amplitude goes to zero at threshold due to
unitarity, as it should be, while the Breit-Wigner approximation does not. At
higher $M_{I}$ the situation is opposite. This discrepancy should be stressed
since in conventional experimental partial-wave analyses it is customary to
use Breit-Wigner forms for the $\sigma$, although in some cases more
sophisticated forms are taken, like in Ref. newbes  based upon the formulation
of Ref. zoubugg . Altogether we find about a factor of two difference in the
evaluation of the width by both methods, the realistic one being the one
obtained in terms of $t_{\pi\pi,\pi\pi}^{I=0}$, of course.

Finally, if we conduct the convolution of the result in Fig. 4 with the
$a_{1}(1260)$ mass distribution according to Eq. (33) we find

$\Gamma=13.4\,\textrm{ MeV}\,.$ (40)

Compared with the $a_{1}(1260)$ width of Eq. (32), it gives a branching ratio
of

${\cal B}(a_{1}(1260)\to\pi\sigma)\simeq(3.3-4.0)\%$ (41)

considering the errors of Eq. (32).

###  III.1 Uncertainty analysis

The branching ratio of Eq. (41) reflects only the uncertainties in the width
of Eq. (33). To estimate uncertainties from the use of other models we resort
to a simple method by looking at the couplings obtained by different groups
which are tabulated in Table 5 of Ref. sigma . We simply scale the result of
Eq. (40) by the squared of the ratio of the couplings shown in that table. We
show the results in Table 1 and find a range of results from about $13$ MeV to
$22$ MeV for the width of $a_{1}(1260)$ to $\sigma\pi$.

Ref. | $g_{\sigma\pi\pi}(\mathrm{GeV})$ | $\Gamma(\mathrm{MeV})$
---|---|---
Caprini  | $3.58\pm 0.03$ | $19.3$
Pelaezrios  | $3.5$ | $18.3$
Masjuan ; Caprini2  | $3.8\pm 0.4$ | $21.6$
This work / npa  | $2.99$ | $13.4$
Table 1: Estimation of the uncertainties in $\Gamma$ from the coupling
$g_{\sigma,\pi\pi}$ with some $g_{\sigma,\pi\pi}$ taken from Ref. sigma .

The branching ratios obtained, $(3.3-6.6)$%, are in the range of some of the
small values quoted in Eqs. (5), but much smaller than the rate quoted in Ref.
Asner , Eq. (1). We should certainly note the large dispersion of experimental
results in Eqs. (1), (5) that might come from treating the $\sigma$ as an
ordinary resonance as discussed before. The formulation presented here should
help perform more accurate analyses of data in future analyses. Actually, with
this work finished, there has been a recent experiment newbes  for the
$D^{+}_{s}\to\pi^{+}\pi^{-}\pi^{+}\eta$ decay with emphasis on different
partial wave contributions, and from where one can induce the branching ratio
that we have evaluated. From Table I of that work one finds a weight for
$a_{1}(1260)^{+}\eta$, with $a_{1}^{+}$ observed in $\rho^{0}\pi^{+}$,
$W_{1}=55.4\pm 3.9\pm 2.0$. Then $a_{1}^{+}\eta$, with $a_{1}^{+}$ observed in
$f_{0}(500)\pi^{+}$, has a weight of $W_{2}=8.1\pm 1.9\pm 2.1$. But this only
accounts for $f_{0}(500)\to\pi^{+}\pi^{-}$. The $\pi^{0}\pi^{0}$ decay has a
strength of $\nicefrac{{1}}{{2}}$ of that of $\pi^{+}\pi^{-}$. Furthermore,
since the weight for $a_{1}(1260)$ decay into $\rho^{+}\pi^{0}$ is the same as
for $\rho^{0}\pi^{+}$, and the total branching ratio of $a_{1}\to\rho\pi$ is
$60\%$ PDG , we have

$\dfrac{\Gamma(a_{1}\to\sigma\pi)}{\Gamma_{a_{1}}}=\dfrac{3W_{2}/2}{2W_{1}/0.6}=(6.6\pm
2.4)\%,$ (42)

where we have summed errors in quadrature. This number is in good agreement
with our results. It is interesting to see that with our $a_{1}\to\rho\pi$
coupling of Eq. (10), we obtain a branching fraction of $a_{1}$ to $\rho\pi$
of 51%, close to the one reported in the PDG PDG . If, instead of the rather
uncertain experimental $a_{1}\to\rho\pi$ 60% branching ratio, we use our
theoretical value of 51%, the numbers in Eq. (42) change to $(5.6\pm 2)\%$.
Hence, a range of $\Gamma(a_{1}\to\sigma\pi)/\Gamma_{a_{1}}\approx(3.6-9)\%$
is a safe estimate from the experiment of Ref. newbes .

##  IV Conclusions

We have evaluated the decay width of the $a_{1}(1260)$ resonance into
$\pi\sigma(f_{0}(500))$ taking into account that the $a_{1}(1260)$ is a state
dynamically generated by the vector-pseudoscalar interaction, concretely with
the channels $\rho\pi$ and $K^{*}\bar{K}$, and the $\sigma$ arises from the
unitarized interaction of pseudoscalar-pseudoscalar mesons, mostly of two
pions. Given this starting point, the mechanism that leads to the
$a_{1}\to\pi\sigma$ decay is given by a triangle mechanism in which the
$a_{1}(1260)$ decays to $\rho\pi$, the $\rho$ decays to $\pi\pi$ and one of
these pions, together with the one coming from the $a_{1}\to\rho\pi$ decay,
fuse into the $\sigma$ meson. The mechanism allows to have the intermediate
particles on shell, but does not develop a triangle singularity, yet it
provides a relatively large strength for that decay. We find decay widths in
this channel of the order of $13-22\,\textrm{ MeV}$, which amount to a
branching fraction of $(3.3-6.6)$%. This is a small fraction and several
experiments, so far largely contradictory, have been devoted to estimate that
width. On the other hand, the recent BESIII experiment newbes  allows to
extract this branching ratio to be in the range $(3.6-9)\%$, which overlaps
with the range that we obtain.

One of the findings is that the width obtained depends much on how the mass
distribution of the $\sigma$ is considered. We found that treating the
$\sigma$ as an ordinary Breit-Wigner resonance, using realistic coupling, mass
and width obtained in theoretical works, gives rise to results which differ
substantially from those obtained by using the imaginary part of the $I=0$,
$\pi\pi$ scattering amplitude as the mass distribution. This should serve as a
warning for future experimental analyses of this partial decay width. The work
of Ref. newbes  uses a more elaborate and realistic picture for the $\sigma$
resonance extracted from Ref. zoubugg .

Any information on the $\pi\sigma$ and other subdominant branching ratios of
the $a_{1}(1260)$ helps assess systematic effects in recent efforts to
calculate the properties of three-body resonances in lattice QCD Mai:2021nul .
In particular, the smallness of that branching ratio found in this study
strengthens the truncation to $\pi\rho$ channels made in the finite-volume
frameworks of Refs. RocaFinite ; Mai:2021nul .

The results of this study are tied to the assumptions about the nature of both
the $a_{1}(1260)$ and the $\sigma$ resonances as being dynamically generated
from the pseudoscalar-vector and pseudoscalar-pseudoscalar interaction. Given
the ongoing discussion on the nature of the axial vector and scalar meson
resonances, the precise determination of this width should provide valid
information concerning this relevant issue, and it is our hope that the
present work stimulates further efforts in this direction.

###### Acknowledgements.

R. M. acknowledges support from the CIDEGENT program with Ref.
CIDEGENT/2019/015 and from the spanish national grants PID2019-106080GB-C21
and PID2020-112777GB-I00. This work is partly supported by the National
Natural Science Foundation of China under Grants No. 11975083 and No.
11947413. This work is also partly supported by the Spanish Ministerio de
Economia y Competitividad and European FEDER funds under Contracts No.
FIS2017-84038-C2-1-P B and by Generalitat Valenciana under contract
PROMETEO/2020/023. This project has received funding from the European Unions
Horizon 2020 research and innovation programme under grant agreement No.
824093 for the “STRONG-2020” project. This material is based upon work
supported by the National Science Foundation under Grant No. PHY-2012289 and
the U.S. Department of Energy under Award Number DE-SC0016582 and DE-
AC05-06OR23177. MD is grateful for the hospitality of the University of
Valencia, where part of this work was done.

##  Appendix A Coupling of $\sigma$ to $\pi\pi$

Figure 6: Amplitudes $\pi\pi\to\pi\pi$ through an isoscalar source.

We consider here the couplings of $\sigma$ to $\pi^{+}\pi^{-},\pi^{0}\pi^{0}$
taking into account the identity of the $\pi^{0}$. An isoscalar coupling to
$\pi\pi$ is given by a Lagrangian of the type,

$\displaystyle{\cal
L}=\alpha\,\vec{\phi}\cdot\vec{\phi}=\alpha\,(-1)^{\mu}\phi_{\mu}\phi_{-\mu},\quad\mu=-1,0,+1$
(43)

with

$\displaystyle\phi_{-1}=\frac{1}{\sqrt{2}}(\phi_{1}+i\phi_{2}),$
$\displaystyle\phi_{0}=\phi_{3},$
$\displaystyle\phi_{+1}=-\frac{1}{\sqrt{2}}(\phi_{1}-i\phi_{2}),$ (44)

where $\phi_{-}$ destroys a $\pi^{-}$, creates a $\pi^{+}$, while $\phi_{+}$
destroys a $\pi^{+}$ and creates a $\pi^{-}$.

Using the Lagrangians of Eq. (43) and the two ways to destroy or create
$\pi^{+}\pi^{-}$, or $\pi^{0}\pi^{0}$, we obtain

$\displaystyle
t^{I=0}_{\pi^{+}\pi^{-},\pi^{+}\pi^{-}}=t^{I=0}_{\pi^{+}\pi^{-},\pi^{0}\pi^{0}}=t^{I=0}_{\pi^{0}\pi^{0},\pi^{0}\pi^{0}}=\frac{(2\alpha)^{2}}{s-m^{2}_{\sigma}},$
(45)

and using the $\pi\pi(I=0)$ states of Eq. (6) we find

$\displaystyle
t^{I=0}_{\pi\pi,\pi\pi}=\frac{1}{3}(4t^{I=0}_{\pi^{+}\pi^{-},\pi^{+}\pi^{-}}+4t^{I=0}_{\pi^{+}\pi^{-},\pi^{0}\pi^{0}}+t^{I=0}_{\pi^{0}\pi^{0},\pi^{0}\pi^{0}})$
$\displaystyle=\frac{9}{3}t^{I=0}_{\pi^{+}\pi^{-},\pi^{0}\pi^{0}}=\frac{9}{3}\frac{(2\alpha)^{2}}{s-m^{2}_{\sigma}}\equiv\frac{g^{2}_{\sigma,\pi\pi}}{s-m^{2}_{\sigma}}.$
(46)

Hence we find:

$\displaystyle g_{\sigma,\pi^{+}\pi^{-}}\equiv
g_{\sigma,\pi^{0}\pi^{0}}=2\alpha=\frac{1}{\sqrt{3}}g_{\sigma,\pi\pi}.$ (47)

##  Appendix B Boosted momenta

The vector $\vec{q}\,^{\prime}$ defined in the rest frame of the $\sigma$ is
given by

$\displaystyle\vec{q}\,^{\prime}=\left[\left(\frac{E_{\sigma}}{M_{I}}-1\right)\frac{\vec{q}\cdot\vec{p}_{\sigma}}{\vec{p}\,^{2}_{\sigma}}-\frac{q^{0}}{M_{I}}\right]\vec{p}_{\sigma}+\vec{q},$
(48)

where

$\displaystyle\vec{p}_{\sigma}=-\vec{p}_{\pi},$ $\displaystyle
q^{0}=\omega(\vec{q}),$ $\displaystyle
E_{\sigma}=M_{a_{1}}-\omega(\vec{p}_{\pi}),$
$\displaystyle|\vec{p}_{\pi}|=\frac{\lambda^{1/2}(M^{2}_{a_{1}},m^{2}_{\pi},M^{2}_{I})}{2M_{a_{1}}},$
(49)

and $M_{I}$ is the invariant mass carried by the $\sigma$ in the diagrams of
Figs. 1(a), 2, and 6.

## References

  * (1) S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985).
  * (2) J. Vijande, F. Fernandez and A. Valcarce, J. Phys. G 31, 481 (2005).
  * (3) M. F. M. Lutz and E. E. Kolomeitsev, Nucl. Phys. A 730, 392 (2004).
  * (4) L. Roca, E. Oset and J. Singh, Phys. Rev. D 72, 014002 (2005).
  * (5) L. S. Geng, E. Oset, L. Roca and J. A. Oller, Phys. Rev. D 75, 014017 (2007).
  * (6) Y. Zhou, X. L. Ren, H. X. Chen and L. S. Geng, Phys. Rev. D 90, no. 1, 014020 (2014).
  * (7) M. C. Birse, Z. Phys. A 355, 231 (1996).
  * (8) H. Nagahiro, L. Roca, A. Hosaka and E. Oset, Phys. Rev. D 79, 014015 (2009).
  * (9) J. M. Dias, G. Toledo, L. Roca and E. Oset, Phys. Rev. D 103, 116019 (2021).
  * (10) M. Bando, T. Kugo and K. Yamawaki, Phys. Rept. 164, 217 (1988).
  * (11) M. Harada and K. Yamawaki, Phys. Rept. 381, 1 (2003).
  * (12) U. G. Meissner, Phys. Rept. 161, 213 (1988).
  * (13) G. Ecker, J. Gasser, H. Leutwyler, A. Pich and E. de Rafael, Phys. Lett. B 223, 425 (1989).
  * (14) S. X. Nakamura, Phys. Rev. D 93, no. 1, 014005 (2016).
  * (15) S. Prelovsek, C. B. Lang, D. Mohler and M. Vidmar, PoS LATTICE2011, 137 (2011).
  * (16) L. Roca and E. Oset, Phys. Rev. D 85, 054507 (2012).
  * (17) C. B. Lang, L. Leskovec, D. Mohler and S. Prelovsek, JHEP 1404, 162 (2014).
  * (18) A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards and D. J. Wilson, Phys. Rev. D 100, no. 5, 054506 (2019).
  * (19) M. Mai, A. Alexandru, R. Brett, C. Culver, M. Döring, F. X. Lee and D. Sadasivan, [arXiv:2107.03973 [hep-lat]].
  * (20) M. Mai, B. Hu, M. Doring, A. Pilloni and A. Szczepaniak, Eur. Phys. J. A 53, 177 (2017) [arXiv:1706.06118 [nucl-th]].
  * (21) M. Mai and M. Döring, Eur. Phys. J. A 53, no.12, 240 (2017) [arXiv:1709.08222 [hep-lat]].
  * (22) F. Aceti, J. M. Dias and E. Oset, Eur. Phys. J. A 51, 48 (2015).
  * (23) F. Aceti, J. J. Xie and E. Oset, Phys. Lett. B 750, 609 (2015).
  * (24) J. J. Wu, X. H. Liu, Q. Zhao and B. S. Zou, Phys. Rev. Lett. 108, 081803 (2012).
  * (25) F. Aceti, W. H. Liang, E. Oset, J. J. Wu and B. S. Zou, Phys. Rev. D 86, 114007 (2012).
  * (26) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. Lett. 108, 182001 (2012).
  * (27) N. N. Achasov, A. A. Kozhevnikov and G. N. Shestakov, Phys. Rev. D 92, 036003 (2015).
  * (28) M. C. Du and Q. Zhao, Phys. Rev. D 100, no. 3, 036005 (2019).
  * (29) Y. Huang, C. j. Xiao, Q. F. L., R. Wang, J. He and L. Geng, Phys. Rev. D 97, no. 9, 094013 (2018).
  * (30) Y. Huang, M. Z. Liu, J. X. Lu, J. J. Xie and L. S. Geng, Phys. Rev. D 98, no. 7, 076012 (2018).
  * (31) J. J. Xie, G. Li and X. H. Liu, Chin. Phys. C 44, no. 11, 114104 (2020).
  * (32) Z. M. Ding, H. Y. Jiang and J. He, Eur. Phys. J. C 80, no. 12, 1179 (2020).
  * (33) Y. Huang, J. X. Lu, J. J. Xie and L. S. Geng, Eur. Phys. J. C 80, no. 10, 973 (2020).
  * (34) Y. K. Hsiao, Y. Yu and B. C. Ke, Eur. Phys. J. C 80, no. 9, 895 (2020).
  * (35) X. Z. Ling, M. Z. Liu, J. X. Lu, L. S. Geng and J. J. Xie, arXiv:2102.05349 [hep-ph].
  * (36) R. Karplus, C. M. Sommerfield and E. H. Wichmann, Phys. Rev. 111, 1187 (1958).
  * (37) L. D. Landau, Nucl. Phys. 13, no. 1, 181 (1959) [Sov. Phys. JETP 10, no. 1, 45 (1960)] [Zh. Eksp. Teor. Fiz. 37, no. 1, 62 (1959)].
  * (38) S. Coleman and R. E. Norton, Nuovo Cim. 38, 438 (1965).
  * (39) M. Bayar, F. Aceti, F. K. Guo and E. Oset, Phys. Rev. D 94, no. 7, 074039 (2016).
  * (40) F. K. Guo, X. H. Liu and S. Sakai, Prog. Part. Nucl. Phys. 112, 103757 (2020).
  * (41) C. Adolph et al. [COMPASS Collaboration], Phys. Rev. Lett. 115, no. 8, 082001 (2015).
  * (42) X. H. Liu, M. Oka and Q. Zhao, Phys. Lett. B 753, 297 (2016).
  * (43) M. Mikhasenko, B. Ketzer and A. Sarantsev, Phys. Rev. D 91, no. 9, 094015 (2015).
  * (44) F. Aceti, L. R. Dai and E. Oset, Phys. Rev. D 94, no. 9, 096015 (2016).
  * (45) M. G. Alexeev et al. [COMPASS Collaboration], arXiv:2006.05342 [hep-ph].
  * (46) N. Ikeno, R. Molina and E. Oset, arXiv:2103.01712 [nucl-th].
  * (47) R. Molina, N. Ikeno and E. Oset, arXiv:2102.05575 [nucl-th].
  * (48) M. Bashkanov et al., Phys. Rev. Lett. 102, 052301 (2009).
  * (49) D. Sadasivan, M. Mai, H. Akdag and M. Döring, Phys. Rev. D 101, no. 9, 094018 (2020); Erratum: [Phys. Rev. D 103, no. 1, 019901 (2021)].
  * (50) X. Zhang and J. J. Xie, Commun. Theor. Phys. 70, no. 1, 060 (2018).
  * (51) L. R. Dai, L. Roca and E. Oset, Eur. Phys. J. C 80, no.7, 673 (2020).
  * (52) P. A. Zyla et al. [Particle Data Group], Prog. Theor. Exp. Phys. 2020, no. 8, 083C01 (2020) and 2021 update.
  * (53) D. M. Asner et al. [CLEO Collaboration], Phys. Rev. D 61, 012002 (2000).
  * (54) P. Salvini et al. [OBELIX Collaboration], Eur. Phys. J. C 35, 21 (2004).
  * (55) R. R. Akhmetshin et al. [CMD-2 Collaboration], Phys. Lett. B 466, 392 (1999).
  * (56) R. S. Longacre, Phys. Rev. D 26, 82 (1982).
  * (57) S. U. Chung et al., Phys. Rev. D 65, 072001 (2002).
  * (58) G. Colangelo, J. Gasser and H. Leutwyler, Nucl. Phys. B 603, 125 (2001).
  * (59) J. R. Pelaez, Phys. Rept. 658, 1 (2016).
  * (60) H. Kamano, S. X. Nakamura, T. S. H. Lee and T. Sato, Phys. Rev. D 84, 114019 (2011) [arXiv:1106.4523 [hep-ph]].
  * (61) G. Janssen, J. W. Durso, K. Holinde, B. C. Pearce and J. Speth, Phys. Rev. Lett. 71, 1978-1981 (1993).
  * (62) M. Wagner and S. Leupold, Phys. Lett. B 670, 22-26 (2008) [arXiv:0708.2223 [hep-ph]].
  * (63) M. Wagner and S. Leupold, Phys. Rev. D 78, 053001 (2008) [arXiv:0801.0814 [hep-ph]].
  * (64) S. Schael et al. [ALEPH], Phys. Rept. 421, 191-284 (2005) [arXiv:hep-ex/0506072 [hep-ex]].
  * (65) J. A. Oller and E. Oset, Nucl. Phys. A 620, 438 (1997); Erratum: [Nucl. Phys. A 652, 407 (1999)].
  * (66) H. Albrecht et al. [ARGUS Collaboration], Z. Phys. C 58, 61 (1993).
  * (67) J. A. Oller and E. Oset, Phys. Rev. D 60, 074023 (1999).
  * (68) S. Sakai, E. Oset and A. Ramos, Eur. Phys. J. A 54, 10 (2018).
  * (69) D. Gamermann, J. Nieves, E. Oset and E. Ruiz Arriola, Phys. Rev. D 81, 014029 (2010).
  * (70) W. H. Liang and E. Oset, Phys. Lett. B 737, 70-74 (2014).
  * (71) M. Alekseev et al. [COMPASS], Phys. Rev. Lett. 104, 241803 (2010).
  * (72) M. Ablikim et al. [BESIII], [arXiv:2106.13536 [hep-ex]].
  * (73) D. V. Bugg, B. S. Zou and A. V. Sarantsev, Nucl. Phys. B 471, 59-89 (1996).
  * (74) I. Caprini, Phys. Rev. D 77, 114019 (2008).
  * (75) J. R. Pelaez and G. Rios, Phys. Rev. D 82, 114002 (2010).
  * (76) P. Masjuan, J. Ruiz de Elvira and J. J. Sanz-Cillero, Phys. Rev. D 90, no.9, 097901 (2014).
  * (77) I. Caprini, P. Masjuan, J. Ruiz de Elvira and J. J. Sanz-Cillero, Phys. Rev. D 93, no.7, 076004 (2016).


# Iterative solution of relativistic Boltzmann equation in curved spacetime
with application to kinetic coefficients

Long Cui1, Xin Hao2, and Liu Zhao1
1School of Physics, Nankai University, Tianjin 300071, China
2School of Physics, Hebei Normal University, Shijiazhuang 050024, China
email<EMAIL_ADDRESS>
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

Under relaxation time approximation, we obtain an iterative solution to the
relativistic Boltzmann equation in generic stationary spacetime. This solution
provides a scheme to study non-equilibrium system order by order. As a
specific example, we analytically calculated the covariant expressions of the
particle flow and the energy momentum tensor up to the first order in
relaxation time. Finally and most importantly, we present all 14 kinetic
coefficients for a neutral system, which are verified to satisfy the Onsager
reciprocal relation and guarantee a non-negative entropy production.

##  1 Introduction

Relativistic kinetic theory set up a powerful framework to study non-
equilibrium phenomena in curved spacetime. In general, at long wavelength and
low frequency, non-equilibrium systems can be excellently described by
hydrodynamics which, from a modern perspective, is an effective theory
consists of a gradient expansion about the local equilibrium state. At linear
level, the property of the fluid system is totally determined by a set of
phenomenological transport coefficients which can be measured by experiment
and can be determined by calculations in the underlying microscopic theory.
Although there are multiple approaches with different domain of validity to
calculate the transport coefficients, in a curved spacetime background the
relativistic kinetic theory turns out to be the most efficient and convenient
option with many important applications ranging from stability theory [1, 2,
3], astrophysics[4, 5] to cosmology [6, 7, 8].

Historically, kinetic theory is developed to study thermodynamic behaviors of
classical gaseous systems. However, from quantum field theoretical point of
view, many thermodynamic properties of weakly coupled systems can be also
obtained using kinetic theory as an effective theory. One of the most well-
known example is the effective kinetic theory for quark-gluon plasma where at
high temperature or large density the gauge coupling constant becomes
sufficiently small to allow for perturbative calculations [9, 10, 11].
Although LHC experiment is the main motivation for interests in quark-gluon
plasma, such a state also occurs in the early universe and the core of some
neutron stars where the effect of spacetime curvature is un-negligible.
Unfortunately, non-equilibrium thermodynamics in curved spacetime, even for
classical system, remains to some extent an open question. In this direction,
relativistic kinetic theory in curved spacetime is expected to be still alive
and well applicable.

The earliest statistical description of the equilibrium of relativistic gases
began with Jüttner [12], who extended the Maxwell-Boltzmann distribution of
equilibrium gases to the relativistic case. Based on the curved spacetime,
Tolman and Ehrenfest [13, 14] were the first to note that relativistic gases
can only be in equilibrium when the effects of the temperature gradient and
the gravitational field cancel out. Tauber and Weinberg [15] had earlier
formulated the kinetic theory under general relativity. Israel [16]
established the non-equilibrium distribution function and transport
coefficient of relativistic gas by using the Chapman-Enskog method. Lindquist
[17] calculated the transport equation of a gaseous system composed of zero
rest mass particles in spherically symmetric spacetime. Kremer [18, 19]
calculated the bulk viscosity and shear viscosity coefficients of gases in
Schwarzschild spacetime under the post-Newton approximation, obtained
Fourier’s laws for a single gas, and Fick’s law for gas mixtures. However,
there is a scarce literature on the study of relativistic gases that deviate
from equilibrium in general curved spacetime in the framework of relativistic
kinetic theory.

Our goal is to calculate the kinetic coefficients in a covariant formalism in
a generic stationary spacetime. In the previous works [20, 21], the solution
of the relativistic Boltzmann equation was constructed through the gravito-
electromagnetism analogy which facilitate the discussion of particle and
energy transport, but not the viscosity phenomena. In the present work, we
construct the solution of the relativistic Boltzmann equation by use of an
iterative procedure in terms of the relaxation time, which allows us to
calculate all the kinetic coefficients to any order in relaxation time. Each
of these kinetic coefficients can be expressed as a function of temperature
and chemical potential. Compared with the previous literature, our formalism
is simpler and more intuitive, and also applies to the case of degenerate
gases. For simplicity, this paper includes only kinetic coefficients up to the
first order in relaxation time, and focuses on the gaseous system composed of
massive neutral particles, which act as a probe system in the background
spacetime. In this framework, the computations are fully analytical and
clearly covariant.

It is worth emphasizing that calculations are not carried out in a specific
hydrodynamic frame [22], since the kinetic coefficients reflects the response
of the “flow” to the thermodynamic force, which should not be related to the
choice of hydrodynamic frame. By imposing appropriate constraints, our results
can be reduced to the well-known Eckart frame or Landau frame.

The paper is structured as follows. In section 2, we review the relativistic
Boltzmann equation and the detailed balance distribution. Section 3 discusses
the deviation from the detailed balance and iteratively solves the Boltzmann
equation under the relaxation time approximation. In section 4, we use the
first order iterative solution to calculate all the transport coefficients up
to the first order in relaxation time. Section 5 verifies that the above
result satisfies the Onsager reciprocal relation and lead to a non-negative
entropy production.

We adopt the metric signature $(-,+,\cdots,+)$, where the dimension of space
is $d$. The Greek letters $\mu,\nu,\cdots=0,1,\cdots,d$ refer to spacetime
indices and the Latin letters with hat
$\hat{a},\hat{b}=\hat{0},\hat{1},\cdots,\hat{d};\enspace\hat{i},\hat{j},\cdots=\hat{1},\hat{2},\cdots,\hat{d}$
refer to basis indices. In addition, quantities with an overbar (such as
$\bar{f},\bar{T},\bar{\mu},\cdots$) indicate the values in detailed balance.

##  2 Relativistic Boltzmann equation and detailed balance

Let us start with a brief review of the kinetic theory description of
relativistic fluids. The background spacetime $M$ is taken as general as
possible, and the tangent bundle of the spacetime manifold $M$ is denoted as
$TM$. In relativistic kinetic theory, the one particle distribution function
(1PDF) $f$ is defined on the future mass shell bundle [23, 24]

$\Gamma_{m}^{+}:=\\{(x,p)\in TM|p^{\mu}p^{\nu}g_{\mu\nu}=-m^{2}c^{2}\\},\quad
x\in M,p\in T_{x}M,$ (1)

where $p^{\mu}$ is the momentum and $m$ is the mass of particles.

The particle flow $N^{\mu}$, energy momentum tensor $T^{\mu\nu}$ and entropy
flow $S^{\mu}$ can be expressed in terms of the 1PDF $f(x,p)$ as

$N^{\mu}=c\int\boldsymbol{\varpi}p^{\mu}f,\quad
T^{\mu\nu}=c\int\boldsymbol{\varpi}p^{\mu}p^{\nu}f,$ (2)
$S^{\mu}=-k_{B}c\int\boldsymbol{\varpi}p^{\mu}f\left[\log\left(\frac{h^{d}f}{f^{\ast}\mathfrak{g}}\right)-\frac{\mathfrak{g}\log
f^{\ast}}{\varsigma h^{d}f}\right],$ (3)

where $k_{B}$ is Boltzmann constant, $h$ is Planck constant, $\mathfrak{g}$ is
degree of degeneracy, $\varsigma=0,1,-1$ denote the non-degenerate, bosonic
and fermionic cases respectively,
$f^{\ast}=1+\varsigma\mathfrak{g}^{-1}h^{d}f$ and
$\displaystyle\boldsymbol{\varpi}=\frac{\sqrt{g}}{|p_{0}|}(\mathrm{d}p)^{d}$
is the invariant volume element in the momentum space (in which
$g=|\det(g_{\mu\nu})|$).

It is customary to decompose $N^{\mu},T^{\mu\nu}$ using the proper velocity
$U^{\mu}$ ($U^{\mu}U_{\mu}=-c^{2}$) of some prescribed observer $\mathcal{O}$,

$\displaystyle N^{\mu}$ $\displaystyle=$ $\displaystyle n\,U^{\mu}+j^{\mu},$
$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\epsilon\,U^{\mu}U^{\nu}+\frac{1}{c^{2}}q^{\mu}U^{\nu}+\frac{1}{c^{2}}q^{\nu}U^{\mu}+(P+\Pi)\Delta^{\mu\nu}+\Pi^{\mu\nu},$
(4)

where $n$ is the particle number density, $\varepsilon$ is the energy density,
$j^{\mu}$ is the particle flux, $q^{\mu}$ is the energy flux, $P$ is the
hydrostatic pressure, $\Pi$ is the dynamic pressure, $\Pi^{\mu\nu}$ is the
deviatoric stress tensor, and
$\displaystyle\vspace{2pt}\Delta^{\mu\nu}=g^{\mu\nu}+\frac{1}{c^{2}}U^{\mu}U^{\nu}$
is the induced metric on the spacelike hypersurface with normal vector field
$U^{\mu}$. All the above mentioned hydrodynamic variables are measured by the
observer $\mathcal{O}$, and they satisfy the orthogonal relations
$j^{\mu}U_{\mu}=q^{\mu}U_{\mu}=\Pi^{\mu\nu}U_{\mu}=0$.

The evolution of the 1PDF is determined by the relativistic Boltzmann equation

$\mathcal{L}_{\mathcal{H}}f=\mathcal{C}(x,p),$ (5)

where

$\mathcal{L}_{\mathcal{H}}=p^{\mu}\cfrac{\partial}{\partial
x^{\mu}}-\Gamma^{\mu}_{\alpha\beta}p^{\alpha}p^{\beta}\cfrac{\partial}{\partial
p^{\mu}}$ (6)

is the Liouville vector field which is tangent to both $TM$ and
$\Gamma_{m}^{+}$, and

$\mathcal{C}(x,p)=\int\boldsymbol{\varpi}_{2}\boldsymbol{\varpi}_{3}\boldsymbol{\varpi}_{4}[W_{x}(p_{3}+p_{4}\mapsto
p+p_{2})f_{3}f_{4}f^{\ast}f_{2}^{\ast}-W_{x}(p+p_{2}\mapsto
p_{3}+p_{4})ff_{2}f_{3}^{\ast}f_{4}^{\ast}]$ (7)

is the collision integral [25] where $W_{x}$ is referred to as the transition
probability. Further, requiring that the collision process satisfies
microscopic reversibility $W_{x}(p_{3}+p_{4}\mapsto
p+p_{2})=W_{x}(p+p_{2}\mapsto p_{3}+p_{4})$, the collision integral can be
reduced to

$\mathcal{C}(x,p)=\int\boldsymbol{\varpi}_{2}\boldsymbol{\varpi}_{3}\boldsymbol{\varpi}_{4}[W_{x}(p_{3}+p_{4}\mapsto
p+p_{2})(f_{3}f_{4}f^{\ast}f_{2}^{\ast}-ff_{2}f_{3}^{\ast}f_{4}^{\ast})].$ (8)

According to the symmetry of collision integral (8), the divergence-free
nature of $N^{\mu}$ and $T^{\mu\nu}$ is guaranteed and the H-theorem follows
[23],

$\nabla_{\mu}N^{\mu}=0,\quad\nabla_{\mu}T^{\mu\nu}=0,$ (9)
$\nabla_{\mu}S^{\mu}\geqslant 0.$ (10)

When the total entropy of the system reaches the maximum, the distribution
function no longer evolves (collision integral $\mathcal{C}(x,p)=0$), the
state of the fluid is referred to as detailed balance and the corresponding
1PDF reads

$\bar{f}=\frac{\mathfrak{g}}{h^{d}}\frac{1}{\mathrm{e}^{\bar{\alpha}-\bar{\mathcal{B}}_{\mu}p^{\mu}}-\varsigma}.$
(11)

For a fluid composed of massive neutral particles in detailed balance state,
$\bar{\alpha}$ is a constant scalar and $\bar{\mathcal{B}}^{\mu}$ is a
timelike Killing vector field. This in turn requires that the underlying
spacetime must be stationary. We can express $\bar{\mathcal{B}}^{\mu}$ as
$\bar{\mathcal{B}}^{\mu}=\bar{\beta}\bar{U}^{\mu}$ where $\bar{U}^{\mu}$
satisfies $\bar{U}^{\mu}\bar{U}_{\mu}=-c^{2}$, and we shall see later that
$\bar{U}^{\mu}$ represents the proper velocity of the fluid in detailed
balance.

Substituting eq.(11) into eqs.(2) and (3), we can obtain the particle flow
$\bar{N}^{\mu}$, the energy momentum tensor $\bar{T}^{\mu\nu}$ and the entropy
flow $\bar{S}^{\mu}$ in detailed balance

$\displaystyle\bar{N}^{\mu}$ $\displaystyle=$
$\displaystyle\bar{n}\,\bar{U}^{\mu}$ (12) $\displaystyle=$
$\displaystyle\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\bar{J}_{d-1,1}\bar{U}^{\mu},$
$\displaystyle\bar{T}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\bar{\epsilon}\,\bar{U}^{\mu}\bar{U}^{\nu}+\bar{P}\bar{\Delta}^{\mu\nu}$
(13) $\displaystyle=$
$\displaystyle\frac{\mathfrak{g}}{h^{d}}(mc)^{d+1}c\mathcal{A}_{d-1}\left(\bar{J}_{d-1,2}\frac{1}{c^{2}}\bar{U}^{\mu}\bar{U}^{\nu}+\frac{1}{d}\bar{J}_{d+1,0}\bar{\Delta}^{\mu\nu}\right),$
$\displaystyle\bar{S}^{\mu}$ $\displaystyle=$
$\displaystyle\bar{s}\,\bar{U}^{\mu}$ (14) $\displaystyle=$ $\displaystyle
k_{B}\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\left(\bar{\alpha}\bar{J}_{d-1,1}+\bar{\beta}mc^{2}\bar{J}_{d-1,2}+\frac{1}{d}\bar{\beta}mc^{2}\bar{J}_{d+1,0}\right)\bar{U}^{\mu},$

where $\bar{J}_{m,n}$ is a function in
$(\bar{\alpha},\bar{\zeta}=\bar{\beta}mc^{2})$ defined as

$\bar{J}_{m,n}(\bar{\alpha},\bar{\zeta})\equiv\int_{0}^{\infty}\frac{\sinh^{m}\vartheta\cosh^{n}\vartheta}{\mathrm{e}^{\bar{\alpha}+\bar{\zeta}\cosh\vartheta}-\varsigma}\mathrm{d}\vartheta,$
(15)

$\mathcal{A}_{d-1}$ is the area of the $(d-1)$-dimensional unit sphere
$S^{d-1}$, and
$\displaystyle\vspace{4pt}\bar{\Delta}^{\mu\nu}=g^{\mu\nu}+\frac{1}{c^{2}}\bar{U}^{\mu}\bar{U}^{\nu}$.
It is easy to see that $\bar{N}^{\mu}$ is parallel to $\bar{U}^{\mu}$, and
also $\bar{U}^{\mu}$ is the unique normalized timelike eigenvector of
$\bar{T}^{\mu\nu}$. This means that the particle transport and energy
transport directions are parallel to $\bar{U}^{\mu}$. For this reason,
$\bar{U}^{\mu}$ is interpreted as the proper velocity of the fluid.

From eqs.(12)-(14) we can extract the particle number density $\bar{n}$,
energy density $\bar{\epsilon}$, pressure $\bar{P}$ and entropy density
$\bar{s}$ measured by the comoving observer $\bar{\mathcal{O}}$ (i.e. that
with proper velocity $\bar{U}^{\mu}$) in detailed balance, which are all
function in $(\bar{\alpha},\bar{\zeta})$:

$\displaystyle\bar{n}$ $\displaystyle=$
$\displaystyle-\frac{1}{c^{2}}\bar{U}_{\mu}\bar{N}^{\mu}=\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\bar{J}_{d-1,1},$
$\displaystyle\bar{\epsilon}$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\bar{U}_{\mu}\bar{U}_{\nu}\bar{T}^{\mu\nu}=\frac{\mathfrak{g}}{h^{d}}m^{d+1}c^{d+2}\mathcal{A}_{d-1}\bar{J}_{d-1,2},$
$\displaystyle\bar{P}$ $\displaystyle=$
$\displaystyle\frac{1}{d}\bar{\Delta}_{\mu\nu}\bar{T}^{\mu\nu}=\frac{\mathfrak{g}}{h^{d}}m^{d+1}c^{d+2}\frac{\mathcal{A}_{d-1}}{d}\bar{J}_{d+1,0},$
$\displaystyle\bar{s}$ $\displaystyle=$
$\displaystyle-\frac{1}{c^{2}}\bar{U}_{\mu}\bar{S}^{\mu}=k_{B}\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\left(\bar{\alpha}\bar{J}_{d-1,1}+\bar{\zeta}\bar{J}_{d-1,2}+\frac{1}{d}\bar{\zeta}\bar{J}_{d+1,0}\right).$

It’s easy to verify that these four scalars satisfy the local Euler relation
and the Gibbs-Duhem relation

$\bar{s}=k_{B}\left(\bar{\alpha}\bar{n}+\bar{\beta}\bar{\epsilon}+\bar{\beta}\bar{P}\right),$
(16)
$-\bar{s}\,\mathrm{d}\left(\frac{1}{k_{B}\bar{\beta}}\right)+\mathrm{d}\bar{P}+\bar{n}\,\mathrm{d}\left(\frac{\bar{\alpha}}{\bar{\beta}}\right)=0,$
(17)

which implies that $\bar{\alpha},\bar{\beta}$ have further thermodynamic
correspondence

$\bar{\alpha}=-\frac{\bar{\mu}}{k_{B}\bar{T}},\quad\bar{\beta}=\frac{1}{k_{B}\bar{T}},$
(18)

where $\bar{\mu}$ is the chemical potential and $\bar{T}$ is the temperature
in detailed balance.

Recall that $\bar{\alpha}$ is constant scalar and
$\bar{\mathcal{B}}^{\mu}=\bar{\beta}\bar{U}^{\mu}$ is Killing, i.e.

$\nabla_{\mu}\bar{\alpha}=0,\quad\nabla_{(\mu}\bar{\mathcal{B}}_{\nu)}=0.$
(19)

Substituting eq.(18) into eq.(19) and decomposing the results into scalar,
vector and tensor parts, we get the following equations. For the scalar part,
we have

$\bar{U}^{\mu}\nabla_{\mu}\bar{T}=0,\quad\bar{U}^{\mu}\nabla_{\mu}\bar{\mu}=0,\quad\nabla_{\mu}\bar{U}^{\mu}=0,$
(20)

which means that the temperature and chemical potential of the fluid in
detailed balance remain constant along the direction of motion, and the fluid
has no expansion. For the vector part, we have

$\nabla_{\mu}\bar{T}+\frac{\bar{T}}{c^{2}}\bar{U}^{\nu}\nabla_{\nu}\bar{U}_{\mu}=0,\quad\nabla_{\mu}\bar{\mu}+\frac{\bar{\mu}}{c^{2}}\bar{U}^{\nu}\nabla_{\nu}\bar{U}_{\mu}=0,$
(21)

these two equations encode the well-known Tolman-Ehrenfest and the Klein
effects. For the tensor part,

$\Delta^{\mu\rho}\Delta^{\nu\sigma}\nabla_{(\rho}\bar{U}_{\sigma)}-\frac{1}{d}\nabla_{\rho}\bar{U}^{\rho}\bar{\Delta}^{\mu\nu}=0,$
(22)

which indicates that the fluid has no shear effect in detailed balance.

In later calculations, we will not replace the variables
$(\bar{\alpha},\bar{\beta})$ by $(\bar{T},\bar{\mu})$ for convenience. In
terms of $(\bar{\alpha},\bar{\beta})$, eqs.(20) and (21) can be written as

$\bar{U}^{\mu}\nabla_{\mu}\bar{\alpha}=0,\quad\bar{U}^{\mu}\nabla_{\mu}\bar{\beta}=0,\quad\nabla_{\mu}\bar{U}^{\mu}=0,$
(23)
$\Delta^{\mu\nu}\nabla_{\nu}\bar{\alpha}=0,\quad\nabla_{\mu}\bar{\beta}-\frac{\bar{\beta}}{c^{2}}\bar{U}^{\nu}\nabla_{\nu}\bar{U}_{\mu}=0.$
(24)

We will see in Section 4 that the deviation from 0 of eqs. (22), (23) and (24)
leads to correction terms for the particle flow and the energy momentum
tensor.

##  3 Deviation from the detailed balance

For the near equilibrium fluid, we can reasonably use the local equilibrium
assumption to approximate the zeroth order 1PDF as

$f^{(0)}=\frac{\mathfrak{g}}{h^{d}}\frac{1}{\mathrm{e}^{\alpha-\mathcal{B}_{\mu}p^{\mu}}-\varsigma},$
(25)

and we still have the hydrodynamics implications

$\alpha=-\cfrac{\mu}{k_{B}T},\quad\mathcal{B}^{\mu}=\beta
U^{\mu}=\dfrac{1}{k_{B}T}U^{\mu},$ (26)

where $\mu$ is the local chemical potential, $T$ is the local temperature and
$U^{\mu}$ is still referred to as the velocity of the fluid. However, $\alpha$
is no longer the constant and $\mathcal{B}^{\mu}$ needs not be Killing.
Therefore $\nabla_{\mu}\alpha$ and $\nabla_{(\mu}\mathcal{B}_{\nu)}$ provide a
measure for how far the state of the fluid deviates from detailed balance. It
should be remarked that the previous works [20, 21] considered only the case
in which only the magnitude of $\mathcal{B}^{\mu}$ differs from the detailed
balance value. However, in this work, we allow both the magnitude and the
direction of $\mathcal{B}^{\mu}$ differ from the detailed balance case.

A subtle point that needs to be clarified is that for non-equilibrium systems
(even if the system is only slightly out of equilibrium), there is no prior
definition for temperature, chemical potential and fluid velocity. One is free
to define temperature, chemical potential and fluid velocity in different
ways, but different definitions are often not equivalent [26]. Choosing a set
of definitions of temperature, chemical potential and fluid velocity is known
as selecting a hydrodynamic frame [22]. Different hydrodynamic frames affect
the division of fluid orders without changing the physics itself. In the
subsequent calculations, we just keep the form of $f^{(0)}$ without
constraints on $\alpha$ and $\mathcal{B}^{\mu}$, which means that our
calculation is not carried out in a specific hydrodynamic frame. When
discussing the transport coefficients in Section 4, we will return to the
specific hydrodynamic frame by imposing constraints on $\alpha,\beta$ and
$U^{\mu}$.

The full Boltzmann equation is an extremely complicated integral-differential
equation, where the collision integral makes it difficult to solve
analytically. In this work, we adopt the standard relaxation time
approximation by replacing the collision integral with an Anderson-Witting-
like collision model

$\mathcal{C}(x,p)\simeq-\frac{\varepsilon}{c^{2}\tau}(f-f^{(0)})=\frac{U_{\mu}p^{\mu}}{c^{2}\tau}(f-f^{(0)}),$
(27)

where $\varepsilon=-U_{\mu}p^{\mu}$ is the energy of a single particle
measured by the observer $\mathcal{O}$, and $\tau$ is referred to as the
relaxation time which represents the time scale for the system to restore
balance.

Next, we shall try to solve the Boltzmann equation using a standard iterative
procedure. To the first order in relaxation time, we have

$\mathcal{L}_{H}f^{(0)}=-\frac{\varepsilon}{c^{2}\tau}(f_{[1]}-f^{(0)}),$ (28)

where $f_{[1]}=f^{(0)}+f^{(1)}$ denotes the first order iterative solution. In
general, the iterative solution of order $n$ satisfy

$\mathcal{L}_{H}f_{[n-1]}=-\frac{\varepsilon}{c^{2}\tau}(f_{[n]}-f^{(0)}),$
(29)

where we can solve the expression for $f_{[n]}$:

$f_{[n]}=\sum_{i=0}^{n}f^{(i)}=\sum_{i=0}^{n}(-1)^{i}\left(\frac{c^{2}\tau}{\varepsilon}\right)^{i}(\mathcal{L}_{H})^{i}f^{(0)}.$
(30)

While there are potential issues regarding the convergence of the iteration
procedure even in the non-relativistic case [27, 28, 29], it is known that a
cut-off at the first or second order yields reasonable results for describing
the behaviors of lower order fluids. Following this spirit, in the rest part
of this work, we will focus on the first order iterative solution.

Substituting eq.(25) into eq.(6), we have

$\mathcal{L}_{H}f^{(0)}=(-p^{\mu}p^{\nu}\nabla_{\mu}\mathcal{B}_{\nu}+p^{\mu}\nabla_{\mu}\alpha)\frac{\partial
f^{(0)}}{\partial\alpha}.$ (31)

Thus the first order iterative solution reads

$f_{[1]}=f^{(0)}+f^{(1)}=f^{(0)}-\frac{c^{2}\tau}{\varepsilon}(-p^{\mu}p^{\nu}\nabla_{\mu}\mathcal{B}_{\nu}+p^{\mu}\nabla_{\mu}\alpha)\frac{\partial
f^{(0)}}{\partial\alpha}.$ (32)

It is clear that $\nabla_{(\mu}\mathcal{B}_{\nu)}$, $\nabla_{\mu}\alpha$
represents the deviation of the system from the detailed balance.

##  4 First order kinetic coefficients

Using eqs.(2) and (32), we can analytically calculate the particle flow and
the energy momentum tensor (see Appendix A for details). At the zeroth order,
we have

$\displaystyle N^{(0)\mu}$
$\displaystyle=c\int\boldsymbol{\varpi}p^{\mu}f^{(0)}=n^{(0)}U^{\mu}$
$\displaystyle=\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}J_{d-1,1}U^{\mu},$
(33) $\displaystyle T^{(0)\mu\nu}$
$\displaystyle=c\int\boldsymbol{\varpi}p^{\mu}p^{\nu}f^{(0)}=\frac{1}{c^{2}}\epsilon^{(0)}U^{\mu}U^{\nu}+P\Delta^{\mu\nu}$
$\displaystyle=\frac{\mathfrak{g}}{h^{d}}m^{d+1}c^{d+2}\mathcal{A}_{d-1}\left(J_{d-1,2}\frac{1}{c^{2}}U^{\mu}U^{\nu}+\frac{1}{d}J_{d+1,0}\Delta^{\mu\nu}\right),$
(34)

where $J_{m,n}$ is a function in $(\alpha,\zeta=\beta mc^{2})$

$J_{m,n}(\alpha,\zeta)\equiv\int_{0}^{\infty}\frac{\sinh^{m}\vartheta\cosh^{n}\vartheta}{\mathrm{e}^{\alpha+\zeta\cosh\vartheta}-\varsigma}\mathrm{d}\vartheta.$
(35)

At the first order, using the relation $\mathcal{B}^{\mu}=\beta U^{\mu}$, we
can simplify the results as

$\displaystyle N^{(1)\mu}$
$\displaystyle=c\int\boldsymbol{\varpi}p^{\mu}f^{(1)}$
$\displaystyle=n^{(1)}U^{\mu}+j^{(1)\mu},$ (36) $\displaystyle T^{(1)\mu\nu}$
$\displaystyle=c\int\boldsymbol{\varpi}p^{\mu}p^{\nu}f^{(1)}$
$\displaystyle=\frac{1}{c^{2}}\epsilon^{(1)}U^{\mu}U^{\nu}+\frac{1}{c^{2}}q^{(1)\mu}U^{\nu}+\frac{1}{c^{2}}q^{(1)\nu}U^{\mu}+\Pi^{(1)}\Delta^{\mu\nu}+\Pi^{(1)\mu\nu},$
(37)

where the first order hydrodynamic variables are

$\displaystyle n^{(1)}$ $\displaystyle=-\frac{1}{c^{2}}U_{\mu}N^{(1)\mu}$
$\displaystyle=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d}\left(\frac{\partial
J_{d-1,1}}{\partial\alpha}U^{\nu}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\nu}\nabla_{\nu}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,0}}{\partial\alpha}\beta\nabla_{\nu}U^{\nu}\right),$ (38)
$\displaystyle\epsilon^{(1)}$
$\displaystyle=\frac{1}{c^{2}}U_{\mu}U_{\nu}T^{(1)\mu\nu}$
$\displaystyle=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}c\left(\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\rho}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d-1,3}}{\partial\alpha}U^{\rho}\nabla_{\rho}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,1}}{\partial\alpha}\beta\nabla_{\rho}U^{\rho}\right),$ (39)
$\displaystyle\Pi^{(1)}$
$\displaystyle=\frac{1}{d}T^{(1)\mu\nu}\Delta_{\mu\nu}$
$\displaystyle=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\frac{c}{d}\left(\frac{\partial
J_{d+1,0}}{\partial\alpha}U^{\rho}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}U^{\rho}\nabla_{\rho}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\beta\nabla_{\rho}U^{\rho}\right),$ (40)
$\displaystyle j^{(1)\mu}$
$\displaystyle=\Delta^{\mu}_{\enspace\nu}N^{(1)\nu}$
$\displaystyle=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d}\frac{c^{2}}{d}\left(\frac{\partial
J_{d+1,-1}}{\partial\alpha}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d+1,0}}{\partial\alpha}\left(\nabla_{\nu}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\nu}\right)\right)\Delta^{\mu\nu},$
(41) $\displaystyle q^{(1)\mu}$
$\displaystyle=-U_{\nu}T^{(1)\nu\sigma}\Delta^{\mu}_{\enspace\sigma}$
$\displaystyle=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\frac{c^{3}}{d}\left(\frac{\partial
J_{d+1,0}}{\partial\alpha}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}\left(\nabla_{\rho}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\sigma}\right)\right)\Delta^{\rho\mu},$
(42) $\displaystyle\Pi^{(1)\mu\nu}$
$\displaystyle=T^{(1)\rho\sigma}\Delta^{\mu}_{\enspace\rho}\Delta^{\nu}_{\enspace\sigma}-\frac{1}{d}T^{(1)\rho\sigma}\Delta_{\rho\sigma}\Delta^{\mu\nu}$
$\displaystyle=\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+2}\frac{2c^{2}}{(d+2)d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\beta\left(\Delta^{\rho\mu}\Delta^{\sigma\nu}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta^{\mu\nu}\right).$
(43)

From the tensor part
$\displaystyle\Pi^{(1)\mu\nu}=-\eta\left(\Delta^{\rho\mu}\Delta^{\sigma\nu}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta^{\mu\nu}\right)$,
we can read off the shear viscosity coefficient

$\eta=-c^{2}\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+2}\frac{2}{(d+2)d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\beta.$ (44)

According to the definition of the special function $J_{m,n}$ (35), we can
easily conclude that $\eta>0$ and $\eta$ increases with temperature.

Now, introducing a new variable $u=\zeta(\cosh\vartheta-1)$, we can rewrite
the integral $J_{m,n}$ as

$J_{m,n}=\int_{0}^{\infty}\frac{\left(\frac{u^{2}}{\zeta^{2}}+2\frac{u}{\zeta}\right)^{\frac{m-1}{2}}\left(\frac{u}{\zeta}+1\right)^{n}}{\zeta(\mathrm{e}^{\alpha+\zeta+u}-\varsigma)}\mathrm{d}u.$

In the high temperature limit, $\beta\rightarrow 0$, $\zeta=\beta
mc^{2}\rightarrow 0$, we have

$\displaystyle J_{m,n}$ $\displaystyle\rightarrow$
$\displaystyle\left(\frac{1}{\zeta}\right)^{m+n}\int_{0}^{\infty}\frac{u^{m+n-1}}{\mathrm{e}^{\alpha+u}-\varsigma}\mathrm{d}u=\left(\frac{1}{\zeta}\right)^{m+n}\frac{1}{\varsigma}\Gamma(m+n)\operatorname{Li}_{m+n}(\mathrm{e}^{-\alpha}\varsigma),$
(45) $\displaystyle\frac{\partial J_{m,n}}{\partial\alpha}$
$\displaystyle\rightarrow$
$\displaystyle-\left(\frac{1}{\zeta}\right)^{m+n}\frac{1}{\varsigma}\Gamma(m+n)\operatorname{Li}_{m+n-1}(\mathrm{e}^{-\alpha}\varsigma).$
(46)

It is easy to see that the shear viscosity coefficient $\eta$ diverges in the
high temperature limit due to the asymptotic behavior $\cfrac{\partial
J_{d+3,-1}}{\partial\alpha}\beta\sim\left(\cfrac{1}{\zeta}\right)^{d+1}$.

The vector part $\\{j^{(1)\mu},q^{(1)\mu}\\}$ has been investigated in the
previous works [20, 21], and it has been shown that the Onsager reciprocal
relation and Wiedemann-Franz law hold, and the heat and electric
conductivities have divergent behavior near the ergosphere. Nevertheless, it
is meaningful to show again the procedure for obtaining the heat conductivity
coefficient because now we are using a different solution of the Boltzmann
equation.

In Eckart frame, $n^{(1)}=0,\epsilon^{(1)}=0,j^{(1)}=0$. The last condition
$j^{(1)}=0$ gives a constraint between $\nabla_{\nu}\alpha$ and
$\displaystyle\nabla_{\nu}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\nu}$:

$\frac{\partial
J_{d+1,-1}}{\partial\alpha}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d+1,0}}{\partial\alpha}\left(\nabla_{\nu}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\nu}\right)=0.$
(47)

In this frame, the energy flow is simply the heat flow. We can eliminate
$\nabla_{\nu}\alpha$ and making use of the relation $T=\dfrac{1}{k_{B}\beta}$
to express the heat flow as

$q^{(1)\mu}=-\kappa\left(\nabla_{\rho}T+\frac{T}{c^{2}}U^{\rho}\nabla_{\rho}U_{\sigma}\right)\Delta^{\rho\mu},$
(48)

where the heat conductivity coefficient reads

$\kappa=-\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+2}\frac{c^{4}}{d}\frac{1}{k_{B}T^{2}}\left(-\frac{\left(\dfrac{\partial
J_{d+1,0}}{\partial\alpha}\right)^{2}}{\dfrac{\partial
J_{d+1,-1}}{\partial\alpha}}+\frac{\partial
J_{d+1,1}}{\partial\alpha}\right).$ (49)

It is east to verify that $\kappa>0$ and $\kappa$ increases with temperature.
In the high temperature limit, $\kappa\sim\left(\cfrac{1}{\zeta}\right)^{d}$.

For the scalar part $\\{n^{(1)},\epsilon^{(1)},\Pi^{(1)}\\}$, the bulk
viscosity coefficient is usually discussed in the Eckart frame or Landau
frame. These two frames require $n^{(1)}=0,\epsilon^{(1)}=0$, which gives us
two constraints

$\displaystyle\frac{\partial
J_{d-1,1}}{\partial\alpha}U^{\nu}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\nu}\nabla_{\nu}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,0}}{\partial\alpha}\beta\nabla_{\nu}U^{\nu}$ $\displaystyle=0,$ (50)
$\displaystyle\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\rho}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d-1,3}}{\partial\alpha}U^{\rho}\nabla_{\rho}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,1}}{\partial\alpha}\beta\nabla_{\rho}U^{\rho}$ $\displaystyle=0.$ (51)

From these two equations, we can express $U^{\nu}\nabla_{\nu}\alpha$ and
$U^{\nu}\nabla_{\nu}\beta$ in terms of $\nabla_{\nu}U^{\nu}$, so that
$\Pi^{(1)}$ can be written as

$\Pi^{(1)}=-\xi\nabla_{\rho}U^{\rho},$ (52)

where the bulk viscosity coefficient $\xi$ reads

$\displaystyle\xi$
$\displaystyle=c^{2}\tau\frac{\mathfrak{g}}{h^{d}}\frac{\mathcal{A}_{d-1}}{d^{2}}(mc)^{d+2}\beta\cdot$
$\displaystyle~{}~{}~{}~{}\left(\frac{\left(\dfrac{\partial
J_{d+1,0}}{\partial\alpha}\right)^{2}\dfrac{\partial
J_{d-1,3}}{\partial\alpha}+\left(\dfrac{\partial
J_{d+1,1}}{\partial\alpha}\right)^{2}\dfrac{\partial
J_{d-1,1}}{\partial\alpha}-2\dfrac{\partial
J_{d-1,2}}{\partial\alpha}\dfrac{\partial
J_{d+1,1}}{\partial\alpha}\dfrac{\partial
J_{d+1,0}}{\partial\alpha}}{-\left(\dfrac{\partial
J_{d-1,2}}{\partial\alpha}\right)^{2}+\dfrac{\partial
J_{d-1,1}}{\partial\alpha}\dfrac{\partial
J_{d-1,3}}{\partial\alpha}}-\frac{\partial
J_{d+3,-1}}{\partial\alpha}\right),$ (53)

which, in the high temperature limit, behaves as
$\xi\sim\left(\cfrac{1}{\zeta}\right)^{d-3}$. The different asymptotic
behaviors for the shear viscosity coefficient $\eta$ and the bulk viscosity
coefficient $\xi$ indicate that these two viscosities are of different order
of magnitude in nature. The above results extends Kremer’s works [18, 19] on
shear and bulk viscosity to generic stationary backgrounds.

##  5 Onsager reciprocal relation and local entropy production

According to the definition of entropy flow $S^{\mu}$ (3), up to the zeroth
order in relaxation time, we can verify the covariant Euler relation and the
covariant Gibbs-Duhem relation

$\displaystyle S^{(0)\mu}=$
$\displaystyle~{}-k_{B}c\int\boldsymbol{\varpi}p^{\mu}f^{(0)}\left[\log\left(\frac{h^{d}f^{(0)}}{{f^{(0)}}^{\ast}\mathfrak{g}}\right)-\frac{{\mathfrak{g}\log
f^{(0)}}^{\ast}}{\varsigma h^{d}f^{(0)}}\right]$ $\displaystyle=$
$\displaystyle~{}k_{B}(\Xi^{\mu}-\mathcal{B}_{\nu}T^{(0)\mu\nu}+\alpha
N^{(0)\mu}),$ (54)
$\mathrm{d}\Xi^{\mu}=T^{(0)\mu\nu}\mathrm{d}\mathcal{B}_{\nu}-N^{(0)\mu}\mathrm{d}\alpha,$
(55)

where

$\Xi^{\mu}=P\mathcal{B}^{\mu}=c\int\boldsymbol{\varpi}p^{\mu}\frac{\mathfrak{g}}{\varsigma
h^{d}}{\log\bigg{(}1+\frac{\varsigma h^{d}}{\mathfrak{g}}f\bigg{)}}.$ (56)

Considering that in the first order, the deviation $f^{(1)}$ from the local
equilibrium is a small correction, then

$\displaystyle S^{(1)\mu}=$
$\displaystyle~{}-k_{B}c\int\boldsymbol{\varpi}p^{\mu}f^{(1)}\log\frac{f^{(0)}}{\mathfrak{g}/h^{d}+\varsigma
f^{(0)}}$ $\displaystyle=$
$\displaystyle~{}k_{B}c\int\boldsymbol{\varpi}p^{\mu}f^{(1)}(\alpha-\mathcal{B}_{\nu}p^{\nu})$
$\displaystyle=$ $\displaystyle~{}k_{B}(\alpha
N^{(1)\mu}-\mathcal{B}_{\nu}T^{(1)\mu\nu}).$ (57)

Therefore, up to the first order, we can conclude that

$S^{\mu}=S^{(0)\mu}+S^{(1)\mu}=k_{B}(P\mathcal{B}^{\mu}-\mathcal{B}_{\nu}T^{\mu\nu}+\alpha
N^{\mu}).$ (58)

Up to the first order, using the conservation equation
$\nabla_{\mu}N^{\mu}=0$, $\nabla_{\mu}T^{\mu\nu}=0$ and the corollary derived
from the Gibbs-Duhem relation

$\nabla_{\mu}\Xi^{\mu}=T^{(0)\mu\nu}\nabla_{\mu}\mathcal{B}_{\nu}-N^{(0)\mu}\nabla_{\mu}\alpha,$
(59)

we can obtain

$\displaystyle\frac{1}{k_{B}}\nabla_{\mu}S^{\mu}=$
$\displaystyle~{}-T^{(1)\mu\nu}\nabla_{\mu}\mathcal{B}_{\nu}+N^{(1)\mu}\nabla_{\mu}\alpha$
$\displaystyle=$
$\displaystyle~{}n^{(1)}U^{\mu}\nabla_{\mu}\alpha+\epsilon^{(1)}U^{\mu}\nabla_{\mu}\beta-\Pi^{(1)}\beta\nabla_{\mu}U^{\mu}$
$\displaystyle~{}+j^{(1)\mu}\Delta_{\mu}^{\enspace\nu}\nabla_{\nu}\alpha+q^{(1)\mu}\Delta_{\mu}^{\enspace\rho}\left(\nabla_{\rho}\beta-\frac{\beta}{c^{2}}U^{\sigma}\nabla_{\sigma}U_{\rho}\right)$
$\displaystyle~{}-\Pi^{(1)\mu\nu}\beta\left(\Delta_{\mu}^{\enspace\rho}\Delta_{\nu}^{\enspace\sigma}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta_{\mu\nu}\right).$
(60)

It is worth noting that in addition to the vector part
$\\{j^{(1)\mu},q^{(1)\mu}\\}$, the entropy production also contains
contributions from the scalar part $\\{n^{(1)},\epsilon^{(1)},\Pi^{(1)}\\}$
and the tensor part $\\{\Pi^{(1)\mu\nu}\\}$.

Recalling the expressions of the first order quantity of hydrodynamics
(38)-(43), we can simply write the constitutive relations of the first order
hydrodynamic variables as

$\displaystyle n^{(1)}$
$\displaystyle=\nu_{1}U^{\nu}\nabla_{\nu}\alpha+\nu_{2}U^{\nu}\nabla_{\nu}\beta-\nu_{3}\beta\nabla_{\nu}U^{\nu},$
$\displaystyle\epsilon^{(1)}$
$\displaystyle=\rho_{1}U^{\nu}\nabla_{\nu}\alpha+\rho_{2}U^{\nu}\nabla_{\nu}\beta-\rho_{3}\beta\nabla_{\nu}U^{\nu},$
$\displaystyle\Pi^{(1)}$
$\displaystyle=\pi_{1}U^{\nu}\nabla_{\nu}\alpha+\pi_{2}U^{\nu}\nabla_{\nu}\beta-\pi_{3}\beta\nabla_{\nu}U^{\nu},$
$\displaystyle j^{(1)\mu}$
$\displaystyle=\sigma_{1}\Delta^{\mu\nu}\nabla_{\nu}\alpha+\sigma_{2}\Delta^{\mu\nu}\left(\nabla_{\nu}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\nu}\right),$
(61) $\displaystyle q^{(1)\mu}$
$\displaystyle=\kappa_{1}\Delta^{\mu\nu}\nabla_{\nu}\alpha+\kappa_{2}\Delta^{\mu\nu}\left(\nabla_{\nu}\beta-\frac{\beta}{c^{2}}U^{\rho}\nabla_{\rho}U_{\nu}\right),$
$\displaystyle\Pi^{(1)\mu\nu}$
$\displaystyle=-\eta\left(\Delta^{\rho\mu}\Delta^{\sigma\nu}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta^{\mu\nu}\right),$

where all the kinetic coefficients
$\\{\nu_{i},\sigma_{i},\rho_{i},\kappa_{i},\pi_{i},\eta\\}$ are calculated in
eqs.(38)-(43). It is easy to see that

$\displaystyle\boldsymbol{F}\equiv$
$\displaystyle~{}\biggl{[}U^{\mu}\nabla_{\mu}\alpha,U^{\mu}\nabla_{\mu}\beta,-\beta\nabla_{\mu}U^{\mu},\Delta^{\mu\nu}\nabla_{\nu}\alpha,\Delta^{\mu\rho}\left(\nabla_{\rho}\beta-\frac{\beta}{c^{2}}U^{\sigma}\nabla_{\sigma}U_{\rho}\right),-\beta\left(\Delta^{\mu\rho}\Delta^{\nu\sigma}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta^{\mu\nu}\right)\biggl{]}$
(62)

are the generalized force which cause the local entropy production and vanish
in the detailed balance state (c.f. eqs.(23),(24) and (22)).

It remains to verify that all the kinetic coefficients in (61) satisfy a more
complete Onsager reciprocity relation. This can be easily achieved by
observing that the coefficient matrix is not only symmetric in the vector part
of the constitutive relation, but also symmetric in the scalar part, i.e.

$\left(\begin{array}[]{cc}\sigma_{1}&\sigma_{2}\\\
\kappa_{1}&\kappa_{2}\end{array}\right),\quad\left(\begin{array}[]{ccc}\nu_{1}&\nu_{2}&\nu_{3}\\\
\rho_{1}&\rho_{2}&\rho_{3}\\\ \pi_{1}&\pi_{2}&\pi_{3}\end{array}\right)$ (63)

are both symmetric.

Another interesting thing is to investigate whether the approximate solution
(32) ensure a non-negative entropy production at the first order in relaxation
time. Substituting eq.(61) into eq.(60) we can obtain a quadratic form

$\frac{1}{k_{B}}\nabla_{\mu}S^{\mu}=\boldsymbol{F}\left(\begin{array}[]{ccc|cc|c}\nu_{1}&(\nu_{2}+\rho_{1})/2&(\nu_{3}+\pi_{1})/2&0&0&0\\\
(\nu_{2}+\rho_{1})/2&\rho_{2}&(\rho_{3}+\pi_{2})/2&0&0&0\\\
(\nu_{3}+\pi_{1})/2&(\rho_{3}+\pi_{2})/2&\pi_{3}&0&0&0\\\ \hline\cr
0&0&0&\sigma_{1}&(\sigma_{2}+\kappa_{1})/2&0\\\
0&0&0&(\sigma_{2}+\kappa_{1})/2&\kappa_{2}&0\\\ \hline\cr
0&0&0&0&0&\eta/\beta\end{array}\right)\boldsymbol{F}^{T},$ (64)

and the non-negative entropy production means that all eigenvalues of the
quadratic form are non-negative.

For the tensor part $\\{\Pi^{(1)\mu\nu}\\}$, the only transport coefficient
$\eta>0$ ensures that the eigenvalue of the tensor part is non-negative.

For the vector part $\\{j^{(1)\mu},q^{(1)\mu}\\}$, the condition that the
quadratic form is non-negative can be written as $\sigma_{1}\geqslant
0,4\sigma_{1}\kappa_{2}-(\sigma_{2}+\kappa_{1})^{2}\geqslant 0$. Substituting
the value of the transport coefficient, it is easy to verify that both
conditions are satisfied

$\displaystyle\sigma_{1}=$
$\displaystyle~{}-c^{2}\tau\frac{\mathfrak{g}}{h^{d}}\frac{\mathcal{A}_{d-1}}{d}(mc)^{d}\frac{\partial
J_{d+1,-1}}{\partial\alpha}>0,$ (65) $\displaystyle
4\sigma_{1}\kappa_{2}-(\sigma_{2}+\kappa_{1})^{2}=$
$\displaystyle~{}4\left(\tau\frac{\mathfrak{g}}{h^{d}}(mc)^{d+1}\mathcal{A}_{d-1}\frac{c^{3}}{d}\right)^{2}\left[\frac{\partial
J_{d+1,-1}}{\partial\alpha}\frac{\partial
J_{d+1,1}}{\partial\alpha}-\left(\frac{\partial
J_{d+1,0}}{\partial\alpha}\right)^{2}\right]>0.$ (66)

Figure 1: The smallest eigenvalue of the dimensionless coefficient matrix for
the scalar part. The left plot depicts the smallest eigenvalue as a function
of $\alpha^{\ast}=\alpha+\zeta$ at fixed $\zeta=\beta mc^{2}=0.01$. The right
plot depicts the smallest eigenvalue as a function of $\zeta$ at fixed
$\alpha^{\ast}=1$. The curve with $\varsigma=1$ in the left plot does not
cross the vertical line at $\alpha^{\ast}=0$ because the chemical potential
$\mu$ of relativistic Bose gas is not greater than $mc^{2}$.

For the scalar part $\\{n^{(1)},\epsilon^{(1)},\Pi^{(1)}\\}$ (which
corresponds to the upper-left block of the matrix in eq.(64)), it is too
difficult to calculate the eigenvalues analytically. However, based on
numerical calculations, no eigenvalue less than zero is found. We plot the
dimensionless coefficient matrix’s smallest eigenvalue for the scalar part,
see Figure 1.

##  6 Conclusions and remarks

In this work, we have iteratively solved the relativistic Boltzmann equation
under the relaxation time approximation in generic stationary spacetime,
calculated the first order hydrodynamic variables using the first order
solution, and analyzed the first order kinetic coefficients. Moreover, given
the important role of Onsager reciprocal relation in conventional non-
equilibrium statistical physics [30, 31], we further analyzed these kinetic
coefficients. Our calculations show that, up to the first order of the
relaxation time the kinetic coefficients satisfy a more generalized Onsager
reciprocal relation in the relativistic context, and lead to a non-negative
entropy production up to the first order in relaxation time.

A point need to be clarified is that the kinetic coefficients
$\\{\nu_{i},\sigma_{i},\rho_{i},\kappa_{i},\pi_{i}\\}$ are different from the
transport coefficients $\\{\kappa,\xi\\}$. The kinetic coefficients are only
relevant to the microscopic theory and not to the choice of hydrodynamic
frame, while the definition of transport coefficients depends on the
particular hydrodynamic frame. By imposing constraints on the variables
$\\{\alpha,\beta,U^{\mu}\\}$, it is possible to restrict our results to some
fixed hydrodynamic frame. In the spirit of BDNK theory [32, 33], a causal,
stable, and locally well-posed first order relativistic fluid can be obtained
in some particular hydrodynamic frame.

For relativistic systems, the choice of observer can greatly affect the
phenomena. Therefore, a question that deserves further consideration is the
effect of the choice of observer on the constitutive relations. In addition,
at the first order in relaxation time, all the kinetic coefficients are
scalars. At the second order of the relaxation time, the kinetic coefficients
will be tensors, which lead to more abundant phenomena. We hope to investigate
these questions in later study.

##  Appendix A Calculation of particle flow and energy momentum tensor

We work in orthonormal basis $\\{(e_{\hat{a}})^{\mu}\\}$ obeying
$\eta_{\hat{a}\hat{b}}=g_{\mu\nu}(e_{\hat{a}})^{\mu}(e_{\hat{b}})^{\nu}$.
Without loss of generality, we require that $U^{\mu}=c(e_{\hat{0}})^{\mu}$,
which implies that the induced metric can be expressed as
$\Delta^{\mu\nu}=\delta^{\hat{i}\hat{j}}(e_{\hat{i}})^{\mu}(e_{\hat{j}})^{\nu}$.

For massive particles, the momentum $p^{\hat{a}}=p^{\mu}(e^{\hat{a}})_{\mu}$
can be parameterized by mass shell conditions

$p^{\hat{a}}=mc(\cosh\vartheta,n^{\hat{i}}\sinh\vartheta),$ (67)

where $n^{\hat{i}}\in S^{d-1}$ is a spacelike unit vector. Then the momentum
space volume element be represented as

$\boldsymbol{\varpi}=\frac{(\mathrm{d}p)^{d}}{|p_{\hat{0}}|}=\frac{|\boldsymbol{p}|^{d-1}\mathrm{d}|\boldsymbol{p}|\mathrm{d}\Omega_{d-1}}{p^{\hat{0}}}=(mc\sinh\vartheta)^{d-1}\mathrm{d}\vartheta\mathrm{d}\Omega_{d-1},$
(68)

where $\mathrm{d}\Omega_{d-1}$ is the volume element of the
$(d-1)$-dimensional unit sphere $S^{d-1}$. In this way, the integration in the
momentum space can be decomposed into integration over
$\vartheta\in(0,\infty)$ and integration over the unit sphere $S^{d-1}$.

Here we list some useful integration formulae,

$\displaystyle\int\mathrm{d}\Omega_{d-1}$ $\displaystyle=\mathcal{A}_{d-1},$
(69) $\displaystyle\int n^{\hat{i}}n^{\hat{j}}\mathrm{d}\Omega_{d-1}$
$\displaystyle=\frac{1}{d}\mathcal{A}_{d-1}\delta^{\hat{i}\hat{j}},$ (70)
$\displaystyle\int
n^{\hat{i}}n^{\hat{j}}n^{\hat{k}}n^{\hat{l}}\mathrm{d}\Omega_{d-1}$
$\displaystyle=\frac{3}{(d+2)d}\mathcal{A}_{d-1}\delta^{(\hat{i}\hat{j}}\delta^{\hat{k}\hat{l})},$
(71) $\int n^{\hat{i}}\mathrm{d}\Omega_{d-1}=\int
n^{\hat{i}}n^{\hat{j}}n^{\hat{k}}\mathrm{d}\Omega_{d-1}=0.$ (72)

Using eqs.(67)-(72), we can further calculate

$\displaystyle\int\boldsymbol{\varpi}p^{\mu}f^{(0)}=$
$\displaystyle~{}\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d}J_{d-1,1}\frac{1}{c}U^{\mu},$
(73) $\displaystyle\int\boldsymbol{\varpi}p^{\mu}p^{\nu}f^{(0)}=$
$\displaystyle~{}\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\left(J_{d-1,2}\frac{1}{c^{2}}U^{\mu}U^{\nu}+\frac{1}{d}J_{d+1,0}\Delta^{\mu\nu}\right),$
(74)
$\displaystyle\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}f^{(0)}=$
$\displaystyle~{}\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d}\left(J_{d-1,1}\frac{1}{c^{2}}U^{\mu}U^{\nu}+\frac{1}{d}J_{d+1,-1}\Delta^{\mu\nu}\right),$
(75)
$\displaystyle\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}p^{\sigma}f^{(0)}=$
$\displaystyle~{}\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\left(J_{d-1,2}\frac{1}{c^{3}}U^{\mu}U^{\nu}U^{\sigma}+\frac{3}{d}J_{d+1,0}\frac{1}{c}U^{(\mu}\Delta^{\nu\sigma)}\right),$
(76)
$\displaystyle\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}p^{\sigma}p^{\rho}f^{(0)}=$
$\displaystyle~{}\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+2}\biggl{(}J_{d-1,3}\frac{1}{c^{4}}U^{\mu}U^{\nu}U^{\sigma}U^{\rho}$
$\displaystyle~{}+\frac{6}{d}J_{d+1,1}\frac{1}{c^{2}}U^{(\mu}U^{\nu}\Delta^{\sigma\rho)}+\frac{3}{(d+2)d}J_{d+3,-1}\Delta^{(\mu\nu}\Delta^{\sigma\rho)}\biggl{)}.$
(77)

The following calculations are then straightforward,

$\displaystyle N^{(1)\mu}=$
$\displaystyle~{}c\int\boldsymbol{\varpi}p^{\mu}f^{(1)}$ $\displaystyle=$
$\displaystyle-c^{3}\tau\int\boldsymbol{\varpi}\frac{1}{\varepsilon}p^{\mu}(-p^{\nu}p^{\sigma}\nabla_{\nu}\mathcal{B}_{\sigma}+p^{\nu}\nabla_{\nu}\alpha)\frac{\partial
f^{(0)}}{\partial\alpha}$ $\displaystyle=$
$\displaystyle-c^{2}\tau\left(-\nabla_{\nu}\mathcal{B}_{\sigma}\frac{\partial}{\partial\alpha}\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}p^{\sigma}f^{(0)}+\nabla_{\nu}\alpha\frac{\partial}{\partial\alpha}\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}f^{(0)}\right)$
$\displaystyle=$
$\displaystyle-c^{2}\tau\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\biggl{[}\left(\frac{\partial
J_{d-1,1}}{\partial\alpha}\frac{1}{c^{2}}U^{\nu}\nabla_{\nu}\alpha-m\frac{\partial
J_{d-1,2}}{\partial\alpha}\frac{1}{c^{2}}U^{\nu}U^{\sigma}\nabla_{\nu}\mathcal{B}_{\sigma}-m\frac{1}{d}\frac{\partial
J_{d+1,0}}{\partial\alpha}\Delta^{\nu\sigma}\nabla_{\nu}\mathcal{B}_{\sigma}\right)U^{\mu}$
$\displaystyle+\frac{1}{d}\left(-2m\frac{\partial
J_{d+1,0}}{\partial\alpha}U^{\sigma}\nabla_{(\sigma}\mathcal{B}_{\nu)}+\frac{\partial
J_{d+1,-1}}{\partial\alpha}\nabla_{\nu}\alpha\right)\Delta^{\mu\nu}\biggl{]}$
$\displaystyle=$
$\displaystyle-\tau\frac{\mathfrak{g}}{h^{d}}(mc)^{d}\mathcal{A}_{d-1}\biggl{[}\left(\frac{\partial
J_{d-1,1}}{\partial\alpha}U^{\nu}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\nu}\nabla_{\nu}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,0}}{\partial\alpha}\beta\nabla_{\nu}U^{\nu}\right)U^{\mu}$
$\displaystyle+\frac{c^{2}}{d}\left(\frac{\partial
J_{d+1,-1}}{\partial\alpha}\nabla_{\nu}\alpha+mc^{2}\frac{\partial
J_{d+1,0}}{\partial\alpha}\left(\nabla_{\nu}\beta-\frac{1}{c^{2}}\beta
U^{\rho}\nabla_{\rho}U_{\nu}\right)\right)\Delta^{\mu\nu}\biggl{]},$ (78)
$\displaystyle T^{(1)\mu\nu}=$
$\displaystyle~{}c\int\boldsymbol{\varpi}p^{\mu}p^{\nu}f^{(1)}$
$\displaystyle=$
$\displaystyle-c^{3}\tau\int\boldsymbol{\varpi}\frac{1}{\varepsilon}p^{\mu}p^{\nu}(-p^{\rho}p^{\sigma}\nabla_{\rho}\mathcal{B}_{\sigma}+p^{\sigma}\nabla_{\sigma}\alpha)\frac{\partial
f^{(0)}}{\partial\alpha}$ $\displaystyle=$
$\displaystyle-c^{2}\tau\left(-\nabla_{\rho}\mathcal{B}_{\sigma}\frac{\partial}{\partial\alpha}\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}p^{\rho}p^{\sigma}f^{(0)}+\nabla_{\sigma}\alpha\frac{\partial}{\partial\alpha}\int\boldsymbol{\varpi}\frac{1}{p^{\hat{0}}}p^{\mu}p^{\nu}p^{\sigma}f^{(0)}\right)$
$\displaystyle=$
$\displaystyle-c\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\biggl{[}\frac{1}{c^{2}}\left(-mc^{2}\frac{\partial
J_{d-1,3}}{\partial\alpha}\frac{1}{c^{2}}U^{\sigma}U^{\rho}\nabla_{\rho}\mathcal{B}_{\sigma}-mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,1}}{\partial\alpha}\nabla_{\rho}\mathcal{B}_{\sigma}\Delta^{\sigma\rho}+\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\sigma}\nabla_{\sigma}\alpha\right)U^{\mu}U^{\nu}$
$\displaystyle+\frac{1}{d}\left(-2mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}\frac{1}{c^{2}}U^{\sigma}\nabla_{(\rho}\mathcal{B}_{\sigma)}+\frac{\partial
J_{d+1,0}}{\partial\alpha}\nabla_{\rho}\alpha\right)2U^{(\mu}\Delta^{\nu)\rho}$
$\displaystyle+\frac{1}{d}(-mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}\frac{1}{c^{2}}\nabla_{\rho}\mathcal{B}_{\sigma}U^{\sigma}U^{\rho}-mc^{2}\frac{1}{d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\nabla_{\rho}\mathcal{B}_{\sigma}\Delta^{\sigma\rho}+\frac{\partial
J_{d+1,0}}{\partial\alpha}U^{\sigma}\nabla_{\sigma}\alpha)\Delta^{\mu\nu}$
$\displaystyle-mc^{2}\frac{2}{(d+2)d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\left(\nabla_{(\rho}\mathcal{B}_{\sigma)}\Delta^{\mu\sigma}\Delta^{\nu\rho}-\frac{1}{d}\nabla_{\rho}\mathcal{B}_{\sigma}\Delta^{\sigma\rho}\Delta^{\mu\nu}\right)\biggl{]}$
$\displaystyle=$
$\displaystyle-c\tau\frac{\mathfrak{g}}{h^{d}}\mathcal{A}_{d-1}(mc)^{d+1}\biggl{[}\frac{1}{c^{2}}\left(\frac{\partial
J_{d-1,2}}{\partial\alpha}U^{\sigma}\nabla_{\sigma}\alpha+mc^{2}\frac{\partial
J_{d-1,3}}{\partial\alpha}U^{\sigma}\nabla_{\sigma}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+1,1}}{\partial\alpha}\beta\nabla_{\rho}U^{\rho}\right)U^{\mu}U^{\nu}$
$\displaystyle+\frac{1}{d}\left(\frac{\partial
J_{d+1,0}}{\partial\alpha}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}\left(\nabla_{\rho}\beta-\frac{1}{c^{2}}\beta
U^{\rho}\nabla_{\rho}U_{\sigma}\right)\right)(\Delta^{\rho\nu}U^{\mu}+\Delta^{\rho\mu}U^{\nu})$
$\displaystyle+\frac{1}{d}\left(\frac{\partial
J_{d+1,0}}{\partial\alpha}U^{\rho}\nabla_{\rho}\alpha+mc^{2}\frac{\partial
J_{d+1,1}}{\partial\alpha}U^{\rho}\nabla_{\rho}\beta-
mc^{2}\frac{1}{d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\beta\nabla_{\rho}U^{\rho}\right)\Delta^{\mu\nu}$
$\displaystyle-mc^{2}\frac{2}{(d+2)d}\frac{\partial
J_{d+3,-1}}{\partial\alpha}\beta\left(\Delta^{\rho\mu}\Delta^{\sigma\nu}\nabla_{(\rho}U_{\sigma)}-\frac{1}{d}\nabla_{\rho}U^{\rho}\Delta^{\mu\nu}\right)\biggl{]}.$
(79)

## Acknowledgement

This work is supported by the National Natural Science Foundation of China
under the grant No. 12275138 and by the Hebei NSF under the Grant No.
A2021205037.

## Data Availability Statement

This work is purely theoretical and contains only analytic analysis. Hence
there is no associated numeric data.

## Declaration of competing interest

The authors declare no competing interest.

## References

  * [1] Z. Banach and S. Piekarski, “Perturbation theory based on the Einstein–Boltzmann system. I. Illustration of the theory for a Robertson–Walker geometry,” J. Math. Phys. 35 no. 9, (1994) 4809–4831.
  * [2] D. Fajman, J. Joudioux, and J. Smulevici, “The stability of the Minkowski space for the Einstein-Vlasov system,” Anal. Part. Diff. Eq. 14 no. 2, (2021) 425–531, [arXiv:1707.06141 [gr-qc]].
  * [3] G. Rein, “Stability and instability results for equilibria of a (relativistic) self-gravitating collisionless gas – A review,” [arXiv:2305.02098 [gr-qc]].
  * [4] P. Rioseco and O. Sarbach, “Accretion of a relativistic, collisionless kinetic gas into a Schwarzschild black hole,” Class. Quant. Grav. 34 no. 9, (2017) 095007, [arXiv:1611.02389 [gr-qc]].
  * [5] A. Cieślik, P. Mach, and A. Odrzywolek, “Accretion of the relativistic Vlasov gas in the equatorial plane of the Kerr black hole,” Phys. Rev. D 106 no. 10, (2022) 104056, [arXiv:2208.04218 [gr-qc]].
  * [6] S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69 (2004) 023503, [arXiv:astro-ph/0306304].
  * [7] C. A. Agon, J. F. Pedraza, and J. Ramos-Caro, “Kinetic theory of collisionless self-gravitating gases: Post-Newtonian polytropes,” Phys. Rev. D 83 (2011) 123007, [arXiv:1104.5262 [gr-qc]].
  * [8] G. V. Vereshchagin and A. G. Aksenov, Relativistic kinetic theory: with applications in astrophysics and cosmology.  Cambridge University Press, 2017.
  * [9] U. W. Heinz, “Kinetic theory for nonabelian plasmas,” Phys. Rev. Lett. 51 (1983) 351.
  * [10] H. T. Elze and U. W. Heinz, “Quark-gluon transport theory,” Phys. Rept. 183 (1989) 81–135.
  * [11] P. B. Arnold, “Quark-gluon plasmas and thermalization,” Int. J. Mod. Phys. E 16 (2007) 2555–2594, [arXiv:0708.0812 [hep-ph]].
  * [12] F. Jüttner, “Das maxwellsche gesetz der geschwindigkeitsverteilung in der relativtheorie,” Annalen der Physik 339 no. 5, (1911) 856–882.
  * [13] R. C. Tolman, “On the weight of heat and thermal equilibrium in general relativity,” Phys. Rev. 35 (1930) 904–924.
  * [14] R. Tolman and P. Ehrenfest, “Temperature equilibrium in a static gravitational field,” Phys. Rev. 36 no. 12, (1930) 1791–1798.
  * [15] G. E. Tauber and J. W. Weinberg, “Internal state of a gravitating gas,” Phys. Rev. 122 no. 4, (1961) 1342–1365.
  * [16] W. Israel, “Relativistic kinetic theory of a simple gas,” J. Math. Phys. 4 no. 9, (1963) 1163–1181.
  * [17] R. W. Lindquist, “Relativistic transport theory,” Annals Phys. 37 no. 3, (1966) 487–518.
  * [18] G. M. Kremer, “Relativistic gas in a Schwarzschild metric,” J. Stat. Mech. 1304 (2013) P04016, [arXiv:1212.5573 [gr-qc]].
  * [19] G. M. Kremer, “Diffusion of relativistic gas mixtures in gravitational fields,” Physica A 393 (2014) 76–85, [arXiv:1303.6463 [gr-qc]].
  * [20] S. Liu, X. Hao, S.-F. Liu, and L. Zhao, “Covariant transport equation and gravito-conductivity in generic stationary spacetimes,” Eur. Phys. J. C 82 no. 12, (2022) 1100, [arXiv:2210.10907 [gr-qc]].
  * [21] X. Hao, S. Liu, and L. Zhao, “Gravito-thermal transports, Onsager reciprocal relation and gravitational Wiedemann-Franz law,” [arXiv:2306.04545 [cond-mat.stat-mech]].
  * [22] P. Kovtun, “Lectures on hydrodynamic fluctuations in relativistic theories,” J. Phys. A 45 (2012) 473001, [arXiv:1205.5040 [hep-th]].
  * [23] O. Sarbach and T. Zannias, “Relativistic kinetic theory: An introduction,” AIP Conf. Proc. 1548 no. 1, (2013) 134–155, [arXiv:1303.2899 [gr-qc]].
  * [24] O. Sarbach and T. Zannias, “The geometry of the tangent bundle and the relativistic kinetic theory of gases,” Class. Quant. Grav. 31 (2014) 085013, [arXiv:1309.2036 [gr-qc]].
  * [25] C. Cercignani and G. M. Kremer, The relativistic Boltzmann equation: Theory and applications.  Springer, 2002.
  * [26] P. Kovtun, “Temperature in relativistic fluids,” Phys. Rev. D 107 no. 8, (2023) 086012, [arXiv:2210.15605 [gr-qc]].
  * [27] A. Santos, J. J. Brey, and J. W. Dufty, “Divergence of the Chapman-Enskog expansion,” Phys. Rev. Lett. 56 (1986) 1571–1574.
  * [28] G. S. Denicol and J. Noronha, “Divergence of the Chapman-Enskog expansion in relativistic kinetic theory,” [arXiv:1608.07869 [nucl-th]].
  * [29] N.-X. Chen and B.-H. Sun, “Note on divergence of the Chapman–Enskog expansion for solving Boltzmann equation,” Chinese Physics Letters 34 no. 2, (2017) 020502, [arXiv:1808.10739 [physics.gen-ph]].
  * [30] L. Onsager, “Reciprocal relations in irreversible processes. I.,” Phys. Rev. 37 no. 4, (1931) 405–426.
  * [31] L. Onsager, “Reciprocal relations in irreversible processes. II.,” Phys. Rev. 38 no. 12, (1931) 2265–2279.
  * [32] F. S. Bemfica, M. M. Disconzi, C. Rodriguez, and Y. Shao, “Local well-posedness in Sobolev spaces for first-order conformal causal relativistic viscous hydrodynamics,” [arXiv:1911.02504 [math.AP]].
  * [33] P. Kovtun, “First-order relativistic hydrodynamics is stable,” JHEP 10 (2019) 034, [arXiv:1907.08191 [hep-th]].