($\quantity[\sqrt{s}]=\mathrm{GeV}$), for the $uu\to uu$, $ud\to ud$, $us\to
us$ and $ss\to ss$ processes for $T=215\leavevmode\nobreak\ \mathrm{MeV}$ and
$T=250\leavevmode\nobreak\ \mathrm{MeV}$. On the left-hand side panel,
LABEL:sigma_uuuu_udud_Klev and 8(c), we use the usual 3-momentum
regularization proposed in Ref. [11] whereas, in the right-hand side panel,
LABEL:sigma_uuuu_udud_Yama and 8(d), we use the 3-momentum sphere intersection
regularization, proposed in this work. The large dots mark the minimum center
of mass energy requirements.

(c) (d)

(a)

(b)

(c) (d) Figure 9: Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as a
function of the center of mass energy ($\left[\sqrt{s}\right]=\mathrm{GeV}$),
for the $u\overline{d}\to u\overline{d}$, $u\overline{u}\to u\overline{u}$,
$u\overline{u}\to d\overline{d}$ and $u\overline{u}\to s\overline{s}$
processes, for $T=215\leavevmode\nobreak\ \mathrm{MeV}$ in the top and
$T=250\leavevmode\nobreak\ \mathrm{MeV}$ in the bottom panels. On the left-
hand side panels, 9(a) and 9(c), the usual 3-momentum regularization proposed
in Ref. [11] is used whereas, in the right-hand side panel, 9(b) and 9(d), the
3-momentum sphere intersection regularization, proposed in this work, is used.
The large dots mark the minimum center of mass energy requirements.

In Fig. V.1, we show the cross section for the $uu\to uu$, $ud\to ud$, $us\to
us$ and $ss\to ss$ quark-quark scattering processes, as a function of the
center of mass energy ($\sqrt{s}$), for the two temperatures mentioned before,
$T=215$ MeV and $T=250$ MeV. On the left-hand side panels,
LABEL:sigma_uuuu_udud_Klev and 8(c), we calculated the cross sections using
the usual 3-momentum regularization proposed in Ref. [11] while, on the right-
hand side panels, LABEL:sigma_uuuu_udud_Yama and 8(d), we used the approach
proposed in this work with the 3-momentum sphere intersection regularization.
From this result we can clearly observe that using the two different
regularization procedures result in very different cross sections for both
$uu\to uu$ and $ud\to ud$ processes. The reason behind this huge difference
can be traced back to the values of the mesons 3-momenta that are exchanged in
each channel, for each scattering process. Since we are considering the limit
of equal up and down quark masses, $M_{u}=M_{d}=M_{l}$, (and also equal
chemical potentials, $\mu_{u}=\mu_{d}=\mu_{s}=0$), for both scatterings, the
zero component of the mesons momentum in the $t$\- and $u$-channels is
$k_{0}=0$ while, the absolute value of the exchanged 3-momentum is
$\absolutevalue{\bm{k}}=\sqrt{-4M_{l}^{2}+s+t}$ (see Eqs.
(LABEL:k0_tChannel)$-$(LABEL:absk_uChannel)). Thus, this calculation is
dominated by evaluations of quark-antiquark polarization functions which have
a finite absolute value of external 3-momentum $\absolutevalue{\bm{k}}$. As
discussed in previous sections (and in the Appendix D), the different
regularization procedures lead to equal results if $\absolutevalue{\bm{k}}=0$
but differ when considering $\absolutevalue{\bm{k}}\geq 0$. Firstly, this is
caused because the quark-antiquark polarization function in the 3-momentum
sphere intersection regularization has contributions coming from the one
fermion, one loop integral, $\mathcal{A}$, regularized differently from the
integral $A$ used in the usual 3-momentum regularization. Secondly, as
explained before, in the usual 3-momentum regularization the two fermion, one
loop integral, $B_{0}$ is evaluated with a change of variables which breaks
certain symmetries of the integral when $\absolutevalue{\bm{k}}\geq 0$. In the
regularization proposed here, that is not the case, resulting in a different
$B_{0}$ function for finite $\absolutevalue{\bm{k}}$. The union between these
features explains the different behavior observed in this figure and others
that will be analyzed later. Focusing on the results with the 3-momentum
sphere intersection regularization (Fig. LABEL:sigma_uuuu_udud_Yama), we can
observe a minimum for both processes, right after the reaction threshold
energy. After this minimum, the cross sections increase, almost linearly with
the center of mass energy, with the cross section for the $ud\to ud$ process
being slightly larger, for the same temperature, when compared with the $uu\to
uu$ process due to the flavor factors that multiply the different meson
contributions, as noted in Ref. [12]. Furthermore, the effect of temperature
is to overall decrease the cross sections for a given value of center of mass
energy. We note that the minimum center of mass energy for each process (large
dots) also decrease with increasing temperature due to the larger restoration
of chiral symmetry. The effect of increasing temperature is known to restore
chiral symmetry and decrease the quark effective masses towards their bare
values and the threshold energies for the reaction to occur, on the center of
mass frame, depends only on the outgoing particles masses. Additionally, in
the $us\to us$ process considering the simple 3-momentum regularization, Fig.
8(c), there is a peak like structure not visible in the case of the new
3-momentum sphere intersection regularization, Fig. 8(d). In Ref. [12], this
structure was also observed and the authors claimed that its origin is a
cutoff artifact. Our calculation seems to more or less agree with such
statement: since we do not observe such peak, we can root the origin of such
structure to regularization artifacts but not specifically to the presence of
a cutoff since the new regularization procedure, proposed in this work, also
employs the use of a 3-momentum cutoff.

Next, we turn or attention to Fig. V.1 where we present the results for the
cross section of quark-antiquark processes with only light quarks in the
initial state. This requirement includes the following quark-antiquark
processes: $u\overline{d}\to u\overline{d}$, $u\overline{u}\to u\overline{u}$,
$u\overline{u}\to d\overline{d}$ and $u\overline{u}\to s\overline{s}$. Again,
we show results for $T=215\leavevmode\nobreak\ \mathrm{MeV}$ in the top
panels, $T=250\leavevmode\nobreak\ \mathrm{MeV}$ in the bottom panels whilst
comparing the results of using both the usual 3-momentum regularization, given
in the left panels (Figs. 9(a) and 9(c)), with the results of using the
regularization proposed in this work, shown in the right panels (Figs. 9(b)
and 9(d)). The most striking feature of these results is that both
regularization give almost exactly the same results, with very tiny
differences observed at larger center of mass energies, for both temperatures.
This can be traced back to the fact the quark-antiquark processes depend on
both the $s$\- and $t$-channels. In the case of meson propagators in the
$s$-channel, they are always evaluated with zero external 3-momentum,
$\absolutevalue{\bm{k}}^{(s)}=0$ and $k_{0}^{(s)}=\sqrt{s}$. Since both
regularization yield the same results at zero external 3-momentum, it implies
that all contributions to the differential cross section coming from the
$s$-channel are identical in both regularization schemes. In the case of meson
propagators in the $t$-channel, such is not the case since they are evaluated
for $k_{0}^{(t)}=0$ and $\absolutevalue{\bm{k}}^{(t)}=\sqrt{-t}$. However, the
integration range over $t$ for processes with $m_{1}=m_{2}$ and $m_{3}=m_{4}$,
are zero exactly at $s=s_{\mathrm{min}}$ and increase with increasing values
of center of mass energy, $\sqrt{s}$. Thus, for small values of $\sqrt{s}$
only small values of external 3-momentum are considered and only in the
$t$-channel, with the contributions from the $s$-channel being all for zero
external 3-momentum. So, this behavior explains why, for these processes, the
cross sections are almost identical in both regularization schemes, becoming
slightly different only with increasing center of mass energies. The large
peaks observed in this figures is due the fact that these cross sections are
being evaluated near the Mott temperatures of the pions and kaons, meaning
their behavior will be dominated by the existence of these resonances at
certain values of center of mass energy [12].

(a)

(b)

(c) (d) Figure 10: Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as a
function of the center of mass energy ($\left[\sqrt{s}\right]=\mathrm{GeV}$),
for the $u\overline{s}\to u\overline{s}$, $s\overline{s}\to u\overline{u}$ and
$s\overline{s}\to s\overline{s}$ processes, for $T=215\leavevmode\nobreak\
\mathrm{MeV}$ in the top and $T=250\leavevmode\nobreak\ \mathrm{MeV}$ in the
bottom panels. On the left-hand side panels, 10(a) and 10(c), the usual
3-momentum regularization proposed in Ref. [11] is used whereas, in the right-
hand side panel, 10(b) and 10(d), the 3-momentum sphere intersection
regularization, proposed in this work is used. The large dots mark the minimum
center of mass energy requirements.

Finally, in Fig. V.1, we present a comparison of cross sections for quark-
antiquark processes involving a strange antiquark in the initial state. The
processes shown in this figure are: $u\overline{s}\to u\overline{s}$,
$s\overline{s}\to u\overline{u}$ and $s\overline{s}\to s\overline{s}$. The
first observation, present in both regularization schemes is the diverging
nature of the $s\overline{s}\to u\overline{u}$ cross section at the energy
threshold. As pointed out in Ref. [12], this is a kinematic singularity caused
by the exothermic nature of this particular quark-antiquark process. Here, the
cross section for process $s\overline{s}\to s\overline{s}$ is quite different
for different regularizations but is extremely similar for the other two
processes. Again, the resonant structure observed is due to the proximity to
the Mott temperatures of the kaon resonances.

####  V.3.2 Results at finite temperature and chemical potential in the
3-momentum sphere intersection regularization scheme

In this section we show the results of evaluating quark-quark, quark-antiquark
and antiquark-antiquark processes at finite (fixed) temperature and
(increasing) quark chemical potential, as a function of the center of mass
energy, $\sqrt{s}$, allowing the study of the influence of finite density
effects on the cross sections. For the temperature we used the previous value
of $T=250\leavevmode\nobreak\ \mathrm{MeV}$ while, for the chemical potential,
we used 6 different values,
$\mu=\quantity{0,100,200,300,400,500,600}\leavevmode\nobreak\ \mathrm{MeV}$.
Additionally, we consider equal quark chemical potentials,
$\mu=\mu_{u}=\mu_{d}=\mu_{s}$. In this case we will have 16 different
processes instead of the 11 different processes at $\mu=0$. The extra
processes that have to be calculated are:
$\overline{u}\,\overline{u}\to\overline{u}\,\overline{u}$,
$\overline{u}\,\overline{d}\to\overline{u}\,\overline{d}$,
$\overline{u}\,\overline{s}\to\overline{u}\,\overline{s}$,
$\overline{s}\,\overline{s}\to\overline{s}\,\overline{s}$ and
$s\overline{u}\to s\overline{u}$. These processes are degenerate to their
particle analogues at zero chemical potential (when switching quarks to
antiquarks and antiquark to quarks) however, when considering finite chemical
potential, they will be different.

(a)

(b)

(c) (d) Figure 11:  Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as a
function of the center of mass energy ($\left[\sqrt{s}\right]=\mathrm{GeV}$),
for the $uu\to uu$, $\overline{u}\,\overline{u}\to\overline{u}\,\overline{u}$,
$ud\to ud$, $\overline{u}\,\overline{d}\to\overline{u}\,\overline{d}$
processes, for $T=250\leavevmode\nobreak\ \mathrm{MeV}$ and increasing
chemical potentials. The large dots mark the minimum center of mass energy
requirements for a given processes to occur. Here, the 3-momentum sphere
intersection regularization is used.

(a)

(b)

(c) (d) Figure 12:  Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as a
function of the center of mass energy ($\left[\sqrt{s}\right]=\mathrm{GeV}$),
for the $us\to us$, $\overline{u}\,\overline{s}\to\overline{u}\,\overline{s}$,
$ss\to ss$, $\overline{s}\,\overline{s}\to\overline{s}\,\overline{s}$
processes, for $T=250\leavevmode\nobreak\ \mathrm{MeV}$ and increasing
chemical potentials. The large dots mark the minimum center of mass energy
requirements for a given processes to occur. Here, the 3-momentum sphere
intersection regularization is used.

We start by presenting all the quark-quark and antiquark-antiquark cross
sections in Figs. V.3.2 and V.3.2. When considering finite chemical potential
and, consequently, finite baryon density, the medium is characterized by an
asymmetry between the number of particles and antiparticles. Such asymmetry is
responsible, for instance, for the difference observed in the meson behavior
of certain meson states at finite chemical potential and the masses of the
charged kaons and their scalar particles are different in the presence of non-
zero baryon density [10]. Since the cross sections are highly dependent on the
meson propagators and, taking into account that finite chemical potential
already breaks the symmetry between some charged mesons, it is natural to
expect it to also break the symmetry between processes which, at zero chemical
potential, are identical. Hence, in these figures, one can observe the effect
of breaking the symmetry between charged mesons at finite density. Of course,
the observed differences between quark-quark and antiquark-antiquark at finite
chemical potential is also connected to the Pauli blocking factor considered
to evaluate the cross section, see Eq. (LABEL:cross_section_definition).

In Fig. 11(a) and 11(b) we show the results for the processes $uu\to uu$ and
$\overline{u}\,\overline{u}\to\overline{u}\,\overline{u}$. As expected, at
zero chemical potential, these processes yield the same result however, at
finite chemical potential, they have different behaviors. As the chemical
potential increases, the cross section for the $uu\to uu$ process consistently
decreases across all values of the center of mass energy. Conversely, for the
$\overline{u}\,\overline{u}\to\overline{u}\,\overline{u}$ process, the cross
section demonstrates an increase with chemical potential for large center of
mass energy values. However, for small energies, the behavior of the cross
section becomes more complex and does not follow a straightforward trend. The
$ud\to ud$ cross section and
$\overline{u}\,\overline{d}\to\overline{u}\,\overline{d}$ cross sections, as
well as the for the $us\to us$ and
$\overline{u}\,\overline{s}\to\overline{u}\,\overline{s}$ processes, shown in
Figs. 11(c), 11(d), 12(a) and 12(b) respectively, share the same overall
behavior that was previously described about processes involving only the $u$
quark.

The minimum center of mass energies for these processes to occur, represented
by the colored dots in Fig. V.3.2, are independent of the considered process
and are equal for a specific chemical potential. This is a consequence of the
isospin symmetry considered in this study which yields equal effective quark
masses for the up and down quarks. Naturally, processes involving the strange
quark have a higher threshold energy due to the higher strange quark effective
mass. As expected from the restoration of chiral symmetry in the medium, the
threshold energies for all processes gets smaller with increasing chemical
potential since the quarks effective masses tend to their respective bare
values. This feature is common to all results at finite chemical potential.

We highlight that, in our calculation, the quark-quark process $us\to us$ and
the equivalent process $su\to su$ are identical. From the physics point of
view this is obvious however, from the practical point of view this is only
possible because the 3-momentum sphere intersection regularization does not
ruin this feature, as explained before. If one uses the usual 3-momentum
regularization, at finite chemical potential, the cross sections for the
processes $us\to us$ and $su\to su$ will be different.

(a)

(b)

(c)

(d) (e) (f) Figure 13:  Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as
a function of the center of mass energy
($\left[\sqrt{s}\right]=\mathrm{GeV}$), for the $u\overline{u}\to
u\overline{u}$, $u\overline{d}\to u\overline{d}$, $u\overline{u}\to
d\overline{d}$, $u\overline{u}\to s\overline{s}$, $s\overline{s}\to
u\overline{u}$ and $s\overline{s}\to s\overline{s}$ processes, for
$T=250\leavevmode\nobreak\ \mathrm{MeV}$ and increasing chemical potentials.
The large dots mark the minimum center of mass energy requirements for a given
processes to occur. Here, the 3-momentum sphere intersection regularization is
used.

(a)

(b)

Figure 14:  Cross sections ($\left[\sigma\right]=\mathrm{mb}$), as a function
of the center of mass energy ($\left[\sqrt{s}\right]=\mathrm{GeV}$), for the
$u\overline{s}\to u\overline{s}$ and $s\overline{u}\to s\overline{u}$
processes, for $T=250\leavevmode\nobreak\ \mathrm{MeV}$ and increasing
chemical potentials. The large dots mark the minimum center of mass energy
requirements for a given processes to occur. Here, the 3-momentum sphere
intersection regularization is used.

We now turn our attention to the quark-antiquark cross sections. In Fig. V.3.2
we have the processes with zero net strangeness,which include:
$u\overline{u}\to u\overline{u}$, $u\overline{d}\to u\overline{d}$,
$u\overline{u}\to d\overline{d}$, $u\overline{u}\to s\overline{s}$,
$s\overline{s}\to u\overline{u}$ and $s\overline{s}\to s\overline{s}$. The
cross sections of quark-antiquark processes with strangeness are presented in
Fig. 14, and include: $u\overline{s}\to u\overline{s}$ and $s\overline{u}\to
s\overline{u}$ processes. As in the case of zero chemical potential, several
of these quantities display a resonance structure due to the proximity of the
values of temperature and chemical potential to the Mott temperature of the
pions and kaons. Of course, when considering finite density, for every
chemical potential one can find the respective Mott temperature. Thus, when
considering in-medium systems, the melting point of the meson bound states is
defined by a point in the space of possible thermodynamic variables, in this
case, there will be a Mott temperature and a Mott chemical potential. The
study of meson behavior in the medium, within the context of the model used
here, is beyond the scope of this work and is left as future work. However, in
the range of chemical potentials considered here, the Mott temperature of
pions and kaons decreases with increasing temperature. Such effect can be
clearly seen in some of the cross sections of Figs. V.3.2 and 14. The peak
caused by the close proximity to the Mott temperatures of the pion and kaons
gets smaller with increasing chemical potential. This can be observed in the
cross sections of the following processes: $u\overline{u}\to u\overline{u}$,
$u\overline{d}\to u\overline{d}$, $u\overline{u}\to d\overline{d}$,
$u\overline{s}\to u\overline{s}$ and $s\overline{u}\to s\overline{u}$.

The kinematic divergence present in the exothermic cross section of the
process $s\overline{s}\to u\overline{u}$ is present for all considered
chemical potentials however, it moves to smaller values of $\sqrt{s}$. The
reduction of the strange quark mass, driven by the mechanism of restoration of
chiral symmetry, is responsible for this behavior since the threshold energy
for this process to occur gets smaller with increasing chemical potential.

The cross section for the $u\overline{s}\to u\overline{s}$ and
$s\overline{u}\to s\overline{u}$ are very similar. As mentioned, they also
display a peak structure related to the proximity to meson melting points.
However, the peaks are not located at the same points. These processes, in the
$s$-channel, are mediated by the charged kaons, see Table
LABEL:s_channel_prop_qqbar. While $u\overline{s}\to u\overline{s}$ involves
the positively charged $K^{+}$ , the process $s\overline{u}\to s\overline{u}$
exchanges a negatively charged $K^{-}$. Again, at finite density these meson
are not degenerate, displaying different propagators, masses, widths and
consequently, different melting points (labeled by a Mott chemical potential
and temperature). Thus, in each process, the peak structure is located near
its respective charged kaon melting point and their difference in essence,
amounts to the charged kaons possessing slightly different melting points as
chemical potential increases.

It is expected that the analyzing the effective meson-quark couplings at
finite temperature and chemical potential [5], could provide additional
information regarding the fact that cross sections for quark-quark processes
gets smaller with chemical potential while, in general, the antiquark-
antiquark cross sections gets larger. However, such study goes beyond the
scope of the present work.

##  VI Conclusions

In this work we derived the one and two fermion line integrals from a general
NJL model. We considered that the NJL Lagrangian contained the usual Dirac
contribution for fermions fields and different scalar-pseudoscalar quark-quark
interactions, contained in $\mathcal{L}_{\mathrm{int}}$. The formalism
developed here is applicable to any number of scalar-pseudoscalar multi-quark
interactions. For instance, one could consider the usual four quark scalar-
pseudoscalar interaction, the ’t Hooft interaction, 8-quark interactions or
explicit chiral symmetry breaking interactions. In order to derive the one and
two fermion line integrals, we performed the linear and quadratic expansions
of the product between quark bilinear operators in the mean field
approximation. From these expansions one can extract thermodynamics quantities
and meson propagators. The latter were then used to evaluate quark-quark and
quark-antiquark cross sections at finite temperature and chemical potential.
Furthermore, the meson propagators obtained in this work are equivalent to the
ones used in other works where they are derived from the Bethe-Salpeter
equation in the random phase approximation.

In order to perform these tasks, we introduced a new approach to control the
diverging nature of the model, the so-called 3-momentum sphere intersection
regularization. This new procedure is essential in order to maintain the
symmetries of the model. For instance, it ensures that the cross sections for
the processes $us\to us$ and $su\to su$ are the same for all values of
temperature and chemical potential, which does not happen in the usual
3-momentum regularization.

As a practical implementation of the developed formalism, and to demonstrate
the differences between the 3-momentum sphere intersection regularization and
the usual 3-momentum regularization, we considered the standard $SU(3)$ NJL
model, with four and six quark interactions [1], to evaluate quark-quark and
quark-antiquark cross sections. We observed quantitative and qualitative
differences when comparing quark-quark cross sections in both schemes. For
quark-antiquark cross sections however, the results are similar. This is
expected since, from the kinematic point of view, the differential cross
sections for quark-antiquark processes do not have major dependencies on the
absolute value of the external 3-momentum, $\absolutevalue{\bm{k}}$, meaning
that both regularization schemes will give very similar results. This is a
feature of the 3-momentum sphere intersection regularization: it depends on
the absolute value of the external 3-momentum, $\absolutevalue{\bm{k}}$. Since
we are regularizing the integral by evaluating it in the region of integration
corresponding to the intersection of two spheres with centers distanced
$\absolutevalue{\bm{k}}$ and radii $\Lambda$, when $\absolutevalue{\bm{k}}=0$,
both spheres completely intersect and we end up with the usual 3-momentum
regularization.

Thus, not only this new approach to the 3-momentum regularization scheme
yields different quantitative and qualitative results for the quark-quark
cross sections, it also provides a way to correctly account for the symmetries
of the system. These points are very important when using the cross sections
evaluated from this model to obtain the quark relaxation time and use it to
predict the in-medium transport properties of the quark plasma. The impact of
the regularization scheme proposed in this work on the evaluation of transport
coefficients like shear/bulk viscosities and thermal/electric conductivities
is currently being analyzed and being prepared for future publications.

As future work, we would like to use the Pauli-Villars regularization in order
to regularize the model. We are specially interested on studying the
differences to the cross sections, between using the 3-momentum sphere
intersection regularization, used in this work, and applying the Pauli-Villars
regularization. As discussed, it is extremely important to use a
regularization procedure which, like the one presented here, is able to
preserve the symmetries between meson and antimesons.

Another extension of this work is to include the effect of magnetic fields in
the calculation of the one and two fermion line integrals. Such extension can
be very important to study systems under the presence of intense magnetic
fields such as the case of heavy ion collision experiments or to study the
matter inside magnetars. In such case, the Pauli-Villars regularization may be
an extremely useful approach in order to study the effects of an external
magnetic field. When considering a finite magnetic field, the 3-momentum
integration is broken down in to an integration over the momenta aligned with
the magnetic field, usually $p_{z}$, and the two dimensional integration over
the transverse momenta, $p_{x}-p_{y}$, which is replaced by a sum over
discrete values, the Landau levels. If considering a 3-momentum cutoff, the
transverse component of the momentum should also be restricted, i.e., cutting
off the sum. However, when dealing with an integration from 0 to infinity, as
is the case in the Pauli-Villars regularization, this sum can be performed
analytically, greatly simplifying the calculations.

We also wish to include in the model a coupling of the quark fields to an
external static background of gluons. In such extended model, one is able to
also study the deconfinement transition by studying the Polyakov loop order
parameter at finite temperature and chemical potential. Such improved models
are named Polyakov$-$Nambu$-$Jona-Lasinio (PNJL) models. In this work we used
the pion Mott temperature as an indication for the system to be in the regime
of deconfined matter. However, using the PNJL model would allow for a better
description of the deconfinement transition and a better qualitative and
quantitative indication of the presence of deconfined matter. Not only that,
the PNJL is known to yield better qualitative results when compared to lattice
QCD calculations, yielding a better qualitative prediction of the quark-quark
and quark-antiquark cross sections. One can even use the so-called entangled-
PNJL model where the coupling constants of the model are also dependent on the
Polyakov loop order parameter.

Finally, the formalism developed here can also be extended to allow for the
inclusion of pseudoscalar condensation at finite temperature and chemical
potential, e.g., pion and kaon condensation [81, 82]. In asymmetric matter,
such as the one created in heavy-ion collisions experiments or present inside
neutron stars, the isospin chemical potential can be sufficiently large for
meson condensation to occur [83, 84, 85, 86, 82].

## Acknowledgments

This work was supported by a research grant under Project No. PTDC/FIS-
NUC/29912/2017, funded by national funds through FCT (Fundação para a Ciência
e a Tecnologia, I.P, Portugal) and co-financed by the European Regional
Development Fund (ERDF) through the Portuguese Operational Program for
Competitiveness and Internationalization, COMPETE 2020, by the project
CERN/FIS-PAR/0040/2019 and by national funds from FCT, within the Projects No.
UIDB/04564/2020 and No. UIDP/04564/2020.

##  Appendix A Linear expansion of the product between $N$ commuting operators

An operator, $\hat{\mathcal{O}}_{i}$, can be written as its own expectation
value, plus a small perturbation $\delta\hat{\mathcal{O}}_{i}$ around it [28]:

$\displaystyle\hat{\mathcal{O}}_{i}$
$\displaystyle=\expectationvalue*{\hat{\mathcal{O}}_{i}}+\quantity(\hat{\mathcal{O}}_{i}-\expectationvalue*{\hat{\mathcal{O}}_{i}})=\expectationvalue*{\hat{\mathcal{O}}_{i}}+\delta\hat{\mathcal{O}}_{i}.$
(198)

The perturbation is defined as:

$\displaystyle\delta\hat{\mathcal{O}}_{i}=\hat{\mathcal{O}}_{i}-\expectationvalue*{\hat{\mathcal{O}}_{i}}.$
(199)

Using these expressions, we can write the product between two commuting
operators, $\hat{\mathcal{O}}_{1}$ and $\hat{\mathcal{O}}_{2}$, as:

$\displaystyle\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}$
$\displaystyle=\quantity(\expectationvalue*{\hat{\mathcal{O}}_{1}}+\delta\hat{\mathcal{O}}_{1})\quantity(\expectationvalue*{\hat{\mathcal{O}}_{2}}+\delta\hat{\mathcal{O}}_{2})$
$\displaystyle=\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\delta\hat{\mathcal{O}}_{2}+\delta\hat{\mathcal{O}}_{1}\expectationvalue*{\hat{\mathcal{O}}_{2}}+\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}.$
(200)

The so-called linear product between these two operators,
$\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]$, is
defined by neglecting the quadratic term
$\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}$ and then, re-writing
the perturbation, $\delta\hat{\mathcal{O}}_{i}$, using its definition (see Eq.
(199)). Thus, one can write:

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

The process to define the linear expansion of the product between three
commuting operators is exactly the same. Using Eq. (198), we can write the
product between three operators as:

$\displaystyle\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}$
$\displaystyle=\quantity(\expectationvalue*{\hat{\mathcal{O}}_{1}}+\delta\hat{\mathcal{O}}_{1})\quantity(\expectationvalue*{\hat{\mathcal{O}}_{2}}+\delta\hat{\mathcal{O}}_{2})\quantity(\expectationvalue*{\hat{\mathcal{O}}_{3}}+\delta\hat{\mathcal{O}}_{3})$
$\displaystyle=\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}+\expectationvalue*{\hat{\mathcal{O}}_{2}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\delta\hat{\mathcal{O}}_{1}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{3}}\delta\hat{\mathcal{O}}_{2}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\expectationvalue*{\hat{\mathcal{O}}_{2}}\delta\hat{\mathcal{O}}_{3}$
$\displaystyle+\expectationvalue*{\hat{\mathcal{O}}_{3}}\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}+\expectationvalue*{\hat{\mathcal{O}}_{2}}\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{3}+\expectationvalue*{\hat{\mathcal{O}}_{1}}\delta\hat{\mathcal{O}}_{2}\delta\hat{\mathcal{O}}_{3}+\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}\delta\hat{\mathcal{O}}_{3}.$
(202)

In this case, one neglects the terms proportional to both
$(\delta\hat{\mathcal{O}})^{2}$ and $(\delta\hat{\mathcal{O}})^{3}$. Then,
using Eq. (199), one gets:

$\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}.$
(203)

This procedure can be performed to products with any number of operators.
Thus, the linear product between $N$ operators can be obtained, by the three
following operations: firstly, writing each operator as in Eq. (198),
secondly, neglecting terms proportional to $(\delta\hat{\mathcal{O}})^{n}$
with $n\geq 2$, thirdly, writing the perturbations using the definition stated
in Eq. (199).

Interestingly, one can obtain the linear product between $N$ operators
iteratively, by elevating the linear product between two operators,
$\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]$
(explicitly given in Eq. (201)), to a definition. For completeness, one must
also consider that the linear expansion of an operator or a constant (products
of mean fields) is the operator or constant itself,
$\mathrm{L}\left[\hat{\mathcal{O}}\right]=\hat{\mathcal{O}}$ and
$\mathrm{L}\left[\mathrm{cte}\right]=\mathrm{cte}$. Using these definitions,
one can build the linear expansion of any number of operators by iteratively
using the following defining property:

$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N}\right]=\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N-1}\right]\hat{\mathcal{O}}_{N}\right].$
(204)

Here, the particular set of $N-1$ operators chosen to be inside the nested
linear operation does not matter. Recursively, one can continue breaking down
the innermost linear product, creating an iterative set of operations until
the defining case of two operators is reached:
$\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]$.
Explicitly, one can write:

$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N}\right]=\mathrm{L}\left[\mathrm{L}\left[\mathrm{L}\left[\ldots\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{4}\ldots\hat{\mathcal{O}}_{N-1}\right]\hat{\mathcal{O}}_{N}\right].$
(205)

For instance, following these definitions, the linear expansion of the product
between three commuting operators is given by:

$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]$
$\displaystyle=\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]\hat{\mathcal{O}}_{3}\right]=\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{2}\right]=\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{1}\right].$
(206)

Here, it is highlighted the fact that the particular combination of two
operators chosen in the inner linear operation does not matter. By
substituting the definition of
$\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]$ in the
above (see Eq. (201)), one recovers the linear product between three operators
obtained in Eq. (203). For instance, for four operators, one can write:

$\displaystyle\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\hat{\mathcal{O}}_{4}\right]$
$\displaystyle=\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{4}\right]$
$\displaystyle=\mathrm{L}\left[\mathrm{L}\left[\mathrm{L}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{4}\right].$
(207)

Here, we only wrote one particular choice for the nested linerization
operations. There are several other possible permutations among the operators
that give the same result.

Due to its recursive nature, applying the definition given in Eq. (204) to get
the linear product between $N$ operators can be quite cumbersome. Fortunately,
the linear product between $N=n+1$ operators, with $n\geq 1$, can also be
obtained via the following formula:

$\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}}.$
(208)

This formula can be proved by mathematical induction. First, one considers the
base case and show that the proposition is valid for $n=1$. Substituting
$n=1$, in Eq. (208), it yields the linear 2-operator product given in Eq.
(201). This establishes the validity of the formula for the $n=1$ case. The
next step is to assume the validity of the formula written in Eq. (208) for a
given $n=m$, and then show that it also holds for the case $n=m+1$. Thus, one
has to show that the following equality holds:

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

The right-hand side of this equation can be written as:

$\displaystyle\mathrm{RHS}$
$\displaystyle=\quantity(\sum_{i=1}^{m+1}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}+\frac{\hat{\mathcal{O}}_{m+2}}{\expectationvalue*{\hat{\mathcal{O}}_{m+2}}}-\quantity(m+1))\prod_{j=1}^{m+2}\expectationvalue*{\hat{\mathcal{O}}_{j}}.$
(210)

Now, consider the left-hand side of Eq. (209). Writing the $m+2$ contribution
in the product explicitly and using the defining property of the linear
product between $N$ operators given in Eq. (204), we can write the left-hand
side as:

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

Here, we assumed that Eq. (208) is valid for some $m$ and used the fact that
$\mathrm{L}\left[\hat{\mathcal{O}}_{m+2}\right]=\hat{\mathcal{O}}_{m+2}$. The
first term can be expanded using the base case, i.e., the result for the
linear product of two operators given in Eq. (201). We write:

$\displaystyle\mathrm{LHS}$
$\displaystyle=\quantity(\sum_{i=1}^{m+1}\frac{1}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}\quantity(\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}+\frac{\hat{\mathcal{O}}_{m+2}}{\expectationvalue*{\hat{\mathcal{O}}_{m+2}}}-1)\expectationvalue*{\hat{\mathcal{O}}_{i}}\expectationvalue*{\hat{\mathcal{O}}_{m+2}}-m\hat{\mathcal{O}}_{m+2})\prod_{j=1}^{m+1}\expectationvalue*{\hat{\mathcal{O}}_{j}}$
$\displaystyle=\quantity(\sum_{i=1}^{m+1}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}+\frac{\hat{\mathcal{O}}_{m+2}}{\expectationvalue*{\hat{\mathcal{O}}_{m+2}}}-\quantity(m+1))\prod_{j=1}^{m+2}\expectationvalue*{\hat{\mathcal{O}}_{j}}$
$\displaystyle=\mathrm{RHS}.$ (212)

Since both sides of Eq. (209) are equal, the formula also holds for $n=m+1$,
completing the inductive step. Hence, by mathematical induction, Eq. (208) is
valid for $n\geq 1$.

##  Appendix B Quadratic expansion of the product between $N$ commuting
operators

Similarly to the linear product between $N$ commuting operators, the quadratic
product between $N$ commuting operators can be be obtained by writing each
operator as in Eq. (198), neglecting terms proportional to cubic perturbations
or higher, $(\delta\hat{\mathcal{O}})^{n}$ for $n\geq 3$, and then re-writing
the perturbations using Eq. (199).

Consider the product between three operators, $\hat{\mathcal{O}}_{1}$,
$\hat{\mathcal{O}}_{2}$ and $\hat{\mathcal{O}}_{3}$. It was explicitly written
out in Eq. (202). Neglecting the term
$\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}\delta\hat{\mathcal{O}}_{3}$
and re-writing the perturbations, yields its quadratic expansion:

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

Following the same steps, the quadratic product between four operators can be
written as:

$\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}}.$
(214)

In order to write this expression, the terms proportional to both
$\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}\delta\hat{\mathcal{O}}_{3}$
and
$\delta\hat{\mathcal{O}}_{1}\delta\hat{\mathcal{O}}_{2}\delta\hat{\mathcal{O}}_{3}\delta\hat{\mathcal{O}}_{4}$
were neglected, after each operator in the product
$\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\hat{\mathcal{O}}_{4}$
was written using Eq. (198).

In a similar way to the linear expansion of the product between $N$ operators
described in the previous section, one can obtain the quadratic product
between $N$ operators iteratively, by taking the quadratic product between
three operators,
$\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]$,
as a definition. Additionally, one must consider the following statements:
$\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\right]=\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}$,
$\mathrm{Q}\left[\hat{\mathcal{O}}\right]=\hat{\mathcal{O}}$ and
$\mathrm{Q}\left[\mathrm{cte}\right]=\mathrm{cte}$. Using these constraints,
the quadratic expansion of the product of any number of operators can be
obtained recursively using the following defining property:

$\displaystyle\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N}\right]=\mathrm{Q}\left[\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N-1}\right]\hat{\mathcal{O}}_{N}\right].$
(215)

The particular set of $N-1$ operators, chosen to be inside the nested
quadratic operation, does not matter. Recursively, one can continue breaking
down the innermost quadratic product, creating an iterative set of operations
until the defining case is reached:
$\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]$.
Explicitly, one can write:

$\displaystyle\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\ldots\hat{\mathcal{O}}_{N}\right]=\mathrm{Q}\left[\mathrm{Q}\left[\mathrm{Q}\left[\ldots\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{4}\right]\hat{\mathcal{O}}_{5}\ldots\hat{\mathcal{O}}_{N-1}\right]\hat{\mathcal{O}}_{N}\right].$
(216)

As an example, the quadratic product between four operators, can be obtained
by writing:

$\displaystyle\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\hat{\mathcal{O}}_{4}\right]$
$\displaystyle=\mathrm{Q}\left[\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]\hat{\mathcal{O}}_{4}\right].$
(217)

Substituting
$\mathrm{Q}\left[\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}\hat{\mathcal{O}}_{3}\right]$
by its definition, given in Eq. (213), one recovers the expression given in
Eq. (214).

Instead of using the recursive approach laid out above, the quadratic product
between $N=n+2$ operators, with $n\geq 1$, can be conveniently obtained using
the following formula:

$\displaystyle\mathrm{Q}\left[\prod_{i=1}^{n+2}\hat{\mathcal{O}}_{i}\right]$
$\displaystyle=\quantity(\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))\prod_{k=1}^{n+2}\expectationvalue*{\hat{\mathcal{O}}_{k}}.$
(218)

As for the linear product expansion, this formula can be proved by
mathematical induction. First we consider the base case and show that the
expression is valid for $n=1$, i.e., for the quadratic expansion of the
product between three operators. Using this formula one can readily obtain the
result presented in Eq. (213), validating the formula for the base case with
$n=1$. Next, assuming that the formula holds for $n=m$, one must demonstrate
that it is also valid for $n=m+1$. Hence, we have to show that the following
relations holds ($n=m+1$ case):

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

Let us focus on the right-hand side of Eq. (219). Separating the $m+3$
contribution in the sums, we can write:

$\displaystyle\mathrm{RHS}=$
$\displaystyle\bigg{(}\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+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})+\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}$
$\displaystyle\qquad-\quantity(m+1)\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}-\quantity(m+1)\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}+\frac{\quantity(m+1)}{2}\quantity(m+2)\bigg{)}\prod_{k=1}^{m+3}\expectationvalue*{\hat{\mathcal{O}}_{k}}.$
(220)

Consider the left-hand side of Eq. (219). One can separate the $i=m+3$
contribution to the product and, using the defining property of Eq. (215),
write:

$\displaystyle\mathrm{LHS}$
$\displaystyle=\mathrm{Q}\left[\quantity(\prod_{i=1}^{m+2}\hat{\mathcal{O}}_{i})\hat{\mathcal{O}}_{m+3}\right]$
$\displaystyle=\mathrm{Q}\left[\mathrm{Q}\left[\prod_{i=1}^{m+2}\hat{\mathcal{O}}_{i}\right]\hat{\mathcal{O}}_{m+3}\right]$
$\displaystyle=\Bigg{(}\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+2}\frac{\mathrm{Q}\left[\hat{\mathcal{O}}_{i}\hat{\mathcal{O}}_{j}\hat{\mathcal{O}}_{m+3}\right]}{\expectationvalue*{\hat{\mathcal{O}}_{i}}\expectationvalue*{\hat{\mathcal{O}}_{j}}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\quantity(1-\delta_{ij})$
$\displaystyle\qquad\qquad\qquad-m\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}+\frac{m}{2}\quantity(m+1)\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\Bigg{)}\prod_{k=1}^{m+3}\expectationvalue*{\hat{\mathcal{O}}_{k}}.$
(221)

Here, we applied the base case (see Eq.(218)). Defining the first term inside
the brackets by $t$, we can write it as:

$\displaystyle t$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+2}\frac{\mathrm{Q}\left[\hat{\mathcal{O}}_{i}\hat{\mathcal{O}}_{j}\hat{\mathcal{O}}_{m+3}\right]}{\expectationvalue*{\hat{\mathcal{O}}_{i}}\expectationvalue*{\hat{\mathcal{O}}_{j}}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\quantity(1-\delta_{ij}).$
(222)

Using Eq. (213), we can write $t$ as:

$\displaystyle
t=\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+2}\frac{1}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}\frac{1}{\expectationvalue*{\hat{\mathcal{O}}_{j}}}\frac{1}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\quantity(1-\delta_{ij})\Bigg{(}$
$\displaystyle\hat{\mathcal{O}}_{i}\hat{\mathcal{O}}_{j}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}+\hat{\mathcal{O}}_{m+3}\hat{\mathcal{O}}_{i}\expectationvalue*{\hat{\mathcal{O}}_{j}}+\hat{\mathcal{O}}_{j}\hat{\mathcal{O}}_{m+3}\expectationvalue*{\hat{\mathcal{O}}_{i}}$
$\displaystyle\quad-\hat{\mathcal{O}}_{i}\expectationvalue*{\hat{\mathcal{O}}_{j}}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}-\hat{\mathcal{O}}_{m+3}\expectationvalue*{\hat{\mathcal{O}}_{i}}\expectationvalue*{\hat{\mathcal{O}}_{j}}$
$\displaystyle\qquad\quad-\hat{\mathcal{O}}_{j}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}\expectationvalue*{\hat{\mathcal{O}}_{i}}+\expectationvalue*{\hat{\mathcal{O}}_{i}}\expectationvalue*{\hat{\mathcal{O}}_{j}}\expectationvalue*{\hat{\mathcal{O}}_{m+3}}\Bigg{)}.$
(223)

After some straightforward algebra, we can write $t$ as:

$\displaystyle
t=\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+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})$
$\displaystyle+\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\quantity(m+1)\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}$
$\displaystyle-\quantity(m+1)\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}$
$\displaystyle-\frac{1}{2}\quantity(m+2)\quantity(m+1)\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}+\frac{1}{2}\quantity(m+2)\quantity(m+1).$
(224)

Where we have used,
$\sum_{i,j=1}^{m+2}\quantity(1-\delta_{ij})=\quantity(m+2)\quantity(m+1)$ and
$\sum_{i,j=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}\quantity(1-\delta_{ij})=\quantity(m+1)\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}$.
Finally, the left-hand side of Eq. (219) can be written as:

$\displaystyle\mathrm{LHS}=$
$\displaystyle\bigg{(}\frac{1}{2}\sum_{i=1}^{m+2}\sum_{j=1}^{m+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})+\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}$
$\displaystyle\qquad-\quantity(m+1)\sum_{i=1}^{m+2}\frac{\hat{\mathcal{O}}_{i}}{\expectationvalue*{\hat{\mathcal{O}}_{i}}}-\quantity(m+1)\frac{\hat{\mathcal{O}}_{m+3}}{\expectationvalue*{\hat{\mathcal{O}}_{m+3}}}+\frac{\quantity(m+1)}{2}\quantity(m+2)\bigg{)}\prod_{k=1}^{m+3}\expectationvalue*{\hat{\mathcal{O}}_{k}}$
$\displaystyle=$ $\displaystyle\mathrm{RHS}.$ (225)

Hence, the left-hand side is equal to the right-hand side given in Eq. (220),
the formula also holds for $n=m+1$, completing the inductive step. Hence, by
mathematical induction, Eq. (218) is valid for $n\geq 1$.

##  Appendix C Defining the region of integration for the $B_{0}$ and
$\mathcal{A}$ functions after changing variables

In order to complete the change of variables performed in the two fermion, one
loop integral (in the $\absolutevalue{\bm{k}}>0$ case) defined in Eqs.
(61)$-$(65), one must find the limits of integration. The case of equal masses
was performed in Ref. [71]. By construction, the following equality holds,
$E[M_{j},\bm{p}-\bm{k}]^{2}-E\quantity[M_{i},\bm{p}]^{2}=\bm{k}^{2}-2\absolutevalue{\bm{p}}\absolutevalue{\bm{k}}\cos\quantity[\theta]-M_{i}^{2}+M_{j}^{2}=2E\varepsilon$,
with $E\quantity[M,\bm{p}]$ defined in Eq. (47) while, $E$ and $\varepsilon$,
are the new integration variables, defined in Eqs. (61) and (62),
respectively. Solving this equation with respect to $\cos\quantity[\theta]$
and squaring both sides, yields:

$\displaystyle\cos\quantity[\theta]^{2}=\frac{\quantity(\bm{k}^{2}-M_{i}^{2}+M_{j}^{2}-2E\varepsilon)^{2}}{4\bm{p}^{2}\bm{k}^{2}}.$
(226)

In order to have a right-hand side which depends only on the new integration
variables and the external 3-momentum, $\absolutevalue{\bm{k}}$, we can solve
$E\quantity[M_{i},\bm{p}]=\sqrt{\bm{p}^{2}+M_{i}^{2}}=\quantity(E-\frac{\varepsilon}{2})$
for $\bm{p}^{2}$ (see Eqs. (61) and (62)) and substitute in the above
equation. Since, in the original integration bounds, the angle
$\theta\in\quantity[0,\pi]$, it implies that $\cos^{2}[\theta]\leq 1$. This
allows us to define the following inequality:

$\displaystyle
P\quantity[\varepsilon,E,M_{i},M_{j},\absolutevalue{\bm{k}}]=\quantity(\bm{k}^{2}-M_{i}^{2}+M_{j}^{2}-2E\varepsilon)^{2}-4\bm{k}^{2}\quantity(\quantity(E-\frac{\varepsilon}{2})^{2}-M_{i}^{2})\leq
0.$ (227)

In the above, for simplicity, we have defined the function
$P\quantity[\varepsilon,E,M_{i},M_{j},\absolutevalue{\bm{k}}]$. It is easier
to obtain the integration region if we write the function
$P\quantity[\varepsilon,E,M_{i},M_{j},\absolutevalue{\bm{k}}]$ as a polynomial
in $\varepsilon$, with coefficients which depend on the remaining variables,
$E$, $M_{i}$, $M_{j}$ and $\bm{k}^{2}$. The resulting inequality is quadratic
in $\varepsilon$:

$\displaystyle
P\quantity[\varepsilon,E,M_{i},M_{j},\absolutevalue{\bm{k}}]=C_{0}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]+C_{1}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]\varepsilon+C_{2}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]\varepsilon^{2}\leq
0,$ (228)

where,

$\displaystyle C_{0}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\bm{k}^{4}+2\bm{k}^{2}\quantity(M_{i}^{2}+M_{j}^{2}-2E^{2})+\quantity(M_{i}^{2}-M_{j}^{2})^{2},$
(229) $\displaystyle C_{1}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=4E\quantity(M_{i}^{2}-M_{j}^{2}),$ (230) $\displaystyle
C_{2}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=4E^{2}-\bm{k}^{2}.$ (231)

The concavity of the inequality can be found by checking the sign of the
coefficient in front of the term proportional to $\varepsilon^{2}$,
$C_{2}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]$. One can show that this
coefficient is always bigger than zero due to the so-called triangle
inequality: in a normed vector space $V$, one can write the triangle
inequality:
$\absolutevalue{\bm{x}+\bm{y}}\leq\absolutevalue{\bm{x}}+\absolutevalue{\bm{y}}$,
for all $\bm{x}$, $\bm{y}$ that belongs to $V$. Trivially, we can write
$\absolutevalue{\bm{k}}=\absolutevalue{\quantity(\bm{k}-\bm{p})+\bm{p}}$ and,
making use of the triangle inequality, it follows that:

$\displaystyle\absolutevalue{\bm{k}}$
$\displaystyle\leq\absolutevalue{\bm{k}-\bm{p}}+\absolutevalue{\bm{p}}$
$\displaystyle\Rightarrow\bm{k}^{2}$
$\displaystyle\leq\quantity(\absolutevalue{\bm{k}-\bm{p}})^{2}+\quantity(\absolutevalue{\bm{p}})^{2}+2\absolutevalue{\bm{k}-\bm{p}}\absolutevalue{\bm{p}}.$
(232)

One can always write
$\absolutevalue{\bm{p}}\leq\sqrt{\bm{p}^{2}+M_{i}^{2}}=E\quantity[M_{i},\bm{p}]$
and
$\absolutevalue{\bm{k}-\bm{p}}\leq\sqrt{\quantity(\bm{k-p})^{2}+M_{j}^{2}}=E[M_{j},\bm{p}-\bm{k}]$,
since $M_{i}$ and $M_{j}$ are real and $M_{i},M_{j}>0$. The above inequality
then becomes:

$\displaystyle\bm{k}^{2}$
$\displaystyle<\quantity(E[M_{j},\bm{p}-\bm{k}]+E\quantity[M_{i},\bm{p}])^{2}=4E^{2}$
$\displaystyle\Rightarrow 0$ $\displaystyle<4E^{2}-\bm{k}^{2}.$ (233)

Thus, the coefficient of the quadratic term,
$C_{2}\quantity[E,M_{i},M_{j},\absolutevalue{\bm{k}}]=4E^{2}-\bm{k}^{2}$, is
always positive, meaning that it corresponds to a parabola with the concavity
turned upwards. Hence, for every fixed value of $E$, the region of integration
is an upwards parabola. Since we are interested in the region delimited by the
values of $\varepsilon$ and $E$ for which
$P\quantity[\varepsilon,E,M_{i},M_{j},\absolutevalue{\bm{k}}]\leq 0$, this
implies that the region of integration corresponds to the interior of the
parabola (the region between the two zeros).

To finish determining the integration region, we need to apply the restriction
coming from the 3-momentum cutoff in the 3-momentum sphere intersection
regularization (see Section II.3). These restrictions can readily be found by
demanding that, $\absolutevalue{\bm{p}}\leq\Lambda$ and
$\absolutevalue{\bm{p}-\bm{k}}\leq\Lambda$. For the first constraint, we get:

$\displaystyle\absolutevalue{\bm{p}}\leq\Lambda\Rightarrow$
$\displaystyle\quantity(E-\frac{\varepsilon}{2})^{2}\leq\Lambda^{2}+M_{i}^{2}\Rightarrow-\sqrt{\Lambda^{2}+M_{i}^{2}}\leq
E-\frac{\varepsilon}{2}\leq\sqrt{\Lambda^{2}+M_{i}^{2}}.$ (234)

For the second constraint, one gets:

$\displaystyle\absolutevalue{\bm{p}-\bm{k}}$
$\displaystyle\leq\Lambda\Rightarrow\quantity(E+\frac{\varepsilon}{2})^{2}\leq\Lambda^{2}+M_{j}^{2}\Rightarrow-\sqrt{\Lambda^{2}+M_{j}^{2}}\leq
E+\frac{\varepsilon}{2}\leq\sqrt{\Lambda^{2}+M_{j}^{2}}.$ (235)

At this point, one specific choice of integration order must be chosen in
order to evaluate the integral. As previously discussed it is quite useful to
utilize both integration orders since, in certain terms, one can analytically
perform the innermost integration when using the order
$\int\differential{\varepsilon}\int\differential{E}$ while, for other terms,
this analytical integration can only be performed for the other integration
order, $\int\differential{E}\int\differential{\varepsilon}$.

It will also be quite useful to evaluate the integral $R_{ij}$, defined by:

$\displaystyle
R_{ij}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]=\frac{\absolutevalue{\bm{k}}}{2\pi}\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential[3]{\bm{p}}}{E\quantity[M_{i},\bm{p}]E[M_{j},\bm{p}-\bm{k}]}.$
(236)

Where the integration is to be made considering the 3-momentum sphere
intersection regularization. Changing to the $E$ and $\varepsilon$ variables
(see Eqs. (61) and (62) and (63)), one obtains:

$\displaystyle
R_{ij}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{E}\differential{\varepsilon}.$
(237)

Hence, one immediately finds that $2\pi R_{ij}$ corresponds to the volume of
the integration region delimited by the 3-momentum sphere intersection
regularization in the new variables $E$ and $\varepsilon$. Evaluating this
quantity will be quite useful in order to: $(1)$ check that the change of
variables is valid; $(2)$ verify that both integration orders,
$\int\differential{\varepsilon}\int\differential{E}$ and
$\int\differential{E}\int\differential{\varepsilon}$, yield the same results.
In a later section we will also show how to evaluate this integral ($R_{ij}$),
without making the change of variables, allowing for a multiple numerical
cross check for the validity of the proposed integration procedure.

###  C.1  Integration order $\differential{\varepsilon}\differential{E}$

In this section we will consider the
$\differential{\varepsilon}\differential{E}$ integration order, hence we will
integrate first over the $E$ variable. In other words, we write the
integration of a function $f(E,\varepsilon)$, as:

$\displaystyle\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{\varepsilon}\differential{E}f(E,\varepsilon)=\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}f(E,\varepsilon)$
(238)

and define the different functions,
$E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$,
$\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$
and
$\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$,
which encompass the integration region and are functions of the masses
($M_{i}$ and $M_{j}$), the absolute value of the external 3-momentum
($\absolutevalue{\bm{k}}$), and the 3-momentum cutoff, $\Lambda$. In the
special case in which $f(E,\varepsilon)=1$, one obviously gets the integral
$R_{ij}$, see Eq. (237).

To achieve this, we start by solving Eq. (227) with respect to $E$. This
yields two solutions ($E_{\mathrm{min}}^{\pm}$):

$E_{\mathrm{min}}^{\pm}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]=\begin{cases}-\frac{(M_{i}^{2}-M_{j}^{2})}{4\varepsilon}-\frac{\varepsilon\quantity(M_{i}^{2}+M_{j}^{2})}{2(M_{i}^{2}-M_{j}^{2})},&\text{if}\;\varepsilon=\pm\absolutevalue{\bm{k}},\\\
\frac{\varepsilon(M_{i}^{2}-M_{j}^{2})\pm\sqrt{\bm{k}^{2}\quantity(\bm{k}^{2}+(M_{i}-M_{j})^{2}-\varepsilon^{2})\quantity(\bm{k}^{2}+(M_{i}+M_{j})^{2}-\varepsilon^{2})}}{2\quantity(\absolutevalue{\bm{k}}-\varepsilon)\quantity(\absolutevalue{\bm{k}}+\varepsilon)},&\text{otherwise}.\end{cases}$

In the special case in which $M_{i}=M_{j}=M$ one gets:

$\displaystyle
E_{\mathrm{min}}^{\pm}\quantity[M,M,\absolutevalue{\bm{k}},\varepsilon]=\pm\frac{\sqrt{\bm{k}^{2}\quantity(\bm{k}^{2}-\varepsilon^{2})\quantity(\bm{k}^{2}+4M^{2}-\varepsilon^{2})}}{2\quantity(\absolutevalue{\bm{k}}-\varepsilon)\quantity(\absolutevalue{\bm{k}}+\varepsilon)}.$
(240)

From the restrictions imposed by the cutoff in Eqs. (234) and (235), one can
further constrain the integration region. One can define the following
function from solving theses equation with respect to $E$:

$\displaystyle E_{\Lambda_{i}}^{\pm}[\Lambda,M_{i},\varepsilon]$
$\displaystyle=+\frac{\varepsilon}{2}\pm\sqrt{\Lambda^{2}+M_{i}^{2}},$ (241)
$\displaystyle E_{\Lambda_{j}}^{\pm}[\Lambda,M_{j},\varepsilon]$
$\displaystyle=-\frac{\varepsilon}{2}\pm\sqrt{\Lambda^{2}+M_{j}^{2}}.$ (242)

By definition, the variable $E$ is positive (see Eqs. (55), (56) and (61)).
Hence, we are interested only on the functions, $E_{\Lambda_{i}}^{+}$ and
$E_{\Lambda_{j}}^{+}$. Furthermore, these functions intersect at a certain
point, $\quantity(\varepsilon_{\Lambda_{s}},E_{\Lambda_{s}})$. Such point is
the solution to the equation, $E_{\Lambda_{i}}^{+}=E_{\Lambda_{j}}^{+}$ and it
is given by:

$\displaystyle\varepsilon_{\Lambda_{s}}\quantity[\Lambda,M_{i},M_{j}]$
$\displaystyle=\sqrt{\Lambda^{2}+M_{j}^{2}}-\sqrt{\Lambda^{2}+M_{i}^{2}},$
(243) $\displaystyle E_{\Lambda_{s}}\quantity[\Lambda,M_{i},M_{j}]$
$\displaystyle=\frac{1}{2}\quantity(\sqrt{\Lambda^{2}+M_{i}^{2}}+\sqrt{\Lambda^{2}+M_{j}^{2}}).$
(244)

Also, at this point, the restriction coming from the cutoff demands one to
switch from one branch, to another. This can be more easily understood after a
graphical analysis. Thus, we can unify the constraint coming from the
3-momentum cutoff in the following equation:

$E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon]=\begin{cases}E_{\Lambda_{j}}^{+}[\Lambda,M_{j},\varepsilon],&\text{if}\;\varepsilon>\varepsilon_{\Lambda_{s}}\quantity[\Lambda,M_{i},M_{j}],\\\
E_{\Lambda_{i}}^{+}[\Lambda,M_{i},\varepsilon],&\text{otherwise}.\end{cases}$

At this point we can define the lower bound of the $\differential{E}$
integration,
$E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$. It is
simply given by:

$\displaystyle
E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]=E_{\mathrm{min}}^{+}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon].$
(246)

The upper bound of the $\differential{E}$ integration,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$, however, is
not so simple. When dealing with equal fermion masses ($M_{i}=M_{j}=M$), the
constraint coming from the cutoff defined previously, completely fixes the
upper bound. However, when dealing with different masses, this is not true and
a switch from $E_{\Lambda}$ to the function $E_{\mathrm{min}}^{-}$, can occur.
This is also better understood after a visual analysis of the interplay
between the functions. In order to define
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$ we will need
to find other quantities of interest where some of the previously mentioned
quantities intersect. Firstly, consider a possible intersection between the
function $E_{\mathrm{min}}^{-}$ and $E_{\Lambda_{i}}^{+}$, i.e., solutions of
the equation,
$E_{\mathrm{min}}^{-}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon_{\gamma_{1}}]=E_{\Lambda_{i}}^{+}[\Lambda,M_{i},\varepsilon_{\gamma_{1}}]$.
This results in four different intersection points, with only one being
meaningful to the integration region (the other ones are outside the region).
Its $\quantity(\varepsilon_{\gamma_{1}},E_{\gamma_{1}})$ values are given by:

$\displaystyle\varepsilon_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\sqrt{\quantity(\absolutevalue{\bm{k}}-\Lambda)^{2}+M_{j}^{2}}-\sqrt{\Lambda^{2}+M_{i}^{2}},$
(247) $\displaystyle
E_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2}\quantity(\sqrt{\quantity(\absolutevalue{\bm{k}}-\Lambda)^{2}+M_{j}^{2}}+\sqrt{\Lambda^{2}+M_{i}^{2}}).$
(248)

Likewise, the solution of interest to the equation
$E_{\mathrm{min}}^{-}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon_{\gamma_{2}}]=E_{\Lambda_{j}}^{+}[\Lambda,M_{j},\varepsilon_{\gamma_{2}}]$
is:

$\displaystyle\varepsilon_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\sqrt{\Lambda^{2}+M_{j}^{2}}-\sqrt{\quantity(\absolutevalue{\bm{k}}-\Lambda)^{2}+M_{i}^{2}},$
(249) $\displaystyle
E_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2}\quantity(\sqrt{\Lambda^{2}+M_{j}^{2}}+\sqrt{\quantity(\absolutevalue{\bm{k}}-\Lambda)^{2}+M_{i}^{2}}).$
(250)

Secondly, the equations $E_{\mathrm{min}}^{+}=E_{\mathrm{min}}^{-}$ are equal
in four different points, with two of them being useful for our purposes,
$\varepsilon_{\alpha_{1}}$ and $\varepsilon_{\alpha_{2}}$. They can be
calculated to yield:

$\displaystyle\varepsilon_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=-\sqrt{\bm{k}^{2}+(M_{i}-M_{j})^{2}}=-\varepsilon_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}].$
(251)

To each of these $\varepsilon$ points there are corresponding $E$ values,
given by:

$\displaystyle E_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2}\frac{\quantity(M_{i}+M_{j})}{\quantity(M_{i}-M_{j})}\sqrt{\bm{k}^{2}+(M_{i}-M_{j})^{2}}=-E_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}].$
(252)

From the above, one can observe that, $E_{\alpha_{1}},E_{\alpha_{2}}\to\infty$
for the degenerate mass case, meaning that the two equations,
$E_{\mathrm{min}}^{+}$ and $E_{\mathrm{min}}^{-}$, are never equal. However,
for different masses, there will exist two switching points, symmetric with
respect to $\varepsilon$ (see Eq. (251)). Thus, as previously aluded to,
increasing the difference between the masses, causes $E_{\alpha_{1}}$ and
$E_{\alpha_{2}}$ to decrease and the upper bound,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$, ceases to
be defined by $E_{\Lambda}$ switching to $E_{\mathrm{min}}^{-}$. In the case
where $M_{i}>M_{j}$, one can write:

$E_{\mathrm{max}}^{(M_{i}>M_{j})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]=\begin{cases}E_{\mathrm{min}}^{-}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon],&\begin{array}[]{l}\text{if}\;\varepsilon<\varepsilon_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]\;\text{and}\\\
\phantom{\text{if}\;}E_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]<E_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon],&\;\text{otherwise}.\end{cases}$

In the opposite mass case, $M_{j}>M_{i}$, the upper bound is:

$E_{\mathrm{max}}^{(M_{j}>M_{i})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]=\begin{cases}E_{\mathrm{min}}^{-}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon],&\begin{array}[]{l}\text{if}\;\varepsilon>\varepsilon_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]\;\text{and}\\\
\phantom{\text{if}\;}E_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]<E_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon],&\;\text{otherwise}.\end{cases}$

As explained before, for the degenerate mass case, the upper bound is
trivially given by the cutoff restriction only, i.e., $E_{\Lambda}$. Gathering
all this information, one can write the upper bound of the $\differential{E}$
integration,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}$, as:

$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]=\begin{cases}E_{\mathrm{max}}^{(M_{i}>M_{j})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon],&\text{if}\;M_{i}>M_{j},\\\
E_{\mathrm{max}}^{(M_{j}>M_{i})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon],&\text{if}\;M_{i}<M_{j},\\\
E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon],&\text{otherwise}.\end{cases}$

Finally, taking into account all the intersections of interest, one can write
the lower bound of the $\differential{\varepsilon}$ integration,
$\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}$,
as:

$\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\begin{cases}\varepsilon_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],&\begin{array}[]{l}\text{if}\;M_{i}>M_{j}\;\text{and}\\\
\phantom{\text{if}\;}E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]]\geq
E_{\alpha_{1}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
\varepsilon_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\;\text{otherwise}.\end{cases}$

The upper bound, on the other hand, is given by:

$\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\begin{cases}\varepsilon_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],&\begin{array}[]{l}\text{if}\;M_{i}<M_{j}\;\text{and}\\\
\phantom{\text{if}\;}E_{\Lambda}[\Lambda,M_{i},M_{j},\varepsilon_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]]\geq
E_{\alpha_{2}}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
\varepsilon_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\;\text{otherwise}.\end{cases}$

One can check that the size of the integration region gets smaller with
increasing $\absolutevalue{\bm{k}}$. This is a consequence of the
regularization used in this work. As explained earlier, increasing the value
of $\absolutevalue{\bm{k}}$ is equivalent to increasing the distance between
the center of two 3-momentum spheres of radii $\Lambda$. For
$\absolutevalue{\bm{k}}=2\Lambda$ the spheres no longer intersect and the
integration region should be zero. Hence, to improve the numerical
implementation of the calculation, one can additionally include a condition to
only calculate the integral if $0<\absolutevalue{\bm{k}}<2\Lambda$.

In Fig. 15 we show the region of integration
$\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}$,
for fixed values of $\Lambda$, $\absolutevalue{\bm{k}}$, $M_{i}$ and $M_{j}$.
We consider the three possible scenarios of $M_{i}=M_{j}$, $M_{i}>M_{j}$ and
$M_{i}<M_{j}$, in the left, middle and right panels, respectively. For
simplicity, we fix the cutoff to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$
and consider different values of $\absolutevalue{\bm{k}}$ and masses. The goal
is to provide different scenarios for easier reproduction of these results. We
also show the values for the different areas resulting from integrating the
region (see Eq. (237)). As discussed before, calculating these areas
corresponds to evaluating the $R_{ij}$ integral first defined in Eq (236).

(a) $M_{i}=M_{j}$

(b) $M_{i}>M_{j}$

(c) $M_{i}<M_{j}$

Figure 15:  Region of the integration
$\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}$,
for fixed values of $\Lambda$, $\absolutevalue{\bm{k}}$ and different values
of $M_{i}$, $M_{j}$, with $M_{i}=M_{j}$ in the left panel, $M_{i}>M_{j}$ in
the middle and $M_{j}>M_{i}$ in the right panel. The full black line is the
upper limit of the integration over the variable $E$,
$E_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]$
(Eq. (C.1)) while, the dashed black line, is the inferior bound,
$E_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},\varepsilon]$
(Eq. (246)). The red dots represent the minimum and maximum values of the
integration over $\varepsilon$,
$\varepsilon_{\mathrm{min}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
and
$\varepsilon_{\mathrm{max}}^{(\differential{\varepsilon}\differential{E})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$,
respectively (Eqs. (C.1) and (C.1)). In each plot we also give the values for
the area of each integration region, which corresponds to evaluating the
$R_{ij}$ integral, defined in Eq. (236).

###  C.2  Integration order $\differential{E}\differential{\varepsilon}$

In this section we will consider the
$\differential{E}\differential{\varepsilon}$ integration order, meaning that
one first integrates over the $\varepsilon$ variable. We write:

$\displaystyle\int_{\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential{\phi}}{2\pi}\differential{E}\differential{\varepsilon}f(E,\varepsilon)=\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}f(E,\varepsilon).$
(258)

Again, if $f(E,\varepsilon)=1$, one recovers the integral $R_{ij}$, see Eq.
(237).

For this case, the starting point is finding the solutions to Eq. (227) with
respect to $\varepsilon$. This yields two solutions
($\varepsilon_{\mathrm{min}}^{\pm}$):

$\varepsilon_{\mathrm{min}}^{\pm}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},E]=\begin{cases}-\frac{(M_{i}^{2}-M_{j}^{2})}{4E}-\frac{2E\quantity(M_{i}^{2}+M_{j}^{2})}{(M_{i}^{2}-M_{j}^{2})},&\text{if}\;E=\pm\nicefrac{{\absolutevalue{\bm{k}}}}{{2}},\\\
\frac{-2E\quantity(M_{i}^{2}-M_{j}^{2})\pm\sqrt{\bm{k}^{2}\quantity(4E^{2}-\bm{k}^{2}-(M_{i}-M_{j})^{2})\quantity(4E^{2}-\bm{k}^{2}-(M_{i}+M_{j})^{2})}}{4E^{2}-\bm{k}^{2}},&\text{otherwise}.\end{cases}$

For equal masses, $M_{i}=M_{j}=M$, one gets:

$\displaystyle\varepsilon_{\mathrm{min}}^{\pm}\quantity[M,M,\absolutevalue{\bm{k}},E]=\pm\frac{\sqrt{\bm{k}^{2}\quantity(4E^{2}-\bm{k}^{2})\quantity(4E^{2}-\bm{k}^{2}-4M^{2})}}{4E^{2}-\bm{k}^{2}}.$
(260)

As before, the cutoff restrictions to the integration region can be obtained
by solving Eqs. (234) and (235) with respect to $\varepsilon$. It yields:

$\displaystyle\varepsilon_{\Lambda_{i}}^{\pm}[\Lambda,M_{i},E]$
$\displaystyle=2E\pm 2\sqrt{\Lambda^{2}+M_{i}^{2}},$ (261)
$\displaystyle\varepsilon_{\Lambda_{j}}^{\pm}[\Lambda,M_{j},E]$
$\displaystyle=-2E\pm 2\sqrt{\Lambda^{2}+M_{j}^{2}}.$ (262)

The upper and lower bounds of the $\differential{\varepsilon}$ integration,
$\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}$
and
$\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}$
can be defined from these equations. Starting with the upper bound, it is
defined by the functions $\varepsilon_{\mathrm{min}}^{+}$ and
$\varepsilon_{\Lambda_{j}}^{+}$ which intersect at the point
$\quantity(E_{\gamma_{2}},\varepsilon_{\gamma_{2}})$, previously defined in
Eqs. (249) and (250). Hence, the upper bound of the
$\differential{\varepsilon}$ integration,
$\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}$,
can be defined as:

$\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},E]=\begin{cases}\varepsilon_{\Lambda_{j}}^{+}[\Lambda,M_{j},E],&\text{if}\;E\geq
E_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],\\\
\varepsilon_{\mathrm{min}}^{+}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},E],&\text{otherwise}.\end{cases}$

Similarly, the lower bound is defined by the functions
$\varepsilon_{\mathrm{min}}^{-}$ and $\varepsilon_{\Lambda_{i}}^{-}$. These
functions intersect at the point
$\quantity(\varepsilon_{\gamma_{1}},E_{\gamma_{1}})$, defined in Eqs. (247)
and (248). Thus, we can write:

$\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},E]=\begin{cases}\varepsilon_{\Lambda_{i}}^{-}[\Lambda,M_{i},E],&\text{if}\;E\geq
E_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],\\\
\varepsilon_{\mathrm{min}}^{-}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}},E],&\text{otherwise}.\end{cases}$

As with the integration order considered in the previous section,
$\differential{\varepsilon}\differential{E}$, depending if one is considering
equal or different fermion masses, the lower bound of the $\differential{E}$
integration,
$\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}$
may be defined by the intersection of the functions
$\varepsilon_{\mathrm{min}}^{-}$ with $\varepsilon_{\mathrm{min}}^{+}$, by the
intersection of $\varepsilon_{\mathrm{min}}^{+}$ with
$\varepsilon_{\Lambda_{i}}^{-}$ or by the intersection of
$\varepsilon_{\mathrm{min}}^{-}$ with $\varepsilon_{\Lambda_{j}}^{+}$. When
the fermion masses are equal, the lower bound is given by the first
intersection mentioned, i.e., when
$\varepsilon_{\mathrm{min}}^{-}=\varepsilon_{\mathrm{min}}^{+}$. This occurs
at the point $\quantity(E_{\alpha},\varepsilon_{\alpha})$ defined by:

$\displaystyle E_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=\frac{1}{2}\sqrt{\bm{k}^{2}+\quantity(M_{i}+M_{j})^{2}},$ (265)
$\displaystyle\varepsilon_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=-\frac{\quantity(M_{i}^{2}-M_{j}^{2})\sqrt{\bm{k}^{2}+(M_{i}+M_{j})^{2}}}{\quantity(M_{i}+M_{j})^{2}}.$
(266)

If $M_{i}>M_{j}$, the lower bound of the $\differential{E}$ integration can
also be located at the above intersection point, $E_{\alpha}$. However, for
certain values of the parameters $\Lambda$ $M_{i}$, $M_{j}$ and
$\absolutevalue{\bm{k}}$, the lower bound of the $\differential{E}$
integration can occur at the intersection between
$\varepsilon_{\mathrm{min}}^{+}$ and $\varepsilon_{\Lambda_{i}}^{-}$. This
intersection also occurs at the point
$\quantity(\varepsilon_{\gamma_{1}},E_{\gamma_{1}})$, defined in Eqs. (247)
and (248). Hence, we can write:

$E_{\mathrm{min}}^{(M_{i}>M_{j})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\begin{cases}E_{\gamma_{1}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\begin{array}[]{l}\text{if}\;\varepsilon_{\Lambda_{i}}^{-}[\Lambda,M_{i},\varepsilon_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]]>\varepsilon_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
E_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],&\;\text{otherwise}.\end{cases}$

Likewise, for the opposite mass case, in which $M_{j}>M_{i}$, one can realize
that the lower bound of the integration can occur at the point $E_{\alpha}$
or, it may occur at the intersection between $\varepsilon_{\mathrm{min}}^{-}$
and $\varepsilon_{\Lambda_{j}}^{+}$, defined by the point
$\quantity(E_{\gamma_{2}},\varepsilon_{\gamma_{2}})$ given in Eqs. (249) and
(250). Thus, we write:

$E_{\mathrm{min}}^{(M_{j}>M_{i})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\begin{cases}E_{\gamma_{2}}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\begin{array}[]{l}\text{if}\;\varepsilon_{\Lambda_{j}}^{+}[\Lambda,M_{i},\varepsilon_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}]]<\varepsilon_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],\end{array}\\\
E_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],&\;\text{otherwise}.\end{cases}$

Collecting these conditions, we can finally define:

$E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\begin{cases}E_{\mathrm{min}}^{(M_{i}>M_{j})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\text{if}\;M_{i}>M_{j},\\\
E_{\mathrm{min}}^{(M_{j}>M_{i})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}],&\text{if}\;M_{i}<M_{j},\\\
E_{\alpha}\quantity[M_{i},M_{j},\absolutevalue{\bm{k}}],&\text{otherwise}.\end{cases}$

The upper bound of the $\differential{E}$ integration is much simpler to
define. It is defined by the intersection between the cutoff restrictions,
which takes place when
$\varepsilon_{\Lambda_{i}}^{-}=\varepsilon_{\Lambda_{j}}^{+}$, in the formerly
defined point $\quantity(E_{\Lambda_{s}},\varepsilon_{\Lambda_{s}})$ given in
Eqs. (243) and (244). In such a way, the upper limit is then simply given by:

$\displaystyle
E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=E_{\Lambda_{s}}\quantity[\Lambda,M_{i},M_{j}].$
(270)

In Fig. 16 we show the region of integration
$\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{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}$,
for fixed values of $\Lambda$, $\absolutevalue{\bm{k}}$, $M_{i}$ and $M_{j}$.
Again, we consider the three possible scenarios of $M_{i}=M_{j}$,
$M_{i}>M_{j}$ and $M_{i}<M_{j}$, in the left, middle and right panels,
respectively. It is now easy to check that one gets exactly the same areas
(or, identically, values for the $R_{ij}$ integral) using this integration
order, for the same set of parameters, when comparing with the other possible
integration order whose areas are shown in Fig. 15.

(a) $M_{i}=M_{j}$

(b) $M_{i}>M_{j}$

(c) $M_{i}<M_{j}$

Figure 16:  Region of the integration
$\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{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}$,
for fixed values of $\Lambda$, $\absolutevalue{\bm{k}}$ and different values
of $M_{i}$, $M_{j}$, with $M_{i}=M_{j}$ in the left panel, $M_{i}>M_{j}$ in
the middle and $M_{j}>M_{i}$ in the right panel. The full black line is the
upper limit of the integration over the variable $\varepsilon$,
$\varepsilon_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},E]$
(Eq. (C.2)) while, the dashed black line, is the inferior bound,
$\varepsilon_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}},E]$
(Eq. (C.2)). The red dots represent the minimum and maximum values of the
integration over $E$,
$E_{\mathrm{min}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
and
$E_{\mathrm{max}}^{(\differential{E}\differential{\varepsilon})}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$,
respectively ((Eq. (C.2)) and (Eq. (270))). In each plot we also give the
values for the area of each integration region, which corresponds to
evaluating the $R_{ij}$ integral, defined in Eq. (236).

###  C.3 Calculating the integral $R_{ij}$ without changing variables

In this section we show how to evaluate the integral $R_{ij}$, defined in Eq.
(236), without changing to the $E$ and $\varepsilon$ variables. This is useful
for numerical cross checks, as well as evaluating the vacuum limit of the one
fermion line integral in the 3-momentum sphere intersection regularization,
$\mathcal{A}$.

As discussed in Section II.3, using the 3-momentum sphere intersection
regularization amounts to performing an integral over the region defined by
the intersection between two spheres,
$\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}$. One of the spheres is
centered at the origin while the other is distanced $\absolutevalue{\bm{k}}$
from the first, in the $p_{z}$ axis,
$\mathbb{I}_{\absolutevalue{\bm{k}}}^{\Lambda}=\mathbb{S}_{0}^{\Lambda}\cap\mathbb{S}_{\absolutevalue{\bm{k}}}^{\Lambda}$.
Such region can be observed in Fig. 3. The major difficulty of evaluating the
integral $R_{ij}$ over this 3-dimensional region lays in the fact that it is
not mirror-symmetric in the $p_{z}$ axis. In order to make the region
symmetrical we are allowed to consider, in Eq. (236), the following change of
variables, $\bm{p}\to\bm{p}+\nicefrac{{\bm{k}}}{{2}}$. With this
transformation, we are allowed to write the $R_{ij}$ integral as:

$\displaystyle
R_{ij}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=\frac{\absolutevalue{\bm{k}}}{2\pi}\int_{{\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}}\frac{\differential[3]{\bm{p}}}{E\quantity[M_{i},\bm{p}+\frac{\bm{k}}{2}]E[M_{j},\bm{p}-\frac{\bm{k}}{2}]}.$
(271)

Where, the new region of integration, is defined by the intersection of the
two spheres distanced $\absolutevalue{\bm{k}}$ from each other, with one
centered at
$\quantity(p_{x},p_{y},p_{z})=\quantity(0,0,\nicefrac{{\absolutevalue{\bm{k}}}}{{2}})$
and the other centered at
$\quantity(p_{x},p_{y},p_{z})=\quantity(0,0,-\nicefrac{{\absolutevalue{\bm{k}}}}{{2}})$
i.e.,
${\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}=\mathbb{S}_{\nicefrac{{\absolutevalue{\bm{k}}}}{{2}}}^{\Lambda}\cap\mathbb{S}_{-\nicefrac{{\absolutevalue{\bm{k}}}}{{2}}}^{\Lambda}$.
This new integration region is mirror-symmetric with respect to all axis,
allowing for a simpler parametrization of the integration bounds.

Figure 17:  Region of integration in the 3-momentum sphere intersection
regularization scheme used to regularize the $\mathcal{A}$ integral.

In Fig. 17 we show the 2-dimensional cross section of the intersection
${\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}$, in the
$p_{y}-p_{z}$ plane. Consider the more general case of an integral, $I$, with
an integrand, $F$, which is independent of the azimuthal angle, $\phi$:

$\displaystyle
I=\int_{{\mathbb{I}^{\prime}}_{\absolutevalue{\bm{k}}}^{\Lambda}}\differential[3]{\bm{p}}F[\absolutevalue{\bm{p}},\theta].$
(272)

This integral can be separated into two parts: one corresponding to the upper
hemisphere and the other to the lower hemisphere of the intersection. Consider
also the functions $B_{U}$ and $B_{L}$ which parametrize the boundary of the
upper and lower hemispheres, respectively (see Fig. 17). Hence, we can write,
$I=I_{U}+I_{L}$, where:

$\displaystyle\frac{I_{U}}{2\pi}$
$\displaystyle=\int_{0}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{0}^{\nicefrac{{\pi}}{{2}}}\differential{\theta}\sin[\theta]F[\absolutevalue{\bm{p}},\theta]-\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{0}^{B_{U}}\differential{\theta}\sin[\theta]F[\absolutevalue{\bm{p}},\theta],$
(273) $\displaystyle\frac{I_{L}}{2\pi}$
$\displaystyle=\int_{0}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{0}^{\nicefrac{{\pi}}{{2}}}\differential{\theta}\sin[\theta]F[\absolutevalue{\bm{p}},\theta]-\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{B_{L}}^{\pi}\differential{\theta}\sin[\theta]F[\absolutevalue{\bm{p}},\theta].$
(274)

In the above, we first consider the integral of a sphere of radius
$p_{\mathrm{Max}}$ (the distance from the origin to the point where both
boundaries intersect, $B_{U}=B_{L}$, which occurs when
$\theta=\nicefrac{{\pi}}{{2}}$) and than subtract the extra volume coming from
the boundaries, $B_{U}$ and $B_{L}$ , up to this maximum sphere. The quantity
$p_{\mathrm{Min}}$ is the distance from the origin to the boundaries when
$\theta=0$. The quantities $B_{U}$, $B_{L}$, $p_{\mathrm{Min}}$ and
$p_{\mathrm{Max}}$, will be defined later. For simplicity, we apply the change
of variables, $\alpha=\cos[\theta]$ and change the integration bounds
accordingly to get:

$\displaystyle-\frac{I}{2\pi}$
$\displaystyle=\int_{0}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{1}^{-1}\differential{\alpha}F[\absolutevalue{\bm{p}},\alpha]$
$\displaystyle-\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{1}^{B_{U}}\differential{\alpha}F[\absolutevalue{\bm{p}},\alpha]$
$\displaystyle-\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{B_{L}}^{-1}\differential{\alpha}F[\absolutevalue{\bm{p}},\alpha].$
(275)

This expression can be further simplified by breaking the
$\differential{\absolutevalue{\bm{p}}}$ integral in the first term at the
point $p_{\mathrm{Min}}$ and by rearranging the integral. Finally, we can
write the integral as:

$\displaystyle I$
$\displaystyle=-2\pi\quantity(\int_{0}^{p_{\mathrm{Min}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{1}^{-1}\differential{\alpha}F[\absolutevalue{\bm{p}},\alpha]+\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{B_{U}}^{B_{L}}\differential{\alpha}F[\absolutevalue{\bm{p}},\alpha]).$
(276)

At this point we define the functions that define the integration bounds. For
the different hemispheres, the equations which define the boundaries are given
by:

$\displaystyle\bm{p}^{2}\sin[\theta]^{2}\cos[\phi]^{2}+\bm{p}^{2}\sin[\theta]^{2}\sin[\phi]^{2}+\quantity(\absolutevalue{\bm{p}}\cos[\theta]\pm\frac{\absolutevalue{\bm{k}}}{2})^{2}=\Lambda^{2}.$
(277)

Where the plus sign is for the upper hemisphere and the minus sign for the
lower hemisphere. Solving these equations with respect to $\cos[\theta]$,
yields:

$\displaystyle\cos[\theta]=\pm\frac{\Lambda^{2}-\nicefrac{{\bm{k}^{2}}}{{4}}-\bm{p}^{2}}{\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}.$
(278)

Again, the plus sign is for the upper region while the minus sign is for the
lower hemisphere. Thus, in the $\alpha=\cos[\theta]$ variable we can find
write $B_{U}$ and $B_{L}$ as:

$\displaystyle
B_{U}\quantity[\Lambda,\absolutevalue{\bm{k}},\absolutevalue{\bm{p}}]=\frac{\Lambda^{2}-\nicefrac{{\bm{k}^{2}}}{{4}}-\bm{p}^{2}}{\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}=-B_{L}\quantity[\Lambda,\absolutevalue{\bm{k}},\absolutevalue{\bm{p}}].$
(279)

The quantities $p_{\mathrm{Min}}$ and $p_{\mathrm{Max}}$ can be derived by
plugging $\theta=0$ and $\theta=\nicefrac{{\pi}}{{2}}$ in Eq. (278), to yield:

$\displaystyle p_{\mathrm{Min}}\quantity[\Lambda,\absolutevalue{\bm{k}}]$
$\displaystyle=\Lambda-\frac{\absolutevalue{\bm{k}}}{2},$ (280) $\displaystyle
p_{\mathrm{Max}}\quantity[\Lambda,\absolutevalue{\bm{k}}]$
$\displaystyle=\sqrt{\Lambda^{2}-\frac{\bm{k}^{2}}{4}}.$ (281)

We are finally ready to evaluate the integral $R_{ij}$. Applying the tools
developed above, we can write
$F[\absolutevalue{\bm{p}},\theta]=1/\quantity(E\quantity[M_{i},\bm{p}+\frac{\bm{k}}{2}]E[M_{j},\bm{p}-\frac{\bm{k}}{2}])$.
Hence, we can explicitly write:

$\displaystyle R_{ij}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]=$
$\displaystyle-\absolutevalue{\bm{k}}\int_{0}^{p_{\mathrm{Min}}}\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}}\cos[\theta]+M_{i}^{2}}\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}-\bm{p}\absolutevalue{\bm{k}}\cos[\theta]+M_{j}^{2}}}$
$\displaystyle-\absolutevalue{\bm{k}}\int_{p_{\mathrm{Min}}}^{p_{\mathrm{Max}}}\differential{\absolutevalue{\bm{p}}}\bm{p}^{2}\int_{B_{U}}^{B_{L}}\frac{\differential{\alpha}}{\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}+\bm{p}\absolutevalue{\bm{k}}\cos[\theta]+M_{i}^{2}}\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}-\bm{p}\absolutevalue{\bm{k}}\cos[\theta]+M_{j}^{2}}}.$
(282)

We can calculate the innermost integrals analytically:

$\displaystyle\int\differential{\alpha}$
$\displaystyle\frac{\absolutevalue{\bm{k}}\bm{p}^{2}}{\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}+\absolutevalue{\bm{p}}\absolutevalue{\bm{k}}\cos[\theta]+M_{i}^{2}}\sqrt{\bm{p}^{2}+\frac{\bm{k}^{2}}{4}-\absolutevalue{\bm{p}}\absolutevalue{\bm{k}}\cos[\theta]+M_{j}^{2}}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad=-2\absolutevalue{\bm{p}}\arcsin\quantity[\frac{\sqrt{\frac{\bm{k}^{2}}{4}+M_{j}^{2}+\bm{p}^{2}-\alpha\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}}{\sqrt{2}\sqrt{\frac{\bm{k}^{2}}{4}+\frac{M_{i}^{2}+M_{j}^{2}}{2}+\bm{p}^{2}}}]+\mathrm{cte}.$
(283)

Using this primitive, we are finally able to write $R_{ij}$ as:

$\displaystyle R_{ij}\quantity[\Lambda,M_{i},M_{j},\absolutevalue{\bm{k}}]$
$\displaystyle=2\sum_{\eta=\pm
1}\eta\int_{0}^{\Lambda-\frac{\absolutevalue{\bm{k}}}{2}}\differential{\absolutevalue{\bm{p}}}\absolutevalue{\bm{p}}\arcsin\quantity[\frac{\sqrt{\frac{\bm{k}^{2}}{4}+M_{j}^{2}+\bm{p}^{2}+\eta\absolutevalue{\bm{k}}\absolutevalue{\bm{p}}}}{\sqrt{2}\sqrt{\frac{\bm{k}^{2}}{4}+\frac{M_{i}^{2}+M_{j}^{2}}{2}+\bm{p}^{2}}}]$
$\displaystyle+2\sum_{\eta=\pm
1}\eta\int_{\Lambda-\frac{\absolutevalue{\bm{k}}}{2}}^{\sqrt{\Lambda^{2}-\frac{\bm{k}^{2}}{4}}}\differential{\absolutevalue{\bm{p}}}\absolutevalue{\bm{p}}\arcsin\quantity[\frac{\sqrt{M_{j}^{2}+\eta\Lambda^{2}+(1-\eta)(\bm{p}^{2}+\frac{\bm{k}^{2}}{4})}}{\sqrt{2}\sqrt{\frac{\bm{k}^{2}}{4}+\frac{M_{i}^{2}+M_{j}^{2}}{2}+\bm{p}^{2}}}].$
(284)

In Table 9 we give some results from solving the integral which can then be
directly compared to the results obtained when using the variables $E$ and
$\varepsilon$, given in Figs. 15 and 16.

$\Lambda\leavevmode\nobreak\ \quantity[\text{GeV}]$ | $M_{i}\leavevmode\nobreak\ \quantity[\text{GeV}]$ | $M_{j}\leavevmode\nobreak\ \quantity[\text{GeV}]$ | $\absolutevalue{\bm{k}}\leavevmode\nobreak\ \quantity[\text{GeV}]$ | $R_{ij}\leavevmode\nobreak\ \quantity[\text{GeV}^{2}]$
---|---|---|---|---
$1.0$ | $0.4$ | $0.4$ | $1.2$ | $0.2348$
$1.0$ | $0.3$ | $0.6$ | $0.8$ | $0.3583$
$1.0$ | $0.6$ | $0.3$ | $0.8$ | $0.3583$
$1.0$ | $0.4$ | $0.7$ | $0.2$ | $0.1476$
$1.0$ | $0.7$ | $0.4$ | $0.2$ | $0.1476$
$1.0$ | $1.2$ | $0.2$ | $0.5$ | $0.2334$
$1.0$ | $0.2$ | $1.2$ | $0.5$ | $0.2334$

Table 9:  Values for the $R_{ij}$ integral, defined in Eq. (236), for
different values of the cutoff, $\Lambda$, effective masses, $M_{i}$ and
$M_{j}$ and absolute value of the 3-momentum, $\absolutevalue{\bm{k}}$. The
$R_{ij}$ integral is identical to the area defined by the integral
$\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}$
and
$\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{max}}^{(\differential{E}\differential{\varepsilon})}}\differential{\varepsilon}$

(see Eq. (237) and Figs. 15 and 16 for the evaluation of the $R_{ij}$ using
the $E$ and $\varepsilon$ variables).

##  Appendix D Numerical results for the two fermion line integral, $B_{0}$

In this section we show numerical results for the two fermion line integral,
$B_{0}$, by choosing a particular value for the cutoff, $\Lambda$, and
constant values for the fermion masses, $M_{i}$ and $M_{j}$. This means we are
not supposing any model for the underlying fermions dynamics, meaning that the
fermion masses are not temperature and/or chemical potential dependent. We
chose to perform such study in order to isolate the temperature and chemical
potential dependence of the $B_{0}$ function. We show several results for this
function, in a diverse set of scenarios, in order to simplify the task of
reproducibility. As already stated, the calculation of this integral is
extremely important to studies involving mesons masses, the phase diagram of
strongly interacting matter and the transport properties of the quark-gluon
plasma in the relaxation time approximation, in the context of the NJL model
and similar models.

Figure 18: $B_{0}$ integral, in the vacuum, as a function of the zero
component of the external momentum, normalized by the cutoff, $k_{0}/\Lambda$.
The cutoff is fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the
fermion masses are equal and fixed to $M_{i}=M_{j}=0.4\Lambda$ and the
magnitude of the external 3-momentum, $\absolutevalue{\bm{k}}=0.0$. The full
and dashed lines correspond to the real and imaginary parts part of the
$B_{0}$ integral, respectively. In the left panel it was considered the mass
shift, $M^{2}\to M^{2}-i\epsilon$ while, in the right panel, the external
momentum shift was considered, $k_{0}\to k_{0}+i\epsilon$.

In Fig. 18, we show the result for the real and imaginary parts of the $B_{0}$
function, in the vacuum, as a function of the zero component of the external
momentum, normalized by the cutoff, $k_{0}/\Lambda$. We fixed the 3-momentum
cutoff, $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion masses are
equal and fixed to $M_{i}=M_{j}=0.4\Lambda$ and the magnitude of the external
3-momentum is zero, $\absolutevalue{\bm{k}}=0.0$.

In Ref. [11], other values for the above parameters are chosen, yielding a
slightly different result. However, when using the same parameter set, we were
able to obtain exactly the same results presented in Ref. [11]. As already
discussed, our approach only differs from the one proposed in that work if one
is considering the case with $\absolutevalue{\bm{k}}>0$, being exactly the
same if $\absolutevalue{\bm{k}}=0$.

As discussed previously, one must consider a pole shift in order to apply the
Sokhotski–Plemelj theorem. The way in which this shift is done, only changes
the behavior of the imaginary parts of the $B_{0}$ integral, as shown in Eqs.
(50) and (52). In each panel of this figure, we show the result for a specific
shift: on the left the mass shift is used, $M^{2}\to M^{2}-i\epsilon$ while,
on the right, the external momentum shift is used, $k_{0}\to k_{0}+i\epsilon$.
As predicted, using the mass shift yields a $\imaginary[B_{0}]$ function that
is symmetric with respect to $k_{0}$ while, using the external momentum shift,
$\imaginary[B_{0}]$ is antisymmetric. By considering the opposite sign in the
shifts, $M^{2}\to M^{2}+i\epsilon$ and $k_{0}\to k_{0}-i\epsilon$, gives an
overall minus factor in the imaginary parts, yielding a total of four possible
different scenarios for the imaginary part of $B_{0}$. In the following
numerical results, we always use the mass shift, $M^{2}\to M^{2}-i\epsilon$.

Although the specific choice for the pole shift changes the properties of the
imaginary part, it does not change the results for the cross sections in the
NJL model, for example. Indeed, in order to study cross sections in the NJL
model, the meson propagators are calculated which are dependent on the quark-
antiquark polarization function, as discussed in earlier sections of this
work. These quantities are calculated using the $B_{0}$ function, and display
the same symmetry properties as the $B_{0}$ integral. However, in the
differential cross section calculation, one only needs products of meson
propagators with the complex conjugate of other meson propagators, rendering
the calculation invariant under a specific overall sign of the imaginary part.
For more details see Section LABEL:cross_sections.

Figure 19: $B_{0}$ integral, in the vacuum, as a function of the zero
component of the external momentum, normalized by the cutoff, $k_{0}/\Lambda$
(left panel) and as a function of the magnitude of the external 3-momentum,
normalized by the cutoff, $\absolutevalue{\bm{k}}/\Lambda$ (right panel). The
full and dashed lines correspond to the real and imaginary parts part of the
$B_{0}$ integral, respectively. The cutoff is fixed to
$\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion masses are equal and
fixed to $M_{i}=M_{j}=0.4\Lambda$.

In Fig. 19, the $B_{0}$ function is shown as a function of the external
momentum, in the vacuum. On the left panel, it is shown as a function of the
zero component of the external momentum, $k_{0}/\Lambda$, while, on the right
panel, it is shown as a function of the magnitude of the external 3-momentum,
$\absolutevalue{\bm{k}}/\Lambda$. In the left panel, four different values for
$\absolutevalue{\bm{k}}$ are considered and analogously, in the right panel,
six values $k_{0}$ are fixed, making the information in both plots complement
each other. In the calculation of $B_{0}$ as a function of $k_{0}$, we only
show results for positive $k_{0}$ since the function is symmetric with respect
to $k_{0}$.

From the results of the $B_{0}$ as a function of $\absolutevalue{\bm{k}}$
(right panel) one can see the effect of the 3-momentum sphere intersection
regularization. The absolute value of the external 3-momentum,
$\absolutevalue{\bm{k}}$, is the distance between two spheres, each with
radii, $\Lambda$. At, $\absolutevalue{\bm{k}}=2\Lambda$ the spheres no longer
intersect and the $B_{0}$ function should be zero. Indeed that is the observed
behavior in the right panel of Fig. 19. In fact this is observed in all the
calculations shown in this work, even at finite temperature and finite
chemical potential, as we will see later. Moreover, we can observe that the
presence of an imaginary part is not constant, certain scenarios do not have a
finite imaginary part.

In the results concerning the dependence on $k_{0}$ (see left panel), we have
fixed $\absolutevalue{\bm{k}}$ to four specific values:
$\absolutevalue{\bm{k}}/\Lambda=\quantity{0.0,0.5,1.0,1.5}$. We do not
consider larger values because, for $\absolutevalue{\bm{k}}\geq 2\Lambda$,
within our approach, the $B_{0}$ function is automatically zero, due to the
regularization proposed in this work. Interestingly, we can also observe from
this plot that the starting $k_{0}$ point for a finite imaginary part,
$\imaginary[B_{0}]$, is $\absolutevalue{\bm{k}}$ dependent while, the ending
point is not. Also, differently from the previous case, at $k_{0}=2\Lambda$
the function remains finite. This is not surprising since the regularization
is made only over the loop’s 3-momentum while, the integration over the zero
component of the momentum inside the loop is made using the Matusbara
formalism, with the sum calculated from negative to positive infinity (see Eq.
(46)).

Figure 20: $B_{0}$ integral, for zero chemical potential, as a function of the
zero component of the external momentum, normalized by the cutoff,
$k_{0}/\Lambda$. The full and dashed lines correspond to the real and
imaginary parts part of the $B_{0}$ integral, respectively. The cutoff is
fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion masses are
equal and fixed to $M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix the
magnitude of the external 3-momentum,
$\absolutevalue{\bm{k}}/\Lambda=\quantity{0.0,0.5,1.0,1.5}$ and the
temperature was fixed to different values, corresponding to lines with
different colors.

Figure 21: $B_{0}$ integral, for zero chemical potential, as a function of the
magnitude of the external 3-momentum, normalized by the cutoff,
$\absolutevalue{\bm{k}}/\Lambda$. The full and dashed lines correspond to the
real and imaginary parts part of the $B_{0}$ integral, respectively. The
cutoff is fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion
masses are equal and fixed to $M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix
the magnitude of the zero component of the external momentum,
$k_{0}/\Lambda=\quantity{0.0,0.5,1.0,2.0}$ and the temperature was fixed to
different values, corresponding to lines with different colors.

In Figs. D and D, we show several results at finite temperature and zero
chemical potential ($\mu_{i}=\mu_{j}=0$). In the first, $B_{0}$ is displayed
as a function of $k_{0}$ while, in the latter, as a function of
$\absolutevalue{\bm{k}}$. Following the previous pattern, the cutoff and
masses are still fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$ and
$M_{i}=M_{j}=0.4\Lambda$, respectively. In Fig. D, specific values for the
ratio $\absolutevalue{\bm{k}}/\Lambda$ are fixed while, in Fig. D, some values
of $k_{0}$ are fixed instead.

For finite $\absolutevalue{\bm{k}}$ (observe Fig. D), the starting and ending
points for a finite imaginary part, are not temperature dependent. Indeed, in
each panel of this figure, the imaginary part stars and ends at the same
values of $k_{0}$. This can be understood from the equations which define the
imaginary parts: they can be non-zero exclusively inside precisely defined
kinematic regions, which are not temperature or chemical potential dependent,
see Eqs. (76) and (77).

When comparing the zero with the finite temperature cases, we can observe that
the there is an extra imaginary contribution, which does not exist at zero
temperature, present at low values of $k_{0}$, for finite values of
$\absolutevalue{\bm{k}}$. This extra contribution to the imaginary part, not
present in the vacuum, is caused by the presence of a finite contribution
coming from the scattering part of the function,
$\imaginary[B_{0,\mathrm{s}}]$. In the vacuum only the
$\imaginary[B_{0,\mathrm{p}}]$ is non-zero, since only pair creation and pair
annihilation processes are possible. However, at finite temperature, the
scattering channel opens. In these figures we can also observe the effect of
the regularization, independent of the temperature, at
$\absolutevalue{\bm{k}}=2\Lambda$, the $B_{0}$ function is zero.

Figure 22: $B_{0}$ integral, for zero temperature, as a function of the zero
component of the external momentum, normalized by the cutoff, $k_{0}/\Lambda$.
The full and dashed lines correspond to the real and imaginary parts part of
the $B_{0}$ integral, respectively. The cutoff is fixed to
$\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion masses are equal and
fixed to $M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix the magnitude of the
external 3-momentum,
$\absolutevalue{\bm{k}}/\Lambda=\quantity{0.0,0.5,1.0,1.5}$ and the chemical
potentials are equal and fixed to different values, corresponding to lines
with different colors.

Figure 23: $B_{0}$ integral, for zero temperature, as a function of the
magnitude of the external 3-momentum, normalized by the cutoff,
$\absolutevalue{\bm{k}}/\Lambda$. The full and dashed lines correspond to the
real and imaginary parts part of the $B_{0}$ integral, respectively. The
cutoff is fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion
masses are equal and fixed to $M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix
the magnitude of the zero component of the external momentum,
$k_{0}/\Lambda=\quantity{0.0,0.5,1.0,2.0}$ and the chemical potentials are
equal and fixed to different values, corresponding to lines with different
colors.

Next, we analyse the complementary physical scenario from the one discussed
above, considering zero temperature and finite chemical potential. For
simplicity, we considered degenerate chemical potentials, i.e.,
$\mu_{i}=\mu_{j}$. However, the formalism presented is perfectly suited to
study cases with different chemical potentials, $\mu_{i}\neq\mu_{j}$. Such
scenario arise in the study of neutron star matter and Heavy-Ion Collision
experiments where asymmetric nuclear/quark matter, is present. Following the
previous analysis, we show the $B_{0}$ integral as a function of
$k_{0}/\Lambda$ and $\absolutevalue{\bm{k}}/\Lambda$, in Figs. D and Fig. D,
respectively.

At zero temperature, only for values of chemical potential higher than the
fermions masses one gets results that differ from the vacuum result. This
characteristic is to be expected since at zero temperature, $T=0$, the Fermi-
Dirac distributions present in the $B_{0}$ integral act as step-functions.
Hence, below a certain energy regime, only the vacuum contributes: only when
the momenta is higher the Fermi’s momenta, for a fixed chemical potential,
does this momentum contributes to the momenta integration. This feature is
sometimes termed by Silver Blaze property [87, 88].

Previously, when discussing Fig. D, we argued that the values of $k_{0}$,
which define the intervals in which the imaginary part is finite, are neither
temperature or chemical potential dependent. However, Fig. D seems to
contradict this statement, since increasing the chemical potential drives the
onset of the imaginary part to larger values of $k_{0}$. However, this is a
feature of the Eqs. (72) and (73) themselves and not caused by the restriction
imposed by the endpoints of the integration. In these figures, as expected, we
continue to observe the general effect of the regularization, with $B_{0}=0$
for $\absolutevalue{\bm{k}}\geq 2\Lambda$.

Figure 24: $B_{0}$ integral, as a function of the zero component of the
external momentum, normalized by the cutoff, $k_{0}/\Lambda$. The full and
dashed lines correspond to the real and imaginary parts part of the $B_{0}$
integral, respectively. The cutoff is fixed to $\Lambda=1\leavevmode\nobreak\
\mathrm{GeV}$, the fermion masses are equal and fixed to
$M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix the magnitude of the external
3-momentum, $\absolutevalue{\bm{k}}/\Lambda=\quantity{0.0,0.5,1.0,1.5}$ and
the temperature and chemical potentials are fixed to different values,
corresponding to lines with different colors.

Figure 25: $B_{0}$ integral, for zero temperature, as a function of the
magnitude of the external 3-momentum, normalized by the cutoff,
$\absolutevalue{\bm{k}}/\Lambda$. The full and dashed lines correspond to the
real and imaginary parts part of the $B_{0}$ integral, respectively. The
cutoff is fixed to $\Lambda=1\leavevmode\nobreak\ \mathrm{GeV}$, the fermion
masses are equal and fixed to $M_{i}=M_{j}=0.4\Lambda$. In each panel, we fix
the magnitude of the zero component of the external momentum,
$k_{0}/\Lambda=\quantity{0.0,0.5,1.0,2.0}$ and the temperature and chemical
potentials are fixed to different values, corresponding to lines with
different colors.

Lastly, in Figs. D and D, we show the completely in-medium scenario, where
finite temperature and finite degenerate chemical potential are considered.
The general behavior is very similar to the ones found in the previously
discussed scenarios, presenting only quantitative differences.

In all scenarios considered in this section, a general pattern can be observed
in the real part of the integral, $\real[B_{0}]$: it changes behavior when
there is a switch on/off of the imaginary part, be it as a function of $k_{0}$
or as a function of $\absolutevalue{\bm{k}}$. This can be understood taking
into account the nature of the integration for a specific set of values of
$\\{k_{0},\absolutevalue{\bm{k}}\\}$. For the imaginary contribution to be
non-zero, there must exist a singularity inside the integration bounds.
Whenever one is dealing with a pair of values for
$\\{k_{0},\absolutevalue{\bm{k}}\\}$, which lays outside the integration
bounds, there should not be any singularity within the integration. When the
opposite happens, there is a singular point inside the integration region and
the now improper integral is assigned a value by the application of the
Cauchy-Principal value.

In the case where there is a model for the fermion masses, which encodes the
temperature and chemical potential dependence of the fermion masses, the
interpretation of the results gets more complicated. For instance, in the case
of the NJL model, one foresees that as temperature increases, at zero chemical
potential, the quark effective mass decreases. A similar pattern is expect at
finite chemical potential. In such cases, the behavior of the $B_{0}$ function
as a function of temperature and chemical potential is different and outside
the scope of this section. Here our aim was to simply analyze the
thermodynamic structure of the two fermion line integral $B_{0}$ by itself.

## References

  * Rehberg  _et al._ [1996a] P. Rehberg, S. Klevansky, and J. Hufner, Hadronization in the SU(3) Nambu-Jona-Lasinio model, Phys.Rev. C53, 410 (1996a), arXiv:hep-ph/9506436 [hep-ph] .
  * Costa  _et al._ [2004] P. Costa, M. C. Ruivo, Y. L. Kalinovsky, and C. A. de Sousa, Pseudoscalar mesons in hot, dense matter, Phys. Rev. C 70, 025204 (2004), arXiv:hep-ph/0304025 .
  * Costa  _et al._ [2005] P. Costa, M. C. Ruivo, C. A. de Sousa, and Y. L. Kalinovsky, Analysis of the U(A)(1) symmetry-breaking and restoration effects on scalar-pseudoscalar spectrum, Phys. Rev. D71, 116002 (2005), arXiv:hep-ph/0503258 .
  * Nagahiro  _et al._ [2006] H. Nagahiro, M. Takizawa, and S. Hirenzaki, eta- and eta-prime-mesic nuclei and U(A)(1) anomaly at finite density, Phys. Rev. C 74, 045203 (2006), arXiv:nucl-th/0606052 .
  * Costa  _et al._ [2009] P. Costa, M. C. Ruivo, C. A. de Sousa, H. Hansen, and W. M. Alberico, Scalar-pseudoscalar meson behavior and restoration of symmetries in SU(3) PNJL model, Phys. Rev. D79, 116003 (2009), arXiv:0807.2134 [hep-ph] .
  * Contrera  _et al._ [2010] G. A. Contrera, D. G. Dumm, and N. N. Scoccola, Meson properties at finite temperature in a three flavor nonlocal chiral quark model with Polyakov loop, Phys. Rev. D 81, 054005 (2010), arXiv:0911.3848 [hep-ph] .
  * Inagaki  _et al._ [2011] T. Inagaki, D. Kimura, H. Kohyama, and A. Kvinikhidze, Nonet meson properties in Nambu–Jona-Lasinio model with dimensional versus cutoff regularization, Phys. Rev. D 83, 034005 (2011), arXiv:1010.0309 [hep-ph] .
  * Blanquier [2014] E. Blanquier, Hadronization within a Nambu–Jona-Lasinio model with a Polyakov loop, Phys. Rev. C 89, 065204 (2014).
  * Blaschke  _et al._ [2017] D. Blaschke, A. Dubinin, A. Radzhabov, and A. Wergieluk, Mott dissociation of pions and kaons in hot, dense quark matter, Phys. Rev. D 96, 094008 (2017), arXiv:1608.05383 [hep-ph] .
  * Costa and Pereira [2019] P. Costa and R. C. Pereira, Phase Diagram, Scalar-Pseudoscalar Meson Behavior and Restoration of Symmetries in (2 + 1) Polyakov-Nambu-Jona-Lasinio Model, Symmetry 11, 507 (2019), arXiv:1904.05805 [hep-ph] .
  * Rehberg and Klevansky [1996] P. Rehberg and S. P. Klevansky, One loop integrals at finite temperature and density, Annals Phys. 252, 422 (1996), arXiv:hep-ph/9510221 .
  * Rehberg  _et al._ [1996b] P. Rehberg, S. P. Klevansky, and J. Hufner, Elastic scattering and transport coefficients for a quark plasma in SU-f(3) at finite temperatures, Nucl. Phys. A 608, 356 (1996b), arXiv:hep-ph/9607263 .
  * Zhuang  _et al._ [1995] P. Zhuang, J. Hufner, S. P. Klevansky, and L. Neise, Transport properties of a quark plasma and critical scattering at the chiral phase transition, Phys. Rev. D 51, 3728 (1995).
  * Sasaki and Redlich [2010] C. Sasaki and K. Redlich, Transport coefficients near chiral phase transition, Nucl. Phys. A 832, 62 (2010), arXiv:0811.4708 [hep-ph] .
  * Marty  _et al._ [2013] R. Marty, E. Bratkovskaya, W. Cassing, J. Aichelin, and H. Berrehrah, Transport coefficients from the Nambu-Jona-Lasinio model for $SU(3)_{f}$, Phys. Rev. C 88, 045204 (2013), arXiv:1305.7180 [hep-ph] .
  * Soloveva  _et al._ [2021] O. Soloveva, D. Fuseau, J. Aichelin, and E. Bratkovskaya, Shear viscosity and electric conductivity of a hot and dense QGP with a chiral phase transition, Phys. Rev. C 103, 054901 (2021), arXiv:2011.03505 [nucl-th] .
  * Vogl and Weise [1991] U. Vogl and W. Weise, The Nambu and Jona Lasinio model: Its implications for hadrons and nuclei, Prog.Part.Nucl.Phys. 27, 195 (1991).
  * Klevansky [1992] S. P. Klevansky, The Nambu-Jona-Lasinio model of quantum chromodynamics, Rev. Mod. Phys. 64, 649 (1992).
  * Hatsuda and Kunihiro [1994] T. Hatsuda and T. Kunihiro, QCD phenomenology based on a chiral effective Lagrangian, Phys. Rept. 247, 221 (1994), arXiv:hep-ph/9401310 .
  * Buballa [2005] M. Buballa, NJL model analysis of quark matter at large density, Phys. Rept. 407, 205 (2005), arXiv:hep-ph/0402234 .
  * [21] G. Ripka, Quarks bound by chiral fields: The quark-structure of the vacuum and of light mesons and baryons, oxford, UK: Clarendon Pr. (1997) 205 p.
  * Osipov  _et al._ [2006] A. A. Osipov, B. Hiller, and J. da Providencia, Multi-quark interactions with a globally stable vacuum, Phys. Lett. B634, 48 (2006), arXiv:hep-ph/0508058 .
  * Osipov  _et al._ [2007] A. A. Osipov, B. Hiller, A. H. Blin, and J. da Providencia, Effects of eight-quark interactions on the hadronic vacuum and mass spectra of light mesons, Annals Phys. 322, 2021 (2007), arXiv:hep-ph/0607066 .
  * Osipov  _et al._ [2013a] A. Osipov, B. Hiller, and A. Blin, Light Quark Masses in Multi-Quark Interactions, Eur.Phys.J. A49, 14 (2013a), arXiv:1206.1920 [hep-ph] .
  * Osipov  _et al._ [2013b] A. Osipov, B. Hiller, and A. Blin, Effective multiquark interactions with explicit breaking of chiral symmetry, Phys.Rev. D88, 054032 (2013b), arXiv:1309.2497 [hep-ph] .
  * Reinhardt and Alkofer [1988] H. Reinhardt and R. Alkofer, Instanton Induced Flavor Mixing in Mesons, Phys.Lett. B207, 482 (1988).
  * Câmara Pereira  _et al._ [2020] R. Câmara Pereira, J. a. Moreira, and P. Costa, The strange critical endpoint and isentropic trajectories in an extended PNJL model with eight Quark interactions, Eur. Phys. J. A 56, 214 (2020), arXiv:2006.02385 [hep-ph] .
  * Pereira [2021] R. C. Pereira,  _Quantum Chromodynamics Phase Diagram Under Extreme Conditions_ , Ph.D. thesis, Coimbra U. (2021).
  * Hubbard [1959] J. Hubbard, Calculation of partition functions, Phys. Rev. Lett. 3, 77 (1959).
  * Stratonovich [1957] R. L. Stratonovich, On a Method of Calculating Quantum Distribution Functions, Soviet Physics Doklady 2, 416 (1957).
  * Ebert [1998] D. Ebert, Bosonization in particle physics, Lect. Notes Phys. 508, 103 (1998), arXiv:hep-ph/9710511 .
  * Ball [1989] R. Ball, Chiral gauge theory, Physics Reports 182, 1 (1989).
  * Novikov  _et al._ [1984] V. A. Novikov, M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Calculations in External Fields in Quantum Chromodynamics. Technical Review, Fortsch. Phys. 32, 585 (1984).
  * Fraser [1985] C. M. Fraser, Calculation of Higher Derivative Terms in the One Loop Effective Lagrangian, Z. Phys. C 28, 101 (1985).
  * Aitchison and Fraser [1985] I. J. R. Aitchison and C. M. Fraser, Derivative Expansions of Fermion Determinants: Anomaly Induced Vertices, Goldstone-Wilczek Currents and Skyrme Terms, Phys. Rev. D 31, 2605 (1985).
  * Oertel  _et al._ [2001] M. Oertel, M. Buballa, and J. Wambach, Meson loop effects in the NJL model at zero and nonzero temperature, Phys. Atom. Nucl. 64, 698 (2001), arXiv:hep-ph/0008131 .
  * Nikolov  _et al._ [1996] E. N. Nikolov, W. Broniowski, C. V. Christov, G. Ripka, and K. Goeke, Meson loops in the Nambu-Jona-Lasinio model, Nucl.Phys. A608, 411 (1996), arXiv:hep-ph/9602274 [hep-ph] .
  * Costa  _et al._ [2003] P. Costa, M. C. Ruivo, and Y. L. Kalinovsky, Pseudoscalar neutral mesons in hot and dense matter, Phys. Lett. B 560, 171 (2003), arXiv:hep-ph/0211203 .
  * Blaschke  _et al._ [2014] D. Blaschke, M. Buballa, A. Dubinin, G. Roepke, and D. Zablocki, Generalized Beth–Uhlenbeck approach to mesons and diquarks in hot, dense quark matter, Annals Phys. 348, 228 (2014), arXiv:1305.3907 [hep-ph] .
  * Blaschke and Ebert [2017] D. Blaschke and D. Ebert, Variational path-integral approach to back-reactions of composite mesons in the Nambu–Jona-Lasinio model, Nucl. Phys. B 921, 753 (2017), arXiv:1703.08964 [hep-ph] .
  * Haber [2023] H. Haber, Notes on the matrix exponential and logarithm (2023), [Online; last downloaded at 01-September-2023].
  * Bernard  _et al._ [1988] V. Bernard, R. L. Jaffe, and U. G. Meissner, Strangeness Mixing and Quenching in the Nambu-Jona-Lasinio Model, Nucl. Phys. B308, 753 (1988).
  * Saionji  _et al._ [2023] T. Saionji, D. Jido, and M. Harada, Inverse mass ordering of light scalar mesons in the Nambu–Jona-Lasinio model, PTEP 2023, 033D01 (2023), arXiv:2210.01993 [hep-ph] .
  * Lutz  _et al._ [1992] M. F. M. Lutz, S. Klimt, and W. Weise, Meson properties at finite temperature and baryon density, Nucl. Phys. A 542, 521 (1992).
  * Ebert and Reinhardt [1986] D. Ebert and H. Reinhardt, Effective chiral meson lagrangian with weinberg and ksfr relations from microscopic four-quark interactions, Physics Letters B 173, 453 (1986).
  * Alkofer and Reinhardt [1995] R. Alkofer and H. Reinhardt,  _Chiral Quark Dynamics_, Lecture Notes in Physics Monographs (Springer Berlin Heidelberg, 1995).
  * Zuk [1986] J. A. Zuk, Expansion of the one-loop effective action in covariant derivatives, Phys. Rev. D 33, 3645 (1986).
  * Osipov and Hiller [2001a] A. A. Osipov and B. Hiller, Generalized proper-time approach for the case of broken isospin symmetry, Phys. Rev. D63, 094009 (2001a), arXiv:hep-ph/0012294 .
  * Osipov and Hiller [2001b] A. A. Osipov and B. Hiller, Large mass invariant asymptotics of the effective action, Phys.Rev. D64, 087701 (2001b), arXiv:hep-th/0106226 [hep-th] .
  * Salcedo [2001] L. L. Salcedo, Generalized heat kernel coefficients, Eur. Phys. J. direct 3, 14 (2001), arXiv:hep-th/0107133 .
  * Osipov and Volkov [2017] A. Osipov and M. Volkov, Chiral transformations of spin-1 mesons in the non-symmetric vacuum, Annals of Physics 382, 50 (2017).
  * Osipov  _et al._ [2004] A. Osipov, H. Hansen, and B. Hiller, Long distance expansion for the njl model with su(3) and ua(1) breaking, Nuclear Physics A 745, 81 (2004).
  * Vassilevich [2003] D. V. Vassilevich, Heat kernel expansion: User’s manual, Phys. Rept. 388, 279 (2003), arXiv:hep-th/0306138 .
  * Dunne [2008] G. V. Dunne, Functional determinants in quantum field theory, J. Phys. A 41, 304006 (2008), arXiv:0711.1178 [hep-th] .
  * Osipov and Hiller [2001c] A. A. Osipov and B. Hiller, Inverse mass expansion of the one-loop effective action, Phys. Lett. B515, 458 (2001c).
  * Moreira  _et al._ [2015] J. Moreira, J. Morais, B. Hiller, A. A. Osipov, and A. H. Blin, Strange quark matter in the presence of explicit symmetry breaking interactions, Phys. Rev. D91, 116003 (2015), arXiv:1409.0336 [hep-ph] .
  * Meissner  _et al._ [1990] T. Meissner, E. R. Arriola, and K. Goeke, Regularization scheme dependence of vacuum observables in the Nambu-Jona-Lasinio model, Zeitschrift fur Physik A Hadrons and Nuclei 336, 91 (1990).
  * Kahana and Lavelle [1993] D. Kahana and M. Lavelle, On the Pauli-Villars regularization scheme in the NJL model, Phys. Lett. B 298, 397 (1993).
  * Inagaki  _et al._ [2012] T. Inagaki, D. Kimura, H. Kohyama, and A. Kvinikhidze, Nonet meson properties in Nambu-Jona-Lasinio model with dimensional regularization at finite temperature and chemical potential, Phys. Rev. D 85, 076002 (2012), arXiv:1110.5898 [hep-ph] .
  * Inagaki  _et al._ [2013] T. Inagaki, D. Kimura, H. Kohyama, and A. Kvinikhidze, Regularization parameter independent analysis in Nambu-Jona-Lasinio model, Int. J. Mod. Phys. A 28, 1350164 (2013), arXiv:1302.5667 [hep-ph] .
  * Kohyama  _et al._ [2015] H. Kohyama, D. Kimura, and T. Inagaki, Regularization dependence on phase diagram in Nambu–Jona-Lasinio model, Nucl. Phys. B 896, 682 (2015), arXiv:1501.00449 [hep-ph] .
  * Torres Bobadilla  _et al._ [2021] W. J. Torres Bobadilla _et al._ , May the four be with you: Novel IR-subtraction methods to tackle NNLO calculations, Eur. Phys. J. C 81, 250 (2021), arXiv:2012.02567 [hep-ph] .
  * Xue  _et al._ [2022] K. Xue, X. Yu, and X. Wang, Do we need to use regularization for the thermal part in the njl model? *, Chinese Physics C 46, 013103 (2022).
  * Costa  _et al._ [2010] P. Costa, H. Hansen, M. C. Ruivo, and C. A. de Sousa, How parameters and regularization affect the PNJL model phase diagram and thermodynamic quantities, Phys. Rev. D 81, 016007 (2010), arXiv:0909.5124 [hep-ph] .
  * Hansen  _et al._ [2020] H. Hansen, R. Stiele, and P. Costa, Quark and Polyakov-loop correlations in effective models at zero and nonvanishing density, Phys. Rev. D 101, 094001 (2020), arXiv:1904.08965 [hep-ph] .
  * Costa  _et al._ [2020] P. Costa, R. Câmara Pereira, and C. Providência, Role of the conserved charges in the chiral symmetry restoration phase transition, Phys. Rev. D 102, 054010 (2020), arXiv:2009.01781 [hep-ph] .
  * Kapusta [1989] J. I. Kapusta,  _Finite-temperature field theory_ (Cambridge University Press, 1989).
  * Fujii and Ohtani [2004] H. Fujii and M. Ohtani, Sigma and hydrodynamic modes along the critical line, Phys. Rev. D 70, 014016 (2004), arXiv:hep-ph/0402263 .