{'fraction_non_alphanumeric': 0.046032424866840324, 'fraction_numerical': 0.020566851988647408, 'mean_word_length': 3.8901140684410644, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
One Loop Integrals at Finite Temperature and Density

The technique of decomposing Feynman diagrams at the one loop level into elementary integrals is generalized to the imaginary time Matsubara formalism. The three lowest integrals, containing one, two and three fermion lines, are provided in a form that separates out the real and imaginary parts of these complex functions, according to the input arguments, in a fashion that is suitable for numerical evaluation. The forms given can be evaluated for arbitrary values of temperature, particle mass, particle momenta and chemical potential.


I. INTRODUCTION
The imaginary time Matsubara formalism for Green functions [1] is a venerated method of dealing with problems of finite temperature. This is so because of the fact that Wick's theorem can easily be shown to hold for the Matsubara operators.
As such, this method of dealing with finite temperature systems has found wide application, primarily in condensed matter physics [2], but also recently in problems relevant for high energy physics. In particular, it has been recently used in studies of effective QCD Lagrangians (see, e. g. the reviews [3] and [4] and references cited therein), through which an understanding of hot and dense matter is sought [5][6][7][8][9]. However, to date most such studies have often been constrained by specific parameter choices to enable tractability -often degenerate particle masses, zero chemical potential or special kinematics are chosen. Furthermore, often only the principal values of these integrals are considered, with the complex nature being completely ignored. An extension of those calculations to physically relevant cases with arbitrary parameters and generalized kinematics, as may be appropriate for example, in calculating cross sections in an SU (2) or SU(3) effective model of QCD [7][8][9], becomes technically extremely complicated. It is thus the purpose of this paper to address the technical problems associated with such integrals, and present them in a form suitable for the numerical evaluation of both their real and imaginary components. This work has come about through a study of hadronization processes within the Nambu-Jona-Lasinio model [9], and is intended to provide the technical knowledge required for such functions for researchers who are active in this field.
The associated computer code conformal to the analysis presented here is freely For calculations at zero temperature and density, it has been shown [10], that all occurring transition amplitudes at the one loop level can be decomposed into a number of elementary integrals, which in turn can be classified by the number of particle lines they contain. Since this decomposition technique can be generalized to the Matsubara formalism, the aim of this paper is to provide suitable forms for the elementary integrals containing one, two and three fermion lines, which we term A, B 0 and C 0 , respectively, which can be used for numerical evaluation. (The calculation can be easily taken over to boson loops on replacing the Fermi distribution by the Bose distribution.) We analyze these integrals, explicitly displaying their complex nature that is a function of the input parameters. The parameters that may then be specified can be of the most general form, i. e. arbitrary temperatures, particle masses, chemical potentials and kinematics. As such, the associated computer code can be treated as a library-like 'black box'.
We note at this point that the calculation of the one loop integral A is straightforward. Its evaluation here serves merely to remind the reader of the usage of the Matsubara formalism and it also renders this work complete. The evaluation of B 0 , on the other hand, serves to demonstrate the complexities involved in a consistent evaluation of this function. Finally we deal with the three line integral C 0 , which then is highly nontrivial and complicated.
It is shown in this paper that each of integrals B 0 and C 0 (and thus implicitly all higher one loop integrals) can be written as a sum over terms that are constructed from a single generic function that is taken at different values of its arguments.
The analytical structure of these generic functions is connected with the number of nonvanishing external momenta. We are thus required to distinguish between several cases in evaluating the generic functions when constructing a physical quantity.
As is usually the case in relativistic quantum mechanics, the one and two fermion line integrals A and B 0 are divergent, which means that they have to be regularized.
We do this using a three momentum cutoff Λ, which means that we integrate over all momenta with | p| < Λ. The three fermion line integral C 0 is however convergent.
We nevertheless perform our analysis of this function including a three momentum cutoff. The reasons for doing this are twofold: (i) the extension to the case Λ → ∞ is simple to obtain and provides a check of our calculations, and (ii) for models such as that of Nambu and Jona-Lasinio extended to QCD, it is regarded as internally consistent within the model to restrict all quark momenta in all loops.
Having dealt with the analytical complexities of the loop integrals, we demonstrate in the last section the usage of the integral decomposition in a simple physical example, viz. that of calculating the kaon masses in the SU(3) Nambu-Jona-Lasinio model. This involves the integrals A and B 0 . Use of C 0 may be found, for example, in Ref. [9]. This paper is organized as follows: in Sec. II, we define the functions A, B 0 and C 0 . These functions are brought into a form suitable for numerical integration in Secs. III, IV and V, respectively. At the end of Sec. IV, we also give some technical hints for the evaluation of our results. In Sec. VI, we give an example for the application of our calculations. Further examples can also be found in Ref. [9].
In Sec. VII we summarize and conclude. Some special formulae are proven in an appendix.
We define the functions A, B 0 and C 0 to be the analytical continuation of the one, two and three fermion line integrals, , and C 0 (k, q, δ k, q , m 1 , µ 1 , m 2 , µ 2 , iν m , m 3 , µ 3 , iα l , β, Λ) = (2.3) , in which the complex frequencies iν m and iα l are to be continued to their values on the real axis after the Matsubara summation on n is carried out. It is to be understood that the limit η → 0 is to be taken after the Matsubara summation. In the above expressions, we have introduced the abbreviations where the symbols m i and µ i denote the masses and chemical potentials of the fermions running on the individual lines of the diagram. The m i are to be understood to include a negative infinitesimal imaginary part, m i → m i −iǫ, ǫ > 0. In Eqs. (2.1) to (2.3), β is the inverse temperature, β = 1/T . The three momentum cutoff is denoted by Λ, and all integrals are considered for | p| ≤ Λ.
The complex four momenta (iν m ; k) and (iα l ; q), are the four momenta that enter the loop. The frequencies iν m and iα l are bosonic in nature, while the frequencies iω n are fermionic, After analytical continuation, iν m and iα l become the zeroth components associated with the three momenta k and q, i. e. iν m → k 0 , iα l → q 0 . One may note at this point that, due to rotational invariance, the functions B 0 and C 0 do not fully depend on both k and q; B 0 depends only on k = | k|, while C 0 depends only on k, q = | q| and the enclosed angle δ k, q , as is implied by the arguments specified on the left hand side of Eqs. (2.2) and (2.3).

III. THE CALCULATION OF A
The calculation of A is simple and straightforward. After evaluating the Matsubara sum by contour integration in the usual fashion [2], one has where we have introduced the Fermi distribution function f (x) = 1 exp(βx) + 1 (3.2) and the energy cutoff where κ = max(m, |µ|) and κ ′ = √ κ 2 − m 2 .

A. Matsubara Summation
The two line integral B 0 (k, m 1 , µ 1 , m 2 , µ 2 , iν m , β, Λ) can be analyzed in a similar fashion to A(m, µ, β, Λ). The denominator of Eq. (2.2) has the four zeros from which we obtain, after the summation over n, the terms (4.2) To eliminate the angular dependence in the arguments of the Fermi functions, we carry out the substitution in the third and fourth terms of Eq. (4.2) and obtain where the meaning of E 1 and E 2 has changed to become We have also introduced Note that λ will be real after the analytical continuation. One should further note that the substitution Eq. (4.3) has been performed under the assumption that Λ ≫ k. Thus no corresponding shift of variables in the integration limits has been carried out at this point.
Without loss of generality, we may choose the coordinate system in such a way that k points in the z direction. Then the integration over the angle φ yields a factor 2π. It is now useful to consider each term occurring in Eq. (4.4) separately. These correspond to the individual poles in the original expression and have a common generic structure. One may therefore decompose B 0 into the four terms in terms of the generic functioñ Numerical computation of B 0 thus relies purely on the evaluation ofB ± 0 . The analytical structure of this integral differs for k = 0 and k > 0, so that these two cases must be considered separately. We examine these in the following two subsections. We may also note that the symmetry relation can be inferred from Eq. (4.4).
B. Calculation ofB ± 0 for k = 0 For k = 0, the integral over x is trivial, and we obtaiñ (4.10) (Here and in the following, we drop the arguments ofB ± 0 for convenience.) The analytical continuation of the Matsubara frequency iν m to real values leads to poles in the integrand at E = E 0 , where To calculate the integral in this case, we recall that the particle where P denotes the Cauchy principal value. From this, we find where p 0 = E 2 0 − m 2 and the Θ function ensures that the imaginary part appears only if E 0 lies within the integration interval [m, Λ E ]. Figure 2 shows B 0 as a function of k 0 /Λ for the parameter set k = 0, µ 1 = µ 2 = 0, One recognizes that at small k 0 , B 0 is a real function, whereas for k 0 ≥ 2m an imaginary part emerges. At k 0 > 2Λ E , the imaginary part vanishes due to the presence of the cutoff.
For k > 0, the generic functionB ± 0 has the form , (4.14) where we have already indicated that the infinitesimal imaginary part in the integrand is determined by the sign of λ. In the limit ǫ → 0, the integrand contains a singularity at On performing the x integration and applying Eq. (4.12), one obtains The integral for the imaginary part can be done analytically: The integrand for the real part in Eq. (4.17) contains logarithmic poles at the energies E = E 1/2 which are defined by These energies also constitute the limits of integration for the imaginary part. Squaring Eq. (4.19) leads to a quadratic equation for E 1/2 , which has the solutions 2. One solution E 1 lies inside the integration interval, the other one, E 2 lies outside. It follows that the integral for the imaginary part has to be performed from E 1 to Λ E in this case.
3. The two solutions E 1 < E 2 lie within the integration interval. In this case, the integral for the imaginary part has to be performed from E 1 to E 2 .
Since we will deal with integration limits of this type very often, we will adopt the convention that the solutions E 1/2 of Eq. (4.20) are tailored as described above.
Eq. (4.17) then takes the form Together with Eq. (4.20), Eq. (4.21) provides a sufficient basis for the numerical calculation ofB ± 0 . It can be shown that in the limit k → 0, Eq. (4.21) continously approaches Eq. (4.13). Thus the generic functionB ± 0 is determined, and consequently B 0 via Eq. (4.7). Figure 3 shows B 0 for the same parameter set as was used in Fig. 2, but with k = 100MeV. For k 0 < k, Eq. (4.20) has one solution, from which one obtains an imaginary part. For k 0 > 2m, Eq. (4.20) has two solutions, which become cut off at high k 0 , so that the imaginary part continously goes to zero. Figure 4 shows another special case. Here we plot B 0 as a function of k for k 0 = 0. In this case, we have no imaginary part at all [7].

D. Programming Considerations
At this point, we deviate from the analytical calculation to discuss some of the basic aspects of programming such integrals. Comments made here are relevant for the B 0 and C 0 integrals, the first of which is the least complicated integral containing multiple singularities. We discuss this now, bearing in mind that similar techniques are required for the more complex C 0 integral.
To calculate B 0 for the most general case via the decomposition (4.7), one has to take several things into consideration. First, the type of singularity determines the type of integration routine that should be utilized. One requires three such routines: (i) one routine which integrates smooth functions, (ii) one which integrates functions containing logarithmic or other integrable singularities, and (iii) one which computes the Cauchy principal value of an integral containing a 1/x singularity. Examples for this can be found in the literature [11,12].
On entry, the program has to decide, which of the forms (4.13), pertinent to the case k = 0, or (4.21), pertinent for k > 0, to use. It is necessary to handle the two cases separately, since (i) Eq. (4.13) and Eq. (4.21) have a different analytical structure, and (ii) Eq. (4.21) contains a factor 1/k. In actual practice, the case k = 0 is needed very often. After that, the position of singularities can be computed according to Eq. (4.11) or Eq. (4.20), respectively. The appropriate integration routine has to be chosen, depending on the result of this computation.
Note that although the decomposition as given Eq. (4.7) is the most general one possible, it does not always provide the optimal way of calculating B 0 . If, for example, one knows a priori that an actual calculation will be confined to m 1 = m 2 , µ 1 = µ 2 and k = 0, the four terms of Eq. (4.7) can be combined, yielding the simplified expression The result of this equation stays finite in the limit k 0 → 0. The functionsB ± 0 , however, individually diverge in this limit, so that a computation of B 0 by decomposition suffers from cancellations at small k 0 . Note that similar simplifications follow for

A. Matsubara Summation
In this section, we examine the analytical structure of the three fermion line integral C 0 , and decompose it into a sum involving a single generic function.
The three fermion line function given in Eq. (2.3) has poles at with E 1 , E 2 and E 3 as specified in Eq. (2.4). In analogy to Eq. (4.6), we define the complex variables We again note that these parameters become real numbers after the analytical continuation is performed.
After summing the Matsubara frequencies, one obtains a sum over six terms that arise from each of the poles, It is useful to make the substitution p → k− p in the third and fourth terms, p → q− p in the fifth and sixth terms, in order to get rid of the angular dependence in the arguments of the Fermi functions. Again this is performed under the assumption that Λ ≫ k, q. After this, we obtain a decomposition of C 0 of the form where we have introduced the three fermion line generic functioñ In analogy to Eq. (4.9), one can derive a symmetry relation for the function C 0 , as being As was the case for the two line generic functionB ± 0 ,C ± 0 is calculated for several special cases. We list these in order of increasing difficulty: 1. The case k = q = 0.
4. The case k > 0, q > 0, p and k not collinear -general case.
In the following sections, we will derive results for each of these special cases.
In this case, Eq. (5.5) reduces tõ where we have again omitted the arguments ofC ± 0 for convenience and have included a small imaginary part in the denominators to do the analytical continuation. The constants E 1 and E 2 are defined by We suppose therefore that if one of them lies inside the integration interval. After taking the limit ǫ → 0, we obtaiñ  Without loss of generality, we assume that q points in the z direction. With this assumption, the scalar product p q occurring in Eq. (5.5) simplifies to pq cos θ. The integration over φ becomes trivial and one obtains where now and E 1 is as defined in Eq. (5.8a). The angular integral is singular if |λ 2 (E − E 2 )| < pq. This equation has exactly the same structure as Eq. (4.16), so Eq. (4.20) and the following discussion can be applied here too. We label the endpoints of the singular interval E 21 and E 22 . The angular integral becomes ForC ± 0 we have accordinglỹ In this expression, we again have to take the limit ǫ → 0. The final expression we find is (5.14) A numerical evaluation of Eq. (5.14) has to take into account Cauchy singularities at E = E 1 and (integrable) logarithmic singularities at E = E 21 and E = E 22 .
D. Calculation ofC ± 0 for k > 0, q > 0, q and k collinear Again, we assume without loss of generality that q and k have the form where η = ±1. With this assumption, we have p q = pq cos θ and p k = ηpk cos θ in Eq. (5.5). The φ integration is again trivial, and one obtains and E 2 from Eq. (5.11). The integral over x becomes (cf. [13], integral No. 12.8) In what follows, we label the energies for which is fulfilled E 1 1/2 , the energies for which is fulfilled E 2 1/2 . These can be calculated explicitly via Eq. (4.20). Note that the discussion following Eq. (4.20) applies here too.
The second and third terms of Eq. (5.18) still contain a singularity at This is, however, not true for the first term, since the logarithmic factor goes to zero at E = E 0 , making the integrand continous here. So we again add a small imaginary part in the denominators where appropriate and obtaiñ For the numerical evaluation of the real part, one has to take into account the logarithmic poles at E = E 1 1/2 and E = E 2 1/2 . For the imaginary part, the Cauchy singularity at E = E 0 must be integrated.
In Fig. 6, we show m 2 C 0 as a function of k 0 /Λ for q = q 0 = k 0 /2 and all other parameters as in Fig. 5. This special case has been considered previously in Ref. [8].
As can be seen from Eq. (5.4), the computation of C 0 requires a computation ofC ± 0 for k = 0, q > 0, as well as for k, q collinear.

E. Calculation ofC ± 0 for the General Case
In this subsection, we handle the general case of arbitrary momenta in the evaluation ofC ± 0 from Eq. (5.5). We choose the coordinate system in such a way that From this choice we obtain p k = pk(cos θ cos δ + sin θ sin δ cos φ) (5.25a) p q = pq cos θ .
It is useful to introduce the abbreviations In terms of these variables,C ± 0 takes the form .
The additional difficulty introduced by the generality of handling arbitrary momenta is evident in the φ integration. As will be seen in the following, this can however still be dealt with analytically.
To perform the φ integration, we use the formula , (5.28) which we prove in Appendix A. With this formula, one immediately obtains the result In the following subsections, we examine the remaining angular integrals that are required.
1. The θ Integration for ∆ 1 < 0 In the next step, we calculate the integral We first note that for x = ±1, ∆ 1 is always nonnegative. ∆ 1 has zeros at Consequently two conditions must be fulfilled in order to obtain ∆ 1 ≤ 0: (i) |c 1 | < 1 and (ii) x 1 ≤ x ≤ x 2 . From this it follows that we can apply the formula where the factor sgn(x + c 2 ) has been taken at x = (x 1 + x 2 )/2 = −c 1 cos δ. ∆ 0 can be obtained from ∆ 1 by substituting −c 2 for x: It remains to compute the integral In doing so, one easily notes that c 1 + x cos δ vanishes only at points where ∆ 1 ≤ 0, so that the sgn term can be treated effectively as a constant. The real part can be computed from [14], integral No. 231.10. Bearing the definition of x 1/2 from Eq. (5.32) in mind, one obtains and F is either one of the functions Eq. (5.37) can be further simplified by Up until this point, the calculation in this section has been purely analytical.
However, although F + and F − are equivalent in specifying the indefinite integral pertinent to Eq. (5.36), they are not equally well behaved numerically if either |c 1 | or |c 2 | is approximately equal to one. It is therefore necessary at this point to devise a numerically stable procedure, and we give our algorithm here. To this end, it is useful to define Ξ as Using the identity one can show that the following algorithm gives the correct result: This is implemented in the numerical procedure.

Putting the Parts together
The final expression forC ± 0 is now given. We may writẽ where Ξ has to be taken from Eq. (5.43). In the spirit of Sec. IV D, in order to evaluate this integral, one has to identify the singular points of the integrand. It is easy to see that the following points are singularities: • Energies, at which |c 1 | = 1 or |c 2 | = 1. These can be calculated from Eq. (4.20).
• Energies, at which ∆ 0 = 0. To compute these energies, one notes that c 1 and c 2 can be written as where a i and b i can be obtained from Eq. (5.26). With these constants, the equation ∆ 0 = 0 can be cast into the form which can be easily solved for E.
• Energies, at which |c 1 | = | cos δ|. These can be computed by All of these points correspond to integrable singularities. On closer inspection, one also obtains that the energies, at which the arguments of the logarithmic function (in F ± ) vanish, or the energies, at which the argument of the arccosine (in G) becomes greater than one, are included in these cases.
To illustrate C 0 for the general case, we plot m 2 C 0 in Fig. 7 as a function of k 0 /Λ with k = 3/4k 0 , q = 3/4q 0 and δ = π/6. All other parameters are as in Fig. 5.

VI. APPLICATION: KAON MASSES IN THE NJL MODEL
In this section, we give a specific example for a calculation using our elementary integrals. It is also intended for the benefit of our reader, who wishes to compare his/her results with our numbers. To this end, we compute the kaon masses in the framework of the SU(3) Nambu-Jona-Lasinio model [3,4,9] as a function of temperature. In this model, the quark masses have to be determined from the coupled gap equations with the current quark masses m 0q , m 0s , the number of colors N c and the coupling constants G and K. S f denotes the finite temperature quark propagator The trace of the propagator can be expressed in terms of the function A: where (k 0 , k) is the total four momentum of the quarks, Π(k 0 , k) the irreducible polarization and K eff an effective coupling constant, which is given by The irreducible polarization Π(k 0 , k) is computed from the diagram in Fig. 8.
where E q = p 2 + m 2 q and E s = ( p − k) 2 + m 2 s . To cast this into a form which contains the functions A and B 0 , one uses the identity After continuing iν m to k 0 , one immediately obtains the result The kaon mass is now computed using the dispersion relation Using the numbers m 0q = 5.5MeV, m 0s = 140.7MeV, µ q = µ s = 0, Λ = 602.3MeV, GΛ 2 = 1.835 and KΛ 5 = 12.36, one obtains the zero temperature results m q = 367.7MeV and m s = 549.5MeV from Eq. (6.4). Solving Eq. (6.10) with these numbers yields m K = 497.7MeV. The temperature dependence of m q , m s and m K is computed using the same formulae and the result is depicted in Fig. 9. Further examples can be found in Ref. [9].

VII. SUMMARY AND CONCLUSION
In this paper, we have presented the technical aspects required for a numerical evaluation of the one, two and three fermion line integrals at finite temperature, structuring each as a sum of terms of a specific generic integral taken at different values of its arguments, that are often required for field theoretic calculations such as within the Nambu-Jona-Lasinio model. Both real and imaginary parts of these functions are explicitly calculated. Concomitantly, we have illustrated all functions graphically for certain parameter sets, in order to enable the user to verify his/her own calculation. A simple example, employing two of the integrals has also been shown. A computer source code is also available that routinizes these integrals. The generalization to bosonic line integrals can also be seen to follow in an analogous fashion.

ACKNOWLEDGMENTS
We would like to thank J. Hüfner for his encouragement and support of this project, and those many people who have requested that we make this work generally available. We are also indebted to the Bundesministerium für Bildung und Our assertion is The real part of this equation can be easily obtained from integral tables. (See e. g. [14], integral No. 331.41.) The singularity in the integral only appears for |b| > |a|, which gives the Θ function for the imaginary part. In this case, singularities appear at φ = arccos(−a/b) and The first of these singularities appears in the interval [0, π], the second one in the interval [π, 2π]. We split the integral into two parts and focus on the first of these two parts. Let φ 0 = arccos(−a/b) and the function f (φ) be defined by It can be easily seen that f (φ) is continous in the interval [0, π] and since f (φ 0 ) is positive. Our standard formula Eq. (4.12) can be applied to this integral to yield The second integral on the right hand side of Eq. (A3) can be treated in the same fashion and is found to give the same contribution. QED.
The Θ function on the right hand side of Eq. (B2) emerges from the fact that the singularity occurs only if The real part of Eq. (B2) can be looked up in integral tables. (See e. g. [13], integral No. 221.7.) For ∆ > 0, one obtains For ∆ < 0, one obtains Combining Eqs. (B2), (B5) and (B7), proves the assertion.