Equilibrium thermodynamic susceptibilities for a dense degenerate Dirac field

Parity preserving relativistic fluids in four spacetime dimensions admit seven independent thermodynamic susceptibilities at the second order in the hydrodynamic derivative expansion. We compute all parity-even second order thermodynamic susceptibilities for a free massive Dirac field at zero temperature and a nonzero chemical potential, based on the Kubo formulas reported in P. Kovtun and A. Shukla, Kubo formulas for thermodynamic transport coefficients, J. High Energy Phys. 10 (2018) 007. We also compute the second order constitutive relations for the energy-momentum tensor and the conserved current in the absence of external gauge fields.


I. INTRODUCTION
Hydrodynamics provides an effective low energy description of physical systems near thermal equilibrium [2][3][4]. The description is valid on length scales much larger than the typical microscopic length scales associated with the system, such as the mean free path or the correlation length. In the hydrodynamic regime, the conserved currents of the system admit a derivative expansion in terms of the hydrodynamic variables such as temperature, chemical potential, fluid velocity etc. The dynamical equations of hydrodynamics express the conservation of the energy-momentum tensor and other conserved currents in the system.
At each order in the hydrodynamic derivative expansion one can have several transport coefficients quantifying the response of the system to external perturbations. Some of these transport coefficients are dissipative in nature, contributing to entropy production out of equilibrium, whereas others are non-dissipative and exist even in the limit of thermal equilibrium [5][6][7]. Interestingly, as was first noted in [8,9], the nondissipative transport coefficients follow from an equilibrium generating functional which itself admits a derivative expansion, with various thermodynamic susceptibilities appearing as coefficients of different terms in the expansion. The variation of the equilibrium generating functional with respect to external sources gives the constitutive relations for the conserved currents in the system, which contain various transport coefficients that do not vanish in the equilibrium limit. Thus, in thermal equilibrium, one can think of the susceptibilities entering the generating functional to be the more fundamental objects, and the transport coefficients appearing in the constitutive relations to be derived entities.
For relativistic fluids with a conserved current in equilibrium in four spacetime dimensions there are no thermodynamic susceptibilities at the first order in derivatives, whereas at the second order there are nine [8]. Out of the nine susceptibilities, seven are parity even. Kubo formulas for all the second order thermodynamic susceptibilities in terms of the correlations functions of the conserved currents were derived in [1]. As it turns out, the parity-even susceptibilities can be computed using zero-frequency flat space two-point correlation functions of the energy-momentum tensor and the conserved current in the absence of sources, whereas the parityodd susceptibilities require three-point functions. * ashish@uvic.ca In this paper we make use of the Kubo formulas of [1] to compute the seven parity-even equilibrium thermodynamic susceptibilities for a free massive Dirac field at zero temperature and non-zero chemical potential. Such a dense and degenerate matter configuration is of interest in astrophysical settings, such as white dwarfs and neutron stars [10][11][12][13], and in QCD [14][15][16]. We also compute the constitutive relations for the energy-momentum tensor and the conserved current, when the theory is not coupled to an external gauge field conjugate to the current.
The paper is organized as follows. Following [1], sections II A and II B review the formalism of equilibrium generating functional, the derivative expansion, and quote the Kubo formulae for the parity-even second order thermodynamic susceptibilities. Section II C then reviews the formalism for Dirac fields, mentioning in particular the energy-momentum tensor and the two-point function which we make use of in our computations. The thermodynamic susceptibilities are then computed in section III, which also includes a discussion of the conformal limit. We compute the second order constitutive relations and end with a discussion in section IV. Appendix A contains a discussion of the effects of non-zero temperature.

A. Thermal equilibrium and the generating functional
We consider a macroscopic system with a conserved current in equilibrium. The system can be coupled to an external nondynamical metric g µν , and an external non-dynamical gauge field A µ corresponding to the conserved current. Equilibrium implies the presence of a timelike Killing vector V µ . Coordinates in which V µ = (1, 0) correspond to the matter at rest. In the thermodynamic frame [9], the fluid four-velocity u µ , the temperature T , and the chemical potential µ can be defined via the Killing vector V µ as Here T 0 is the equilibrium temperature in the matter rest frame, and Λ V is a gauge function introduced to ensure the gauge independence of µ. The condition for the system to be in equilibrium translates to where L V denotes the Lie derivative with respect to V µ . As has been shown in [8,9], physical properties of the fluid in equilibrium, in particular the zero-frequency correlation functions over length scales much bigger than the length scales associated with the microscopic dynamics of the fluid, are encoded in a generating functional W. The generating functional is a functional of the external sources. For the system under consideration, the generating functional is a functional of the background metric g µν and the background gauge field A µ , which act as sources for the energy-momentum tensor and the conserved current, respectively. The generating functional can be expressed as where F [g µν , A µ ] is a local function of the sources. The variation of the generating functional gives with T µν , J µ respectively denoting the energy-momentum tensor and the conserved current. We assume that the system is free of any quantum anomalies, in which case the generating functional is both diffeomorphism as well as gauge invariant. This translates to where F µν = ∂ µ A ν − ∂ ν A µ is the field strength tensor for the background field A µ . In the hydrodynamic regime, when the external sources vary over length scales much larger than the microscopic length scales associated with the fluid such as the mean free path, the density F [g µν , A µ ] admits a derivative expansion. The expansion is in terms of the derivatives of the sources as well as the fluid variables defined in eq. (1). Thus, to compute the generating functional one is interested in finding out all the diffeomorphism and gauge invariant objects that can be constructed out of the sources and the fluid variables upto a given order in derivatives.
For instance, to the zeroth order in the derivatives we have only two invariants, T and µ. The generating functional then takes the form where p(µ, T ) is the equilibrium pressure of the system 1 , and the dots denote terms with higher derivatives. There are no invariants at the first order in derivatives in equilibrium, as the equilibrium conditions eq. (2) guarantee that the scalars u λ ∂ λ T, u λ ∂ λ µ, ∇.u all vanish.
At the second order in derivatives, interestingly, there are nine invariants [8]. In terms of these the generating functional can be represented as where s (2) n denotes the n th second order invariant, and f n (µ, T ) denotes the corresponding thermodynamic susceptibility [1]. In terms of the acceleration a µ = u ν ∇ ν u µ , the vorticity Ω µ = ǫ µναβ u ν ∇ α u β , 2 the electric field E µ = F µν u ν , and the magnetic field B µ = 1 2 ǫ µναβ u ν F αβ , the generating functional written in terms of the nine invariants at the second order is 3 where R denotes the Ricci scalar of the background geometry. The susceptibilities f n (µ, T ) are to be determined from the microscopic theory governing the system, just like the pressure. The seven susceptibilities f 1 , · · · f 7 are parity even, whereas f 8 , f 9 are parity odd. Our primary focus in the discussion that follows will be on the parity-even susceptibilities.

B. Constitutive relations and the Kubo formulae
In terms of the fluid four-velocity u µ , the energy-momentum tensor for the fluid can be decomposed as where ∆ µν = g µν + u µ u ν is a projector that projects orthogonal to the fluid four-velocity u µ . In the decomposition eq. (10) for T µν , the energy density E = u µ u ν T µν , the pressure eq. (7) take the form which have the form of a perfect fluid, with ǫ = −p+T ∂p/∂T +µ ∂p/∂µ being the equilibrium energy density, and n = ∂p/∂µ being the equilibrium charge density, provided p(µ, T ) is identified as the equilibrium pressure. 2 We have ǫ µναβ = ε µναβ / √ −g, with ε 0123 = 1. 3 We count the metric and the gauge field as O(1) quantities. Consequently the Ricci scalar is O(∂ 2 ), and the electric and magnetic fields are O(∂).
The derivative counting may differ depending upon the nature of the fluid. For instance, if the fluid is insulating and has no free charges, than the electric field inside it is not screened and can be strong, in which case it is to be counted as O(1) rather than O(∂). Similarly, in magnetohydrodynamics the magnetic field can be strong and has to be counted as O (1). See [17,18] for additional discussion. P = 1 3 ∆ µν T µν , the energy flux Q µ = −∆ µα T αβ u β , and τ µν = T µν , where the angle brackets on a pair of indices denote the symmetric transverse-traceless projection, Clearly, the energy flux is orthogonal to the fluid four-velocity, u µ Q µ = 0, and τ µν is a symmetric tensor both orthogonal to the fluid four-velocity as well as traceless, u µ τ µν = 0, g µν τ µν = 0.
In the hydrodynamic regime, the quantities E, P, Q µ and τ µν admit derivative expansions in terms of the derivatives of the sources and the fluid variables, called the constitutive relations. In the limit when the fluid is in equilibrium, we also have the generating functional eq. (3) at our disposal, whose variation with respect to the metric gives the energy-momentum tensor, eq. (4). Consider for instance a fluid with a conserved charge in equilibrium, with the associated chemical potential being non-zero but the background gauge field set to zero. Upto second order in derivatives the generating functional of such a fluid then has the form of eq. (9), with the background fields E, B set to zero. Its variation with respect to the metric gives the constitutive relations [1] where we have used the notation with a comma denoting a partial derivative with respect to the argument following. The constitutive relations eq. (11) allow for a straight forward computation of the equilibrium correlation functions of the energy-momentum tensor. For instance, the equilibrium two-point function is defined as the variation of the equilibrium one-point function with respect to the source. For matter at rest on a flat background, we have where we have decomposed the background metric as g µν = η µν + δg µν (k) e ik·x . Using eq. (12), the constitutive relations eq. (11) imply Kubo formulae for the second order susceptibilities in terms of zero-frequency two-point correlation functions of the energy-momentum tensor, given by [1] The Kubo formulas above are not the only ones we obtain by varying the energy-momentum tensor; for instance, f 1 can also be computed by using 4f 1 = lim k→0 ∂ 2 G T xx T yy /∂k 2 z . However, the formulas above are simple enough to be used for the computation of the susceptibilities as we shall see later. Note also that G T µν T αβ may include contact term contributions; see [19,20] for a discussion.
Similar to eq. (10), the U (1) current can be expressed as where the charge density N = −u µ J µ , and the spatial current J µ = ∆ µ ν J ν , with u µ J µ = 0. In equilibrium, in the spirit of eq. (12), the two-point functions of the current are defined in terms of the variation of the one-point functions, where the variation in the background field is δA λ (x) ≡ δA λ (k)e ik·x . From the second order generating functional eq. (9) one can compute the current by varying it with respect to A µ , eq. (4), and then vary the current with respect to the metric and the gauge field, eqs. (17), (18), giving Kubo formulae for the susceptibilities f 4 , · · · f 7 , which are (see [1] for details) The constitutive relations for a fluid with a conserved current J µ not coupled to the corresponding background field are where n = ∂p/∂µ is the charge density at zeroth-order.

C. Review of the Dirac formalism
The focus of the present article is on computing the parityeven second order thermodynamic susceptibilities f 1 , · · · f 7 for a free massive Dirac field at zero temperature and a nonzero chemical potential. In the present subsection we briefly review the formalism for Dirac fields and present the essential formulas we will need for computing the susceptibilities from two-point functions of the energy-momentum tensor and current, sec. II B. The action for the free massive Dirac field on a curved spacetime background is given by where Γ µ (x) are the spacetime dependent γ-matrices defined via Γ µ (x) = e µ a (x)γ a , with γ a being the spacetime independent γ-matrices and e µ a (x) being the vierbein fields. The Clifford algebras satisfied by the γ-matrices are The covariant derivative of the spinor field appearing in the action eq.(24) is given by with σ ab = 1 4 [γ a , γ b ], and the spin connection ω ab µ given by For computing the susceptibilities we will need the energymomentum tensor and the current following from the action eq.(24). The energy-momentum tensor for the Dirac field can be computed by using T µν = In flat space, the covariant derivatives in the above equations can be replaced by ordinary derivatives, and the spacetime dependent γ-matrices become spacetime independent. We will find it convenient to work in Euclidean signature rather then Lorentzian. The Euclidean Dirac action in flat space is given by whereγ µ are the Euclidean gamma matrices such thatγ 0 = γ 0 ,γ i = −iγ i . They satisfy the Clifford algebra The trace identities of interest are The Euclidean energy-momentum tensor has the form The conserved current for the global U (1) symmetry of the Dirac action eq. (24) is J µ =Ψγ µ Ψ, which when written in Euclidean signature implies J 0 Finally, the Euclidean space two-point function has the form (35) Notice that the above formulas were written for a Dirac field with a vanishing chemical potential. At a non-zero chemical potential µ for the global U (1) symmetry of the Dirac action, the Euclidean action becomes The effect of the non-zero chemical potential can easily be accommodated into our computations by shifting the zeroth component of the Euclidean four-vector k µ in the two-point function eq. (35) to k 0 → k 0 + iµ. As a warm-up exercise, we can compute the zeroth-order equilibrium energy density, pressure and charge density for the free massive Dirac field at T = 0, µ = 0 using the above results. The energy density is given by the one-point function ǫ = T 00 (ω = k = 0) ′ = − T 00 E (ω = k = 0) ′ . 4,5 This can be thought of as evaluating the one-loop diagram of figure 1. Using the expression for T µν E from eq. (34) and the two-point function eq.(35), one finds that ǫ = 0 for |µ| < m, whereas for |µ| > m we have Similarly, the equilibrium pressure can be computed by using 4 Components of the energy-momentum tensor in Lorentzian and Euclidean signatures are related via T 00 ↔ −T 00 E , T 0i ↔ −iT 0i E , T ij ↔ T ij E . 5 A ′ denotes the removal of the factor of (2π) 4 δ 4 (0) from the correlator.

III. THERMODYNAMIC SUSCEPTIBILITIES
With the aid of the Kubo formulas discussed in sec. II B we can now compute the equilibrium thermodynamic susceptibilities f 1 , · · · f 7 for the free massive Dirac field at zero temperature and non-zero chemical potential. The particle number density vanishes for |µ| < m at T = 0, as can be seen from the Fermi distribution. Thus for |µ| < m the system is essentially in its vacuum state. As will become clear from the computations below, all the susceptibilities vanish for |µ| < m.
The Kubo formulas indicate that the susceptibilities can be computed from one-loop diagrams of the form shown in figure  3. Consider for instance the computation of f 1 , given by the Kubo formula eq. (13). We need to first compute the quantity G T xy T xy , same as G T xy E T xy E , given in terms of the flat space zero-frequency two-point function via (40) From the expression for the energy-momentum tensor eq. (34) it is straight forward to compute T xy E (ω = 0, k), which turns out to be where p µ is a Euclidean four-vector. Making use of eq. (35) we can now compute G T xy E T xy E , which turns out to be where, without any loss of generality, we have assumed that k = (0, 0, k z ). Notice that we have shifted q 0 → q 0 + iµ to take into account the non-zero value of the chemical potential µ. To compute f 1 we now act G T xy E T xy E with the operator ∂ 2 /∂k 2 z and take the limit k z → 0. This leaves out the integration over the Euclidean four-vector q µ , where we have moved to spherical coordinates for the q integral, with q ≡ |q|. The θ, φ integrals are easy to perform and can be carried out first. The q 0 integral can then be performed using contour integration, noting that Thus, whenever the denominator of the integrand for the q 0 integral involves powers of (q 0 + iµ) 2 + q 2 + m 2 greater than 1, we can rewrite the integrand with appropriate number of derivatives 1 m ∂ ∂m extracted out. The left-over q 0 integrand thus only has simple poles at q 0 = −i(µ ± q 2 + m 2 ), and the integral can be performed easily. Both the poles for the q 0 integral lie on the imaginary axis. For µ > 0, to have a nonzero contribution we must have µ − q 2 + m 2 < 0. On the other hand, for µ < 0 we must have µ + q 2 + m 2 > 0 to get a non-vanishing result. Both the conditions essentially restrict the domain of integration for q to be µ 2 − m 2 ≤ q ≤ ∞, along with |µ| > m for a non-zero result.
We next perform the q integral, which requires a UV cutoff Λ to regulate the divergences coming from short distance physics. The cutoff dependence can be removed systematically using the standard renormalization procedure of zerotemperature zero-density quantum field theory. Here we are only interested in µ dependent terms, which appear independently from cutoff dependent terms. Keeping therefore only the µ dependent terms and acting upon the result of the q integration with appropriate factors of 1 m ∂ ∂m extracted out earlier, we get (42) for |µ| > m, and f 1 = 0 for |µ| < m. Similar computations using the Kubo formulae eqs. (14), (15) and eqs. (19) -(22) give the remaining second order susceptibilities for |µ| > m to be Like f 1 , the susceptibilities f 2 , · · · f 7 vanish for |µ| < m. Eqs.
(42) -(48) are some of the main results of this paper. Note that f 1 , · · · f 4 and f 6 are even under charge conjugation µ → −µ, whereas f 5 , f 7 are odd. An interesting limit to consider for the above results is the conformal limit, m → 0. 6 In the conformal limit, the generating functional eq. (9) has to be invariant under a Weyl rescaling g µν →g µν = e −2ξ g µν of the background metric to ensure the tracelessness of the energy-momentum tensor. A quantity χ that transforms under the above Weyl rescaling as χ →χ = e wξ χ is said to have a Weyl weight w. Clearly the The nodes correspond to stress tensor and current insertions, depending upon the susceptibility being computed. k denotes the external momentum, while q is the momentum running in the loop. measure in eq. (9) has w = −4, as √ −g → e −4ξ √ −g under the Weyl rescaling. Thus, for the Weyl invariance of the generating functional eq. (9), each of the terms f n s (2) n must have a Weyl weight of 4. The chemical potential µ and the temperature T , defined through eq. (1), each have w = 1. The parity even quantities Ω 2 and B · Ω have well defined Weyl weights of w = 2, 3 respectively. In the conformal limit, the Weyl weights of f 3 and f 5 from eqs. (44) and (46) are w = 2, 1 respectively, ensuring that the contribution of these two terms to the generating functional is Weyl invariant.
The parity even terms R, a 2 and E · a do not have well defined Weyl weights. However, the combination Thus, for the Weyl invariance of the generating functional we expect that in the conformal limit f 1 = Cµ 2 , f 2 = 6f 1 and f 7 = − 6 ∂f 1 /∂µ. From eqs. (42), (43) and (48) we see that in the conformal limit m → 0 we have f 1 = −µ 2 /48π 2 , f 2 = −µ 2 /8π 2 and f 7 = µ/4π 2 , thereby satisfying the criteria for Weyl invariance, with C = −1/48π 2 . Interestingly, when a conformal field theory is coupled to a background metric and/or a background gauge field, the quantum corrections destroy the tracelessness of the theory. The ensuing trace anomaly is given by [23,24] Here a, c are numbers that depend upon the field content of the theory; for the theory of free massless Dirac fermions we have a = 11/360, c = 1/20. However, as the terms multiplying a, c in the trace anomaly eq. (49) are fourth order in derivatives, they are not relevant for our present discussion as we are working only upto the second order. The last term in eq. (49) is however second order in derivatives, and must be encoded in the generating functional. Here b 0 is the coefficient of the one loop β-function of the coupling e used to minimally couple the theory to external electromagnetic fields, where M is the renormalization scale, and the action for the external gauge fields is − 1 7 For a Dirac fermion we Consider now the terms proportional to E 2 , B 2 in the generating functional, d 4 x √ −g f 4 B 2 + f 6 E 2 , which on variation with respect to the metric give the following contribution to the trace of the energy-momentum tensor, where f ′ n = µ ∂ µ f n + T ∂ T f n . From eq. (49) this should be equal to − 1 , which implies that we should have This is indeed borne out by our results in eqs. (45) and (47). Thus the susceptibilities are well behaved in the conformal limit, including the contribution to the trace anomaly.

IV. DISCUSSION
We have computed the seven parity even equilibrium thermodynamic susceptibilities appearing at the second order in the hydrodynamic derivative expansion for a free massive Dirac field at zero-temperature and non-zero chemical potential. The resulting susceptibilities, eqs. (42) -(48), can be inserted back into eqs. (11) and (23) to derive the constitutive relations for the Dirac field on a curved background. For the energy-momentum tensor the constitutive relations, with λ ≡ µ/m > 1, are 8 Here ǫ, p are zeroth-order energy density and pressure, eqs.
(37) and (38). The constitutive relations for the current are with n denoting the zeroth-order charge density, eq. (39). The coefficients of the two-derivative terms appearing in the con- 7 In the limit e → 0 the gauge field becomes non-dynamical. See section 2.1 of [25] for an interesting discussion. 8 We write the constitutive relations for µ > 0. The results for µ < 0 can be obtained similarly using the appropriate signs in eqs. (42) -(48). stitutive relations above are the transport coefficients of the theory. The susceptibilities f 2 , f 6 and f 7 diverge when the magnitude of the chemical potential approaches m from above, i.e. in the limit |µ| → m + . This behaviour is carried on to some of the transport coefficients appearing in the constitutive relations eqs. (52) and (53). As discussed in appendix A, the apparent discontinuity across |µ|/m = 1 is a consequence of working strictly at zero temperature. At non-zero temperatures the susceptibilities are continuous functions of µ/m.
Second order thermodynamic susceptibilities for a free Dirac field were also computed in [26]. However, the results are for a massless field at finite temperature, rather than the massive case we have considered. Ref. [27] computes the second order constitutive relations using equilibrium three-point functions, which are considerably harder to compute as opposed to the two-point functions we have used. As already mentioned, in equilibrium, the susceptibilities appearing in the generating functional are more fundamental than the transport coefficients in the constitutive relations, which are linear functions of the susceptibilities and their derivatives. Thus it is better to evaluate the susceptibilities first and then compute the transport coefficients from them, a route we have chosen. Besides this, the constitutive relations eqs. (52) and (53) exhibit explicit dependence on the background curvature, something not shown in [27].
It would be interesting to make use of the thermodynamic susceptibilities computed here in a physical situation where second order derivative corrections may become relevant. Ref. [28] for instance discusses the situation where a non-zero value for the susceptibilities can lead to a shift in the effective Newton's constant in the presence of matter.

Appendix A: Effects of non-zero temperature
In this paper, we have computed parity-even second order equilibrium thermodynamic susceptibilities for a massive Dirac field at zero temperature, eqs. (42)-(48). It seems from the results that the susceptibilities f 2 , f 6 and f 7 have a divergent discontinuity across |µ|/m = 1, figure 4. However, as we argue below, this discontinuity is an artifact of strictly imposing the T → 0 limit. At any non-zero value of the temperature the susceptibilities are continuous and well-behaved. To illustrate the point, let us consider the computation of f 7 at a non-zero temperature. The Kubo formula for f 7 is given by eq. (22). Using the expressions for the energy-momentum tensor and the current presented in section II C, the zero temperature expression for f 7 turns out to be whereq 0 is the shorthand for q 0 +iµ. When T = 0, the continuous integral over q 0 in the expression above has to be replaced by a discrete Matsubara sum, where β ≡ 1/T is the inverse temperature. Inserting the above in eq. (A1), performing the derivatives and integrating over the angular coordinates in q-space gives where q ≡ |q| as before, and z = iω n − µ.
The next step in the evaluation of f 7 at non-zero temperatures involves evaluating the Matsubara sum in eq. (A2). This can be done by converting the sum into a contour integral [29,30]. Consider the function β 2 tanh β(z+µ)

2
. This function has poles at z = iω n − µ, each with a unit residue, and is otherwise well behaved on the entire complex plane, including being bounded for the limit z → ∞. This allows us to write the infinite sum in eq. (A2) as a contour integral, where the contour C is shown in figure 5. This contour can be continuously deformed into the contour C 1 ∪ C 2 , figure 6. Now, due to the fact that the integrand in eq. (A3) falls off faster than 1/z for z → ∞, we can close C 1 with a semicircle C ′ 1 to its right and C 2 with a semicircle C ′ 2 to its left, without affecting the value of S. From the residue theorem it is then clear that the contributions to S will only come from the poles at z = ± q 2 + m 2 ≡ z ± , given by with the minus sign coming because we have closed the contours in a clockwise manner. 6. The contour C1 ∪ C2. The separation 2ǫ between C1 and C2 tends to zero. C ′ 1 and C ′ 2 are oriented semicircular arcs of ∞ radius which are used to close the contours C1 and C2, respectively.
The result in eq. (A4) above can be further simplified by noting that where F µ (z) = 1/ e β(z−µ) + 1 is the Fermi-Dirac distribution function. Inserting this in eq. (A4) gives The first term in the expression above vanishes, leaving behind the second term, which when substituted into eq. (A2) gives (A6) Although the residues in eq. (A6) above can be readily computed, the leftover integration on q is difficult to perform. The resulting integrand however is a well-behaved function of q, and for the sake of illustration, can be integrated numerically for different choices of the temperature T . The resulting behaviour of f 7 as a function of µ/m is plotted in figure 7 for µ > 0 and different choices of the temperature. As is evident from the figure, the smaller the temperature, the stronger the jump in f 7 around µ/m = 1, with the limiting case occurring at T = 0 as shown in figure 4.
The apparent discontinuities in f 2 and f 6 also smoothen out once one deviates from the strict zero temperature condition, just like the case for f 7 . One has to be careful while performing the numerical integration step to obtain the plots for f 2 and f 6 at a non-zero temperature due to UV divergences, which have to be carefully subtracted off by adding the appropriate counter terms. The results for f 2 and f 6 then are also continuous and well behaved. Plots depicting the behaviour of f7 as a function of µ/m for µ > 0 at different non-zero temperatures. As the temperature becomes smaller, the jump in f7 as well as the peak height around µ/m = 1 become more prominent, with the limiting case T = 0 plotted in figure 4. In an identical manner the non-zero temperature also removes the discontinuity at µ = −m.

ACKNOWLEDGMENTS
The author would like to thank Pavel Kovtun for insightful discussions, and for comments on a draft of this paper which lead to several improvements. The work is supported by the NSERC of Canada.