#  A new approach to the 3-momentum regularization of the in-medium one and
two fermion line integrals with applications to cross sections in the
Nambu–Jona-Lasinio model

Renan Câmara Pereira<EMAIL_ADDRESS>CFisUC - Center for Physics
of the University of Coimbra, Department of Physics, University of Coimbra,
3004-516 Coimbra, Portugal  João Moreira<EMAIL_ADDRESS>CFisUC - Center for
Physics of the University of Coimbra, Department of Physics, University of
Coimbra, 3004-516 Coimbra, Portugal  Departamento de Física, Universidade da
Beira Interior, 6201-001 Covilhã, Portugal  Pedro Costa<EMAIL_ADDRESS>CFisUC -
Center for Physics of the University of Coimbra, Department of Physics,
University of Coimbra, 3004-516 Coimbra, Portugal  Constança Providência
<EMAIL_ADDRESS>CFisUC - Center for Physics of the University of Coimbra, Department
of Physics, University of Coimbra, 3004-516 Coimbra, Portugal

###### Abstract

We propose the 3-momentum sphere intersection regularization applied to the
one and two fermion line integrals at finite temperature and chemical
potential. The quark-antiquark polarization function in this new
regularization approach is equivalent to the usual 3-momentum regularization,
when the absolute value of the external 3-momentum of the polarization is
zero. Additionally, it respects the particle-antiparticle symmetry of meson
states in the Nambu$-$Jona-Lasinio (NJL) model for all values of temperature
and chemical potential. Without this symmetry, in-medium cross sections
calculated in the 3-momentum regularized NJL model are not consistent. In
order to demonstrate the difference between the usual 3-momentum
regularization with the one proposed in this work, we study the quark-quark
and quark-antiquark cross sections in both regularization schemes. To this end
we use the standard $SU(3)$ NJL model, with four and six quark interactions.
We observe major quantitative and qualitative differences when comparing
quark-quark cross sections in both schemes. The quark-antiquark cross
sections, on the other hand, are very similar in both regularizations, owning
to the equivalence between the regularizations when the absolute value of the
external 3-momentum is zero.

##  I Introduction

The calculation of one loop integrals, at finite temperature and chemical
potential, is an essential part of the study of in-medium quantum field
theory. These integrals arise naturally in calculations performed under the
framework of effective theories of quantum chromodynamics (QCD), like
Nambu$-$Jona-Lasinio (NJL) type models [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. In the
NJL model, the one and two fermion line integrals occupy a central role in the
calculation of meson propagators, meson masses, cross sections, or
thermodynamic quantities, like pressure, energy density and entropy density.
Moreover, quark-quark and quark-antiquark cross sections, derived from the NJL
model, can be used to estimate the quark relaxation times, a vital part of
several theoretical predictions about the in-medium behavior of transport
coefficients like the shear and bulk viscosities, as well as, the electrical
and thermal conductivities of strongly interacting matter [11, 12, 13, 14, 15,
16].

Due to the non-perturbative characteristics of QCD at low energies,
phenomenological models, like the NJL model, become an essential tool to make
qualitative predictions about QCD at finite temperature and density, specially
when the goal is to study the behavior of the quark gluon plasma at
temperatures and chemical potentials for which the ab initio approaches of
QCD, like lattice gauge theory or perturbative based calculations, are not
applicable [17, 18, 19, 20]. Hence, making sure that the calculation of one
loop integrals, within the scope of NJL-type models, is made in a symmetry
preserving manner, is extremely important in order to make better qualitative
and quantitative predictions while also assuring that one is not introducing
non-physical features in the phenomenological calculation.

In this work we follow the terminology given in Ref. [11], and denote the one
fermion line integral by $A$ and the two fermion line integral by $B_{0}$.
Here, we will not study the practical implementation of integrals with a
fermion line counting higher then two. As it is well known, both of these
integrals are divergent and some regularization technique must be employed,
becoming an essential part of the calculations [11, 21]. When studying finite
temperature and chemical potential systems, the 3-momentum regularization is
widely applied [11]. In this framework the integration over the internal
3-momentum of quarks, $\bm{p}$, is carried out up to a maximum value of its
absolute value, $\absolutevalue{\bm{p}}<\Lambda$, with $\Lambda$ the
3-momentum cutoff of the model [21]. The usual way in which the 3-momentum
regularization is performed breaks a symmetry of the system: the quark-
antiquark polarization loop, which mixes light and strange quarks, is not
symmetric with respect to the interchange of quarks [12]. This originates from
a shift of variables used in order to simplify the evaluation of the $B_{0}$
integral [11, 12].

One of the main issues with this feature is that physical scattering processes
which ought to be identical for all values of temperature and chemical
potential, become different. A clear example is the cross sections of the
quark-quark processes, $us\to us$ and $su\to su$ (the same occurs for $ds\to
ds$ and $sd\to sd$). Considering the $us\to us$ scattering in the $u$-channel,
the exchanged mesons are the pseudoscalar $K^{+}$ and the scalar $\kappa^{+}$
while, in the $t$-channel, the exchanged mesons are the flavorless, neutral
mesons $\eta$, $\eta^{\prime}$, $\pi_{0}$, $\sigma$, $f_{0}$ and $a_{0}^{0}$.
For the $su\to su$ process the corresponding antimesons are involved in each
channel. Using a theoretical setup in which the quark-quark processes, $us\to
us$ and $su\to su$, yield different results is an indication that the set of
meson propagators used in the calculations do not display the expected
particle-antiparticle symmetry of meson states. One way to solve this problem
is to follow the recipe introduced in Ref. [12] and, when appropriate, replace
the two fermion line integral present in the quark-antiquark polarization loop
function,
$B_{0}=B_{0}\quantity[M_{u},M_{s},T,\mu_{u},\mu_{s},k_{0},\absolutevalue{\bm{k}}]$,
with an average,
$B_{0}=\frac{1}{2}\quantity(B_{0}\quantity[M_{u},M_{s},T,\mu_{u},\mu_{s},k_{0},\absolutevalue{\bm{k}}]+B_{0}\quantity[M_{s},M_{u},T,\mu_{s},\mu_{u},k_{0},\absolutevalue{\bm{k}}])$.
Here $M_{f}$ and $\mu_{f}$ are the effective mass and effective chemical
potential of the quark of flavor $f$, $T$ is the temperature and $k_{0}$,
$\absolutevalue{\bm{k}}$ refer to the temporal-like and spacial-like parts of
the external momentum, respectively. In the physical scenario where isospin
symmetry does not occur (up and down quarks have different bare quark masses)
or if there is an asymmetry caused by non-degenerate up and down quark
chemical potentials (for example, when considering matter in beta
equilibrium), several other scattering processes are also affected, such as
the $ud\to ud$ and $du\to du$ processes.

In this work we propose a novel approach to the 3-momentum regularization,
namely the 3-momentum sphere intersection regularization scheme. In this new
regularization scheme, we regularize the one loop integrals based on the
number of quark propagators present in the integrand. Each quark propagator
has a 3-momentum which must be smaller than the 3-momentum cutoff,
$\absolutevalue{\bm{p}}<\Lambda$. In the 3-momentum sphere intersection
regularization, the individual momentum of each quark propagators must be
smaller than the cutoff, even when external momenta are considered. As we will
discuss in this paper, this can be achieved by considering that the
integration over the momentum in the $N$ fermion line integrals must be
evaluated inside the region defined by the intersection of $N$ spheres.
Furthermore, in the case in which there are no external momenta applied to the
$N$ fermion line integrals, the spheres completely intersect and the
regularization is equivalent to the usual 3-momentum regularization.
Interestingly, if the spheres do not intersect, the integral is automatically
zero. As we will discuss, this happens for the meson propagators: if the
momentum of the meson is larger than $2\Lambda$, the spheres corresponding to
each quark propagator, inside the quark-antiquark polarization function, no
longer intersect and the polarization function is automatically zero. With
this approach to the 3-momentum regularization, all the symmetries of the
system are maintained and the previously discussed problems with the
scattering processes do not occur.

This paper is organized as follows. In Section II, we derive the one and two
fermion line integrals from the NJL model for any number of scalar and
pseudoscalar quark-quark interactions. Such is accomplished in the mean field
approximation by writing the linear and quadratic expansions of the NJL
lagrangian. We also show how the expansion used to derive the meson
propagators from the quadratic NJL model yields the same meson propagators as
the ones obtained via the commonly employed Bethe-Salpeter equation in the
random phase approximation. Additionally, the 3-momentum sphere intersection
is introduced and discussed. The formalism developed here can be used to
incorporate multi-quark interactions, such as the eight quark proposed in [22,
23], or explicitly chiral symmetry breaking interactions, as those in [24,
25]. In Sections III and IV, we evaluate the two fermion line integral,
$B_{0}$, and the one fermion line integral, $A$, within the new regularization
scheme. In Section V the formalism to evaluate quark-quark and quark-antiquark
cross sections is laid out and a specific version of the $SU(3)$ NJL model,
with four and six quark interactions is introduced in order to study the
effect of using both regularizations on the different scattering processes.
The numerical results regarding quark-quark and quark-antiquark scattering
processes for different scenarios of temperature and chemical potential are
presented and discussed. Finally, in Section VI, conclusions are formulated
and additional research is suggested.

##  II The one and two fermion line integrals from NJL-type models

The one fermion line integral, $A$, and the two fermion line integral,
$B_{0}$, at one loop level, arise when studying NJL-type models. In this
section we show one way to derive the definitions of these integrals from the
NJL model and discuss regularization strategies, namely a new way to
regularize the quark-antiquark polarization function.

Consider a general model describing $N_{f}$ flavors of fermions including
different types of scalar and pseudoscalar interactions and the standard free
Dirac Lagrangian. This can be written as:

$\displaystyle\mathcal{L}[\overline{\psi},\psi]=\overline{\psi}\quantity(i\not{\partial}-m)\psi+\mathcal{L}_{\mathrm{int}}\quantity[s_{a},p_{a}].$
(1)

Here, $m=\operatorname{diag}\quantity{m_{0},m_{1},\ldots,m_{N_{f}-1}}$ is the
quark current mass matrix, the term $\mathcal{L}_{\mathrm{int}}$ contains
several types of quark-quark interactions and is written as a function of
quark bilinear operators, $s_{a}=\overline{\psi}\lambda_{a}\psi$ and
$p_{a}=\overline{\psi}i\gamma^{5}\lambda_{a}\psi$, where $\lambda_{a}$ are the
$N_{f}^{2}$ matrices that span the U$(N_{f})$ algebra. The usual 4-quark
scalar-pseudoscalar interaction, for example, can be written as,
$\mathcal{L}_{4q}\propto(\overline{\psi}\lambda_{a}\psi)^{2}+(\overline{\psi}i\gamma_{5}\lambda_{a}\psi)^{2}=s_{a}^{2}+p_{a}^{2}$
(Einstein summation convention is used).

Using this Lagrangian density, one can obtain the generating functional of the
theory and derive all quantities of interest like the quark propagator, the
thermodynamic potential and (anti)quark-(anti)quark cross sections (here and
henceforth we are referring to quark-quark, quark-antiquark and antiquark-
antiquark processes). However, the presence of more than two quark-quark
interactions at the Lagrangian level renders impossible the exact integration
over the quark fields [26]. This is a natural difficulty when dealing with
interacting field theories and some approximation must be employed in order to
deal with the functional integrations and calculate the previously mentioned
quantities.

In this work, we will use the linear and quadratic expansions of the
Lagrangian, allied to the mean field approximation. The one fermion line
integral, $A$, arises in the linear expansion of the NJL Lagrangian, namely in
the calculation of the quark propagator and thermodynamic quantities. The two
fermion line integral, $B_{0}$, emerges in the quadratic expansion and is
essential in the calculation of the quark-antiquark polarization function,
which is a fundamental step in the calculation of meson mass behavior and the
cross sections. In the following we show how these integrals arise from the
linear and quadratic expansions of the Lagrangian defined in Eq. (1).

###  II.1 The linear Lagrangian

As already discussed, in order to obtain the quark propagator and
thermodynamics of the model defined in Eq. (1), one can consider the linear
expansion of the Lagrangian, allied to the mean field approximation. In this
approximation one transforms all quark-quark interactions into quadratic quark
interactions (linear in bilinear quark operators) by the introduction of
auxiliary fields. Terms which have a larger than linear dependency in the
quantum fluctuations are neglected. Using this technique, one can write a
generating functional integrand that is quadratic in the quark fields.

Consider a quark bilinear operator,
$\hat{\mathcal{O}}=\overline{\psi}\hat{\Gamma}\psi$, where $\hat{\Gamma}$ is
an operator that acts on some internal quark index, like color or flavor.
Without loss of generality, the bilinear can be written as a sum between its
mean field value and a small fluctuation,
$\hat{\mathcal{O}}=\expectationvalue*{\hat{\mathcal{O}}}+\delta\hat{\mathcal{O}}$.
In order to obtain the linear expansion of the product between $N-$bilinear
operators, one only keeps terms that are linear in the fluctuations,
$\delta\hat{\mathcal{O}}$, i.e., terms equal or superior to
$(\delta\hat{\mathcal{O}})^{2}$ are neglected. As a matter of fact, one can
prove that the linear product between $N=n+1$ commuting operators can be
written using the following formula (see Appendix A for the proof) [27, 28]:

$\displaystyle\mathrm{L}\left[\prod_{i=1}^{n+1}\hat{\mathcal{O}}_{i}\right]$
$\displaystyle=\quantity(\sum_{i=1}^{n+1}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}-n)\prod_{j=1}^{n+1}\expectationvalue*{\hat{\mathcal{O}}_{j}}.$
(2)

As an example, the linear product between two and three commuting operators,
can directly be obtained from this equation:

$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]$
$\displaystyle=-\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}+\expectationvalue*{\hat{\mathcal{O}}_{2}}\hat{\mathcal{O}}_{1}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\hat{\mathcal{O}}_{2},$
(3)
$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]$
$\displaystyle=-2\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}+\hat{\mathcal{O}}_{1}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\hat{\mathcal{O}}_{2}\expectationvalue*{\hat{\mathcal{O}}_{3}}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\hat{\mathcal{O}}_{3}.$
(4)

Finally, using this approximation, one can obtain a Lagrangian density,
$\mathcal{L}_{\mathrm{L}}$, that is linear in the quark bilinear operators,

$\displaystyle\mathcal{L}_{\mathrm{L}}\quantity[\overline{\psi},\psi,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]=\overline{\psi}\quantity(i\not{\partial}-M\quantity[m,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}])\psi-\mathcal{U}\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}].$
(5)

Here,
$\expectationvalue*{\bm{s}}=\\{\expectationvalue*{s_{0}},\ldots,\expectationvalue*{s_{N_{f}^{2}-1}}\\}$
and
$\expectationvalue*{\bm{p}}=\\{\expectationvalue*{p_{0}},\ldots,\expectationvalue*{p_{N_{f}^{2}-1}}\\}$
are the sets of scalar and pseudoscalar mean fields obtained after performing
the linearization procedure. These fields correspond to the expectation values
of the $N_{f}^{2}$ scalar and pseudoscalar quark bilinear operators,
$s_{a}=\overline{\psi}\lambda_{a}\psi$ and
$p_{a}=\overline{\psi}i\gamma^{5}\lambda_{a}\psi$, respectively
($a=0,\ldots,N_{f}^{2}-1$). The effective mass matrix, $M$ and the mean field
potential $U$ are functions of these sets. Using this quantity, one can find
the generating functional of the model, within this approximation, to be given
by:

$\displaystyle\mathcal{Z}_{\mathrm{L}}\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}},\overline{\eta},\eta]$
$\displaystyle\propto\int\mathcal{D}\overline{\psi}\mathcal{D}\psi\exp\quantity[i\int\differential[4]{x}\quantity(\mathcal{L}_{\mathrm{L}}\quantity[\overline{\psi},\psi,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]+\overline{\psi}\eta+\overline{\eta}\psi)]$
$\displaystyle\propto\exp\quantity[-i\mathcal{U}\,V]\int\mathcal{D}\overline{\psi}\mathcal{D}\psi\exp\quantity[i\int\differential[4]{x}\quantity(\int\differential[4]{y}\overline{\psi}[x]\,\mathcal{G}_{q}^{-1}\quantity[x,y]\,\psi[y]+\overline{\psi}[x]\eta+\overline{\eta}\psi[x])],$
(6)

where $\eta$ and $\overline{\eta}$ are the quark and antiquark sources, $V$ is
the four-dimensional volume and the inverse quark propagator,
$\mathcal{G}_{q}^{-1}$, was defined as:

$\displaystyle\mathcal{G}_{q}^{-1}\quantity[x,y]$
$\displaystyle=\quantity(i\not{\partial}-M)\delta^{4}\quantity(x-y).$ (7)

For simplicity, we omitted the mean field dependencies present on the quark
effective mass, $M=M[m,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$.

The quark fields can then be integrated out, and one can define the effective
action of the model from the generating functional (with vanishing sources),
$\Gamma\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]=-i\ln\quantity[\mathcal{Z}_{\mathrm{L}}\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]]$,
which only depends on the sets of mean fields, $\expectationvalue{\bm{s}}$ and
$\expectationvalue{\bm{p}}$. Such dependence must be removed by requiring the
effective action to be an extreme with respect to these mean fields i.e.,
$\functionalderivative*{\Gamma\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]}{\expectationvalue{s_{a}}}=\functionalderivative*{\Gamma\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]}{\expectationvalue{p_{a}}}=0$.
At this point one has access to the thermodynamic potential from which several
thermodynamic variables of interest, such as pressure, energy density, entropy
and speed of sound, can be evaluated.

Finally, the one fermion, one loop integral, $A$, can be defined by
considering the quark condensate of flavor $i$ [18]:

$\displaystyle\expectationvalue*{\overline{\psi}_{i}\psi_{i}}$
$\displaystyle=-i\tr\mathcal{G}_{q}^{i}\quantity[x,x]$
$\displaystyle=-4iN_{c}M_{i}\int\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\frac{1}{p^{2}-M_{i}^{2}}$
$\displaystyle=\frac{N_{c}}{4\pi^{2}}M_{i}A\quantity[M_{i}].$ (8)

The one fermion, one loop integral in Minkowski spacetime,
$A\quantity[M_{f}]$, is given by:

$\displaystyle
A\quantity[M_{i}]=-16\pi^{2}i\int_{\mathrm{Reg}}\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\frac{1}{p^{2}-M_{i}^{2}}.$
(9)

Here, $M_{i}$ is the $i$-quark effective mass. The inclusion of temperature
and chemical potential will be done in Section IV. Since this integral is
divergent and some regularization scheme must be applied, we explicitly wrote
the label $\mathrm{Reg}$.

###  II.2 The quadratic Lagrangian

The two fermion line integral, $B_{0}$, can be found when considering the
quadratic expansion of the Lagrangian given in Eq. (1). Such expansion allows
one to calculate several interesting properties, namely: meson propagators,
meson masses, decay rates and (anti)quark-(anti)quark cross sections.

The approach to write a Lagrangian that is quadratic in the quark bilinear
operators (fourth order in the quark fields) is similar to the one used in the
previous section to obtain the linear Lagrangian. Firstly, as before, the
operator is written as the sum,
$\hat{\mathcal{O}}=\expectationvalue*{\hat{\mathcal{O}}}+\delta\hat{\mathcal{O}}$.
However, to get the quadratic expansion of the product between $N$-bilinear
operators, only up to quadratic terms in the fluctuations are kept and higher
order terms are neglected. As before, instead of doing such expansion by hand,
there is a formula from which one can readily obtain the quadratic product
between $N=n+2$ commuting operators, with $n\geq 1$ (see Appendix B for the
proof). It is given by [27, 28]:

$\displaystyle\mathrm{Q}\left[\prod_{i=1}^{n+2}\hat{\mathcal{O}}_{i}\right]$
$\displaystyle=\bigg{(}\frac{1}{2}\sum_{i=1}^{n+2}\sum_{j=1}^{n+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}\frac{\hat{\mathcal{O}}_{j}}{\expectationvalue*{\hat{\mathcal{O}}_{j}}}\quantity(1-\delta_{ij})-n\sum_{i=1}^{n+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}+\frac{n}{2}\quantity(n+1)\bigg{)}\prod_{k=1}^{n+2}\expectationvalue*{\hat{\mathcal{O}}_{k}}.$
(10)

As an example, the quadratic expansion of the product between three and four
operators, are given by:

$\displaystyle\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]=$
$\displaystyle-\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{1}-\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{2}-\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\hat{\mathcal{O}}_{3}$
$\displaystyle+\expectationvalue*{\hat{\mathcal{O}}_{1}}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}+\expectationvalue*{\hat{\mathcal{O}}_{2}}\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{3}+\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}},$
(11)
$\displaystyle\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\hat{\mathcal{O}}_{4}\right]=$
$\displaystyle-2\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{1}-2\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{2}$
$\displaystyle-2\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{3}-2\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{4}$
$\displaystyle+\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{3}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}+\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{4}$
$\displaystyle+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{4}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\hat{\mathcal{O}}_{3}\hat{\mathcal{O}}_{4}+\expectationvalue*{\hat{\mathcal{O}}_{3}}\expectationvalue*{\hat{\mathcal{O}}_{4}}\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}$
$\displaystyle+3\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\expectationvalue*{\hat{\mathcal{O}}_{4}}.$
(12)

Using the quadratic expansion above and neglecting second order terms which
mix scalar and pseudoscalar quark bilinear operators ($s_{a}p_{b}$ and
$p_{a}s_{b}$)111This follows from the consideration that only the scalar
condensates get a non-vanishing ground state and pseudoscalar condensates are
zero., the quadratic Lagrangian density, $\mathcal{L}_{\mathrm{Q}}$, (fourth
order in the quark fields) can be written as:

$\displaystyle\mathcal{L}_{\mathrm{Q}}\quantity[\overline{\psi},\psi,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]=\overline{\psi}\quantity(i\not{\partial}-m)\psi+s_{a}S_{ab}s_{b}+p_{a}P_{ab}p_{b}+s_{a}S_{a}+p_{a}P_{a}+R.$
(13)

Here, $S_{ab}=S_{ab}[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$ and
$P_{ab}=P_{ab}[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$ gather
the coefficients that arise after the quadratic expansion of the operators,
which are quadratic in the variables $s_{a}=\overline{\psi}\lambda_{a}\psi$
and $p_{a}=\overline{\psi}i\gamma^{5}\lambda_{a}\psi$, respectively. They are
the so-called meson projectors for the scalar modes $(S_{ab})$ and the
pseudoscalar modes $(P_{ab})$. The first and zeroth order contributions are
gathered in the coefficients
$S_{a}=S_{a}[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$,
$P_{a}=P_{a}[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$ and
$R=R[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$, respectively.

The generating functional for the system, can be written as:

$\displaystyle\mathcal{Z}_{\mathrm{Q}}\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}},\eta,\overline{\eta}]\propto\int\mathcal{D}\overline{\psi}\mathcal{D}\psi\exp\quantity{i\int\differential[4]{x}\quantity(\mathcal{L}_{\mathrm{Q}}\quantity[\overline{\psi},\psi,\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]+\overline{\psi}\eta+\overline{\eta}\psi)}.$
(14)

In order to integrate out the quark fields, one can bosonize the model by
introducing auxiliary meson fields, $\sigma_{a}$ and $\pi_{a}$, through
Hubbard–Stratonovich transformations [29, 30, 31]. For the $s_{a}$
contributions, consider the following Gaussian integral:

$\displaystyle\exp\quantity{i\int\differential[4]{x}\quantity(s_{a}S_{ab}s_{b})}$
$\displaystyle\propto\int\prod_{c}\mathcal{D}\sigma_{c}\exp\quantity{i\int\differential[4]{x}\quantity(-\frac{1}{4}\sigma_{a}\quantity(S^{-1})_{ab}\sigma_{b}-s_{a}\sigma_{a})}.$
(15)

Considering a change of variables on the functional integral as:
$\sigma_{a}\to\sigma_{a}+\quantity(S_{a}+\Delta_{a})$ where
$\Delta_{a}=M_{a}-m_{a}$ i.e., $\lambda_{a}\Delta_{a}=M-m$, and multiplying
both sides by $\exp\quantity{i\int\differential[4]{x}\quantity(s_{a}S_{a})}$,
we can write222Here we used the fact that the matrix $S$ is symmetric.:

$\displaystyle\exp\\{i\int\differential[4]{x}\quantity(s_{a}S_{ab}s_{b}+s_{a}S_{a})\bigg{\\}}$
$\displaystyle\propto\int\prod_{c}\mathcal{D}\sigma_{c}\exp\\{$ $\displaystyle
i\int\differential[4]{x}\bigg{(}-\frac{1}{4}\sigma_{a}\quantity(S^{-1})_{ab}\sigma_{b}-\frac{1}{2}\sigma_{a}\quantity(S^{-1})_{ab}\quantity(S_{b}+\Delta_{b})-s_{a}\quantity(\sigma_{a}+\Delta_{a})\bigg{)}\bigg{\\}}.$
(16)

With this particular change of variables we are considering the possibility of
scalar condensation in the vacuum. In this work we do not consider the
possibility of pseudoscalar condensation. In this case, for the $p_{a}$
contributions, the approach is similar except that the change of variables is
simply $\pi_{a}\to\pi_{a}+P_{a}$. Using these identities, the action becomes
quadratic in the fermion fields and one can integrate out the quark fields to
yield:

$\displaystyle\mathcal{Z}_{\mathrm{Q}}\quantity[\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$
$\displaystyle\propto\int\prod_{a}\mathcal{D}\sigma_{a}\mathcal{D}\pi_{a}\exp\quantity{i\mathcal{S}_{\mathrm{Q}}\quantity[\bm{\sigma},\bm{\pi},\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]}.$
(17)

Here,
$\mathcal{S}_{\mathrm{Q}}\quantity[\bm{\sigma},\bm{\pi},\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]$
is the action of the system and
$\bm{\sigma}=\\{\sigma_{0},\ldots,\sigma_{N_{f}^{2}-1}\\}$ and
$\bm{\pi}=\\{\pi_{0},\ldots,\pi_{N_{f}^{2}-1}\\}$ are the sets of scalar and
pseudoscalar auxiliary fields (respectively) introduced with the
Hubbard–Stratonovich transformation. The action of the system is explicitly
given by:

$\displaystyle\mathcal{S}_{\mathrm{Q}}\quantity[\bm{\sigma},\bm{\pi},\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]=-i\Tr\ln[i\not{\partial}-M$
$\displaystyle-\lambda_{a}\sigma_{a}-i\gamma_{5}\lambda_{a}\pi_{a}\big{]}$
$\displaystyle+\int\differential[4]{x}\bigg{(}$
$\displaystyle-\frac{1}{4}\sigma_{a}\quantity(S^{-1})_{ab}\sigma_{b}-\frac{1}{4}\pi_{a}\quantity(P^{-1})_{ab}\pi_{b}$
$\displaystyle-\frac{1}{2}\sigma_{a}\quantity(S^{-1})_{ab}\quantity(S_{b}+\Delta_{b})-\frac{1}{2}\pi_{a}\quantity(P^{-1})_{ab}P_{b}+R\bigg{)}.$
(18)

In the above expressions, $\Tr$ is the usual notation for the total trace over
continuous spacetime indexes and discrete internal indexes such as color,
flavor and spinor,
$\Tr\equiv\operatorname{{\tr}_{\mathnormal{x}}}\operatorname{{\tr}_{\mathnormal{c}}}\operatorname{{\tr}_{\mathnormal{f}}}\operatorname{{\tr}_{\mathrm{D}}}$
and the functional trace of an operator, $\hat{O}$, is defined by:
$\operatorname{{\tr}_{\mathnormal{x}}}\hat{O}=\int\differential[d]{x}\bra{x}\hat{O}\ket{x}$
[32, 18]. Here, we will use the functional formalism introduced in Refs. [33,
34, 35], where $\hat{O}$ is considered an operator acting directly on an
fictitious Hilbert space, $\quantity{\ket{x}}$, with the position
($\hat{x}_{\mu}$) and momentum operators ($\hat{p}_{\mu}$) obeying the usual
commutation relation,
$\quantity[\hat{x}_{\mu},\hat{p}_{\nu}]=i\delta_{\mu\nu}$ [33, 34, 35, 32].

The gap equations of the model can be obtained by evaluating the path integral
over the auxiliary meson fields, $\sigma_{a}$ and $\pi_{a}$ in the stationary
phase (or saddle point) approximation [17, 19, 31, 21, 36], which is
equivalent to the mean field or Hartree approximation. It amounts to requiring
the action of the model to be stationary with respect to the auxiliary meson
fields,
$\left.\functionalderivative*{\mathcal{S}_{\mathrm{Q}}}{\sigma_{a}}\right|_{\quantity{\sigma^{\mathrm{st}}_{a},\pi^{\mathrm{st}}_{a}}}=\left.\functionalderivative*{\mathcal{S}_{\mathrm{Q}}}{\pi_{a}}\right|_{\quantity{\sigma^{\mathrm{st}}_{a},\pi^{\mathrm{st}}_{a}}}=0$.
Here, we note that we already allowed for a non-vanishing expectation value of
the field $\sigma_{a}$ in the change of variables,
$\sigma_{a}\to\sigma_{a}+\quantity(S_{a}+\Delta_{a})$, and thus are requiring
that the stationary field configurations must vanish:
$\sigma^{\mathrm{st}}_{a}=\pi^{\mathrm{st}}_{a}=0$. Using these requirements,
one can derive the self-consistent gap equations of the model which relates
the quarks effective mass, $M$, and the expectation values of the quark
bilinear operators, $\expectationvalue{\bm{s}}$ and
$\expectationvalue{\bm{p}}$. This approach leads exactly to the same gap
equations as the ones obtained using the linear Lagrangian in the mean field
approximation presented in the previous section. Indeed, for a specific NJL
Lagrangian, one can check this equivalence explicitly. This correspondence
might be related the exact Hubbard–Stratonovich transformation, used to
transform the quadratic action (in terms of quark bilinear operators) into a
linear action. A proper demonstration, between the equivalence of the gap
equations derived in the linear and quadratic Lagrangians (alongside the
saddle point approximation for the auxiliary meson fields), for any
interactions terms, is left as future work.

The quark-antiquark polarization function, at finite chemical potential, can
be derived from this formalism by considering the expansion of the fermionic
determinant to second order in the mesonic auxiliary fields, $\sigma_{a}$ and
$\pi_{a}$. Finite chemical potential can be included in the calculation by
making the following substitution, $\partial_{\alpha}\to
D_{\alpha}=\partial_{\alpha}-i\mu\delta_{\alpha 0}$, with $\mu$ the quark
chemical potential matrix in flavor space [18]. This matrix can be considered
to be diagonal with components,
$\mu=\operatorname{diag}\quantity(\mu_{0},\mu_{1},\ldots,\mu_{N_{f}-1})$.
Next, we consider the Taylor series expansion of the logarithm of a matrix in
Eq. (18). We write [18, 37, 38, 2, 39, 40, 41]:

$\displaystyle\Tr\ln[i\not{D}-M-\lambda_{a}\sigma_{a}-i\gamma_{5}\lambda_{a}\pi_{a}\big{]}$
$\displaystyle=\Tr\ln[i\not{D}-M\big{]}+\Tr\ln\quantity[1+\quantity[i\not{D}-M]^{-1}\bigg{(}-\lambda_{a}\sigma_{a}-i\gamma_{5}\lambda_{a}\pi_{a}\bigg{)}]$
$\displaystyle=\Tr\ln[i\not{D}-M\big{]}-\sum_{n=1}^{\infty}\Tr\quantity[T_{n}].$
(19)

In the above, the logarithm of a matrix was written using its defining Taylor
series expansion (assuming that the series converges) [41]:

$\displaystyle\ln\quantity[X]$
$\displaystyle=\sum_{n=1}^{\infty}\frac{\quantity(-1)^{n+1}}{n}\quantity(X-I)^{n}.$
(20)

A convergence criterion can be written as: $\absolutevalue{X-I}<1$, with
$\absolutevalue{Y}$ a proper way to evaluate the norm of the matrix $Y$ [41].
Additionally, the operator $T_{n}$ is defined as:

$\displaystyle
T_{n}=\frac{1}{n}\quantity(\quantity[i\not{D}-M]^{-1}\quantity(\lambda_{a}\sigma_{a}+i\gamma_{5}\lambda_{a}\pi_{a}))^{n}.$
(21)

We are interested in evaluating the second order term in the expansion,
$T_{2}$, which is an essential part in the definition of the inverse meson
propagators. To this end, we define the operators,
$\mathcal{G}=\quantity[i\not{D}-M]^{-1}$, $\sigma=\lambda_{a}\sigma_{a}$ and
$\pi=\lambda_{a}\pi_{a}$. We highlight that $\mathcal{G}$ is a matrix in
flavor space due to the structure of the effective quark mass and effective
chemical potential matrices,
$M=\operatorname{diag}\quantity(M_{0},M_{1},\dots,M_{N_{f}-1})$, and
$\mu=\operatorname{diag}\quantity(\mu_{0},\mu_{1},\ldots,\mu_{N_{f}-1})$. With
this notation, we can write:

$\displaystyle T_{2}$
$\displaystyle=\frac{1}{2}\mathcal{G}\quantity(\sigma+i\gamma_{5}\pi)\mathcal{G}\quantity(\sigma+i\gamma_{5}\pi)$
$\displaystyle=T_{2}^{\sigma\sigma}+T_{2}^{\sigma\pi}+T_{2}^{\pi\sigma}+T_{2}^{\pi\pi}.$
(22)

Only the traces of the first and fourth operators are nonzero,
$T_{2}^{\sigma\sigma}$ and $T_{2}^{\pi\pi}$. The trace over spinor indexes of
the “mixed” contributions involving the identity matrix $I_{4}$ and the
$\gamma^{5}$ matrix vanish. Without loss of generality, we can deal with the
quantity,
$T_{2}^{\phi\phi}\quantity[\Gamma_{1},\Gamma_{2}]=\frac{1}{2}\mathcal{G}\Gamma_{1}\phi\mathcal{G}\Gamma_{2}\phi$,
with $\phi$ standing in for $\sigma$ or $\pi$. It yields
$T_{2}^{\sigma\sigma}$ when $\Gamma_{1}=\Gamma_{2}=I_{4}$ ($I_{4}$ is the
four-dimensional unity matrix) and $T_{2}^{\pi\pi}$ in the case in which
$\Gamma_{1}=\Gamma_{2}=i\gamma_{5}$. The operator $\Gamma$ only has discrete
Dirac indexes. The fields $\phi$ and the quantity $\mathcal{G}$ both have
discrete flavor indexes as well as continuous spacetime indexes but are
structureless in color space. Additionally, the operator $\mathcal{G}$ also
carries Dirac indexes. As a result of this internal structure, the evaluation
of the traces over these different vector spaces has to be done taking into
account the noncommutative nature of these operators, i.e., one has to take
into account that, $\quantity[\mathcal{G},\Gamma]\neq 0$ in spinor space,
$\quantity[\mathcal{G},\phi]\neq 0$ both for flavor and continuous spacetime
indexes, while the remaining commutation relations are zero (for instance,
$\quantity[\Gamma,\phi]=0$ in flavor and spinor spaces). The way in which one
evaluates the traces can lead to different mathematical expansions. As we will
explain more thoroughly later, we are interested in obtaining an expansion of
the fermionic determinant which yields the same result for the meson
propagators as the ones obtained via Bethe-Salpeter equation in the random
phase approximation (or, equivalently, the ring approximation [19]). We will
make further comments about deriving the meson propagators from the fermionic
determinant in Section II.2.2.

Lets start by evaluating the trace over flavor degrees of freedom,
$\operatorname{{\tr}_{\mathnormal{f}}}\quantity[\mathcal{G}\Gamma_{1}\phi\mathcal{G}\Gamma_{2}\phi]$.
To this end we introduce the matrices $\hat{e}_{AB}$, with indexes in the
range $A,B=\quantity{0,1,\ldots,N_{f}-1}$ and matrix elements given by:
$\quantity(\hat{e}_{AB})_{ij}=\delta_{Ai}\delta_{Bj}$. Explicitly, they are
given by:

Any $N_{f}\times N_{f}$ complex matrix, $C$, can be written as a linear
combination of the $\hat{e}_{AB}$ matrices alongside $N_{f}^{2}$ complex
coefficients, $c_{AB}$, as: $C=c_{AB}\hat{e}_{AB}$. These matrices are
linearly independent since the zero matrix can only be obtained if all
coefficients vanish. Additionally, one can check that the trace between two
$\hat{e}_{AB}$ matrices is,
$\tr\quantity[\hat{e}_{AB}\hat{e}_{CD}]=\delta_{AD}\delta_{BC}$. The
$N_{f}^{2}$ generators of the $U(N_{f})$ group, $\lambda_{a}$, can be written
in terms of this new basis. Using the completeness relation we can write:
$\lambda_{a}=R_{AB}^{a}\hat{e}_{AB}$. Multiplying both sides by $\hat{e}_{CD}$
and taking the trace, one can get the expression for the $R_{AB}^{a}$
coefficients:

$\displaystyle R_{AB}^{a}=\tr\quantity[\lambda_{a}\hat{e}_{BA}].$ (24)

Using this new set of matrices, the trace over flavor space is:

$\displaystyle\operatorname{{\tr}_{\mathnormal{f}}}\quantity[\mathcal{G}\Gamma_{1}\phi\mathcal{G}\Gamma_{2}\phi]$
$\displaystyle=R_{ij}^{a}R_{ji}^{b}\quantity(\mathcal{G})_{ii}\Gamma_{1}\phi_{a}\quantity(\mathcal{G})_{jj}\Gamma_{2}\phi_{b}.$
(25)

Here, $a,b,i,j=\\{0,1,\ldots,N_{f}-1\\}$ and the coefficients $R^{a}_{AB}$ are
given by Eq. (24). Also, we used the fact that the operator $\mathcal{G}$ is
diagonal in flavor space,
$\mathcal{G}=\operatorname{diag}\quantity(\mathcal{G}_{0},\mathcal{G}_{1},\ldots,\mathcal{G}_{N_{f}-1})$.
For the $N_{f}=3$ case, which will be studied later in this work, the nine
generators of the $U(3)$ group are constituted by the scaled unit matrix plus
the eight Gell-Mann matrices and one can identify the
$\quantity(\mathcal{G})_{ii}$ elements with the up, down and strange quark
propagators, i.e., $i=0=u$ (up), $i=1=d$ (down) and $i=2=s$ (strange) and
$\mathcal{G}=\operatorname{diag}\quantity(\mathcal{G}_{u},\mathcal{G}_{d},\mathcal{G}_{s})$.

Due to the lack of structure in color space, the trace over such degrees of
freedom yields a simple $N_{c}$ factor. The full trace over spacetime and
discrete indexes is given by:

$\displaystyle\Tr\quantity[T_{2}^{\phi\phi}\quantity[\Gamma_{1},\Gamma_{2}]]$
$\displaystyle=R_{ij}^{a}R_{ji}^{b}\frac{N_{c}}{2}\operatorname{{\tr}_{\mathrm{D}}}\int\differential[4]{x}\bra{x}\quantity(\mathcal{G})_{ii}\Gamma_{1}\phi_{a}\quantity(\mathcal{G})_{jj}\Gamma_{2}\phi_{b}\ket{x}.$
(26)

Evaluating it, yields the following result:

$\displaystyle\Tr\quantity[T_{2}^{\phi\phi}\quantity[\Gamma_{1},\Gamma_{2}]]$
$\displaystyle=R_{ij}^{a}R_{ji}^{b}\frac{N_{c}}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\phi_{a}\quantity[k]\phi_{b}\quantity[-k]\int_{\mathrm{Reg}[k]}\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\operatorname{{\tr}_{\mathrm{D}}}\quantity[\frac{1}{\not{p}_{i}-M_{i}}\Gamma_{1}\frac{1}{(\not{p}_{j}-\not{k})-M_{j}}\Gamma_{2}]$
$\displaystyle=-\frac{i}{2}R_{ij}^{a}R_{ji}^{b}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\phi_{a}\quantity[k]\phi_{b}\quantity[-k]\Pi_{ij}\quantity[\Gamma_{1},\Gamma_{2},k]$
$\displaystyle=-\frac{i}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\phi_{a}\quantity[k]\,\Pi_{ab}\quantity[\Gamma_{1},\Gamma_{2},k]\,\phi_{b}\quantity[-k].$
(27)

Here, we defined
$\not{p}_{f}=\gamma^{0}(p_{0}+\mu_{f})-\bm{\gamma}\cdot\bm{p}$, with $\mu_{f}$
the chemical potential of the quark with flavor $f$ [18]. Furthermore, we
applied regularization to both momentum integrations, with a small caveat: we
allowed the regularization of the innermost momentum integration ($p$) to be a
function of the outermost momentum ($k$). The reason for this will be
discussed later. The quark-antiquark polarization function with flavor
indexes, $\Pi_{ij}\quantity[\Gamma_{1},\Gamma_{2},k]$ and its analogue,
$\Pi_{ab}\quantity[\Gamma_{1},\Gamma_{2},k]$, with indexes coming from the
generators of the $U(N_{f})$ group, were defined in this expression as:

$\displaystyle\Pi_{ij}\quantity[\Gamma_{1},\Gamma_{2},k]$
$\displaystyle=iN_{c}\int_{\mathrm{Reg}[k]}\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\frac{\operatorname{{\tr}_{\mathrm{D}}}\quantity[(\not{p}_{i}+M_{i})\Gamma_{1}((\not{p}_{j}-\not{k})+M_{j})\Gamma_{2}]}{(p_{i}^{2}-M_{i}^{2})((p_{j}-k)^{2}-M_{j}^{2})},$
(28) $\displaystyle\Pi_{ab}\quantity[\Gamma_{1},\Gamma_{2},k]$
$\displaystyle=R_{ij}^{a}R_{ji}^{b}\Pi_{ij}\quantity[k].$ (29)

The only trace left to calculate is over Dirac indexes. This can be achieved
by choosing $\Gamma_{1}$ and $\Gamma_{2}$. As previously stated, only the
cases $\Gamma_{1}=\Gamma_{2}=I_{4}$ and $\Gamma_{1}=\Gamma_{2}=i\gamma^{5}$
are nonzero. Thus, we define the scalar quark-antiquark polarization function
as
$\Pi_{ij}^{S}[M_{i},M_{j},\mu_{i},\mu_{j},k]=\Pi_{ij}\quantity[I_{4},I_{4},k]$,
and the pseudoscalar quark-antiquark polarization function by
$\Pi_{ij}^{P}[M_{i},M_{j},\mu_{i},\mu_{j},k]=\Pi_{ij}\quantity[i\gamma^{5},i\gamma^{5},k]$.
Evaluating the trace over spinor space yields:

$\displaystyle\Pi_{ij}^{S|P}[M_{i},M_{j},\mu_{i},\mu_{j},k]=-\frac{N_{c}}{8\pi^{2}}\bigg{\\{}\mathcal{A}$
$\displaystyle[M_{i},\mu_{i},k]+\mathcal{A}[M_{j},\mu_{j},k]$
$\displaystyle+\quantity(\quantity(M_{i}\pm
M_{j})^{2}-\quantity(\mu_{i}-\mu_{j}+k_{0})^{2}+\bm{k}^{2})B_{0}[M_{i},M_{j},\mu_{i},\mu_{j},k]\bigg{\\}}.$
(30)

In the right-hand side of this equation, the positive sign $(+)$ is used for
the scalar function, $\Pi_{ij}^{S}$, while, the negative sign $(-)$, is used
for the pseudoscalar quark-antiquark polarization function, $\Pi_{ij}^{P}$.
Here, we also defined the two fermion line integral, $B_{0}$, in Minkowski
spacetime:

$\displaystyle B_{0}[M_{i},M_{j},\mu_{i},\mu_{j},k]$
$\displaystyle=-16\pi^{2}i\int_{\mathrm{Reg}[k]}\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\frac{1}{(\quantity(p_{0}+\mu_{i})^{2}-\bm{p}^{2}-M_{i}^{2})((p_{0}+\mu_{j}-k_{0})^{2}-(\bm{p}-\bm{k})^{2}-M_{j}^{2})}.$
(31)

The function $\mathcal{A}$ is defined by:

$\displaystyle\mathcal{A}[M_{i},\mu_{i},k]$
$\displaystyle=-16\pi^{2}i\int_{\mathrm{Reg}[k]}\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\frac{1}{\quantity(p_{0}+\mu_{i})^{2}-\bm{p}^{2}-M_{i}^{2}}.$
(32)

The integrand is exactly the same as the one present in the definition of the
one fermion line integral $A$ (see Eq. (9)). However, the function
$\mathcal{A}$ carries the same momentum dependent regularization as the
integral $B_{0}$. In the special case where the same regularization is applied
to both $A$ and $\mathcal{A}$, then $A=\mathcal{A}$. Hence, the integral $A$,
defined in Eq. (9) and the integral $\mathcal{A}$, defined above, are exactly
the same quantity but, within different regularizations.

Finally, we can define the meson propagators of the model. Using the series
expansion of Eq. (19) alongside the trace calculation of Eq. (27). The second
order term, proportional to $\sigma_{a}\sigma_{b}$ and $\pi_{a}\pi_{b}$, of
the quadratic action, $\mathcal{S}_{\mathrm{Q}}$, is given by:

$\displaystyle\mathcal{S}_{\mathrm{Q}}^{(2)}\quantity[\bm{\sigma},\bm{\pi},\expectationvalue{\bm{s}},\expectationvalue{\bm{p}}]=$
$\displaystyle-\frac{1}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\,\sigma_{a}\quantity[k]\quantity(\frac{1}{2}\quantity(S^{-1})_{ab}-\Pi_{ab}^{S}\quantity[k])\sigma_{b}\quantity[-k]$
$\displaystyle-\frac{1}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\,\pi_{a}\quantity[k]\quantity(\frac{1}{2}\quantity(P^{-1})_{ab}-\Pi_{ab}^{P}\quantity[k])\pi_{b}\quantity[-k].$
(33)

Using this result we can define the inverse scalar and pseudoscalar meson
propagators, $D_{S}^{-1}$ and $D_{P}^{-1}$, respectively. They are defined by:

$\displaystyle\quantity(D_{S}^{-1})_{ab}\quantity[k]$
$\displaystyle=\frac{1}{2}\quantity(S^{-1})_{ab}-\Pi_{ab}^{S}\quantity[k],$
(34) $\displaystyle\quantity(D_{P}^{-1})_{ab}\quantity[k]$
$\displaystyle=\frac{1}{2}\quantity(P^{-1})_{ab}-\Pi_{ab}^{P}\quantity[k].$
(35)

As we will see later, the meson propagators are essential in the calculation
of quark-quark and quark-antiquark cross sections.

####  II.2.1 Equivalence with the Bethe-Salpeter equation in the random phase
approximation

Figure 1: Bethe-Salpeter equation in the random phase approximation (RPA) [17,
18, 36, 20]. In the first line the transfer matrix, $\mathcal{T}$, is written
as an infinite sum over quark-antiquark polarization loops, $\Pi$. In the
second line, the equation is expressed in a self-consistent way. The double
lines represents the meson propagator and the solid lines the quark
propagators. The wiggly line between two black dots is used to signal the
tensor nature of the NJL scattering kernel, $\mathcal{K}$, due to its
structure in flavor, color and Dirac spaces.

Another commonly employed method to obtain the meson propagators involves the
Bethe-Salpeter equation, which describes the bound states of two body quantum
systems. In the so-called random phase approximation, this equation is given
by (see Fig. 1 for its diagrammatic representation) [17, 18]:

$\displaystyle i\mathcal{T}\quantity[k]$
$\displaystyle=i\mathcal{K}+i\mathcal{K}\Big{(}-i\Pi\quantity[k]\Big{)}i\mathcal{K}+i\mathcal{K}\Big{(}-i\Pi\quantity[k]\Big{)}i\mathcal{K}\Big{(}-i\Pi\quantity[k]\Big{)}i\mathcal{K}+\ldots$
$\displaystyle=i\mathcal{K}+i\mathcal{K}\Big{(}-i\Pi\quantity[k]\Big{)}i\mathcal{T}\quantity[k].$
(36)

Here, $\mathcal{T}$ is the so-called quark-antiquark scattering matrix,
$\mathcal{K}$ is the scattering kernel and $\Pi\quantity[k]$ is the quark-
antiquark polarization loop (see Fig. 2). These quantities contain several
indexes including color, flavor and spinor.

Figure 2: Diagram representing the scalar/pseudoscalar quark-antiquark
polarization function [18]. The solid lines represent quark propagators and
the wiggly line with a black dot represents $\Gamma_{a}^{X}$, with
$\Gamma_{a}^{S}=\lambda_{a}$, for scalar modes and
$\Gamma_{a}^{P}=i\gamma^{5}\lambda_{a}$, for pseudoscalar modes.

The Bethe-Salpeter equation in the random-phase approximation, applied to NJL-
type models in the mean field approximation (or Hartree approximation [17, 18,
19, 36, 20]), has been extensively applied in the literature to study in-
medium meson behavior, and to evaluate quark-quark and quark-antiquark cross
sections [42, 17, 18, 19, 1, 13, 11, 12, 36, 20, 43]. Such standard approach
[36], is known to be a symmetry preserving one, being compatible with chiral
symmetry and respecting the Goldstone theorem and the Gell-
Mann$-$Oakes$-$Renner relation [17, 18, 19, 37, 36, 20]. In this section we
highlight the equivalence between this method and the functional expansion
derived previously. To this end, one can use Eq. (36), to derive a self-
consistent expression for the scalar and pseudoscalar meson propagator,
$D_{ab}^{S|P}$. Consider the transfer matrix written in terms of the meson
propagator:
$\mathcal{T}=D^{XY}_{ab}\quantity(\Gamma_{a}^{X}\otimes\Gamma_{b}^{Y})$, with
$\Gamma_{a}^{X}$ carrying the tensor structure of the vertex,
$\Gamma_{a}^{S}=\lambda_{a}$ and $\Gamma_{a}^{P}=i\gamma^{5}\lambda_{a}$ [17,
18, 36]. In a theory with only four-Fermi interactions (like the one
considered in this work, see the quadratic Lagrangian in Eq. (13)), the
scattering kernel is given by:
$\mathcal{K}=2g_{ab}^{XY}\quantity(\Gamma_{a}^{X}\otimes\Gamma_{b}^{Y})$, with
$g_{ab}^{XY}=g_{ab}^{X}\,\delta_{XY}$ and $g_{ab}^{S}=S_{ab}$ and
$g_{ab}^{P}=P_{ab}$ for scalar and pseudoscalar degrees of freedom,
respectively [17, 18, 36]. The factor of $2$ comes from the Feynman rules and
is due to the application of Wick’s theorem [17, 18]. The polarization loop,
$\Pi_{ab}^{XY}$ is the quark-antiquark polarization loop defined by [17, 18,
36, 20]:

$\displaystyle\Pi_{ab}^{XY}\quantity[k]$
$\displaystyle=i\int\frac{\differential[4]{p}}{\quantity(2\pi)^{4}}\tr\quantity[\frac{\not{p}+M}{p^{2}-M^{2}}\Gamma_{a}^{X}\frac{(\not{p}-\not{k})+M}{({p}-{k})^{2}-M^{2}}\Gamma_{b}^{Y}].$
(37)

Here, the trace is made over flavor, color and Dirac discrete indexes, i.e.,
$\tr\equiv\operatorname{{\tr}_{\mathrm{D}}}\operatorname{{\tr}_{\mathnormal{c}}}\operatorname{{\tr}_{\mathnormal{f}}}$.
After evaluating the trace, this expression agrees with the polarization loop
obtained via functional methods, defined in Eqs. (28) and (29). A diagrammatic
representation of the quark-antiquark polarization function is shown in Fig.
2. Using these definitions, one can write an equation for the meson propagator
as:

$\displaystyle
D_{ab}^{XY}\quantity(\Gamma_{a}^{X}\otimes\Gamma_{b}^{Y})=2g_{ab}^{XY}\quantity(\Gamma_{a}^{X}\otimes\Gamma_{b}^{Y})+2g_{ac}^{XK}\Pi_{cd}^{KL}D_{db}^{LY}\quantity(\Gamma_{a}^{X}\otimes\Gamma_{b}^{Y}).$
(38)

Since the matrix structure is completely contained in the $\Gamma_{a}^{X}$
quantities, we can write an equation for the scalar function, $D_{ab}^{XY}$:

$\displaystyle D_{ab}^{XY}=2g_{ab}^{XY}+2g_{ac}^{XK}\Pi_{cd}^{KL}D_{db}^{LY}.$
(39)

Considering $D_{ab}^{XY}$ as matrix elements, we can write the above equation
in matrix form [17, 18, 43], $D=2g+2g\Pi D$, which can be solved for the meson
propagator, $D$, to yield [17, 18, 20, 43]:

$\displaystyle D$ $\displaystyle=\quantity[1-2g\Pi]^{-1}2g.$ (40)

Writing the flavor indexes explicitly, we can write the scalar and
pseudoscalar meson propagators, $D^{S|P}_{ab}$ as:

$\displaystyle
D^{S|P}_{ab}=\quantity[\frac{1}{2}\quantity(S|P^{-1})_{ab}-\Pi^{S|P}_{ab}]^{-1}.$
(41)

Hence, the meson propagators obtained using the Bethe-Salpeter equation, in
the random phase approximation, applied to the four-fermion Lagrangian density
of Eq. (13), are identical to the ones obtained from the functional expansion
presented earlier in this work, see Eqs. (34) and (35).

####  II.2.2 Comments on the expansion of the fermionic determinant

In the aforementioned calculations, the meson propagators were acquired
through a specific expansion method applied to the fermionic determinant. Our
approach involved utilizing the series representation of the logarithm of a
matrix to derive the expansion. However, the expansion obtained is not unique
and our particular choice is rooted on the fact that we wanted to derive the
same expression for the meson propagators as the one obtained using the Bethe-
Salpeter equation in the random phase approximation [18, 12, 11, 19, 17, 44,
36].

In a different expansion, one could have introduced in the calculation kinetic
terms for the meson fields. Such step could have been made when evaluating the
functional trace over the continuous spacetime indexes. Since the operator
$\quantity(\mathcal{G})_{jj}$ and the fields $\phi_{a}$ do not commute,
instead of directly evaluating the trace in Eq. (26), one could evaluate the
trace over these degrees of freedom independently, by moving all the functions
of momentum operators to the left and all functions of the position operator
to the right, as introduced in Refs. [33, 34, 35, 32]. This can be
accomplished by using the following identity [34, 35]:

$\displaystyle\phi\mathcal{G}=\mathcal{G}\phi+\mathcal{G}^{2}\quantity[\mathcal{G}^{-1},\phi]+\mathcal{G}^{3}\quantity[\mathcal{G}^{-1},\quantity[\mathcal{G}^{-1},\phi]]+\ldots$
(42)

As a side note, the first term in this identity diverges, requiring
regularization, while the following terms are not divergent [21]. By
substituting $\mathcal{G}=\quantity(\mathcal{G})_{jj}$ and $\phi=\phi_{a}$ one
can understand that the tower of iterated commutators are proportional to
field derivatives like $\partial_{\mu}\phi_{a}$, $p\cdot\partial\phi_{a}$,
$\Box\phi_{a}$ , etc. Following the use of this identity, one can write Eq.
(26) as:

$\displaystyle\Tr\quantity[T_{2}^{\phi\phi}\quantity[\Gamma_{1},\Gamma_{2}]]$
$\displaystyle=R_{ij}^{a}R_{ji}^{b}\frac{N_{c}}{2}\operatorname{{\tr}_{\mathrm{D}}}\int\differential[4]{x}\bra{x}\quantity(\mathcal{G})_{ii}\Gamma_{1}\quantity(\mathcal{G})_{jj}\Gamma_{2}\ket{x}\phi_{a}(x)\phi_{b}(x)+\text{$\phi$
derivative terms}.$ (43)

Hence, in such expansion, the action of the model would contain kinetic terms
for the meson fields, $\sigma_{a}$ and $\pi_{a}$, not present in the our
previous calculation. Additionally, only one fermion line quark loops would
contribute to the calculation without any contribution coming from the two
fermion line integral, $B_{0}$.

We would also like to highlight the fact that such result is only obtained if
one starts the evaluation of the total trace with the trace over flavor space.
Starting the trace evaluation by applying directly the identity of Eq. (42),
as a first step in the trace calculation leads to an infinite tower of terms
proportional to iterated commutators between the mass matrix $M$ and the
fields, $\quantity[M,\phi]$, $\quantity[M,\quantity[M,\phi]]$, …Of course, in
the case in which the mass matrix is diagonal and degenerate in flavor space,
$M=M_{q}\operatorname{diag}\quantity(1,1,\dots,1)$, all these contribution
vanish (with $M_{q}$ the effective quark mass). However in the case in which
this matrix is not diagonal and degenerate, this cancellation does not occur.
This is the case when considering the physical scenario of $N_{f}=3$ in which
the strange quark mass is considerably higher than the up and down quark
masses. Furthermore, when considering finite chemical potential, there is an
extra matrix contribution to the $\mathcal{G}^{-1}$ operator, in the form of
an effective quark chemical potential matrix,
$\mu=\operatorname{diag}\quantity(\mu_{0},\mu_{1},\ldots,\mu_{N_{f}-1})$. Just
like the mass case, in physical scenarios in which this quantity is not
diagonal and degenerate an infinite sequence of iterated commutators between
$\mu$ and the fields appears, $\quantity[\mu,\phi]$,
$\quantity[\mu,\quantity[\mu,\phi]]$, …Finally, if one wishes to allow for the
possibility of pseudoscalar condensation in the medium, the operator
$\mathcal{G}^{-1}$ will certainly contain contributions which are not diagonal
and degenerate.

To end this section we would like to point out another technique that can be
used to extract the meson propagators from the fermionic determinant, the so-
called heat kernel expansion [32, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54].
Such method has been widely studied and developed in the context of two and
three flavor NJL-type models, including different types of quark interactions
at the Lagrangian level [49, 55, 48, 52, 56]. When applying this expansion to
such type of models, one considers a Wick rotation of the operator inside the
fermionic determinant in order to make it possible to obtain its Schwinger
proper-time representation. In general, one can write [52]:

$\displaystyle\ln\absolutevalue{\det\quantity[D]}=-\frac{1}{2}\int_{0}^{\infty}\frac{\differential{t}}{t}\rho\quantity[t,\Lambda^{2}]\Tr\quantity[\mathrm{e}^{-tD_{\mathrm{E}}^{\dagger}D_{\mathrm{E}}}].$
(44)

Here, $\rho\quantity[t,\Lambda^{2}]$ is a function, often called regulator or
kernel, which incorporates the necessary regularization [48],
$\Tr\quantity[\ldots]$ stands for the total trace and the operator
$D_{\mathrm{E}}$, is the Wick rotated Dirac operator in the presence of
background fields $\sigma=\sigma_{a}\lambda_{a}$ and $\pi=\pi_{a}\lambda_{a}$.
Additionally, $D_{\mathrm{E}}^{\dagger}D_{\mathrm{E}}=M^{2}-\partial^{2}+Y$,
with $M=M_{a}\lambda_{a}$ the effective quark mass matrix and $Y$ contains
only contributions coming from the background fields [52, 49]. In the case in
which the quark mass matrix is not diagonal and degenerate, it does not
commute with the rest of the operator and a resummation technique was
developed in order to solve this difficulty [49, 55, 48].

Applying the heat kernel expansion to Eq. (17), would lead to a different
model, with different meson propagators when compared to our approach, which
are equivalent to the ones obtained via the Bethe-Salpeter equation in the
random phase approximation. Thus, it could be insightful to evaluate the NJL
quark-quark and quark-antiquark cross sections from meson propagators derived
from a heat kernel expansion. The results of such study could then be directly
compared to the ones obtained in this work. This could be used to assess the
impact of considering different approximations to extract meson propagators
from a purely fermionic theory as is the case of the NJL model.

###  II.3 The 3-momentum sphere intersection regularization

As previously established, in a 4-dimensional spacetime, the one and the two
fermion lines integrals ($A$ and $B_{0}$) are divergent quantities. Other
integrals with higher fermion line count, like the three fermion line integral
($C_{0}$) are, on the other hand, convergent. When faced with divergent
integrals, some regularization procedure must be applied. Then, the
regularization procedure becomes part of the model in question [21] and should
be carried out in a symmetry preserving manner. The primary objective is to
maintain the expected physical properties of the model unchanged after
applying the regularization.

To attain such goal, a plethora of techniques was developed over the years:
the 4-momentum regularization, the 3-momentum regularization, dimensional
regularization, Pauli-Villars regularization, implicit regularizations and the
proper time regularization [57, 58, 37, 21, 59, 60, 61, 62]. The choice of
regularization is of extreme importance since it introduces a scale, which
informs about the range of applicability of a particular model [21]. Such
choice must also take into account the physical scenario under study. For
instance, in quantum field theory at finite temperature and density, a
compatible regularization scheme must be chosen, e.g., if a full sum over
Matsubara frequencies is performed, one is implicitly considering the full
range of the integration over the temporal component of the momentum.

The fact that the $A$ and $B_{0}$ are divergent, poses the following question:
in a model calculation, does one simply regularize the diverging integrals, or
is the regularization scheme applied to all integrals? For instance, in the
NJL model, thermodynamic quantities like pressure and energy density contain
both the well known divergent zero energy mode term (which corresponds to the
vacuum contribution), as well as, convergent contributions (which are
temperature and chemical potential dependent and can be regarded as in-medium
contributions). The entropy density, in the same model, is defined by
convergent contributions. In cases like this, one must choose where to apply
the regularization (see for instance [63]). Indeed, in some works which only
deal with the one fermion line integral, only the vacuum contribution of the
$A$ integral is regularized [64]. One can also find studies where the
convergent contribution of the $A$ integral is not regularized, but the
convergent contribution of the $B_{0}$ integral is, effectively mixing
regularization strategies [65, 66, 16, 63]. In this work, as it was performed
in Ref. [11], we will use the Matsubara frequency summation to introduce
temperature and apply the 3-momentum regularization to the remaining
integrals. Furthermore, we choose to consider the regularization as part of
the model, affecting all contributions, divergent and convergent alike.

At this point we introduce and discuss the 3-momentum sphere intersection
regularization which will be applied to the quark-antiquark polarization
function $\Pi_{ij}^{S|P}$. In the 3-momentum regularization, the quark
propagator is defined up to the maximum allowed momentum, the cutoff,
$\Lambda$. In such scheme, a quark can have any 3-momentum inside a sphere of
radius $\Lambda$. In the definition of the scalar and pseudoscalar
polarization functions written in Eqs. (28) and (29), there are two quark
propagators which must be integrated over a loop, one with momenta $p$ and
another with $p-k$, which carries an insertion from an external momentum $k$.
In order to truly restrict the 3-momentum of the quarks, both momenta must be
smaller than the cutoff, $\Lambda$. This means that the absolute value of the
quark loop 3-momentum, $\absolutevalue{\bm{p}}$, present in the integral of
the polarization loop functions, must obey:
$0\leq\absolutevalue{\bm{p}}\leq\Lambda\wedge
0\leq\absolutevalue{\bm{p}-\bm{k}}\leq\Lambda$. Such restriction can be
visualized as the intersection between two spheres. Let
$\mathbb{S}^{\Lambda}_{0}$ and $\mathbb{S}_{\absolutevalue{\bm{k}}}^{\Lambda}$
be spheres of radius $\Lambda$, where the first is centered at the origin
($0$) and the second at a point distanced $\absolutevalue{\bm{k}}$ from the
origin. We denote the intersection of these spheres333The volume of the
intersection between two spheres of radius $\Lambda$, distanced k between each
other, can be calculated to yield,
$V_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}=\frac{1}{12}\pi\quantity(2\Lambda-k)^{2}\quantity(4\Lambda+k)$.
In the $k\to 0$ limit one trivially recovers the volume of a sphere of radius
$\Lambda$, $V_{\mathbb{I}_{0}^{\Lambda}}=\frac{4}{3}\pi\Lambda^{3}$.  by
$\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}=\mathbb{S}_{0}^{\Lambda}\cap\mathbb{S}_{\absolutevalue{\bm{k}}}^{\Lambda}$.
Hence, the integrals present in the quark-antiquark polarization functions
will be regularized by the intersection of these two spheres of radii
$\Lambda$, distanced by $\absolutevalue{\bm{k}}$. An interesting property of
the 3-momentum sphere intersection regularization is that it becomes the usual
3-momentum regularization for zero external momentum: in such case both
spheres completely intersect each other. This means that the one fermion line
integral present in the quark-antiquark polarization function
$\mathcal{A}[M_{i},\absolutevalue{\bm{k}}]$, regularized with the sphere
intersection scheme, degenerates to the usual one fermion line integral
$A[M_{i}]$ in the $\absolutevalue{\bm{k}}\to 0$ limit, i.e.,
$A[M_{i}]=\mathcal{A}[M_{i},0]$. In Fig. 3, we show the region of integration
which regularizes the quark-antiquark polarization function, $\Pi_{ij}^{S|P}$.
In the left panel of this figure, we show a particular 3-dimensional view of
this integration region which corresponds to a “convex lens shaped” volume,
defined as $\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}$. In the right panel
of the same figure, we show the two-dimensional cross section of the region in
the $p_{y}-p_{z}$ plane.

Figure 3: Region of integration in the 3-momentum sphere intersection
regularization scheme, used to regularize the quark-antiquark polarization
function, $\Pi_{ij}^{S|P}$. In the left panel: a particular 3-dimensional view
of the intersection of two spheres. The red sphere is centered at the origin
(defined by the equation $\bm{p}^{2}\leq\Lambda^{2}$) and the green sphere is
distanced $\absolutevalue{\bm{k}}$ from it, in the $p_{z}$ direction (defined
by the equation $\quantity(\bm{p}-\bm{k})^{2}\leq\Lambda^{2}$). The blue
region corresponds to the intersection between these two spheres and is
designated by $\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}$. In the right
panel, we show the $p_{y}-p_{z}$ cross section of the same integration region.

Another interesting property of this regularization scheme is its behavior for
$\absolutevalue{\bm{k}}\to\infty$. As a matter of fact, for
$\absolutevalue{\bm{k}}\geq 2\Lambda$, the spheres no longer intersect, the
integration region disappears, the integrals $\mathcal{A}$ and $B_{0}$ are
automatically zero and the quark-antiquark polarization functions vanish. So,
the asymptotic behavior of the meson propagators are given by:

$\displaystyle D^{S|P}$
$\displaystyle=\quantity[\frac{1}{2}\quantity(S|P)^{-1}-\Pi^{S|P}]^{-1}\underset{\absolutevalue{\bm{k}}\to\infty}{\Rightarrow}D_{ab}=2\quantity(S|P)_{ab}.$
(45)

Thus, for large momentum, the meson propagators converge to twice the meson
projectors ($2S_{ab}$ and $2P_{ab}$), which are temperature and chemical
potential dependent quantities.

This novel approach to the 3-momentum regularization, can also be extended to
integrals with three or more quark propagators present in the integrand: in
such cases the region of integration is the intersection between the
appropriate number of spheres. For example, for the three fermion line
integral, $C_{0}$, three quark propagators would be present, meaning that the
integral would be calculated over the intersection between three spheres of
radii $\Lambda$. One sphere located at the origin, other distanced
$\absolutevalue{\bm{k}}$ from the origin and another at distance
$\absolutevalue{\bm{q}}$ from the origin, with an angle
$\delta_{\bm{k}\bm{q}}$ between the 3-vectors $\bm{k}$ and $\bm{q}$. In Fig.
4, we show the $p_{y}-p_{z}$ cross section of the integration region for the
three fermion line integral, within the sphere intersection regularization
scheme. We leave the details of the explicit calculation of the $C_{0}$
integral, within this new regularization, to future work.

Figure 4:  Cross section in the $p_{y}-p_{z}$ plane of the region of
integration for the three fermion line integral, $C_{0}$, in the 3-momentum
sphere intersection regularization scheme.

##  III The in-medium two fermion, one loop integral

In order to introduce finite temperature, one first considers a Wick rotation
to Euclidean spacetime, and then applies the Matsubara formalism (see for
instance Ref. [67]). Through this process Eq. (31) yields:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},i\nu_{m},\absolutevalue{\bm{k}}]=16\pi^{2}T\sum^{+\infty}_{n=-\infty}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}$
$\displaystyle\frac{1}{\quantity(i\omega_{n}+\mu_{i})^{2}-E\quantity[M_{i},\bm{p}]^{2}}\times$
$\displaystyle\qquad\frac{1}{\quantity(i\omega_{n}+\mu_{j}-i\nu_{m})^{2}-E[M_{j},\bm{p}-\bm{k}]^{2}},$
(46)

with $\omega_{n}=(2n+1)\pi T$ denoting the Matsubara frequency for fermions
and $\nu_{m}$ the Matsubara frequency for the external particle. Here, there
are two distinct chemical potentials, $\mu_{i}$ and $\mu_{j}$, each
corresponding to different fermion masses, $M_{i}$ and $M_{j}$ and
$E\quantity[M,\bm{p}]$ is the usual dispersion relation:

$\displaystyle E\quantity[M,\bm{p}]=\sqrt{\bm{p}^{2}+M^{2}}.$ (47)

Applying a Wick rotation to the external momentum, one can perform an
analytical continuation as, $i\nu_{m}\to k_{0}$.

One very important symmetry of this function arises when considering a shift
in the Matsubara frequency as: $\omega_{n}\to-\omega_{n}+\nu_{m}$, which is
allowed in the case of this infinite sum. Applying this substitution and
performing the aforementioned Wick rotation, one can write:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=$
$\displaystyle
B_{0}\quantity[M_{j},M_{i},T,-\mu_{j},-\mu_{i},k_{0},\absolutevalue{\bm{k}}].$
(48)

Hence, the two fermion, one loop integral $B_{0}$, is invariant with respect
to the interchange of the masses, $M_{i}\leftrightarrow M_{j}$, and chemical
potentials, albeit with an extra minus sign, $\mu_{i}\rightarrow-\mu_{j}$ and
$\mu_{j}\rightarrow-\mu_{i}$. Such symmetry can also be understood as
exchanging the fermion in the upper line with the antifermion in the lower
line in the loop diagram shown in Fig. 2, where time flows from left to right.

When dealing with the integration over the 3-momentum, $\bm{p}$, it will
become necessary to introduce an infinitesimal shift in the complex plane,
$\pm i\epsilon$ with $\epsilon>0$, in order to deform the contour around
singularities. As discussed in Ref. [11], one can regard the masses as having
an infinitesimal complex contribution, $M^{2}\to M^{2}-i\epsilon$, or consider
other recipes, such as a shift in the zero component of the external momentum,
$k_{0}\to k_{0}\pm i\epsilon$ [68]. The way in which the $\pm i\epsilon$ term
is introduced in the calculation has an impact in a very important symmetry
property of the $B_{0}$ function. In practice, choosing different recipes
amounts to different time orderings of the single particle Green function.
However, for the purposes of calculating meson masses in the NJL model and
(differential) cross sections, the chosen $\pm i\epsilon$ recipe will not
change the results. The reason for this “recipe invariance” will be discussed
later.

Considering the usual shift in the masses, $M^{2}\to M^{2}-i\epsilon$ (see
Ref. [11]), the $B_{0}$ function takes the form:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=16\pi^{2}T\sum^{+\infty}_{n=-\infty}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}$
$\displaystyle\frac{1}{\quantity(i\omega_{n}+\mu_{i})^{2}-E\quantity[M_{i},\bm{p}]^{2}+i\epsilon}\times$
$\displaystyle\qquad\frac{1}{\quantity(i\omega_{n}+\mu_{j}-k_{0})^{2}-E[M_{j},\bm{p}-\bm{k}]^{2}+i\epsilon}.$
(49)

Applying the simple change of variable, $\omega_{n}\rightarrow-\omega_{n}$, in
the Matsubara sum, one can establish the following connection:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=B_{0}\quantity[M_{i},M_{j},T,-\mu_{i},-\mu_{j},-k_{0},\absolutevalue{\bm{k}}].$
(50)

Hence, for this choice of shift, the $B_{0}$ function is symmetric with
respect to the zero component of the external momentum, $k_{0}$. As discussed
above, another way to introduce the shift is to consider it in the external
momentum, $k_{0}\to k_{0}\pm i\epsilon$. In such case one can write:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=16\pi^{2}T\sum^{+\infty}_{n=-\infty}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}$
$\displaystyle\frac{1}{\quantity(i\omega_{n}+\mu_{i})^{2}-E\quantity[M_{i},\bm{p}]^{2}}\times$
$\displaystyle\qquad\frac{1}{\quantity(i\omega_{n}+\mu_{j}-k_{0}-i\epsilon)^{2}-E[M_{j},\bm{p}-\bm{k}]^{2}}.$
(51)

In this case, a different symmetry relation is obtained, involving the complex
conjugate of the $B_{0}$ function:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=B_{0}\quantity[M_{i},M_{j},T,-\mu_{i},-\mu_{j},-k_{0},\absolutevalue{\bm{k}}]^{\ast}.$
(52)

We highlight that, differently from the previous symmetry, the masses and
chemical potentials do not switch positions, one simply changes,
$\mu_{i}\to-\mu_{i}$, $\mu_{j}\to-\mu_{j}$ and $k_{0}\to-k_{0}$. In this work
we will use the previously introduced mass shift however, to simplify the
notation, we will only write the contribution $-i\epsilon$ explicitly, when it
is necessary. In later sections we will discuss why this symmetry is essential
in order to maintain the symmetry between meson and anti-meson pairs and,
consequently, provide a proper framework to evaluate the (differential) cross
sections of different quark-quark and quark-antiquark processes.

The infinite sum over the Matsubara frequencies can be performed analytically,
for example, by converting it into a contour integration on the complex plane.
In such approach, the Matsubara sum is written as a sum of residues which, in
turn, can be written as a sum of contour integrations where each contour
surrounds a simple pole located at a particular Matsubara frequency. The
contours can be deformed and the Cauchy residue theorem applied in order to
obtain a result for the original Matsubara sum [69]. One gets:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=16\pi^{2}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\Bigg{\\{}$
$\displaystyle-\frac{1}{4E_{1}}\frac{1-2f_{\mathrm{F}}\quantity[E_{1}-\mu_{i},T]}{\quantity(E_{1}-\lambda)^{2}-E_{2}^{2}}$
$\displaystyle\qquad-\frac{1}{4E_{1}}\frac{1-2f_{\mathrm{F}}\quantity[E_{1}+\mu_{i},T]}{\quantity(E_{1}+\lambda)^{2}-E_{2}^{2}}$
$\displaystyle\qquad\qquad-\frac{1}{4E_{2}}\frac{1-2f_{\mathrm{F}}\quantity[E_{2}-\mu_{j},T]}{\quantity(E_{2}+\lambda)^{2}-E_{1}^{2}}$
$\displaystyle\qquad\qquad\qquad-\frac{1}{4E_{2}}\frac{1-2f_{\mathrm{F}}\quantity[E_{2}+\mu_{j},T]}{\quantity(E_{2}-\lambda)^{2}-E_{1}^{2}}\Bigg{\\}},$
(53)

where, $f_{\mathrm{F}}$ is the Fermi-Dirac distribution function, given by:

$\displaystyle f_{\mathrm{F}}\quantity[E,T]$
$\displaystyle=\frac{1}{\mathrm{e}^{E/T}+1},$ (54)

and, for simplicity, we have defined:

$\displaystyle E_{1}$ $\displaystyle=E\quantity[M_{i},\bm{p}],$ (55)
$\displaystyle E_{2}$ $\displaystyle=E[M_{j},\bm{p}-\bm{k}],$ (56)
$\displaystyle\lambda$ $\displaystyle=k_{0}+\mu_{i}-\mu_{j}.$ (57)

Figure 5: Pair creation/annihilation processes [68]. Figure 6: Scattering
processes: absorption and emission [68, 70].

As introduced in [68, 71] it is quite useful to write this function using the
so-called pair creation/annihilation, $B_{0,\mathrm{p}}$, and scattering
contributions, $B_{0,\mathrm{s}}$. In Figs. 5 and 6, we show schematic
diagrams representing the pair creation/annihilation and scattering
contributions, respectively. For a deeper explanation about all these
different processes see Ref. [70]. Using this separation, we can write:

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]+B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}],$
(58)

where, each component is given explicitly by:

$\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=16\pi^{2}\sum_{\eta=\pm
1}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{4E_{1}E_{2}}\frac{1-f_{\mathrm{F}}\quantity[E_{1}+\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E_{2}-\eta\mu_{j},T]}{E_{2}+E_{1}+\eta\lambda-i\epsilon},$
(59) $\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=-16\pi^{2}\sum_{\eta=\pm
1}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{4E_{1}E_{2}}\frac{f_{\mathrm{F}}\quantity[E_{1}-\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E_{2}-\eta\mu_{j},T]}{E_{2}-E_{1}+\eta\lambda-\operatorname{sgn}\quantity[E_{2}-E_{1}]i\epsilon}.$
(60)

Here, we introduced the mass shift, $M^{2}\to M^{2}-i\epsilon$ [11], which
will be necessary in order to use the Sokhotski–Plemelj theorem444Introducing
the mass shift and expanding the difference $E_{2}[\epsilon]-E_{1}[\epsilon]$
around zero yields:
$E_{2}[\epsilon]-E_{1}[\epsilon]=\quantity(E_{2}-E_{1})-i\epsilon\frac{\quantity(E_{2}-E_{1})}{2E_{1}E_{2}}+O\quantity[\epsilon^{2}].$
Since $E_{1},E_{2}>0$, the sign of the first order term will depend on the
difference of these energies and the denominator is absorbed by the
infinitesimal $\epsilon$. . The function $\operatorname{sgn}[x]$ returns the
sign of the quantity $x$, i.e., $\operatorname{sgn}[x]=+1$ if $x>0$ and
$\operatorname{sgn}[x]=-1$ if $x<0$.

In the vacuum, the scattering contribution, $B_{0,\mathrm{s}}$, vanishes and
only $B_{0,\mathrm{p}}$ will contribute: in the vacuum there is no medium for
scattering to occur, only the creation and annihilation processes contribute.
This vacuum contribution is the only one which diverges. For cases with non-
zero temperature and/or density the scattering term will contribute. In the
later case, for $T=0$ and finite chemical potential, one can introduce a Fermi
momentum and easily perform two of the three integrations.

###  III.1 Practical implementation

####  III.1.1  The $\absolutevalue{\bm{k}}>0$ case

Lets discuss the implementation of the $B_{0}$ integral for the case where
$\absolutevalue{\bm{k}}>0$. Since the integrand does not depend on the
azimuthal angle, the integration over such variables is trivial, yielding a
factor of $2\pi$. The remaining 2-dimensional integration is more complicated.
One way to simplify the calculations and perform such integration, is to
consider a change of variables. Here, we extend the approach made for the
degenerate mass case ($M_{i}=M_{j}$) that was introduced in [71] and also
applied in [72], by taking into account different masses. Consider the
following change of variables:

$\displaystyle E$ $\displaystyle=\frac{1}{2}\quantity(E_{2}+E_{1}),$ (61)
$\displaystyle\varepsilon$ $\displaystyle=E_{2}-E_{1}.$ (62)

From these definitions one has the following relations:
$E_{1}=E-\nicefrac{{\varepsilon}}{{2}}$ and
$E_{2}=E+\nicefrac{{\varepsilon}}{{2}}$. Looking at Eq. (62), the quantity
$\varepsilon$ can be both positive or negative while, checking Eq. (61), $E$
is always a positive quantity. We also point out that the quantities $E_{1}$
and $E_{2}$ are always positive since they are energies.

The integration measure will be changed as
$\differential{E}\differential{\varepsilon}=\absolutevalue{\det[J]}\differential{\absolutevalue{\bm{p}}}\differential{\theta}$,
where
$J=J\quantity[\absolutevalue{\bm{p}},\theta]=\frac{\partial(E,\varepsilon)}{\partial(\absolutevalue{\bm{p}},\theta)}$,
is the Jacobian matrix. Performing the derivatives, one can write:

$\displaystyle\det[J]=\absolutevalue{\bm{k}}\frac{\bm{p}^{2}\sin[\theta]}{E_{1}E_{2}}.$
(63)

From the above one can realize that this change of variables is only valid for
the $\absolutevalue{\bm{k}}>0$. We treat the case $\absolutevalue{\bm{k}}=0$
as a special case since it is simpler and a different change of variables must
be used.

Considering the 3-momentum sphere intersection regularization,
$\mathrm{Reg}[\absolutevalue{\bm{k}}]\to\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}$
and performing the proposed change of variables, one can write the integrals
as:

$\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{E}\differential{\varepsilon}\frac{1-f_{\mathrm{F}}\quantity[E-\nicefrac{{\varepsilon}}{{2}}+\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E+\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{j},T]}{E+\nicefrac{{\eta\lambda}}{{2}}-i\epsilon},$
(64) $\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{1}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{E}\differential{\varepsilon}\frac{f_{\mathrm{F}}\quantity[E-\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E+\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{j},T]}{\varepsilon+\eta\lambda-\operatorname{sgn}\quantity[\varepsilon]i\epsilon}.$
(65)

The integration over the azimuthal angle yields $2\pi$ and the integration
over the remaining variables corresponds to evaluating the integral over the
2-dimensional cross section of the intersection between two spheres (see Fig.
3). We have written the 2-dimensional integration without referring the
specific order of integration. In fact it is extremely useful to consider both
order of integrations. For the pair creation/annihilation contribution,
$B_{0,\mathrm{p}}$, we will use the order,
$\int_{E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{E}\int_{\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}$,
while for the scattering contribution, $B_{0,\mathrm{s}}$, the order
$\int_{\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{\varepsilon}\int_{E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{E}$
will be used. In principle, if the integrand is a well behaved function,
switching the order of integration yields the same results [73]. The reason
for using the two possible integration orders is twofold. First, it allows to
perform one of the integrations analytically, resulting in only one numerical
integration in order to evaluate the $B_{0}$ function. Second, it allows for a
straightforward isolation of the simple poles present at
$E=\nicefrac{{-\eta\lambda}}{{2}}$ and $\varepsilon=-\eta\lambda$, allowing
the use of the appropriate numerical methods to evaluate such type of
integrals. For convenience, we define the functions $g_{\mathrm{p}}$ and
$g_{\mathrm{s}}$ as:

$\displaystyle g_{\mathrm{p}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon]$
$\displaystyle=1-f_{\mathrm{F}}\quantity[E-\nicefrac{{\varepsilon}}{{2}}+\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E+\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{j},T],$
(66) $\displaystyle
g_{\mathrm{s}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon]$
$\displaystyle=f_{\mathrm{F}}\quantity[E-\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{i},T]-f_{\mathrm{F}}\quantity[E+\nicefrac{{\varepsilon}}{{2}}-\eta\mu_{j},T].$
(67)

For the pair creation/annihilation contribution, we write:

$\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{E}\frac{1}{E+\nicefrac{{\eta\lambda}}{{2}}-i\epsilon}\int_{\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}g_{\mathrm{p}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon]$
$\displaystyle=\frac{1}{2\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{E}\frac{G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},E]}{E+\nicefrac{{\eta\lambda}}{{2}}-i\epsilon}.$
(68)

While, for the scattering contribution, one gets:

$\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{1}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{\varepsilon}\frac{1}{\varepsilon+\eta\lambda-\operatorname{sgn}\quantity[\varepsilon]i\epsilon}\int_{E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{E}g_{\mathrm{s}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon]$
$\displaystyle=-\frac{1}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{\varepsilon}\frac{G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\varepsilon]}{\varepsilon+\eta\lambda-\operatorname{sgn}\quantity[\varepsilon]i\epsilon}.$
(69)

Both integration regions and all the functions that define them,
$E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}$,
$E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}$,
$\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}$,
$\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}$,
$E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$,
$\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$,
$\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$,
are explicitly derived in Appendix C. We point out that the integration
endpoints will be functions of several variables, such as $M_{i}$, $M_{j}$,
$\absolutevalue{\bm{k}}$ and $\Lambda$. However, here, such dependencies are
omitted for simplicity. The functions $G_{\mathrm{p}}$ and $G_{\mathrm{s}}$
have been defined as the innermost integration for each contribution,
separately, i.e.:

$\displaystyle
G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},E]$
$\displaystyle=\int_{\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}g_{\mathrm{p}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon],$
(70) $\displaystyle
G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\varepsilon]$
$\displaystyle=\int_{E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{E}g_{\mathrm{s}}\quantity[T,\mu_{i},\mu_{j},\eta,E,\varepsilon].$
(71)

The reason for this separations is that these integrals can be computed
analytically, simplifying the calculation of the $B_{0}$ function. For the
pair creation/annihilation contribution, $G_{\mathrm{p}}$, one
gets555Numerically, for low temperatures, it can be useful to make use of the
Puiseux series expansion of $\ln[1+x]$ at $x\to\infty$:
$\ln[1+x]=\ln[x]+1/x-1/(2x^{2})+1/(3x^{3})-1/(4x^{4})+O[1/x^{5}]$.:

$\displaystyle
G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},E]=\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}-\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}$
$\displaystyle+2T\ln\quantity[1+\mathrm{e}^{\quantity(-E+\frac{1}{2}\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}-\eta\mu_{i})/T}]$
$\displaystyle-2T\ln\quantity[1+\mathrm{e}^{\quantity(-E-\frac{1}{2}\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}+\eta\mu_{j})/T}]$
$\displaystyle-2T\ln\quantity[1+\mathrm{e}^{\quantity(-E+\frac{1}{2}\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}-\eta\mu_{i})/T}]$
$\displaystyle+2T\ln\quantity[1+\mathrm{e}^{\quantity(-E-\frac{1}{2}\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}+\eta\mu_{j})/T}].$
(72)

For the scattering contribution, one gets:

$\displaystyle
G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\varepsilon]=$
$\displaystyle-T\ln\quantity[1+\mathrm{e}^{\quantity(-E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}-\frac{1}{2}\varepsilon+\eta\mu_{j})/T}]$
$\displaystyle+T\ln\quantity[1+\mathrm{e}^{\quantity(-E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}+\frac{1}{2}\varepsilon+\eta\mu_{i})/T}]$
$\displaystyle-T\ln\quantity[1+\mathrm{e}^{\quantity(-E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}+\frac{1}{2}\varepsilon+\eta\mu_{i})/T}]$
$\displaystyle+T\ln\quantity[1+\mathrm{e}^{\quantity(-E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}-\frac{1}{2}\varepsilon+\eta\mu_{j})/T}].$
(73)

In the case where a singularity resides within the integration bounds, the
integration is made using the Sokhotski–Plemelj theorem and the $B_{0}$
function presents both real and imaginary components. The real part of the
pair creation/annihilation contribution, $B_{0,\mathrm{p}}$, and the real part
of the scattering contribution, $B_{0,\mathrm{s}}$, are given by:

$\displaystyle\real[B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\Big{]}$
$\displaystyle=\frac{1}{2\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\mathscr{P}\int_{E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{E}\frac{G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},E]}{E+\nicefrac{{\eta\lambda}}{{2}}},$
(74)
$\displaystyle\real[B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\Big{]}$
$\displaystyle=-\frac{1}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\mathscr{P}\int_{\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{\varepsilon}\frac{G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\varepsilon]}{\varepsilon+\eta\lambda}.$
(75)

The symbol $\mathscr{P}$ stands for Cauchy principal value. The respective
imaginary contributions are:

$\displaystyle\imaginary[B_{0,\mathrm{p}}$
$\displaystyle\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\Big{]}=$
$\displaystyle\frac{\pi}{2\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\left\\{\begin{array}[]{lr}G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\nicefrac{{-\eta\lambda}}{{2}}],&\text{if
}\nicefrac{{-\eta\lambda}}{{2}}\in\quantity[E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})},E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}],\\\
0,&\text{otherwise },\end{array}\right.$ (76)
$\displaystyle\imaginary[B_{0,\mathrm{s}}$
$\displaystyle\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\Big{]}=$
$\displaystyle-$
$\displaystyle\frac{\pi}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\left\\{\begin{array}[]{lr}\operatorname{sgn}\quantity[-\eta\lambda]G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},-\eta\lambda],&\text{if
}-\eta\lambda\in\quantity[\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})},\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}],\\\
0,&\text{otherwise }.\end{array}\right.$ (77)

If one considers the external momentum shift, $k_{0}\to k_{0}\pm i\epsilon$,
the real part of both contributions, $\real\quantity[B_{0,\mathrm{p}}]$ and
$\real\quantity[B_{0,\mathrm{s}}]$, remains unchanged, however, the overall
sign in front of the imaginary contributions is altered. The imaginary
contributions are changed as:

$\displaystyle
G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\nicefrac{{-\eta\lambda}}{{2}}]$
$\displaystyle\to\mp\eta
G_{\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},\nicefrac{{-\eta\lambda}}{{2}}],$
(78)
$\displaystyle\operatorname{sgn}\quantity[-\eta\lambda]G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},-\eta\lambda]$
$\displaystyle\to\mp\eta
G_{\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},\eta,\absolutevalue{\bm{k}},-\eta\lambda].$
(79)

Finally, the $B_{0}$ function can be written using the real and imaginary
parts of the pair creation/annihilation and scattering contributions as (see
Eq. (58)):

$\displaystyle
B_{0}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]=$
$\displaystyle\real[B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]+B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\bigg{]}$
$\displaystyle+i\imaginary[B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]+B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},\absolutevalue{\bm{k}}]\bigg{]}.$
(80)

####  III.1.2  The $\absolutevalue{\bm{k}}=0$ case

The $\absolutevalue{\bm{k}}=0$ limiting case cannot be calculated from the
case $\absolutevalue{\bm{k}}>0$. As it was explained earlier (see Section
III.1.1), the Jacobian determinant used when making the change of variables is
zero for $\absolutevalue{\bm{k}}=0$ (see Eq. (63)). Hence, the case
$\absolutevalue{\bm{k}}=0$ must be treated separately. Differently from the
previous case, our calculation for the case $\absolutevalue{\bm{k}}=0$
completely agrees with the one found in [11]. In this case there is no
separation between the spheres and the region of integration is simply the one
inside a sphere of radius $\absolutevalue{\bm{p}}\leq\Lambda$, the usual
3-momentum regularization. Nonetheless, for completeness, we will briefly show
the explicit calculation.

It is easier to consider a particular change of variables for the integral
$B_{0,\mathrm{p}}$ and a different one for the integral $B_{0,\mathrm{s}}$.
Furthermore, in this special case the integrand does not have angular
dependencies, meaning that the integral over the azimuthal angle, $\varphi$,
and polar angle, $\theta$, yield an overall factor of $4\pi$ and an
integration over the norm of the 3-momentum, $\absolutevalue{\bm{p}}$.

Lets start by discussing the $\absolutevalue{\bm{k}}=0$ case of the pair
creation/annihilation contribution, $B_{0,\mathrm{p}}$ by considering the
$\absolutevalue{\bm{k}}=0$ limit of Eq. (59). Solving the equation,
$E=\quantity(E_{2}+E_{1})/2=\quantity(E\quantity[M_{i},{\bm{p}}]+E\quantity[M_{j},{\bm{p}}])/2$,
for the variable $\absolutevalue{\bm{p}}$, one can consider a change of
variables from $\absolutevalue{\bm{p}}\to E$. Using the positive solution, one
gets:

$\displaystyle\absolutevalue{\bm{p}}=\frac{1}{4}\sqrt{\frac{\quantity(2E-M_{i}+M_{j})\quantity(2E-M_{i}-M_{j})\quantity(2E+M_{i}-M_{j})\quantity(2E+M_{i}+M_{j})}{E^{2}}}=P\quantity[M_{i},M_{j},E].$
(81)

The integration measure, $\differential{\absolutevalue{\bm{p}}}$, changes as,
$\differential{\absolutevalue{\bm{p}}}\frac{\bm{p}^{2}}{E_{1}E_{2}}=\differential{E}\frac{P\quantity[M_{i},M_{j},E]}{E}$
and the integral is written as:

$\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]$
$\displaystyle=\sum_{\eta=\pm
1}\int_{E_{0}}^{E_{\Lambda}}\differential{E}\frac{P\quantity[M_{i},M_{j},E]}{E}\frac{1-f_{\mathrm{F}}[\overline{E_{1}}+\eta\mu_{i},T]-f_{\mathrm{F}}[\overline{E_{2}}-\eta\mu_{j},T]}{E+\nicefrac{{\eta\lambda}}{{2}}-i\epsilon}.$
(82)

Here, $E_{0}$ and $E_{\Lambda}$ are the endpoints obtained after the change of
variables and are given by:

$\displaystyle E_{0}$
$\displaystyle=\frac{1}{2}\quantity(E\quantity[M_{i},0]+E\quantity[M_{j},0])=\frac{1}{2}\quantity(M_{i}+M_{j}),$
(83) $\displaystyle E_{\Lambda}$
$\displaystyle=\frac{1}{2}\quantity(E\quantity[M_{i},\Lambda]+E\quantity[M_{j},\Lambda])=\frac{1}{2}\quantity(\sqrt{\Lambda^{2}+M_{i}^{2}}+\sqrt{\Lambda^{2}+M_{j}^{2}}).$
(84)

For simplicity, we introduced the functions $\overline{E_{1}}$ and
$\overline{E_{2}}$, which are the functions defined in Eqs. (55) and (56),
written in terms of $E$ instead of $\bm{p}$. They are explicitly given by:

$\displaystyle\overline{E_{1}}$
$\displaystyle=\overline{E_{1}}\quantity[M_{i},M_{j},E]=2E-\frac{1}{4}\sqrt{\frac{\quantity(4E^{2}-M_{i}^{2}+M_{j}^{2})^{2}}{E^{2}}},$
(85) $\displaystyle\overline{E_{2}}$
$\displaystyle=\overline{E_{2}}\quantity[M_{i},M_{j},E]=2E-\frac{1}{4}\sqrt{\frac{\quantity(4E^{2}+M_{i}^{2}-M_{j}^{2})^{2}}{E^{2}}}.$
(86)

As for the case $\absolutevalue{\bm{k}}>0$, one can use the Sokhotski–Plemelj
theorem to write the real and imaginary parts, of the
$\absolutevalue{\bm{k}}=0$ limit of the $B_{0}$ function. We write:

$\displaystyle\real[$ $\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]\Big{]}=\sum_{\eta=\pm
1}\mathscr{P}\int_{E_{0}}^{E_{\Lambda}}\differential{E}\frac{P\quantity[M_{i},M_{j},E]}{E}\frac{1-f_{\mathrm{F}}[\overline{E_{1}}+\eta\mu_{i},T]-f_{\mathrm{F}}[\overline{E_{2}}-\eta\mu_{j},T]}{E+\nicefrac{{\eta\lambda}}{{2}}},$
(87) $\displaystyle\imaginary[$ $\displaystyle
B_{0,\mathrm{p}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]\Big{]}=$
$\displaystyle\qquad\qquad\pi\sum_{\eta=\pm
1}\left\\{\begin{array}[]{lr}\frac{P\quantity[M_{i},M_{j},\nicefrac{{-\eta\lambda}}{{2}}]}{\nicefrac{{-\eta\lambda}}{{2}}}\quantity(1-f_{\mathrm{F}}[\overline{E_{1}}+\eta\mu_{i},T]-f_{\mathrm{F}}[\overline{E_{2}}-\eta\mu_{j},T]),&\text{if
}\nicefrac{{-\eta\lambda}}{{2}}\in\quantity[E_{0},E_{\Lambda}],\\\
0,&\text{otherwise }.\end{array}\right.$

Again, by considering the external momentum shift, $k_{0}\to k_{0}\pm
i\epsilon$, instead of the mass shift, only the imaginary part is changed with
the following substitution:
$P\quantity[M_{i},M_{j},\nicefrac{{-\eta\lambda}}{{2}}]\to-\eta
P\quantity[M_{i},M_{j},\nicefrac{{-\eta\lambda}}{{2}}]$.

Finally, we deal with the scattering contribution, $B_{0,\mathrm{s}}$ applying
the $\absolutevalue{\bm{k}}=0$ limit to Eq. (60). For this case we consider a
change of variables based on solving the equation,
$\varepsilon=E_{2}-E_{1}=E\quantity[M_{j},{\bm{p}}]-E\quantity[M_{i},{\bm{p}}]$
for the absolute value of the momentum, $\absolutevalue{\bm{p}}$ and get:

$\displaystyle\absolutevalue{\bm{p}}=\frac{1}{2}\sqrt{\frac{\quantity(\varepsilon-
M_{i}+M_{j})\quantity(\varepsilon-
M_{i}-M_{j})\quantity(\varepsilon+M_{i}-M_{j})\quantity(\varepsilon+M_{i}+M_{j})}{\varepsilon^{2}}}=P\quantity[M_{i},M_{j},\varepsilon].$
(89)

In this case, the integration measure after the change of variables is,
$\frac{\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}}{E_{1}E_{2}}=-\frac{\differential{\varepsilon}P\quantity[M_{i},M_{j},\varepsilon]}{\varepsilon}$
and the integral can then be written as:

$\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]$
$\displaystyle=2\sum_{\eta=\pm
1}\int_{\varepsilon_{0}}^{\varepsilon_{\Lambda}}\differential{\varepsilon}\frac{P\quantity[M_{i},M_{j},\varepsilon]}{\varepsilon}\frac{f_{\mathrm{F}}[\hat{E_{1}}-\eta\mu_{i},T]-f_{\mathrm{F}}[\hat{E_{2}}-\eta\mu_{j},T]}{\varepsilon+\eta\lambda-\operatorname{sgn}\quantity[\varepsilon]i\epsilon}.$
(90)

As before, the functions, $\hat{E_{1}}$ and $\hat{E_{2}}$ are the functions
defined in Eqs. (55) and (56), written in terms of $\varepsilon$ instead of
$\bm{p}$ and are given by:

$\displaystyle\hat{E_{1}}$
$\displaystyle=\hat{E_{1}}\quantity[M_{i},M_{j},\varepsilon]=\frac{1}{2}\sqrt{\frac{\quantity(\varepsilon^{2}-M_{i}^{2}+M_{j}^{2})^{2}}{\varepsilon^{2}}}-\varepsilon,$
(91) $\displaystyle\hat{E_{2}}$
$\displaystyle=\hat{E_{2}}\quantity[M_{i},M_{j},\varepsilon]=\frac{1}{2}\sqrt{\frac{\quantity(\varepsilon^{2}+M_{i}^{2}-M_{j}^{2})^{2}}{\varepsilon^{2}}}+\varepsilon.$
(92)

The integration endpoints are defined as:

$\displaystyle\varepsilon_{0}$
$\displaystyle=E\quantity[M_{j},0]-E\quantity[M_{i},0]=M_{j}-M_{i},$ (93)
$\displaystyle\varepsilon_{\Lambda}$
$\displaystyle=E\quantity[M_{j},\Lambda]-E\quantity[M_{i},\Lambda]=\sqrt{\Lambda^{2}+M_{j}^{2}}-\sqrt{\Lambda^{2}+M_{i}^{2}}.$
(94)

Following the previous steps, the real and imaginary parts can be written as:

$\displaystyle\real[$ $\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]\Big{]}=2\sum_{\eta=\pm
1}\mathscr{P}\int_{\varepsilon_{0}}^{\varepsilon_{\Lambda}}\differential{\varepsilon}\frac{P\quantity[M_{i},M_{j},\varepsilon]}{\varepsilon}\frac{f_{\mathrm{F}}[\hat{E_{1}}-\eta\mu_{i},T]-f_{\mathrm{F}}[\hat{E_{2}}-\eta\mu_{j},T]}{\varepsilon+\eta\lambda},$
(95) $\displaystyle\imaginary[$ $\displaystyle
B_{0,\mathrm{s}}\quantity[M_{i},M_{j},T,\mu_{i},\mu_{j},k_{0},0]\Big{]}=$
$\displaystyle\qquad\qquad 2\pi\sum_{\eta=\pm
1}\left\\{\begin{array}[]{lr}\operatorname{sgn}\quantity[-\eta\lambda]\frac{P\quantity[M_{i},M_{j},-\eta\lambda]}{-\eta\lambda}(f_{\mathrm{F}}[\hat{E_{1}}-\eta\mu_{i},T]-f_{\mathrm{F}}[\hat{E_{2}}-\eta\mu_{j},T]),&\text{if
}-\eta\lambda\in\quantity[\varepsilon_{0},\varepsilon_{\Lambda}],\\\
0,&\text{otherwise }.\end{array}\right.$

In the momentum shift case, the imaginary part is modified as:
$\operatorname{sgn}\quantity[-\eta\lambda]P\quantity[M_{i},M_{j},-\eta\lambda]\to-\eta
P\quantity[M_{i},M_{j},-\eta\lambda]$. The complete function, $B_{0}$, for the
case with $\absolutevalue{\bm{k}}=0$, can then be obtained as before, using
Eq. (80) by using the real and imaginary parts calculated above.

In the Appendix D, we show and discuss some numerical results from evaluating
the two fermion line integral, $B_{0}$ for different values of temperature,
chemical potential and external momentum by fixing a cutoff $\Lambda$ and the
effective masses, $M_{i}$ and $M_{j}$.

##  IV The in-medium one fermion, one loop integral

As discussed in earlier sections, different versions of the one fermion line
integral arise when making calculations in the context of the NJL model. On
the one hand, it appears in the calculation of the quark condensate (integral
$A$ defined in Eq. (9)) while, on the other hand, it contributes to the quark-
antiquark polarization function (integral $\mathcal{A}$ defined in Eq. (32)).
In this work, we use the 3-momentum sphere intersection regularization,
implying that the integrals $A$ and $\mathcal{A}$ are different. The first,
$A$, within this regularization, is integrated over the region of the
intersection of one sphere centered in the origin and radius $\Lambda$, with
itself, making it identical to the usual 3-momentum regularization. The
second, $\mathcal{A}$, is a contribution to the quark-antiquark polarization
function, $\Pi_{ij}^{S|P}$, and its integral is calculated inside the region
of intersection between two spheres, distanced $\absolutevalue{\bm{k}}$ from
each other, just like the $B_{0}$ integral from the previous section.

Finite temperature can be introduced in the one fermion line integral, $A$
(and $\mathcal{A}$), in the same manner as for the two fermion, one loop
integral, $B_{0}$, see Section III. Applying these steps to Eq. (98) yields:

$\displaystyle A[M_{i},T,\mu_{i}]$
$\displaystyle=16\pi^{2}T\sum^{+\infty}_{n=-\infty}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{\quantity(i\omega_{n}+\mu_{i})^{2}-E\quantity[M_{i},\bm{p}]^{2}}.$
(97)

Here, $\mu_{i}$ is the chemical potential associated with the fermion of mass
$M_{i}$. Performing the infinite sum over the allowed Matsubara frequencies
yields the well known result:

$\displaystyle A[M_{i},T,\mu_{i}]$
$\displaystyle=-16\pi^{2}\int_{\mathrm{Reg}[\absolutevalue{\bm{k}}]}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{2E\quantity[M_{i},\bm{p}]}\quantity(1-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}]-\mu_{i},T]-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}]+\mu_{i},T])$
$\displaystyle=-4\int_{0}^{\Lambda}\frac{\differential{\absolutevalue{\bm{p}}}}{\sqrt{\bm{p}^{2}+M_{i}^{2}}}\quantity(1-f_{\mathrm{F}}\quantity[\sqrt{\bm{p}^{2}+M_{i}^{2}}-\mu_{i},T]-f_{\mathrm{F}}\quantity[\sqrt{\bm{p}^{2}+M_{i}^{2}}+\mu_{i},T]).$
(98)

Since the integrand only depends on $\absolutevalue{\bm{p}}$, the angular
integrations are trivial, leaving only one integration to be performed
numerically. From the above results, the divergence of the vacuum contribution
to the integral is evident and the necessity to regularize the integral is
explicit. In the vacuum, one gets the well known result:

$\displaystyle A[M_{i},0,0]$
$\displaystyle=-2\Lambda\sqrt{\Lambda^{2}+M_{i}^{2}}-M_{i}^{2}\ln\quantity[\frac{\quantity(\Lambda-\sqrt{\Lambda^{2}+M_{i}^{2}})^{2}}{M_{i}^{2}}].$
(99)

Now, we consider the one fermion line loop integral which arises in the second
order contribution, $\mathcal{A}$ (see Eq. (32)). In the 3-momentum sphere
intersection regularization scheme, we can write it as:

$\displaystyle\mathcal{A}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=-16\pi^{2}\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{2E\quantity[M_{i},\bm{p}]}\quantity(1-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}]-\mu_{i},T]-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}]+\mu_{i},T]).$
(100)

As already mentioned, the region of integration
${\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}$ is formed by the
intersection of two spheres, dislocated on the $p_{z}$ axis (see Fig. 3). In
order to simplify the calculation, one can shift the region of integration in
order to have the center of the intersection coinciding with the origin. This
can be done by performing the following change of variables:
$\bm{p}\to\bm{p}+\nicefrac{{\bm{k}}}{{2}}$. Using this transformation, the
integral becomes:

$\displaystyle\mathcal{A}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=-16\pi^{2}\int_{{\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential[3]{\bm{p}}}{\quantity(2\pi)^{3}}\frac{1}{2E\quantity[M_{i},\bm{p}+\frac{\bm{k}}{2}]}$
$\displaystyle\qquad\qquad\times\quantity(1-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}+\frac{\bm{k}}{2}]-\mu_{i},T]-f_{\mathrm{F}}\quantity[E\quantity[M_{i},\bm{p}+\frac{\bm{k}}{2}]+\mu_{i},T]).$
(101)

Here, ${{\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}}$ is the new
integration region constituted by the intersection of two spheres, one
centered at
$\quantity(p_{x},p_{y},p_{z})=\quantity(0,0,\nicefrac{{\absolutevalue{\bm{k}}}}{{2}})$
and other centered at
$\quantity(p_{x},p_{y},p_{z})=\quantity(0,0,\nicefrac{{-\absolutevalue{\bm{k}}}}{{2}})$,
with the respective centers distanced by $\absolutevalue{\bm{k}}$.
Interestingly, after performing the change of variables to center the region
of integration and integrating over the azimutal angle, we are left with an
extra dependence on the polar angle, $\theta$, making this integral slightly
more complicated to calculate when compared to its analogue with only one
sphere, $A$, given in Eq. (98). Using the steps shown in Appendix C.3 (see
Eqs. (272)-(281)), we can write this integral, in the vacuum limit, as:

$\displaystyle\mathcal{A}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=2\int_{0}^{\Lambda-\frac{\absolutevalue{\bm{k}}}{2}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{1}^{-1}\frac{\differential{\alpha}}{\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}+\bm{p}\absolutevalue{\bm{k}}\alpha+M_{i}^{2}}}$
$\displaystyle\qquad\qquad+2\int_{\Lambda-\frac{\absolutevalue{\bm{k}}}{2}}^{\sqrt{\Lambda^{2}-\frac{\bm{k}^{2}}{4}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{\frac{\Lambda^{2}-\frac{\bm{k}^{2}}{4}-\bm{p}^{2}}{\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}}^{-\frac{\Lambda^{2}-\frac{\bm{k}^{2}}{4}-\bm{p}^{2}}{\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}}\frac{\differential{\alpha}}{\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}+\bm{p}\absolutevalue{\bm{k}}\alpha+M_{i}^{2}}}$
$\displaystyle=-\frac{1}{3\absolutevalue{\bm{k}}}\Bigg{\\{}\bm{k}^{2}\quantity(\sqrt{(\Lambda-\absolutevalue{\bm{k}})^{2}+M_{i}^{2}}-3\sqrt{\Lambda^{2}+M_{i}^{2}})$
$\displaystyle\qquad\qquad+\absolutevalue{\bm{k}}\Lambda\quantity(\sqrt{(\Lambda-\absolutevalue{\bm{k}})^{2}+M_{i}^{2}}+3\sqrt{\Lambda^{2}+M_{i}^{2}})$
$\displaystyle\qquad\qquad+2\quantity(\Lambda^{2}+M_{i}^{2})\quantity(\sqrt{\Lambda^{2}+M_{i}^{2}}-\sqrt{(\Lambda-\absolutevalue{\bm{k}})^{2}+M_{i}^{2}})$
$\displaystyle\qquad\qquad+3\absolutevalue{\bm{k}}M_{i}^{2}\ln\quantity[\frac{1}{M_{i}^{2}}\quantity(\Lambda-\sqrt{\Lambda^{2}+M_{i}^{2}})\quantity(\Lambda-\absolutevalue{\bm{k}}-\sqrt{(\Lambda-\absolutevalue{\bm{k}})^{2}+M_{i}^{2}})]\Bigg{\\}}.$
(102)

In the $\absolutevalue{\bm{k}}\to 0$ case, the two spheres in the
regularization completely intersect and one gets the expected result given in
Eq. (99), i.e., $\mathcal{A}[M_{i},0,0,0]=A[M_{i},0,0]$. Of course, this
equivalence also holds at finite temperature and chemical potential as long as
$\absolutevalue{\bm{k}}=0$.

At finite temperature and chemical potential it is not possible to write this
quantity as an 1-dimensional integral using this approach. However, from the
two fermion line integral $B_{0}$, we know that it is possible to further
simplify this quantity, making an appropriate change of variables. For
simplicity, we will use exactly the same change of variables performed in the
case of the $B_{0}$ integral, with $\absolutevalue{\bm{k}}>0$. Thus,
considering the original form of the $\mathcal{A}$ integral (see Eq. (100)),
we multiply the integrand by one written as
$1=E\quantity[M_{j},\bm{p}-\bm{k}]/E\quantity[M_{j},\bm{p}-\bm{k}]=E_{2}/E_{2}$
(see Eq. (56)) and employ the change of variables used in Section III.1.1 (see
Eqs. (61)$-$(63)) in order to write:

$\displaystyle\mathcal{A}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{2}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{E}\differential{\varepsilon}\quantity(E+\frac{\varepsilon}{2})\quantity(\frac{1}{2}+f_{\mathrm{F}}\quantity[E-\frac{\varepsilon}{2}+\eta\mu_{i},T]).$
(103)

We highlight that this integral does not depend on the value of $M_{j}$
introduced when one multiplies the integrand by $E_{2}/E_{2}$. Indeed, any
value for $M_{j}$ can be used in the calculation and this apparent extra
parameter only arises in the calculation to allow us to use the formalism
developed earlier in order to simplify the calculations. By analogy with the
$B_{0}$ integral, we separate this integral into two contributions:

$\displaystyle\mathcal{A}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=\mathcal{A}_{\mathrm{p}}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]+\mathcal{A}_{\mathrm{s}}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
(104)

where each contribution is given by:

$\displaystyle\mathcal{A}_{\mathrm{p}}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{2}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{E}E\,G_{E}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},E],$
(105)
$\displaystyle\mathcal{A}_{\mathrm{s}}[M_{i},T,\mu_{i},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{2}{\absolutevalue{\bm{k}}}\sum_{\eta=\pm
1}\int_{\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{\varepsilon}\varepsilon\,G_{\varepsilon}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},\varepsilon].$
(106)

Here, we defined:

$\displaystyle G_{E}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},E]$
$\displaystyle=\int_{\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}}^{\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}\quantity(\frac{1}{2}+f_{\mathrm{F}}\quantity[E-\frac{\varepsilon}{2}+\eta\mu_{i},T]),$
(107) $\displaystyle
G_{\varepsilon}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},\varepsilon]$
$\displaystyle=\int_{E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}}^{E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}}\differential{E}\quantity(\frac{1}{2}+f_{\mathrm{F}}\quantity[E-\frac{\varepsilon}{2}+\eta\mu_{i},T]).$
(108)

As a result of making the appropriate change of variables, these integrals can
be evaluated analytically to yield:

$\displaystyle
G_{E}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},E]=\frac{1}{2}\quantity(\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}-\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})})$
$\displaystyle-2T\ln\quantity[\mathrm{e}^{\quantity(\frac{1}{2}\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}-E-\eta\mu_{i})/T}+1]$
$\displaystyle+2T\ln\quantity[\mathrm{e}^{\quantity(\frac{1}{2}\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}-E-\eta\mu_{i})/T}+1],$
(109) $\displaystyle
G_{\varepsilon}\quantity[M_{i},T,\mu_{i},\eta,\absolutevalue{\bm{k}},\varepsilon]=\frac{1}{2}\quantity(E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}-E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})})$
$\displaystyle+T\ln\quantity[\mathrm{e}^{\quantity(\frac{1}{2}\varepsilon-
E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}-\eta\mu_{i})/T}+1]$
$\displaystyle-T\ln\quantity[\mathrm{e}^{\quantity(\frac{1}{2}\varepsilon-
E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}-\eta\mu_{i})/T}+1].$
(110)

In this way, the one fermion line integral, in the 3-momentum sphere
intersection regularization, can be written as a sum of one dimensional
integrals. Taking the vacuum limit, one can also perform the remaining
integrations in order to obtain Eq. (102).

##  V Cross sections for quark-quark and quark-antiquark processes in the NJL
model

The goal of this section is to evaluate the cross sections of the NJL model
and show the influence of the regularization scheme proposed in this work
versus the 3-momentum regularization used in other works. We chose the cross
sections to show the differences between regularization schemes because they
constitute a physical quantity which is dependent on the absolute value of
external momentum, $\absolutevalue{\bm{k}}$, present in the quark-antiquark
polarization function $\Pi_{ij}^{S|P}$. As already discussed, the 3-momentum
sphere intersection regularization scheme is dependent on the value of
$\absolutevalue{\bm{k}}$ and, for $\absolutevalue{\bm{k}}=0$, one gets the
usual 3-momentum regularization. Evaluating thermodynamic quantities (at mean
field level, there are no external momenta), or meson masses (usually
calculated at rest with $\absolutevalue{\bm{k}}=0$) would not yield any
differences when compared to the usual 3-momentum regularization with a cutoff
in all integrals.

###  V.1 The model

In order to calculate the cross sections for quark-quark and quark-antiquark
scatterings, we will use the following $SU(3)$ NJL Lagrangian:

$\displaystyle\mathcal{L}=\overline{\psi}\quantity(i\not{\partial}-m)\psi+\frac{G}{2}\Big{(}\quantity(\overline{\psi}\lambda_{a}\psi)^{2}+(\overline{\psi}i\gamma^{5}\lambda_{a}\psi)^{2}\Big{)}+8\kappa\Big{(}\det\quantity[\overline{\psi}P_{R}\psi]+\det\quantity[\overline{\psi}P_{L}\psi]\Big{)}.$
(111)

Here, $\psi^{T}=\quantity(\psi_{u},\psi_{d},\psi_{s})$ and
$m=\operatorname{diag}\quantity{m_{u},m_{d},m_{s}}$ is the quark current mass
matrix. This popular model includes, besides the free Dirac Lagrangian, a
chirally symmetric set of interactions with a spontaneous chiral symmetry
broken ground state in the vacuum. The choice of model is rooted in the goal
to compare our results with other results found in the literature, where the
regularization of the one fermion and and two fermion line integrals are
performed in a different way than the one proposed here [11]. Hence, we only
include in the Lagrangian the necessary four quark-quark scalar-pseudoscalar
interaction and the ’t Hooft interaction [11, 1]. However, the formalism
developed earlier allows for a straightforward extension to include, at the
Lagrangian level, multi-scalar and multi-pseudoscalar quark interactions, such
as eight quark interactions or explicit chiral symmetry breaking interactions
[25, 56, 74, 27]. Indeed, in the future, we will consider an extended version
of the NJL model which includes eight-quark interactions in order to study the
influence of considering these higher quark interactions in the transport
properties of the NJL model.

The central motivation behind the NJL model is the incorporation of a
mechanism which induces spontaneous chiral symmetry breaking. In its original
form this is achieved through the inclusion of a local four-fermion contact
interaction term. One of the drawbacks of the model in this formulation is the
fact that it exhibits a axial $U_{A}(1)$ symmetry which is not observed in
Nature (there are no degenerate mirror hadrons with opposite parity). The fact
that no corresponding Nambu$-$Goldstone boson is observed led to the so called
“$U_{A}(1)$ problem” which can be solved by the anomalous breaking of this
symmetry through instanton effects [75, 76]. This mechanism can be
incorporated in the model with the inclusion of the ’t Hooft determinant
$2N_{f}-$fermion interaction term [77, 42, 26]. It should be noted that in the
three flavor version of the model this has been shown to lead to vacuum
stability problems which, in turn, can be solved with the inclusion of eight
quark interactions [22, 23]. In the present work, however, we will not
consider these in order to enable a straightforward comparison of our results
with the ones found in the literature.

Using the linear expansion of the Lagrangian, written in Eq. (111), as
explained in previous sections, one can readily obtain the gap equations for
this model:

$\displaystyle M_{i}$
$\displaystyle=m_{i}-2G\expectationvalue*{\overline{\psi}_{i}\psi_{i}}-2\kappa\expectationvalue*{\overline{\psi}_{j}\psi_{j}}\expectationvalue*{\overline{\psi}_{k}\psi_{k}}.$
(112)

Here, $i\neq j\neq k\in\\{u,d,s\\}$. The value of the quark effective masses,
at finite temperature and chemical potential, can be obtained by solving these
equations self-consistently and the quark condensates,
$\expectationvalue*{\overline{\psi}_{i}\psi_{i}}$, can be calculated from Eq.
(8). These gap equations can be obtained by calculating the thermodynamic
potential of the model, $\Omega$, and requiring it to be thermodynamically
consistent. This amounts to requiring the potential to be stationary with
respect to the mean fields, $\expectationvalue*{s_{a}}$, i.e.,
$\functionalderivative*{\Omega}{\expectationvalue*{s_{a}}}=0$.

As pointed out earlier, and has it will become clear in following sections, to
calculate the cross sections for different processes, one will also need to
evaluate different meson propagators at finite temperature and density. The
scalar and pseudoscalar meson propagators, $D_{S}$ and $D_{P}$, were defined
in Eqs. (34) and (35), respectively. The meson propagators are written in
terms of the interaction-independent quark-antiquark polarization function,
$\Pi_{ij}^{S|P}$, and of the interaction-dependent meson projectors $S_{ab}$
and $P_{ab}$. As discussed in previous sections, the polarization function is
independent of the type of interactions included at the Lagrangian level. As a
matter of fact we derived it in Section II.2, without making any choice
regarding the quark interaction Lagrangian. The meson projectors, on the other
hand, are completely dependent on the type of interactions considered at the
Lagrangian (as explained in Section II.2) and are defined as the coefficients
in front of the quadratic products,
$s_{a}s_{b}=(\overline{\psi}\lambda_{a}\psi)(\overline{\psi}\lambda_{b}\psi)$
and
$p_{a}p_{b}=(\overline{\psi}i\gamma^{5}\lambda_{a}\psi)(\overline{\psi}i\gamma^{5}\lambda_{b}\psi)$.
As discussed, these coefficients arise when one makes the quadratic expansion
of the product between quark bilinear operators (see Eq. (13)). The
interaction-dependent meson projectors, $S_{ab}$ and $P_{ab}$, for the
Lagrangian defined in Eq. (111), are given by:

$\displaystyle S_{ab}$
$\displaystyle=\frac{1}{2}\quantity(G\delta_{ab}+\frac{3\kappa}{2}A_{abc}\expectationvalue{s_{c}}),$
(113) $\displaystyle P_{ab}$
$\displaystyle=\frac{1}{2}\quantity(G\delta_{ab}-\frac{3\kappa}{2}A_{abc}\expectationvalue{s_{c}}).$
(114)

Here, $A_{abc}$ is given by:

$\displaystyle
A_{abc}=\frac{2}{3}d_{abc}+\sqrt{\frac{2}{3}}\quantity(\delta_{a0}\delta_{b0}\delta_{c0}-\delta_{a0}\delta_{bc}-\delta_{b0}\delta_{ca}-\delta_{c0}\delta_{ab}).$
(115)

The quantities $d_{abc}$ are the symmetric structure constants of the unitary
group $U(3)$. As previously introduced, the expectation value
$\expectationvalue*{s_{a}}=\expectationvalue*{\overline{\psi}\lambda_{a}\psi}$
and we consider that there is no pseudoscalar condensation, i.e.,
$\expectationvalue*{p_{a}}=\expectationvalue*{\overline{\psi}i\gamma^{5}\lambda_{a}\psi}=0$.
The non-zero scalar projectors, $S_{ab}$, are explicitly given by:

$\displaystyle{S}_{00}$
$\displaystyle=\frac{1}{2}\quantity(G+\frac{2}{3}\kappa\quantity(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+\expectationvalue*{\overline{\psi}_{d}\psi_{d}}+\expectationvalue*{\overline{\psi}_{s}\psi_{s}})),$
$\displaystyle S_{11}$
$\displaystyle=S_{22}=S_{33}=\frac{1}{2}\quantity(G-\kappa\expectationvalue*{\overline{\psi}_{s}\psi_{s}}),$
$\displaystyle S_{44}$
$\displaystyle=S_{55}=S_{66}=S_{77}=\frac{1}{2}\quantity(G-\kappa\expectationvalue*{\overline{\psi}_{d}\psi_{d}}),$
$\displaystyle S_{88}$
$\displaystyle=\frac{1}{2}\quantity(G-\frac{\kappa}{3}\quantity(2\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+2\expectationvalue*{\overline{\psi}_{d}\psi_{d}}-\expectationvalue*{\overline{\psi}_{s}\psi_{s}})),$
$\displaystyle S_{03}$
$\displaystyle=S_{30}=-\frac{\kappa}{2\sqrt{6}}\quantity(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}-\expectationvalue*{\overline{\psi}_{d}\psi_{d}}),$
$\displaystyle S_{08}$
$\displaystyle=S_{80}=-\frac{\sqrt{2}\kappa}{12}\quantity(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+\expectationvalue*{\overline{\psi}_{d}\psi_{d}}-2\expectationvalue*{\overline{\psi}_{s}\psi_{s}}).$
(116)

The non-zero pseudoscalar coefficients, $P_{ab}$, are:

$\displaystyle P_{00}$
$\displaystyle=\frac{1}{2}\quantity(G-\frac{2}{3}\kappa\quantity(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+\expectationvalue*{\overline{\psi}_{d}\psi_{d}}+\expectationvalue*{\overline{\psi}_{s}\psi_{s}})),$
$\displaystyle P_{11}$
$\displaystyle=P_{22}=P_{33}=\frac{1}{2}\quantity(G+\kappa\expectationvalue*{\overline{\psi}_{s}\psi_{s}}),$
$\displaystyle P_{44}$
$\displaystyle=P_{55}=P_{66}=P_{77}=\frac{1}{2}\quantity(G+\kappa\expectationvalue*{\overline{\psi}_{d}\psi_{d}}),$
$\displaystyle P_{88}$
$\displaystyle=\frac{1}{2}\quantity(G+\frac{\kappa}{3}\quantity(2\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+2\expectationvalue*{\overline{\psi}_{d}\psi_{d}}-\expectationvalue*{\overline{\psi}_{s}\psi_{s}})),$
$\displaystyle P_{03}$
$\displaystyle=P_{30}=\frac{\kappa}{2\sqrt{6}}(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}-\expectationvalue*{\overline{\psi}_{d}\psi_{d}}),$
$\displaystyle P_{08}$
$\displaystyle=P_{80}=\frac{\sqrt{2}\kappa}{12}\quantity(\expectationvalue*{\overline{\psi}_{u}\psi_{u}}+\expectationvalue*{\overline{\psi}_{d}\psi_{d}}-2\expectationvalue*{\overline{\psi}_{s}\psi_{s}}).$
(117)

These are exactly the same as the one obtained in Ref. [12] with a small
caveat: in the NJL Lagrangian density defined in Eq. (111), the four quark
scalar pseudoscalar interaction is proportional to $G/2$ while, in Ref. [12],
it is proportional to $G$.

At this point, one can explicitly write down the expressions for the meson
propagators, defined in Eqs. (34) and (35). The indexes of the inverse meson
propagators, $\quantity(D_{S}^{-1})_{ab}$ and $\quantity(D_{P}^{-1})_{ab}$,
run from 0 to 8, the number of generators of the $U(3)$ group. Thus, these
objects can be considered as two matrices of dimension $9\times 9$. It is
quite useful to write the scalar, $\Pi_{ab}^{S}$ and pseudoscalar,
$\Pi_{ab}^{P}$, polarization functions in the quark flavor basis, i.e., switch
from $a,b=\quantity{0,\ldots,8}$ to $\Pi_{ij}^{S|P}$, with
$i,j=\quantity{u,d,s}$. Using Eq. (29), we can write the scalar and
pseudoscalar polarization functions in the quark flavor basis. Explicitly, all
non-zero elements are given by (for the sake of simplicity, we omit the scalar
and pseudoscalar labels, ${S|P}$ ):

$\displaystyle\Pi_{00}$
$\displaystyle=\nicefrac{{2}}{{3}}\quantity(\Pi_{uu}+\Pi_{dd}+\Pi_{ss}),$
$\displaystyle\Pi_{03}$
$\displaystyle=\Pi_{30}=\sqrt{\nicefrac{{2}}{{3}}}\quantity(\Pi_{uu}-\Pi_{dd}),$
$\displaystyle\Pi_{08}$
$\displaystyle=\Pi_{80}=\nicefrac{{\sqrt{2}}}{{3}}\quantity(\Pi_{uu}+\Pi_{dd}-2\Pi_{ss}),$
$\displaystyle\Pi_{11}$ $\displaystyle=\Pi_{22}=\Pi_{ud}+\Pi_{du},$
$\displaystyle\Pi_{12}$
$\displaystyle=-\Pi_{21}=i\quantity(\Pi_{ud}-\Pi_{du}),$
$\displaystyle\Pi_{33}$ $\displaystyle=\Pi_{uu}+\Pi_{dd},$
$\displaystyle\Pi_{38}$
$\displaystyle=\Pi_{83}=\nicefrac{{1}}{{\sqrt{3}}}\quantity(\Pi_{uu}-\Pi_{dd}),$
$\displaystyle\Pi_{44}$ $\displaystyle=\Pi_{55}=\Pi_{us}+\Pi_{su},$
$\displaystyle\Pi_{45}$
$\displaystyle=-\Pi_{54}=i\quantity(\Pi_{us}-\Pi_{su}),$
$\displaystyle\Pi_{66}$ $\displaystyle=\Pi_{77}=\Pi_{ds}+\Pi_{sd},$
$\displaystyle\Pi_{67}$
$\displaystyle=-\Pi_{76}=i\quantity(\Pi_{ds}-\Pi_{sd}),$
$\displaystyle\Pi_{88}$
$\displaystyle=\nicefrac{{1}}{{3}}\quantity(\Pi_{uu}+\Pi_{dd}+4\Pi_{ss}).$
(118)

All the other elements are zero. The polarization functions written in terms
of the flavor indexes are not invariant under index exchange, i.e., $\Pi_{ud}$
and $\Pi_{du}$, $\Pi_{us}$ and $\Pi_{su}$ or $\Pi_{sd}$ and $\Pi_{ds}$ are not
necessarily equal at finite temperature and chemical potential. Indeed, it is
this difference which breaks the particle/antiparticle symmetry in the in-
medium meson behavior [10]. For instance the charged kaons behave differently
at finite chemical potential due to the mass difference between the quark and
anti-quark particles which constitute the kaons: at finite chemical potential,
the bath is more populated by particles than anti-particles.

Additionally, in this model, the inverse meson propagators are block matrices.
Thus, we can write the second order term of the quadratic expansion of the
effective action in Eq. (33) in different contributions. There will be a
sector related to the indexes $a,b=1,2$, $a,b=4,5$, $a,b=6,7$ and $a,b=0,3,8$.
In other models, with different physical considerations (like considering
pion/kaon condensation) or interaction terms at the Lagrangian level, this may
not be possible. The sector with indexes $a,b=1,2$ are related to the charged
pion propagators, $D_{\pi^{+}|\pi^{-}}$. Consider the matrices
$\lambda_{\pm}=\nicefrac{{1}}{{\sqrt{2}}}\quantity(\lambda_{1}\pm
i\lambda_{2})$ with $\lambda_{\pm}$ creating a quark bilinear with the quantum
numbers of the $\pi^{\mp}$, i.e.,
$\pi^{\pm}=\overline{\psi}\lambda_{\mp}\psi$. Using these identities one can
write the $\pi_{1|2}$ fields in terms of $\pi_{\pm}$. Thus, one can write the
$a,b=1,2$ sector of the quadratic effective action of the model (defined in
Eq. (33)) as:

$\displaystyle\mathcal{S}_{\mathrm{Q}}^{(2)}$
$\displaystyle=-\frac{1}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\bigg{\\{}\pi^{+}[k]\quantity(\frac{\quantity(D_{P}^{-1})_{11}}{2}+\frac{i\quantity(D_{P}^{-1})_{12}}{2}-\frac{i\quantity(D_{P}^{-1})_{21}}{2}+\frac{\quantity(D_{P}^{-1})_{22}}{2})\pi^{-}[-k]$
$\displaystyle\qquad\qquad\qquad\quad\;+\pi^{-}[k]\quantity(\frac{\quantity(D_{P}^{-1})_{11}}{2}-\frac{i\quantity(D_{P}^{-1})_{12}}{2}+\frac{i\quantity(D_{P}^{-1})_{21}}{2}+\frac{\quantity(D_{P}^{-1})_{22}}{2})\pi^{+}[-k]+\ldots\bigg{\\}}$
$\displaystyle=-\frac{1}{2}\int_{\mathrm{Reg}}\frac{\differential[4]{k}}{\quantity(2\pi)^{4}}\bigg{\\{}\pi^{+}[k]\quantity(D_{\pi^{-}}^{-1})\pi^{-}[-k]+\pi^{-}[k]\quantity(D_{\pi^{+}}^{-1})\pi^{+}[-k]+\ldots\bigg{\\}}.$
(119)

Here, we have defined the inverse propagators of the charged pions as:

$\displaystyle\quantity(D_{\pi^{+}|\pi^{-}})^{-1}$
$\displaystyle=\frac{\quantity(D_{P}^{-1})_{11}}{2}\mp\frac{i\quantity(D_{P}^{-1})_{12}}{2}\pm\frac{i\quantity(D_{P}^{-1})_{21}}{2}+\frac{\quantity(D_{P}^{-1})_{22}}{2}.$
(120)

In this model, the inverse model propagator with indexes $a,b=1,2$ is given,
in matrix form, by: Figure 8: Cross sections
($\quantity[\sigma]=\mathrm{mb}$), as a function of the center of mass energy