Raman response and shear viscosity in the non-Fermi liquid phase of Luttinger semimetals

Luttinger semimetals refer to a class of strongly spin–orbit coupled fermions with pseudospin-3/2 in three dimensions, whose normal state exhibits quadratic touching of Kramers-degenerate valence and conduction bands at the Brillouin zone center (i.e., the Γ-point). Due to the growing relevance of these strongly correlated compounds nowadays in the field of topological quantum materials, there has been a recent surge of theoretical works investigating various aspects of the physics of these systems [1–15]. Important examples of these compounds include HeTe [16–18], some pyrochlore iridates [19, 20], grey-tin [21], and half-Heusler compounds [22, 23].

Since the Luttinger semimetal harbors an isolated Fermi node at the Brillouin zone center, the electron–electron interactions are not effectively screened in these systems. As a consequence, a minimal model to describe these compounds must include strong (long-range) Coulomb interactions [1, 2]. It is by now well-established that the interactions mediated by the Coulomb forces in Luttinger semimetals stabilize a new non-Fermi liquid (NFL) state — the so-called Luttinger–Abrikosov–Beneslavskii (LAB) phase [24, 25]. The effective field theory for this phase was first studied by Abrikosov and Beneslavskii in the 1970s in a controlled approximation, by using a large-$N$ expansion [24, 25]. This NFL phase was later revisited and reformulated using dimensional regularization and renormalization group (RG) techniques by Moon et al [1], with many interesting new predictions. An important distinction of the LAB phase from other well-known NFLs arising for critical Fermi surfaces [26–46] is that the former represents an NFL fixed point at a Fermi node, rather than for a Fermi surface. From recent analytical works, we do have some other examples of nodal NFLs [47, 48] as well.

Recently, we have calculated the optical conductivity, the dc conductivity, the thermal conductivity, and the thermoelectric response of the LAB phase, using the Kubo formula and the memory matrix approaches [14, 15].Our theoretical results (when applicable) agree qualitatively with recent experimental data obtained for Luttinger semimetal compounds like the pyrochlore iridates [(Y_1−x Pr_x )₂Ir₂O₇] [49]. Therefore, as a further step to find other distinctive experimental consequences of the LAB phase, we calculate the Raman response in this paper. Nowadays, the Raman spectroscopy has become an indispensable experimental technique to study strongly correlated materials (for a comprehensive review, see reference [50]), since it is well-suited to yield valuable information of the non-trivial dynamics of the electronic excitations in these systems.

In view of the NFL nature of the LAB fixed point in d = 3 for the Luttinger model, it is physically reasonable to assume that the effects of electron–phonon interaction and the electron-impurity coupling are negligible in these systems. Therefore, an effective hydrodynamic regime is expected to emerge at low temperatures. One key transport property that characterizes such a regime is the shear viscosity, which measures the dissipative effect related to the internal friction in the electron fluid. In a breakthrough paper in 2005 [51], Kovtun–Son–Starinets proposed that the ratio of the shear viscosity (η) and the entropy density (s) for hydrodynamic systems has a universal lower bound of the form η/s ⩾ 1/(4π) (ℏ and k_B are set to unity here), using the anti-de Sitter/conformal field theory correspondence. If η/s turns out to be close to the lower bound of this inequality, then the system is classified as displaying 'minimal-viscosity', and can be viewed as a measure of a strongly interacting system. Important systems that fit into the 'minimal-viscosity' scenario include graphene at charge neutrality [52], quark–gluon plasma [53], and ultracold quantum gases tuned to the unitarity limit [54].

Not all strongly coupled quantum field theories conform to the lower bound of the Kovtun–Son–Starinets ratio though. Recently, Patel et al [55] has shown that a diverging η/s can appear at low temperatures, in the context of an NFL appearing at the Ising-nematic quantum critical point [35, 38, 39]. This behavior can be traced to the violation of the hyperscaling property of the corresponding field theory model [41]. Technically speaking, the hyperscaling is a property in which the entropy scales as if the system were defined effectively in (d − θ) dimensions, where θ is the so-called hyperscaling violation exponent. Indeed, there have been many studies of quantum critical points, which display varying degrees of hyperscaling violation [55–57]. There are essentially two ways through which this can happen: (1) due to a violation of the naive scaling of the conductivity σ(ω = 0, T) (where T denotes the temperature), which comes together with the naive scaling violation of s; (2) due to a mismatch between the scaling in ω of σ(ω, T = 0) and the T-dependence of σ(ω = 0, T). The second scenario usually appears when there is a dangerously irrelevant deformation in the corresponding low-energy effective field theory, breaking the charge-conjugation symmetry of the (otherwise charge-conjugation symmetric) infrared fixed point [58, 59].

Recently, we have shown that the LAB phase violates the hyperscaling property at low temperatures [14], by calculating the optical conductivity of the model. In this paper, we continue to investigate this NFL phase further by calculating the η/s ratio, using the memory matrix method. From our calculations, we find that this ratio scales in d = 3 as η/s ∼ T^λ−49/38, where 0 < λ < 1. Therefore, we find that η/s diverges in three dimensions as T → 0, instead of saturating to a finite value.

The paper is organized as follows. In section 2, we introduce the model for the Luttinger semimetal, augmented by long-ranged Coulomb interactions. We also describe the low-energy NFL fixed point characterizing the LAB phase, obtained from a dimensional regularization scheme in $d=4-\varepsilon$ spatial dimensions. In section 3, we compute the Raman response at T = 0 using the Kubo formalism, and also at T > 0 using the memory matrix approach. We then calculate the free energy, specific heat, and the entropy density in section 4. Finally, we derive the expressions for the optical (T = 0) and dc (T > 0) viscosities in section 5. This allows us to work out the η/s ratio in section 5.4. We end with some concluding remarks in section 6.

We consider the isotropic Luttinger Hamiltonian [60] captured by

$$\begin{equation}{\mathcal{H}}_{0}=\frac{\vert \mathbf{k}{\vert }^{2}}{2{m}^{\prime }}-\frac{\frac{5\vert \mathbf{k}{\vert }^{2}}{4}-{(\mathbf{k}\cdot \mathbf{J})}^{2}}{2m},\end{equation} \tag{ 2.1 }$$

where $\mathbf{J}$ is the vector angular momentum operator for j = 3/2 states. This effective Hamiltonian emerges from the three-dimensional band-structure of certain spin–orbit coupled systems, and harbors quadratic band crossings at the Brillouin zone center. The energy eigenvalues are

$$\begin{equation}{\varepsilon }_{\pm }(\mathbf{k})=\frac{\vert \mathbf{k}{\vert }^{2}}{2{m}^{\prime }}\pm \frac{\vert \mathbf{k}{\vert }^{2}}{2m},\end{equation} \tag{ 2.2 }$$

where the '+' and '−' signs refer to the conduction and valence bands, respectively. Each of these bands is doubly degenerate.

The model can also be expressed in an equivalent form [1, 2, 14]

$$\begin{equation}{\mathcal{H}}_{0}=\sum\limits _{a=1}^{5}{d}_{a}(\mathbf{k}){{\Gamma}}_{a}+\frac{\vert \mathbf{k}{\vert }^{2}}{2{m}^{\prime }},\quad {d}_{a}(\mathbf{k})=\frac{{\tilde{d}}_{a}\,(\mathbf{k})}{2m},\end{equation} \tag{ 2.3 }$$

using m and m' as the mass parameters that define the quadratic bands. The Γ_a matrices provide a rank-four irreducible representation of the Euclidean Clifford algebra, obeying the anticommutation relation {Γ_a , Γ_b } = 2δ_ab . Furthermore, ${\tilde{d}}_{a}(\mathbf{k})$'s represent the angular momentum l = 2 spherical harmonics [61], which have the explicit forms as follows:

$$\begin{align}\hfill & {\tilde{d}}_{1}(\mathbf{k})=\sqrt{3}{k}_{y}{k}_{z},\qquad {\tilde{d}}_{2}(\mathbf{k})=\sqrt{3}{k}_{x}{k}_{z},\qquad {\tilde{d}}_{3}(\mathbf{k})=\sqrt{3}{k}_{x}{k}_{y},\qquad \hfill \\ \hfill & {\tilde{d}}_{4}(\mathbf{k})=\frac{\sqrt{3}({k}_{x}^{2}-{k}_{y}^{2})}{2},\qquad {\tilde{d}}_{5}(\mathbf{k})=\frac{2{k}_{z}^{2}-{k}_{x}^{2}-{k}_{y}^{2}}{2}.\hfill \end{align} \tag{ 2.4 }$$

We will sometimes use the notation d_k instead of d(k) to avoid cluttering in long equations.

The term $\frac{\vert \mathbf{k}{\vert }^{2}}{2{m}^{\prime }}$ multiplies an identity matrix, causing the band masses of the conduction and valence bands to become unequal. For simplicity, we will consider the case with m' set to infinity (i.e. with equal band-masses) for all calculations ³ , except for the Raman response.

Adding the Coulomb interactions via a non-dynamical scalar boson field φ, the action of this system ⁴ can be straightforwardly written as

$$\begin{equation}{S}_{0}=\int \mathrm{d}\tau \,{\mathrm{d}}^{3}\mathbf{x}\left[\sum\limits _{i=1}^{{N}_{\mathrm{f}}}{\psi }_{i}^{{\dagger}}\left({\partial }_{\tau }+{\mathcal{H}}_{0}+\mathrm{i}e\,\varphi \right){\psi }_{i}+\frac{c}{2}{\left(\nabla \varphi \right)}^{2}\right].\end{equation} \tag{ 2.5 }$$

Here, ψ_i denotes the fermionic spinor with flavor index $i$, and c is a constant equal to 1/(4π).

The bare Green's function for each fermionic flavor is given by

$$\begin{equation}{G}_{0}({k}_{0},\mathbf{k})=\frac{\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }}{-{\left(\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}},\end{equation} \tag{ 2.6 }$$

where $\vert \mathbf{d}(\mathbf{k}){\vert }^{2}=\frac{{\mathbf{k}}^{4}}{4{m}^{2}}$.

We now integrate out φ to express the Coulomb interaction as an effective four-fermion interaction vertex, resulting in the final action

$$\begin{align}\hfill S& =\sum\limits _{i=1}^{{N}_{\mathrm{f}}}\int \frac{\mathrm{d}\tau \,{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}{\psi }_{i}^{{\dagger}}(\tau ,\mathbf{k})\left({\partial }_{\tau }+{\mathcal{H}}_{0}\right){\psi }_{i}(\tau ,\mathbf{k})\hfill \\ \hfill & \quad +\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{2c}\sum\limits _{i,j=1}^{{N}_{\mathrm{f}}}\int \frac{\mathrm{d}\tau \,{\mathrm{d}}^{d}\mathbf{k}\,{\mathrm{d}}^{d}{\mathbf{k}}^{\prime }\,{\mathrm{d}}^{d}\mathbf{q}}{{(2\pi )}^{3d}}V(\mathbf{q}){\psi }_{i}^{{\dagger}}(\tau ,\mathbf{k}+\mathbf{q})\hfill \\ \hfill & \quad \times {\psi }_{i}(\tau ,\mathbf{k}){\psi }_{j}^{{\dagger}}(\tau ,{\mathbf{k}}^{\prime }-\mathbf{q}){\psi }_{j}(\tau ,{\mathbf{k}}^{\prime }),\hfill \end{align} \tag{ 2.7 }$$

in d spatial dimensions, where where V(q) = 1/q². In the momentum space, the Coulomb interaction vertex is given by $\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{2c}V(\vert \mathbf{q}\vert )$, as shown in figure 1. We have also scaled $\frac{{e}^{2}}{c}$ by using the floating mass scale Λ (of the renormalization group flow) to make it dimensionless for d = 4 − ɛ spatial dimensions ⁵ .

Zoom In Zoom Out Reset image size

This action has an NFL phase, which can be accessed in a controlled way using dimensional regularization, by considering the theory in d spatial dimensions. The renormalization flow equation for the coupling constant e gives a stable interacting fixed point [1, 25], with the value

$$\begin{equation}{e}^{\ast }=\sqrt{\frac{60{\pi }^{2}c\varepsilon }{m\left.4+15{N}_{\mathrm{f}}\right.}}.\end{equation} \tag{ 2.8 }$$

This is an NFL with no long-lived quasiparticles at low energies, describing the LAB phase, with the dynamical critical exponent z at the fixed point given by ${z}^{\ast }=2-\frac{4\varepsilon }{4+15{N}_{\mathrm{f}}}$.

For performing the integrals without a finite UV cutoff, we will employ the scheme developed by Moon et al [1], where the radial momentum integrals are performed with respect to a d = 4 − ɛ dimensional measure $\int \frac{\vert \mathbf{k}{\vert }^{3-\varepsilon }\mathrm{d}\vert \mathbf{k}\vert }{{(2\pi )}^{4-\varepsilon }}$, but the Γ matrix structure of d = 3 is retained. Although the angular integrals are performed only over the three-dimensional sphere parameterized by the polar and azimuthal angles (θ, ϕ), the overall angular integral of an isotropic function ${\int }_{\hat{{\Omega}}}\cdot 1$ is taken to be 2π² (which is appropriate for the total solid angle in d = 4), and the angular integrals are normalized accordingly. Therefore, the angular integrations are performed with respect to the following measure:

$$\begin{equation}\int \mathrm{d}S(\dots )\equiv \frac{\pi }{2}{\int }_{0}^{\pi }\mathrm{d}\theta {\int }_{0}^{2\pi }\mathrm{d}\phi \,\mathrm{sin}\,\theta (\dots ),\end{equation} \tag{ 2.9 }$$

where the π/2 is inserted for the sake of normalization. To perform the full loop integrals, we will use the relations shown in appendix A. For evaluating various traces involving the gamma matrices, we will use the identity for the trace of the product of four gamma matrices, i.e.,

$$\begin{equation}\mathrm{Tr}\left[{{\Gamma}}_{a}{{\Gamma}}_{{a}^{\prime }}{{\Gamma}}_{b}{{\Gamma}}_{{b}^{\prime }}\right]=4\left({\delta }_{a{a}^{\prime }}{\delta }_{b{b}^{\prime }}-{\delta }_{ab}{\delta }_{{a}^{\prime }{b}^{\prime }}+{\delta }_{a{b}^{\prime }}{\delta }_{{a}^{\prime }b}\right),\end{equation} \tag{ 2.10 }$$

and also the standard identity for the trace of six gamma matrices [cf equation (A4)]. Lastly, in the rest of the paper, we will set N_f = 1.

In this section, we calculate the Raman response of the LAB phase. Raman response involves coupling the electrons to an electromagnetic field representing the incoming and outgoing photons. Hence, we introduce a gauge coupling in the model via the Peierls substitution k → k + e A/c, where A is the vector potential (this is also equivalent to introducing a scalar field ϕ, such that A = ∇ϕ). As a result, the vector potential dependent terms are given by:

$$\begin{equation}{H}_{\mathbf{A}}=\frac{e}{c}{A}_{\alpha }{J}_{\alpha }(\mathbf{k})+\frac{{e}^{2}}{2{c}^{2}}{A}_{\alpha }{A}_{\beta }{\gamma }_{\alpha \beta }(\mathbf{k}),\end{equation} \tag{ 3.1 }$$

where ${J}_{\alpha }=\frac{{k}_{\alpha }}{{m}^{{\prime}}}+{\partial }_{{k}_{\alpha }}{d}_{a}(\mathbf{k}){{\Gamma}}_{a}$ represents the components of the current operator, and

$$\begin{align}\hfill & {\gamma }_{xx}=\frac{{{\Gamma}}_{0}}{{m}^{\prime }}+\frac{\sqrt{3}{{\Gamma}}_{4}-{{\Gamma}}_{5}}{2m},\qquad {\gamma }_{yy}=\frac{{{\Gamma}}_{0}}{{m}^{\prime }}-\frac{\sqrt{3}{{\Gamma}}_{4}+{{\Gamma}}_{5}}{2m},\qquad \hfill \\ \hfill & {\gamma }_{zz}=\frac{{{\Gamma}}_{0}}{{m}^{\prime }}+\frac{{{\Gamma}}_{5}}{m},\hfill \\ \hfill & {\gamma }_{xy}={\gamma }_{yx}=\frac{\sqrt{3}{{\Gamma}}_{3}}{m},\qquad {\gamma }_{yz}={\gamma }_{zy}=\frac{\sqrt{3}{{\Gamma}}_{1}}{m},\qquad \hfill \\ \hfill & {\gamma }_{xz}={\gamma }_{zx}=\frac{\sqrt{3}{{\Gamma}}_{2}}{m}.\hfill \end{align} \tag{ 3.2 }$$

There are two types of electron–photon couplings, as illustrated in figure 2. We note that in evaluating the loop diagrams for Raman response, the two vertices can only be either (Γ₀, Γ₀), or (Γ_a , Γ_b ) for a, b = 1, ..., 5, as the cross-terms vanish identically.

Zoom In Zoom Out Reset image size

The vector potential can be quantized as:

$$\begin{equation}\mathbf{A}(\mathbf{q})=\sqrt{\frac{h{c}^{2}}{{\omega }_{\mathbf{q}}\mathcal{V}}}\left({\hat{\mathbf{e}}}_{\mathbf{q}}{a}_{-\mathbf{q}}+{\hat{\mathbf{e}}}_{\mathbf{q}}^{\ast }{a}_{\mathbf{q}}^{{\dagger}}\right),\end{equation} \tag{ 3.3 }$$

where $\mathcal{V}$ the volume, and ${a}_{\mathbf{q}}^{{\dagger}}$ (a_q ) is the creation (annihilation) operator of the photons with energy ℏω_q = ℏc|q| having a polarization direction denoted by ${\hat{\mathbf{e}}}_{\mathbf{q}}$. Within the Born approximation, the Raman scattering cross-section is given by [50]:

$$\begin{equation}\frac{{\partial }^{2}\sigma }{\partial {\Omega}\partial {\omega }_{i}}\propto \sum\limits _{F,I}\frac{\mathrm{exp}(-\beta {E}_{I})}{Z}\vert {M}_{FI}{\vert }^{2}\delta ({E}_{F}+{\omega }_{s}-{E}_{I}-{\omega }_{i}),\end{equation} \tag{ 3.4 }$$

where (I, F) denote initial and final states of the Luttinger semimetal, Z is the partition function, and M_FI = ⟨F|M|I⟩ (with M being the effective light-scattering operator). The summation represents a thermodynamic average over all possible initial and final states of the system having energies E_I and E_F , respectively, with the momentum vectors in the solid angle element dΩ. Furthermore, ω = ω_i − ω_s is the net frequency, and q is the net momentum transferred by the photons.

Denoting the photon polarization vectors for incident and scattered light by ${\hat{\mathbf{e}}}_{i}$ and ${\hat{\mathbf{e}}}_{s}$, respectively, the Raman scattering amplitude is given by

$$\begin{equation}{M}_{FI}(t)={\left\langle {\hat{\mathbf{e}}}_{s},F\right\vert }_{t}U(t,0){\left\vert {\hat{\mathbf{e}}}_{i},I\right\rangle }_{0},\end{equation} \tag{ 3.5 }$$

where in the interaction picture, we have

$$\begin{equation}U(t,0)=\mathcal{T}\,{\mathrm{e}}^{-\mathrm{i}{\int }_{0}^{t}\mathrm{d}t\,{H}_{\mathbf{A}}(t)},\qquad {H}_{\mathbf{A}}(t)={\mathrm{e}}^{\text{i}\,{\mathcal{H}}_{0}t}{H}_{\mathbf{A}}\,{\mathrm{e}}^{-\mathrm{i}\,{\mathcal{H}}_{0}t}.\end{equation} \tag{ 3.6 }$$

Let us define the operators

$$\begin{equation}{\rho }_{0}={\psi }^{{\dagger}}\psi ,\qquad {\rho }_{a}={\psi }^{{\dagger}}{{\Gamma}}_{a}\psi ,\end{equation} \tag{ 3.7 }$$

whose two-point correlators will contribute to the Raman response |M_FI |² (see below). Considering Raman scattering in the visible range, the zero momentum limit for the response is a good approximation [50]. Hence, we will compute the correlators $\left\langle {\rho }_{0}{\rho }_{0}\right\rangle (\mathrm{i}\omega )$ and $\left\langle {\rho }_{a}{\rho }_{b}\right\rangle (\mathrm{i}\omega )$ up to two-loop order.

Feynman diagrams for computing |M_FI |² involve vertices of two types, as depicted in figure 2. However, only diagrams consisting solely of green vertices involve non-resonant scatterings, while the others give rise to resonant and mixed scatterings which can be neglected in the low-energy limit [50]. Hence, the leading-order Feynman diagrams for |M_FI |² in the non-resonant scattering limit are only considered here.

3.1. Raman response at T = 0

We now proceed with a perturbative calculation of the Raman response at T = 0.

3.1.1. One-loop contribution

At one-loop order, from figure 3, we obtain

$$\begin{align}{\left\langle {\rho }_{0}{\rho }_{0}\right\rangle }_{1\;\text{loop}}(\mathrm{i}\omega )& =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \,\,\frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\,\mathrm{Tr}\left[{G}_{0}(k+q){G}_{0}(k)\right]\\\ & =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \,\,\frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\,\mathrm{Tr}\left[\frac{\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }}{-{\left.\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}\frac{\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }}{-{\left.\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}\right]\\ & =-4\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \,\,\frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\frac{\left.\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.\left.\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}{\left[-{\left.\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]\left[-{\left.\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{{\prime}}}\right.}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]}=0\quad \text{for}\enspace {m}^{\prime } > m,\end{align} \tag{ 3.8 }$$

where we define henceforth q = (ω, 0).

Zoom In Zoom Out Reset image size

As for a, b ≠ 0, we have:

$$\begin{align}\hfill & {\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{1\;\text{loop}}(\mathrm{i}\omega )=-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \enspace \enspace \frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\,\mathrm{Tr}\left[{{\Gamma}}_{a}{G}_{0}(k+q){{\Gamma}}_{b}{G}_{0}(k)\right]\hfill \\ \hfill & \qquad =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \enspace \enspace \frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\,\mathrm{Tr}\left[{{\Gamma}}_{a}\frac{\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }}{-{\left(\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}{{\Gamma}}_{b}\frac{\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }}{-{\left(\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}\right]\hfill \\ \hfill & \qquad =-4\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \enspace \enspace \frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{d}}\frac{{\delta }_{ab}\left[\left(\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)\left(\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)-\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]+2{d}_{a}(\mathbf{k}){d}_{b}(\mathbf{k})}{\left[-{\left(\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]\left[-{\left(\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }}\right)}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]}\hfill \\ \hfill & \qquad =-\frac{{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}{\delta }_{ab}}{10\pi }.\hfill \end{align} \tag{ 3.9 }$$

Consequently, at zeroth order, the Raman response is proportional to ${\omega }^{1-\frac{\varepsilon }{2}}$ for ω ≫ T.

3.1.2. Two-loop contributions

The first two-loop-order correction involves inserting the one-loop fermion self-energy (Σ₁) corrections into the correlator [cf figures 4(a) and (b)], which takes the form:

$$\begin{align}\hfill {\left\langle {\rho }_{a}\,{\rho }_{b}\right\rangle }_{2\text{loop}}^{(1)}(\mathrm{i}\,\omega )& =-\int \frac{d{k}_{0}}{2\,\pi }\int \frac{{d}^{d}\mathbf{k}}{{(2\,\pi )}^{d}}\mathrm{Tr}\left[{{\Gamma}}_{a}\,{G}_{0}(k+q)\,{{\Sigma}}_{1}(\mathbf{k})\right.\hfill \\ \hfill & \quad \left.\times \,{G}_{0}(k+q)\,{{\Gamma}}_{b}\,\,{G}_{0}(k)\right]-\int \frac{d{k}_{0}}{2\,\pi }\int \frac{{d}^{d}\mathbf{k}}{{(2\,\pi )}^{d}}\hfill \\ \hfill & \quad \times \mathrm{Tr}\left[{{\Gamma}}_{a}\,{G}_{0}(k+q)\,{{\Gamma}}_{b}\,{G}_{0}(k)\,{{\Sigma}}_{1}(\mathbf{k})\,{G}_{0}(k)\right],\hfill \end{align} \tag{ 3.10 }$$

where ${{\Sigma}}_{1}(\boldsymbol{\ell })=-\frac{m{e}^{2}}{15{\pi }^{2}c}{\left(\frac{{\Lambda}}{\vert \boldsymbol{\ell }\vert }\right)}^{\varepsilon }\frac{\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }}{\varepsilon }$ (see reference [2]). Evaluating the integrals for a, b > 0 (for more details, see appendix B), we obtain:

$$\begin{equation}{\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )=\frac{{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{75{\pi }^{3}c}\left[\frac{1}{\varepsilon }-\frac{1}{2}\,\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)\right]{\delta }_{ab}.\end{equation} \tag{ 3.11 }$$

Zoom In Zoom Out Reset image size

The second two-loop-order correction involves the vertex correction [cf figure 4(c)], which can be expressed as:

$$\begin{align}\hfill & {\left\langle {\rho }_{a}\,{\rho }_{b}\right\rangle }_{2\text{loop}}^{(2)}(\mathrm{i}\,\omega )\hfill \\ \hfill & \qquad =\frac{{e}^{2}\,{{\Lambda}}^{\varepsilon }}{c}\int \frac{d{k}_{0}\,d{\ell }_{0}}{{(2\,\pi )}^{2}}\int \frac{{d}^{d}\mathbf{k}\,{d}^{d}\boldsymbol{\ell }}{{(2\,\pi )}^{2d}}\mathrm{Tr}\left[{{\Gamma}}_{a}\,{G}_{0}(k+q)\right.\hfill \\ \hfill & \qquad \quad \times \left.\,\frac{1}{{\boldsymbol{\ell }}^{2}}\,{G}_{0}(k+q+\ell )\,{{\Gamma}}_{b}{G}_{0}(k+\ell )\,{G}_{0}(k)\right]\hfill \\ \hfill & \qquad =\frac{{e}^{2}\,{{\Lambda}}^{\varepsilon }}{c}\int \frac{d{k}_{0}\,d{\ell }_{0}}{{(2\,\pi )}^{2}}\int \frac{{d}^{d}\mathbf{k}\,{d}^{d}\boldsymbol{\ell }}{{(2\,\pi )}^{2d}}\mathrm{Tr}\left[{{\Gamma}}_{a}\,{G}_{0}({k}_{0}+\omega ,\mathbf{k})\right.\hfill \\ \hfill & \qquad \quad \times \left.\,\frac{1}{{\boldsymbol{\ell }}^{2}}\,{G}_{0}({\ell }_{0}+\omega ,\mathbf{k}+\boldsymbol{\ell })\,{{\Gamma}}_{b}\,{G}_{0}({\ell }_{0},\mathbf{k}+\boldsymbol{\ell })\,{G}_{0}({k}_{0},\mathbf{k})\right],\hfill \end{align} \tag{ 3.12 }$$

with q = (ω, 0) and ℓ = (ℓ₀, ℓ ). Performing the above integrals for a, b > 0 (for more details, see appendix B), we finally get:

$$\begin{equation}{\langle {\rho }_{a}{\rho }_{b}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )=-\frac{7{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{600{\pi }^{3}c}\left[\frac{1}{\varepsilon }-\frac{1}{2}\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)\right]{\delta }_{ab}.\end{equation} \tag{ 3.13 }$$

3.1.3. Scaling of the Raman response at T = 0

In order to obtain any renormalized physical observable in the effective field theory formalism, we need to use the fact that $\frac{1}{{\varepsilon }^{2}}$ terms are canceled out by the corresponding counterterms of the renormalized action [62]. We also use the value $m{{e}^{\ast }}^{2}/\left({\pi }^{2}c\right)=60\varepsilon /19$ at the NFL fixed point. Putting together all the terms, the final expression for ⟨ρ_a ρ_b ⟩(iω) up to two-loop order takes the form:

$$\begin{align}\hfill & \langle {\rho }_{a}\,{\rho }_{b}\rangle (\mathrm{i}\,\omega )\hfill \\ \hfill & \qquad ={\langle {\rho }_{a}\,{\rho }_{b}\rangle }_{1\text{loop}}(\mathrm{i}\,\omega )+{\langle {\rho }_{a}\,{\rho }_{b}\rangle }_{2\text{loop}}^{(1)}(\mathrm{i}\,\omega )+{\langle {\rho }_{a}\,{\rho }_{b}\rangle }_{2\text{loop}}^{(2)}(\mathrm{i}\,\omega )\hfill \\ \hfill & \qquad \quad +{\langle {\rho }_{a}\,{\rho }_{b}\rangle }_{\text{counterterms}}(\mathrm{i}\,\omega )\hfill \\ \hfill & \qquad =-\frac{{m}^{2-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}\,{\delta }_{ab}}{10\,\pi }-\frac{{{e}^{\ast }}^{2}\,{m}^{3-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{150\,{\pi }^{3}\,c}\mathrm{ln}\hfill \\ \hfill & \qquad \quad \times \left(\frac{m\,\vert \omega \vert }{{{\Lambda}}^{2}}\right){\delta }_{ab}+\frac{7\,{{e}^{\ast }}^{2}\,{m}^{3-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{1200\,{\pi }^{3}\,c}\mathrm{ln}\left(\frac{m\,\vert \omega \vert }{{{\Lambda}}^{2}}\right){\delta }_{ab}\hfill \\ \hfill & \qquad =-\frac{{m}^{2-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}\,{\delta }_{ab}}{10\,\pi }-\frac{\varepsilon \,{m}^{2-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}\mathrm{ln}\left(\frac{m\,\vert \omega \vert }{{{\Lambda}}^{2}}\right)\,{\delta }_{ab}}{380\,\pi }\hfill \\ \hfill & \qquad =-\frac{{m}^{2-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}\,{\delta }_{ab}}{10\,\pi }\left[1+\frac{\varepsilon \,\mathrm{ln}\left(\frac{m\,\vert \omega \vert }{{{\Lambda}}^{2}}\right)}{38}\right]\hfill \\ \hfill & \qquad \simeq -\frac{{m}^{2-\frac{\varepsilon }{2}}\,\vert \omega {\vert }^{1-\frac{\varepsilon }{2}+\frac{\varepsilon }{38}}\,{\delta }_{ab}}{10\,\pi }{\left(\frac{m}{{{\Lambda}}^{2}}\right)}^{\frac{\varepsilon }{38}},\hfill \end{align} \tag{ 3.14 }$$

after re-exponentiating the correction term coming from the two-loop diagrams. Therefore, the corrected Raman response scales as $\vert \omega {\vert }^{1-\varepsilon /2+\varepsilon /38}\;\stackrel{\varepsilon \,=\,1}{\,=\,}\;\vert \omega {\vert }^{1/2+1/38}$, after including the leading-order corrections.

3.2. Raman response at T > 0

In this subsection, we compute the Raman response at finite temperatures, using the memory matrix approach. This approach is appropriate to compute finite temperature responses in systems where quasiparticles do not exist, and has been explained in great detail in our previous works (e.g., references [14, 15, 63–65]). Just as in references [14, 15], we consider here the coupling to weak short-ranged disorder as the main mechanism for momentum relaxation in the LAB phase. Hence, we add the impurity term ${H}_{\mathrm{i}\mathrm{m}\mathrm{p}}={\sum }_{i}\int {\mathrm{d}}^{3}\mathbf{x}\mathbf{}W(\mathbf{x}){\psi }_{i}^{{\dagger}}(\mathbf{x}){\psi }_{i}(\mathbf{x})$ to the system. We assume the disorder to have a Gaussian distribution, such that $\left\langle \left\langle W(\mathbf{x})\right\rangle \right\rangle =0$ and $\left\langle \left\langle W(\mathbf{x})W({\mathbf{x}}^{\prime })\right\rangle \right\rangle ={W}_{0}^{2}\delta (\mathbf{x}-{\mathbf{x}}^{\prime })$, where ${W}_{0}^{2}$ represents the average of the magnitude-squared of the random potential experienced by the fermions. We point out that, unlike the T = 0 calculation (where we performed a systematic ɛ-expansion), we will now work directly in d = 3 in order to minimize the technical complexity involved. In this formalism, we will then need to implement a hard ultraviolet (UV) cutoff Λ for the momentum integrals, rather than using the dimensional regularization scheme (where the integrals are evaluated by sending the upper limit to ∞).

The T > 0 Raman response can be conventionally defined as

$$\begin{equation}{D}_{\text{Raman}}(\omega ,T)={\int }_{0}^{\infty }\mathrm{d}t\,{\text{e}}^{\text{i}\omega t}\left\langle \left[{\rho }_{a}(t),{\rho }_{b}(0)\right]\right\rangle .\end{equation} \tag{ 3.15 }$$

This correlation function can be approximately written within the memory matrix formalism as (cf references [66, 67])

$$\begin{equation}{D}_{\text{Raman}}(\omega ,T)\approx {\chi }_{{\rho }_{a}{\rho }_{b}}\left(\frac{\mathrm{i}\omega }{{M}_{{\rho }_{b}{\rho }_{b}}(\omega )-\mathrm{i}\omega {\chi }_{{\rho }_{b}{\rho }_{b}}}\right){\chi }_{{\rho }_{b}{\rho }_{a}},\end{equation} \tag{ 3.16 }$$

${\chi }_{{\rho }_{a}{\rho }_{b}}$ is the static thermodynamic susceptibility given by ${\chi }_{{\rho }_{a}{\rho }_{b}}(T)={\int }_{0}^{\beta }\mathrm{d}\tau \langle {\rho }_{a}(\tau ){\rho }_{b}(0)\rangle$, and ${M}_{{\rho }_{a}{\rho }_{b}}(\omega )$ is the memory matrix. The latter is defined by

$$\begin{equation}{M}_{{\rho }_{a}{\rho }_{b}}(\omega ,T)={\int }_{0}^{\beta }\mathrm{d}\tau \left\langle {\dot{\rho }}_{a}(0)\mathcal{Q}\frac{\mathrm{i}}{\omega -\mathcal{Q}\text{L}\mathcal{Q}}\mathcal{Q}{\dot{\rho }}_{b}(\mathrm{i}\tau )\right\rangle ,\end{equation} \tag{ 3.17 }$$

where L is the Liouville operator, and $\mathcal{Q}$ is an operator that projects out of a subspace of nearly-conserved operators.

We note that an important step in the memory matrix formalism is to identify a set of nearly-conserved operators in the system and, subsequently, project the dynamics onto these operators. As explained in references [14, 15], the momentum operator P is a nearly-conserved operator in the system under consideration, that plays a key role in the computation of many transport coefficients in the system within the hydrodynamic regime. However, the momentum operator turns out to have no influence on the Raman response of the system, since it has no overlap with the operator ρ_a . Fortunately, we observe that ρ_a is also a nearly-conserved operator ⁶ , and hence we can use it in place of P. Consequently, to leading order, the memory matrix is given by ${M}_{{\rho }_{a}{\rho }_{b}}(\omega ,T)=\text{Im}\,{G}_{{\dot{\rho }}_{a}{\dot{\rho }}_{b}}^{R}(\omega ,T)/\omega$, where ${G}_{{\dot{\rho }}_{a}{\dot{\rho }}_{b}}^{R}(\omega ,T)={\langle {\dot{\rho }}_{a}(\omega ){\dot{\rho }}_{b}(-\omega )\rangle }_{0}$ is the retarded correlation function for ${\dot{\rho }}_{a}=\mathrm{i}\left[H+{H}_{\mathrm{i}\mathrm{m}\mathrm{p}},{\rho }_{a}\right]$. Here, H represents the total Hamiltonian, and the notation ⟨⋯⟩₀ implies that after the addition of the weak disorder term, the average in the grand-canonical ensemble is taken by employing the non-interacting Hamiltonian of the system, to leading order. This gives the memory matrix as

$$\begin{equation}{M}_{{\rho }_{a}{\rho }_{b}}(\omega ,T)={W}_{0}^{2}\int \frac{{\mathrm{d}}^{3}\mathbf{q}}{{(2\pi )}^{3}}\frac{\text{Im}\,{{\Pi}}^{R}(\mathbf{q},\omega )}{\omega },\end{equation} \tag{ 3.19 }$$

where Π^R (q, ω) = Π(q, iω → ω + i0⁺), and ${\Pi}(\mathbf{q},\mathrm{i}\omega )=-T{N}_{\mathrm{f}}{\sum }_{{k}_{0}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}}{{(2\pi )}^{3}}\,\mathrm{Tr}[{{\Gamma}}_{a}{G}_{0}(k+q){{\Gamma}}_{b}{G}_{0}(k)]$. Using equation (3.16) and ${\chi }_{{\rho }_{a}{\rho }_{b}}\propto {\delta }_{ab}$, we finally obtain

$$\begin{equation}\text{Im}\,{D}_{\text{Raman}}(\omega ,T)=\frac{\omega {M}_{{\rho }_{a}{\rho }_{a}}{\chi }_{{\rho }_{a}{\rho }_{a}}^{2}}{{\omega }^{2}{\chi }_{{\rho }_{a}{\rho }_{a}}^{2}+{M}_{{\rho }_{a}{\rho }_{a}}^{2}}\equiv \frac{\omega \tilde{{\Gamma}}{\chi }_{{\rho }_{a}{\rho }_{a}}}{{\omega }^{2}+{\tilde{{\Gamma}}}^{2}},\end{equation} \tag{ 3.20 }$$

where $\tilde{{\Gamma}}={M}_{{\rho }_{a}{\rho }_{a}}{\chi }_{{\rho }_{a}{\rho }_{a}}^{-1}$. Since analytical expressions of ${M}_{{\rho }_{a}{\rho }_{a}}(T)$ and ${\chi }_{{\rho }_{a}{\rho }_{a}}(T)$ at finite temperatures cannot be obtained in closed-forms, we find their values using standard numerical integration techniques. We then fit the corresponding data in order to extract the temperature dependence of these quantities. This allows us to obtain Im D_Raman(ω, T), and a representative plot is shown in figure 5. From this figure, we observe that the Raman response at T > 0 displays a quasi-elastic peak (QEP) at ${\omega }_{\mathrm{max}}(T)\approx \text{Re}\,{M}_{{\rho }_{a}{\rho }_{a}}(\omega =0,T)/{\chi }_{{\rho }_{a}{\rho }_{a}}(T)$, with the corresponding peak height given by ${D}_{\text{Raman}}({\omega }_{\mathrm{m}\mathrm{a}\mathrm{x}})={\chi }_{{\rho }_{a}{\rho }_{a}}(T)/2$. To leading order, the static thermodynamic susceptibility is given by ${\chi }_{{\rho }_{a}{\rho }_{a}}(T)\sim a+b{T}^{1/2}+c/T$ (where a, b and c are temperature independent constants). Furthermore, our numerical results show that ${M}_{{\rho }_{a}{\rho }_{a}}(\omega =0,T)$ is either a constant, or has an extremely weak T-dependence (i.e., not clearly observed within the numerical accuracy).

Zoom In Zoom Out Reset image size

In this section, we determine the scaling behavior of the free energy at a finite temperature T > 0, first from general scaling arguments, and subsequently using the fixed point theory for the LAB phase, obtained using dimensional regularization [1].

4.1. General arguments

We derive the scaling forms for the case when hyperscaling property is obeyed. In a system with spatial dimension d, dynamical critical exponent z, the free energy F has the scaling dimension

$$\begin{equation}[F]=d+z,\end{equation} \tag{ 4.1 }$$

and we expect the temperature-dependence to be

$$\begin{equation}F(T)\sim {T}^{d/z+1}.\end{equation} \tag{ 4.2 }$$

From this scaling form, we conclude that the specific heat should scale as C(T) ∼ T^d/z . Since the current operator is given by ${J}_{\alpha }(T)=\frac{\delta F}{\delta {A}_{\alpha }}$, where ${A}_{\alpha }$ is the vector potential with scaling dimension equal to unity, the Kubo formula gives the scaling dimension of the optical conductivity σ(ω) as

$$\begin{align}\hfill [\sigma ]& =2[J]-z-[\text{Volume}\;\text{in}\;\text{energy}-\text{momentum}\;\text{space}]\hfill \\ \hfill & =2\left(d+z-1\right)-z-d-z=d-2.\hfill \end{align} \tag{ 4.3 }$$

This implies the scaling form

$$\begin{equation}\sigma (\omega )\sim {\omega }^{(d-2)/z}.\end{equation} \tag{ 4.4 }$$

We showed in reference [14] that this hyperscaling relation for σ is violated, when we include the effects of Coulomb interactions at the LAB fixed point.

4.2. Scaling of the free energy

In this subsection, we compute the free energy density at a finite temperature T > 0 to leading order ɛ. The free energy density receives contributions from two parts, namely, the free fermions and the corrections due to the Coulomb interactions. Since the bosons φ have no dynamics, there is no contribution from any free bosons, unlike other cases studied earlier by one of the authors [41, 43].

The contribution of the free fermions is given by ΔF⁽⁰⁾(T) = F⁽⁰⁾(T) − F⁽⁰⁾(0), such that

$$\begin{equation}{\Delta}{F}^{(0)}(T)=-2\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{4-\varepsilon }}\left[T\sum\limits _{p=\pm }\mathrm{ln}\left(1+{e}^{-\frac{p\vert \mathbf{k}{\vert }^{2}}{2mT}}\right)-\frac{\vert \mathbf{k}{\vert }^{2}}{2m}\right],\end{equation} \tag{ 4.5 }$$

where we have subtracted the infinite contribution from the temperature-independent ground state energy, in order to make the expression finite. We have also included a factor of two to account for the doubly-degenerate bands. Rescaling $\mathbf{k}\to \mathbf{k}\sqrt{2mT}$, we get

$$\begin{align}\hfill {\Delta}{F}^{(0)}(T)& =-\frac{2{\pi }^{\frac{4-\varepsilon }{2}}{\left(2mT\right)}^{3-\frac{\varepsilon }{2}}}{m{\Gamma}\left(\frac{4-\varepsilon }{2}\right)}{\int }_{0}^{\infty }\frac{\mathrm{d}\vert \mathbf{k}\vert \vert \mathbf{k}{\vert }^{3-\varepsilon }}{{(2\pi )}^{4-\varepsilon }}\hfill \\ \hfill & \quad \times \left[\sum\limits _{p=\pm }\mathrm{ln}\left(1+{e}^{-p\vert \mathbf{k}{\vert }^{2}}\right)-\vert \mathbf{k}{\vert }^{2}\right]\hfill \\ \hfill & =-\frac{{\pi }^{\frac{4-\varepsilon }{2}}{\left(2mT\right)}^{3-\frac{\varepsilon }{2}}}{m{\Gamma}\left(\frac{4-\varepsilon }{2}\right)}{\int }_{0}^{\infty }\frac{\mathrm{d}u\,{u}^{1-\varepsilon /2}}{{(2\pi )}^{4-\varepsilon }}\hfill \\ \hfill & \quad \times \left[\mathrm{ln}\left(\frac{\left(1+{e}^{u}\right)\left(1+{e}^{-u}\right)}{{e}^{u}}\right)\right]\hfill \\ \hfill & \simeq -\frac{\eta (3){\left(mT\right)}^{3-\frac{\varepsilon }{2}}}{{\pi }^{2}m}=-\frac{3\zeta (3){\left(mT\right)}^{3-\frac{\varepsilon }{2}}}{4{\pi }^{2}m},\hfill \end{align} \tag{ 4.6 }$$

where $\eta (u)=\frac{{\int }_{0}^{\infty }\mathrm{d}t\,{t}^{u-2}\,\mathrm{ln}(1+{e}^{-u})}{{\Gamma}(u-1)}$ is the Dirichlet eta function, and we have set ɛ = 0 in the numerical prefactor. We have also used the identity $\eta (3)=\frac{3\zeta (3)}{4}$, where ζ(u) is the Riemann zeta function.

The first order interaction correction to the free energy, as shown in figure 6, is given by

$$\begin{align}\hfill {F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}& =\frac{{e}^{2}{{\Lambda}}^{\varepsilon }{T}^{2}}{c}\sum\limits _{{{\Omega}}_{p},{\omega }_{{p}^{\prime }}}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{q}\,{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{8-2\varepsilon }}\hfill \\ \hfill & \quad \times \frac{\mathrm{Tr}\left[G({\omega }_{{p}^{\prime }}+{{\Omega}}_{p},\mathbf{k}+\mathbf{q})G({\omega }_{{p}^{\prime }},\mathbf{k})\right]}{{\mathbf{q}}^{2}},\hfill \end{align} \tag{ 4.7 }$$

where Ω_p and ω_p' are bosonic and fermionic Matsubara frequencies, respectively. From this expression, we need to isolate the pole (in ɛ) contributions. To lowest order in ɛ, these are obtained by evaluating one frequency sum as an integral in the limit T → 0, and the other one at finite temperature (T > 0) (cf reference [68]). To evaluate the frequency summations, we use the following zeta-function regularization identities:

$$\begin{equation}T\sum\limits _{{{\Omega}}_{p}}\frac{1}{\vert {{\Omega}}_{p}{\vert }^{s}}=\frac{2{T}^{1-s}\zeta (s)}{{(2\pi )}^{s}},\qquad T\sum\limits _{{\omega }_{p}}\frac{1}{\vert {\omega }_{p}{\vert }^{s}}=\frac{2{T}^{1-s}\zeta (s,1/2)}{{(2\pi )}^{s}}.\end{equation} \tag{ 4.8 }$$

Zoom In Zoom Out Reset image size

In case Ω_p becomes the continuous frequency for the zero temperature part, we can rewrite the diagram as a fermionic loop with an insertion of the fermionic self-energy at T = 0, which is given by ${{\Sigma}}_{1}(\mathbf{k})\equiv {{\Sigma}}_{1}(\mathbf{k},T=0)=-\frac{m{e}^{2}}{15{\pi }^{2}c}{\left(\frac{{\Lambda}}{\vert \mathbf{k}\vert }\right)}^{\varepsilon }\frac{{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }}{\varepsilon }$ (see reference [2]). We denote this as ${F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}^{(1)}$, which can be evaluated to

$$\begin{align}\hfill & {F}_{\mathit{coul}}^{(1)}=-T\,\sum\limits _{{\omega }_{{p}^{\prime }}}\int \frac{{d}^{4-\varepsilon }\mathbf{k}}{{(2\,\pi )}^{4-\varepsilon }}\mathrm{Tr}\left[{{\Sigma}}_{1}(\mathbf{k})\,G({\omega }_{{p}^{\prime }},\mathbf{k})\right]\hfill \\ \hfill & \quad =\frac{T\,m\,{e}^{2}\,{{\Lambda}}^{\varepsilon }}{15\,{\pi }^{2}\,c\,\varepsilon }\sum\limits _{{\omega }_{{p}^{\prime }}}\int \frac{{d}^{4-\varepsilon }\mathbf{k}}{{(2\,\pi )}^{4-\varepsilon }}\frac{\mathrm{Tr}\left[\left\{{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}G({\omega }_{{p}^{\prime }},\mathbf{k})\right]}{\vert \mathbf{k}{\vert }^{\varepsilon }}\hfill \\ \hfill & \quad =\frac{4\,T\,m\,{e}^{2}\,{{\Lambda}}^{\varepsilon }}{15\,{\pi }^{2}\,c\,\varepsilon }\sum\limits _{{\omega }_{{p}^{\prime }}}\int \frac{{d}^{4-\varepsilon }\mathbf{k}}{{(2\,\pi )}^{4-\varepsilon }}\frac{{\left\vert {\mathbf{d}}_{\mathbf{k}}\right\vert }^{2}}{\vert \mathbf{k}{\vert }^{\varepsilon }\left[{\omega }_{{p}^{\prime }}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right]}\hfill \\ \hfill & \quad =-\frac{{e}^{2}\,{{\Lambda}}^{\varepsilon }\,{m}^{3-\varepsilon }\,{\pi }^{\varepsilon -3}\mathrm{csc}\left(\frac{\pi \varepsilon }{2}\right)}{30\,c\,\varepsilon }\,T\sum\limits _{{\omega }_{{p}^{\prime }}}\frac{1}{\vert {\omega }_{{p}^{\prime }}{\vert }^{\varepsilon -2}}\hfill \\ \hfill & \quad \simeq -\frac{{e}^{2}\,{{\Lambda}}^{\varepsilon }\,\zeta (3){\left(m\,T\right)}^{3-\varepsilon }}{10\,{\pi }^{4}\,c\,\varepsilon }\hfill \\ \hfill & \qquad \left[\text{using}\;\zeta \left(\varepsilon -2,\frac{1}{2}\right)=\frac{-3\,\varepsilon \,{\zeta }^{\prime }(-2)}{4}\;\text{and}\;{\zeta }^{\prime }(-2)=-\frac{\zeta (3)}{4\,{\pi }^{2}}\right]\hfill \\ \hfill & \quad =-\frac{{e}^{2}\,\zeta (3){\left(m\,T\right)}^{3-\frac{\varepsilon }{2}}}{10\,{\pi }^{4}\,c\,\varepsilon }+\frac{{e}^{2}\,\zeta (3){\left(m\,T\right)}^{3-\frac{\varepsilon }{2}}\mathrm{ln}\left(\frac{m\,T}{{{\Lambda}}^{2}}\right)}{20\,{\pi }^{4}\,c},\hfill \end{align} \tag{ 4.9 }$$

where we have set ɛ = 0 in the numerical prefactor.

The second case is when ω_p' becomes the continuous frequency. Then the contribution is

$$\begin{align}\hfill {F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}^{(2)}& =\frac{4{e}^{2}{{\Lambda}}^{\varepsilon }T}{c}\sum\limits _{{{\Omega}}_{p}}\int \frac{{\mathrm{d}\omega }_{{p}^{\prime }}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{q}\,{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{8-2\varepsilon }}\hfill \\ \hfill & \quad \times \frac{\left[-\left({\omega }_{{p}^{\prime }}+{{\Omega}}_{p}\right){\omega }_{{p}^{\prime }}+{\mathbf{d}}_{\mathbf{k}+\mathbf{q}}\cdot {\mathbf{d}}_{\mathbf{k}}\right]}{{\mathbf{q}}^{2}\left[{\left({\omega }_{{p}^{\prime }}+{{\Omega}}_{p}\right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}+\mathbf{q}}{\vert }^{2}\right]\left({\omega }_{{p}^{\prime }}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}\right)}\hfill \\ \hfill & =\frac{4{e}^{2}{{\Lambda}}^{\varepsilon }T}{c}\sum\limits _{{{\Omega}}_{p}}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{q}\,{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{8-2\varepsilon }}\hfill \\ \hfill & \quad \times \frac{\left[-{\mathbf{d}}_{\mathbf{k}}+\frac{d{\left\{\mathbf{k}\cdot \left(\mathbf{k}+\mathbf{q}\right)\right\}}^{2}-{\mathbf{k}}^{2}{\left(\mathbf{k}+\mathbf{q}\right)}^{2}}{4{m}^{2}\left(d-1\right){\mathbf{d}}_{\mathbf{k}}}\right]}{{\mathbf{q}}^{2}\left[{({\mathbf{d}}_{\mathbf{k}}+{\mathbf{d}}_{\mathbf{k}+\mathbf{q}})}^{2}+{{\Omega}}^{2}\right]},\hfill \end{align} \tag{ 4.10 }$$

using ${\mathbf{d}}_{\mathbf{k}}\cdot {\mathbf{d}}_{\mathbf{p}}=\frac{d{\left(\mathbf{k}\cdot \mathbf{p}\right)}^{2}-{\mathbf{k}}^{2}{\mathbf{p}}^{2}}{4{m}^{2}\left(d-1\right)}$. Let θ be the angle between k and q. Making a change of variables as |q| = l|k| (where 0 < l < ∞), and defining $\text{ang}=\sqrt{1+{l}^{2}+2l\,\mathrm{c}\mathrm{o}\mathrm{s}\,\theta }$, we get

$$\begin{align}\hfill {F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}^{(2)}& =\frac{4{e}^{2}{{\Lambda}}^{\varepsilon }{m}^{3-\varepsilon }T}{c}\sum\limits _{{{\Omega}}_{p}}\int \frac{\mathrm{d}l\,{l}^{3-\varepsilon }\,\mathrm{d}\theta }{{(2\pi )}^{4-\varepsilon }}\hfill \\ \hfill & \quad \times \frac{\mathrm{sin}\,\theta {l}^{1-\varepsilon }\left[2l\left(\varepsilon -4\right)\mathrm{cos}\,\theta \left(l\mathrm{cos}\,\theta +2\right)+\left({\text{ang}}^{2}-1\right)(5-\varepsilon )\right]}{{2}^{8-\varepsilon }\left(3-\varepsilon \right){\pi }^{3-2\varepsilon }{\left({\text{ang}}^{2}+1\right)}^{4-\varepsilon }\mathrm{sin}\left(\frac{\pi \varepsilon }{2}\right)\vert {\Omega}{\vert }^{\varepsilon -2}}\hfill \\ \hfill & \simeq -\frac{19{e}^{2}{{\Lambda}}^{\varepsilon }\zeta (3){\left(mT\right)}^{3-\varepsilon }}{576{\pi }^{4}c}+\mathcal{O}\left(\varepsilon \right),\hfill \end{align} \tag{ 4.11 }$$

where we have set ɛ = 0 in the numerical prefactor. To this order of approximation, this part of the interaction correction to the free energy does not contain a pole in ɛ, and hence ${F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}^{(2)}$ does not contribute to the renormalization of the free energy.

Using ${{e}^{\ast }}^{2}=\frac{60{\pi }^{2}c\varepsilon }{19m}$, the final expression for the free energy up to two-loop order takes the form

$$\begin{align}\hfill & {F}^{(0)}(T)-{F}^{(0)}(0)+{F}_{\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{l}}^{(1)}+{F}_{\text{counterterms}}^{(1)}\hfill \\ \hfill & \qquad =-\frac{3\zeta (3){\left(mT\right)}^{3-\frac{\varepsilon }{2}}}{4{\pi }^{2}m}+\frac{{{e}^{\ast }}^{2}\zeta (3){\left(mT\right)}^{3-\frac{\varepsilon }{2}}\mathrm{ln}\left(\frac{mT}{{{\Lambda}}^{2}}\right)}{20{\pi }^{4}c}\hfill \\ \hfill & \qquad =-\frac{3\zeta (3){\left(mT\right)}^{3-\frac{\varepsilon }{2}}}{4{\pi }^{2}m}\left[1-\frac{4\varepsilon }{19}\,\mathrm{ln}\left(\frac{mT}{{{\Lambda}}^{2}}\right)\right]\hfill \\ \hfill & \qquad \simeq -\frac{3\zeta (3){m}^{2-\frac{\varepsilon }{2}}{T}^{3-\frac{\varepsilon }{2}-\frac{4\varepsilon }{19}}}{4{\pi }^{2}}{\left(\frac{{{\Lambda}}^{2}}{m}\right)}^{\frac{4\varepsilon }{19}},\hfill \end{align} \tag{ 4.12 }$$

after re-exponentiating the correction term coming from the two-loop diagrams and using equation (2.8). Note that the term ${F}_{\text{counterterms}}^{(1)}$ indicates that we have included the counterterms obtained from one-loop corrections in order to cancel out the singular $\frac{1}{\varepsilon }$ terms in the renormalized action [62]. Therefore, after including the leading order corrections, the free energy scales as

$$\begin{equation}F(T)\sim {T}^{3-\frac{\varepsilon }{2}-\frac{4\varepsilon }{19}}\;\stackrel{\varepsilon =1}{=}\;{T}^{\frac{5}{2}-\frac{4}{19}}.\end{equation} \tag{ 4.13 }$$

This implies that the specific heat scales as

$$\begin{equation}C(T)\sim {T}^{2-\frac{\varepsilon }{2}-\frac{4\varepsilon }{19}}\;\stackrel{\varepsilon =1}{=}\;{T}^{\frac{3}{2}-\frac{4}{19}}.\end{equation} \tag{ 4.14 }$$

Comparing with equation (4.2), and taking into account the fixed point value ${z}^{\ast }=2-\frac{4\varepsilon }{19}$ for the dynamical critical exponent, we find that there is a hyperscaling violation, which is proportional to ɛ. Finally, the entropy density s is, by definition, the derivative of the free energy with respect to T. Hence, using equation (4.13), we infer that

$$\begin{equation}s(T)\sim {T}^{2-\frac{\varepsilon }{2}-\frac{4\varepsilon }{19}}\;\stackrel{\varepsilon =1}{=}\;{T}^{2-\frac{27}{38}}.\end{equation} \tag{ 4.15 }$$

The momentum flux density tensor, also called the stress tensor, is given by

$$\begin{align}\hfill {T}_{\mu \nu }& =\sum\limits _{M}\left[\frac{\delta \mathcal{L}}{\delta ({\partial }_{\mu }{\zeta }_{M})}{\partial }_{\nu }{\zeta }_{M}-{\partial }_{\mu }\left\{\frac{\delta \mathcal{L}}{\delta \left({\partial }_{\alpha }{\partial }^{\alpha }{\zeta }_{M}\right)}\right\}{\partial }_{\nu }{\zeta }_{M}\right]\hfill \\ \hfill & \quad -{\delta }_{\mu \nu }\mathcal{L},\hfill \end{align} \tag{ 5.1 }$$

where $\mathcal{L}$ is the Lagrangian density, and ζ_M stands for all the quasiparticle fields in the theory. Evaluating this explicitly for the Luttinger semimetal, we get

$$\begin{align}\hfill {T}_{\mu \nu }({q}_{0},\mathbf{q})& =\sum\limits _{i}\int \frac{{\mathrm{d}}^{d}\mathbf{k}}{{(2\pi )}^{3}}\left({k}_{\nu }+{q}_{\nu }/2\right){\tilde{\psi }}_{i}^{{\dagger}}({k}_{0}+{q}_{0},\mathbf{k}+\mathbf{q})\hfill \\ \hfill & \quad \times \left[{\nabla }_{{k}_{\mu }}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right]{\tilde{\psi }}_{i}({k}_{0},\mathbf{k}).\hfill \end{align} \tag{ 5.2 }$$

The shear viscosity is the transport coefficient, which characterizes the relaxation of a transverse momentum gradient back to local equilibrium. Using the Kubo formula [69, 70], the optical shear viscosity is given by

$$\begin{align}\hfill \eta (\omega )& =\underset{\mathbf{q}\to 0}{\mathrm{lim}}\frac{1}{\omega }{\chi }_{{T}_{xy}{T}_{xy}}(\omega ,\mathbf{q})=\underset{\mathbf{q}\to 0}{\mathrm{lim}}\frac{1}{\omega }{\chi }_{{T}_{yz}{T}_{yz}}(\omega ,\mathbf{q})\hfill \\ \hfill & =\underset{\mathbf{q}\to 0}{\mathrm{lim}}\frac{1}{\omega }{\chi }_{{T}_{zx}{T}_{zx}}(\omega ,\mathbf{q}),\hfill \end{align} \tag{ 5.3 }$$

where

$$\begin{equation}{\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(\omega ,\mathbf{q})=\langle {T}_{\alpha \beta }{T}_{\alpha \beta }\rangle (\omega ,\mathbf{q})\end{equation} \tag{ 5.4 }$$

is the autocorrelation function of the component T_αβ of the stress tensor.

5.1. General arguments

The spatial components of T_αβ have the same scaling dimension as the Lagrangian density, and hence

$$\begin{equation}[{T}_{zx}]=d+z,\end{equation} \tag{ 5.5 }$$

which immediately leads to

$$\begin{align}\hfill [\eta ]& =2[{T}_{zx}]-z-[\text{Volume}\;\text{in}\;\text{energy}-\text{momentum}\;\text{space}]\hfill \\ \hfill & =2\left(d+z\right)-z-d-z=d.\hfill \end{align} \tag{ 5.6 }$$

This shows that [η] has the same scaling dimension as the entropy density s, and hence the ratio η/s is dimensionless, as expected. If hyperscaling is not violated, equation (5.6) also implies the scaling form

$$\begin{equation}\eta (\omega )\sim {\omega }^{d/z}\end{equation} \tag{ 5.7 }$$

for optical viscosity (i.e., for ω ≫ T), and

$$\begin{equation}\eta (T)\sim {T}^{d/z}\end{equation} \tag{ 5.8 }$$

for dc viscosity (i.e., for T ≫ ω).

5.2. Optical viscosity

In this subsection, we will find the expression of the optical viscosity at T = 0, using the fixed point theory.

5.2.1. One-loop contribution

The autocorrelation function ${\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(\omega ,\mathbf{0})$ at one-loop level is given by a simple fermionic loop with two insertions of the T_zx operator, as shown in figure 7. For the Luttinger semimetal, it evaluates to

$$\begin{align}\hfill & {\langle {T}_{zx}{T}_{zx}\rangle }_{1\;\text{loop}}(\mathrm{i}\omega )\hfill \\ \hfill & \qquad =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int {\mathrm{d}}^{4-\varepsilon }\mathbf{k}{k}_{x}^{2}\,\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k+q)\right.\hfill \\ \hfill & \qquad \quad \left.\times \left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k)\right]\hfill \\ \hfill & \qquad =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int {\mathrm{d}}^{4-\varepsilon }\mathbf{k}{k}_{x}^{2}\,\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \qquad \left.\times \frac{\mathrm{i}{k}_{0}+\mathrm{i}\omega +{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }}{-{\left(\mathrm{i}{k}_{0}+\mathrm{i}\omega \right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}}\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}\frac{\mathrm{i}{k}_{0}+{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }}{-{\left(\mathrm{i}{k}_{0}\right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}}\right]\hfill \\ \hfill & \qquad =4\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int {\mathrm{d}}^{4-\varepsilon }\mathbf{k}{k}_{x}^{2}\hfill \\ \hfill & \quad \qquad \times \frac{{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\right\}}^{2}({k}_{0}+\omega ){k}_{0}-\frac{1}{2}{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}^{2}\right\}}^{2}+{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\right\}}^{2}{\mathbf{d}}_{\mathbf{k}}^{2}}{\left[-{\left(\mathrm{i}{k}_{0}+\mathrm{i}\omega \right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}\right]\left[-{\left(\mathrm{i}{k}_{0}\right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}\right]}\hfill \\ \hfill & \qquad =-\frac{{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}}{40\pi },\hfill \end{align} \tag{ 5.9 }$$

where q = (ω, 0). Consequently, at zeroth order, the optical viscosity η(ω) is proportional to ${\omega }^{2-\frac{\varepsilon }{2}}$ for ω ≫ T. For d = 4 − ɛ, this scaling then obeys the hyperscaling ω^d/z , as obtained from general arguments.

Zoom In Zoom Out Reset image size

5.2.2. Two-loop contributions

At two-loop order, we need to consider three Feynman diagrams, as shown in figures 8(a)–(c). The first two diagrams (figures 8(a) and (b)) include the fermion self-energy corrections. Evaluating them explicitly, we get

$$\begin{align}\hfill & {\langle {T}_{zx}{T}_{zx}\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )\hfill \\ \hfill & \qquad =-2\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{d}}{k}_{x}^{2}\,\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k+q)\right.\hfill \\ \hfill & \quad \qquad \left.\times {{\Sigma}}_{1}(k+q){G}_{0}(k+q)\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k)\right]\hfill \\ \hfill & \qquad =-\frac{{e}^{2}{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}}{675{\pi }^{3}c\varepsilon }{\left(\frac{{{\Lambda}}^{2}}{m\vert \omega \vert }\right)}^{\varepsilon /2}\hfill \\ \hfill & \quad \qquad +\frac{{e}^{2}{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}\,\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)}{1350{\pi }^{3}c\varepsilon }.\hfill \end{align} \tag{ 5.10 }$$

Note that we have included a factor of 2, since the figures 8(a) and (b) give equal contributions. The calculation of the loop integrals follows the same steps as the computation of the corresponding two-loop current–current correlator in reference [14].

Zoom In Zoom Out Reset image size

The diagram in figure 8(c) is the one with a vertex correction, and it evaluates to

$$\begin{align}\hfill & {\langle {T}_{zx}{T}_{zx}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )=\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{c}\int \frac{\mathrm{d}{k}_{0}\,\mathrm{d}{\ell }_{0}}{{(2\pi )}^{2}}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}\,{\mathrm{d}}^{4-\varepsilon }\boldsymbol{\ell }}{{(2\pi )}^{2(4-\varepsilon )}}{k}_{x}\left({k}_{x}+{\ell }_{x}\right)\hfill \\ \hfill & \qquad \quad \times \frac{\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}({k}_{0}+\omega ,\mathbf{k}){G}_{0}({\ell }_{0}+\omega ,\mathbf{k}+\boldsymbol{\ell })\left\{{\partial }_{{k}_{z}+{\ell }_{z}}{\mathbf{d}}_{\mathbf{k}+\boldsymbol{\ell }}\cdot \mathbf{\Gamma }\right\}{G}_{0}({\ell }_{0},\mathbf{k}+\boldsymbol{\ell }){G}_{0}({k}_{0},\mathbf{k})\right]}{{\boldsymbol{\ell }}^{2}}.\hfill \end{align} \tag{ 5.11 }$$

Again, the details of the evaluation of the loop integrals follow the same steps as the computation of the corresponding two-loop current–current correlator in reference [14], finally leading to

$$\begin{align}\hfill {\langle {T}_{zx}{T}_{zx}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )& =-\frac{37{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}}{57600{\pi }^{3}c\varepsilon }{\left(\frac{{{\Lambda}}^{2}}{m\vert \omega \vert }\right)}^{\varepsilon /2}\hfill \\ \hfill & \quad +\frac{37{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}\,\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)}{115200{\pi }^{3}c}.\hfill \end{align} \tag{ 5.12 }$$

5.2.3. Scaling of optical viscosity at T = 0

Including the counterterms of the renormalized action [62], using equation (2.8) at the LAB fixed point, the final expression for the autocorrelator of the stress tensor T_zx takes the form:

$$\begin{align}\hfill \left\langle {T}_{zx}{T}_{zx}\right\rangle (\mathrm{i}\omega )& ={\left\langle {T}_{zx}{T}_{zx}\right\rangle }_{1\;\text{loop}}(\mathrm{i}\omega )+{\left\langle {T}_{zx}{T}_{zx}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )\hfill \\ \hfill & \quad +{\left\langle {T}_{zx}{T}_{zx}\right\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )+{\left\langle {T}_{zx}{T}_{zx}\right\rangle }_{\text{counterterms}}(\mathrm{i}\omega )\hfill \\ \hfill & =-\frac{{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}}{40\pi }+\frac{{{e}^{\ast }}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)}{1350{\pi }^{3}c\varepsilon }\hfill \\ \hfill & \quad +\frac{37{{e}^{\ast }}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)}{115200{\pi }^{3}c}\hfill \\ \hfill & =-\frac{{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}}}{40\pi }\left[1-\frac{367\varepsilon }{2736}\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)\right]\hfill \\ \hfill & =-\frac{{m}^{2-\frac{\varepsilon }{2}}\vert \omega {\vert }^{3-\frac{\varepsilon }{2}-\frac{367\varepsilon }{2736}}}{40\pi }{\left(\frac{{{\Lambda}}^{2}}{m}\right)}^{\frac{367\varepsilon }{2736}},\hfill \end{align} \tag{ 5.13 }$$

up to two-loop order. As usual, we have re-exponentiated the correction term coming from the two-loop diagrams. Therefore, the corrected optical viscosity in the clean limit scales as

$$\begin{equation}\eta (\omega )\sim {\omega }^{2-\frac{\varepsilon }{2}-\frac{367\varepsilon }{2736}},\end{equation} \tag{ 5.14 }$$

after including the leading order corrections due to the presence of the Coulomb interactions. Comparing with equation (5.6), we find that there is indeed a small hyperscaling violation proportional to ɛ.

5.3. DC viscosity

We now focus on studying the shear viscosity η(ω = 0, T) at T > 0. Following the same strategy explained previously in the context of the Raman response at T > 0, we will use the memory matrix method, and work directly in d = 3 spatial dimensions. Moreover, we will also assume here the presence of short-ranged weak disorder (that breaks translational invariance) as the source of momentum relaxation. In the memory matrix formalism, the shear viscosity is given by

$$\begin{align}\hfill \eta (\omega ,T)& ={\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T)\left[\frac{1}{{M}_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T)-\mathrm{i}\omega {\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T)}\right]\hfill \\ \hfill & \quad \times {\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T),\hfill \end{align} \tag{ 5.15 }$$

where ${\chi }_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T)={\int }_{0}^{\beta }\mathrm{d}\tau \left\langle {T}_{\alpha \beta }(\tau ,\mathbf{q}=0){T}_{\alpha \beta }(0,\mathbf{q}=0)\right\rangle$ is the static susceptibility of the component T_αβ of the stress tensor, and ${M}_{{T}_{\alpha \beta }{T}_{\alpha \beta }}(T)$ is the corresponding memory matrix element (this is defined in a way similar to the one defined in equation (3.17)). We want to emphasize that the stress tensor is also a nearly-conserved operator for the Luttinger semimetal, because the equation of motion of T_αβ due to the effects of weak disorder is

$$\begin{align}\hfill {\dot{T}}_{\alpha \beta }& =\mathrm{i}\left[H+{H}_{\mathrm{i}\mathrm{m}\mathrm{p}},{T}_{\alpha \beta }\right]\hfill \\ \hfill {\Rightarrow}{\dot{T}}_{\alpha \beta }& =\mathrm{i}\sum\limits _{i=1}^{{N}_{\mathrm{f}}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}}{{(2\pi )}^{3}}\int \frac{{\mathrm{d}}^{3}{\mathbf{k}}^{\prime }}{{(2\pi )}^{3}}{k}_{\beta }\tilde{W}({\mathbf{k}}^{\prime })\left[{\psi }_{i}^{{\dagger}}({\mathbf{k}+\mathbf{k}}^{\prime })\right.\hfill \\ \hfill & \quad \left.\times ({\partial }_{{k}_{\alpha }}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }){\psi }_{i}(\mathbf{k})-{\psi }_{i}^{{\dagger}}(\mathbf{k})({\partial }_{{k}_{\alpha }}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }){\psi }_{i}({\mathbf{k}-\mathbf{k}}^{\prime })\right]\hfill \end{align} \tag{ 5.16 }$$

for a generic value of N_f, where $\tilde{W}({\mathbf{k}}^{\prime })$ is the Fourier transform of the impurity coupling strength W(x) (defined in section 3.2).

The susceptibility ${\chi }_{{T}_{zx}{T}_{zx}}^{1\;\text{loop}}(T)$ at one-loop level is given by a simple fermionic loop with two insertions of the T_zx operator, as shown in figure 7, and hence is captured by

$$\begin{align}\hfill {\chi }_{{T}_{zx}{T}_{zx}}^{1\;\text{loop}}(T)& =-T\sum\limits _{{k}_{0}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}}{{(2\pi )}^{3}}{k}_{x}^{2}\,\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times {G}_{0}(k)\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k)\right]\hfill \\ \hfill & =4T\sum\limits _{{k}_{0}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}}{{(2\pi )}^{3}}{k}_{x}^{2}\hfill \\ \hfill & \quad \times \frac{{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\right\}}^{2}{k}_{0}^{2}-\frac{1}{2}{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}^{2}\right\}}^{2}+{\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\right\}}^{2}{\mathbf{d}}_{\mathbf{k}}^{2}}{{\left[-{\left(\mathrm{i}{k}_{0}\right)}^{2}+\vert {\mathbf{d}}_{\mathbf{k}}{\vert }^{2}\right]}^{2}}.\hfill \end{align} \tag{ 5.17 }$$

At two-loop order, we need to consider three Feynman diagrams, as shown in figures 8(a)–(c). We will denote the sum of these contributions by ${\chi }_{{T}_{zx}{T}_{zx}}^{2\;\text{loop}}(T)={\chi }_{{T}_{zx}{T}_{zx}}^{(2,1)}(T)+{\chi }_{{T}_{zx}{T}_{zx}}^{(2,2)}(T)$, such that ${\chi }_{{T}_{zx}{T}_{zx}}^{(2,1)}(T)$ is the sum of first two diagrams [cf figures 8(a) and (b)] which include the fermion self-energy corrections.

We find that

$$\begin{align}\hfill {\chi }_{{T}_{zx}{T}_{zx}}^{(2,1)}(T)& =-2T\sum\limits _{{k}_{0}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}}{{(2\pi )}^{3}}{k}_{x}^{2}\,\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k+q)\right.\hfill \\ \hfill & \quad \left.\times {{\Sigma}}_{T}(k+q){G}_{0}(k+q)\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}(k)\right],\hfill \end{align} \tag{ 5.18 }$$

where the temperature-dependent contribution of the self-energy is ${{\Sigma}}_{T}(\mathbf{k})=-\frac{{e}^{2}{\Lambda}}{15{\pi }^{2}c}\frac{({\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma })}{T}$. Note that we have included a factor of 2 in the expression above, since the figures 8(a) and (b) give equal contributions. This calculation involves steps similar to those involved in the computation of the corresponding two-loop diagram for current–momentum susceptibility in reference [14].

The diagram in figure 8(c) is the one with a vertex correction, and it evaluates to

$$\begin{align}\hfill & {\chi }_{{T}_{zx}{T}_{zx}}^{(2,2)}(T)=\frac{{e}^{2}{T}^{2}}{c}\sum\limits _{{k}_{0},{\ell }_{0}}\int \frac{{\mathrm{d}}^{3}\mathbf{k}\,{\mathrm{d}}^{3}\boldsymbol{\ell }}{{(2\pi )}^{6}}{k}_{x}\left({k}_{x}+{\ell }_{x}\right)\hfill \\ \hfill & \qquad \times \frac{\mathrm{Tr}\left[\left\{{\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }\right\}{G}_{0}({k}_{0},\mathbf{k}){G}_{0}({\ell }_{0},\mathbf{k}+\boldsymbol{\ell })\left\{{\partial }_{{k}_{z}+{\ell }_{z}}{\mathbf{d}}_{\mathbf{k}+\boldsymbol{\ell }}\cdot \mathbf{\Gamma }\right\}{G}_{0}({\ell }_{0},\mathbf{k}+\boldsymbol{\ell }){G}_{0}({k}_{0},\mathbf{k})\right]}{{\boldsymbol{\ell }}^{2}}.\hfill \end{align} \tag{ 5.19 }$$

Again, the details of the above calculation follow steps similar to those in reference [14]. We would like to emphasize that although an analytical expression of ${\chi }_{{T}_{zx}{T}_{zx}}^{2\;\text{loop}}(T)$ cannot be obtained in a closed-form, its value is estimated by using standard numerical integration methods. We use the numerical data into a fitting function in order to extract the temperature dependence, which gives us ${\chi }_{{T}_{zx}{T}_{zx}}(T)=\mathcal{A}+\mathcal{B}{T}^{5/2}+\mathcal{C}/T$, where $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ are temperature-independent constants. In fact, $\mathcal{A}$ and $\mathcal{C}$ are non-universal, since they depend on the hard UV momentum cutoff Λ. At low temperatures, we get ${\chi }_{{T}_{zx}{T}_{zx}}(T)\approx \mathcal{A}+\mathcal{C}/T$.

At leading order, the memory matrix element is given by

$$\begin{align}\hfill & {M}_{{T}_{zx}{T}_{zx}}(T)=-{W}_{0}^{2}\,\underset{\omega \to 0}{\mathrm{lim}}\,\frac{\text{Im}{\left.\left[\int \frac{{\mathrm{d}}^{3}\mathbf{k}\,{\mathrm{d}}^{3}\mathbf{q}}{{(2\pi )}^{6}}{k}_{x}^{2}T\sum\limits _{{k}_{0}}\mathrm{Tr}\left[({\partial }_{{k}_{z}+{q}_{z}}{\mathbf{d}}_{\mathbf{k+q}}\cdot \mathbf{\Gamma }){G}_{0}(\omega +{k}_{0},\mathbf{k}+\mathbf{q})({\partial }_{{k}_{z}}{\mathbf{d}}_{\mathbf{k}}\cdot \mathbf{\Gamma }){G}_{0}({k}_{0},\mathbf{k})\right]\right]\right\vert }_{\text{i}\omega \to \omega +\mathrm{i}\delta }}{\omega }.\hfill \end{align} \tag{ 5.20 }$$

After computing the trace, we perform the analytical continuation, and the resulting integral is evaluated numerically. This yields ${M}_{{T}_{zx}{T}_{zx}}(T)\sim {\mathcal{A}}^{\prime }+{\mathcal{B}}^{\prime }/T$, with ${\mathcal{A}}^{\prime }$ and ${\mathcal{B}}^{\prime }$ being (non-universal) temperature-independent constants that depend only on the UV cutoff Λ. At low temperatures, we obtain ${M}_{{T}_{zx}{T}_{zx}}(T)\approx {\mathcal{B}}^{\prime }/T$.

Using equation (5.15), in the ω → 0 limit, we get $\eta (T)={\chi }_{{T}_{zx}{T}_{zx}}^{2}(T)/{M}_{{T}_{zx}{T}_{zx}}(T)$. Therefore, at low temperatures,

$$\begin{equation}\eta (T)\sim {T}^{\lambda },\quad \text{where}\ 0< \lambda < 1.\end{equation} \tag{ 5.21 }$$

5.4. Scaling of the ratio η/s at T > 0

Equations (5.21) and (4.15) give us the ratio

$$\begin{equation}\eta /s\sim {T}^{\lambda -\frac{49}{38}}\end{equation} \tag{ 5.22 }$$

in d = 3 spatial dimensions. The scaling implies that η/s diverges at low temperatures, instead of saturating the universal bound. This divergent behavior is similar to what was found in NFL metals near the Ising-nematic critical point [55], where the quantum critical phase has a critical Fermi surface, and the divergence is a consequence of the hyperscaling violation due to the presence of this Fermi surface [41]. We also note that this singular dependence of η/s on T is at variance with the results found by earlier works in references [9, 71], where the authors used a quantum Boltzmann equation approach.

In this paper, we have performed a two-loop calculation to determine the Raman scattering response at intermediate frequencies in the LAB phase. We have found that at frequencies larger than the temperature, the Raman response exhibits a power-law behavior, which can be verified experimentally. At low frequencies, the Raman response displays instead a clear quasi-elastic peak. In the second part, we have computed the ratio of the shear viscosity and the entropy density η/s. We have demonstrated that in the physically relevant case (i.e., N_f = 1), this ratio diverges at low temperatures, instead of saturating the Kovtun–Son–Starinets universal bound. This divergent behavior has some similarities to what was found in the quantum critical phase near an Ising-nematic critical point [55]. Consequently, it can be traced to the violation of the hyperscaling property that emerges in the LAB phase.

From our results, we can find the scaling of the thermal diffusivity ${D}_{\mathrm{t}\mathrm{h}}=\kappa /C$, where κ is the thermal conductivity at zero current. Since the scalings of the specific heat and the thermal conductivity are described by $C\sim {T}^{2-\frac{\varepsilon }{2}-\frac{4\varepsilon }{19}}$ [cf equation (4.14)] and κ ∼ uT⁻¹ [15] (where u is a parametrically small prefactor), respectively, we get ${D}_{\mathrm{t}\mathrm{h}}\sim u {T}^{-3+\frac{\varepsilon }{2}+\frac{4\varepsilon }{19}}$.

We thank Xiaoyu Wang for participating in the initial stages of the project. HF acknowledges funding from CNPq under grant numbers 310710/2018-9 and 311428/2021-5.

All data that support the findings of this study are included within the article (and any supplementary files).

We derive a set of useful relations [3, 4] for the vector functions d(k) (whose components d_a (k) are the l = 2 spherical harmonics in d spatial dimensions) and the generalized real d × d Gell–Mann matrices Λ_a (a = 1, 2, ..., N) in d-dimensions. The matrices Λ_a are symmetric, traceless, and orthogonal, satisfying

$$\begin{align}\hfill & \mathrm{Tr}[{{\Lambda}}^{a}{{\Lambda}}^{b}]=2{\delta }_{ab},\hfill \\ \hfill & \sum\limits _{a=1}^{N}{\left({{\Lambda}}^{a}\right)}_{ij}\left({{\Lambda}}_{l{j}^{\prime }}^{a}\right)={\delta }_{il}{\delta }_{j{j}^{\prime }}+{\delta }_{i{j}^{\prime }}{\delta }_{jl}-\frac{2}{d}{\delta }_{ij}{\delta }_{l{j}^{\prime }}.\hfill \end{align} \tag{ A1 }$$

As a result, the index a (or b) runs from 1 to $N=\frac{\left(d-1\right)\left(d+2\right)}{2}$. We define the components of d(k) by

$$\begin{equation}{d}_{a}(\mathbf{k})=\sqrt{\frac{d}{2\left(d-1\right)}}\sum\limits _{i,j=1}^{d}\frac{{k}_{i}{\left({{\Lambda}}^{a}\right)}_{ij}{k}_{j}}{2m}.\end{equation} \tag{ A2 }$$

This gives

$$\begin{equation}\sum\limits _{a=1}^{N}{d}_{a}(\mathbf{k}){d}_{a}(\mathbf{p})=\frac{d{\left(\mathbf{k}\cdot \mathbf{p}\right)}^{2}-{\mathbf{k}}^{2}{\mathbf{p}}^{2}}{4{m}^{2}\left(d-1\right)}.\end{equation} \tag{ A3 }$$

The trace of the six gamma matrices, used for the loop-calculations, is given by

$$\begin{align}\hfill \mathrm{Tr}\left({{\Gamma}}_{a}{{\Gamma}}_{b}{{\Gamma}}_{c}{{\Gamma}}_{d}{{\Gamma}}_{e}{{\Gamma}}_{f}\right)& =4\left[{\delta }_{ab}\left({\delta }_{cd}{\delta }_{ef}-{\delta }_{ce}{\delta }_{df}+{\delta }_{cf}{\delta }_{de}\right)\right.\hfill \\ \hfill & \quad \left.-{\delta }_{ac}\left({\delta }_{bd}{\delta }_{ef}-{\delta }_{be}{\delta }_{df}+{\delta }_{bf}{\delta }_{de}\right)\right.\hfill \\ \hfill & \quad \left.+{\delta }_{ad}\left({\delta }_{bc}{\delta }_{ef}-{\delta }_{be}{\delta }_{cf}+{\delta }_{bf}{\delta }_{ce}\right)\right.\hfill \\ \hfill & \quad \left.-{\delta }_{ae}\left({\delta }_{bc}{\delta }_{df}-{\delta }_{bd}{\delta }_{cf}+{\delta }_{bf}{\delta }_{cd}\right)\right.\hfill \\ \hfill & \quad \left.+{\delta }_{af}\left({\delta }_{bc}{\delta }_{de}-{\delta }_{bd}{\delta }_{ce}+{\delta }_{be}{\delta }_{cd}\right)\right].\hfill \end{align} \tag{ A4 }$$

In this appendix, we provide the details of the two-loop calculations for the T = 0 Raman response.

As explained in the main text, the first two-loop-order correction involves inserting the one-loop fermion self-energy (Σ₁) corrections into the correlator, which takes the form:

$$\begin{align}\hfill {\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )& =-\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{4-\varepsilon }}\,\mathrm{Tr}\left[{{\Gamma}}_{a}{G}_{0}(k+q)\right.\hfill \\ \hfill & \quad \left.\times {{\Sigma}}_{1}(\mathbf{k}){G}_{0}(k+q){{\Gamma}}_{b}{G}_{0}(k)\right]\hfill \\ \hfill & \quad -\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{4-\varepsilon }}\,\mathrm{Tr}\left[{{\Gamma}}_{a}{G}_{0}(k+q)\right.\hfill \\ \hfill & \quad \left.\times {{\Gamma}}_{b}{G}_{0}(k){{\Sigma}}_{1}(\mathbf{k}){G}_{0}(k)\right],\hfill \end{align} \tag{ B1 }$$

where ${{\Sigma}}_{1}(\boldsymbol{\ell })=-\frac{m{e}^{2}}{15{\pi }^{2}c}{\left(\frac{{\Lambda}}{\vert \boldsymbol{\ell }\vert }\right)}^{\varepsilon }\frac{\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }}{\varepsilon }$ (see reference [2]). This gives us

$$\begin{align}\hfill {\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )& =\frac{m{e}^{2}}{15{\pi }^{2}c\varepsilon }\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{4-\varepsilon }}{\left(\frac{{\Lambda}}{\vert \mathbf{k}\vert }\right)}^{\varepsilon }\hfill \\ \hfill & \quad \times \frac{{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{1}}{\left\{-{\mathcal{B}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right\}{\left\{-{\mathcal{A}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right\}}^{2}}\hfill \\ \hfill & \quad +\frac{m{e}^{2}}{15{\pi }^{2}c\varepsilon }\int \frac{\mathrm{d}{k}_{0}}{2\pi }\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}}{{(2\pi )}^{4-\varepsilon }}{\left(\frac{{\Lambda}}{\vert \mathbf{k}\vert }\right)}^{\varepsilon }\hfill \\ \hfill & \quad \times \frac{{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{2}}{{\left\{-{\mathcal{B}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right\}}^{2}\left\{-{\mathcal{A}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right\}},\hfill \end{align} \tag{ B2 }$$

where

$$\begin{equation}\mathcal{A}=\mathrm{i}{k}_{0}+\mathrm{i}\omega -\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }},\qquad \mathcal{B}=\mathrm{i}{k}_{0}-\frac{{\mathbf{k}}^{2}}{2{m}^{\prime }},\end{equation} \tag{ B3 }$$

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{1}& =\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right],\hfill \\ \hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{2}& =\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right].\hfill \end{align} \tag{ B4 }$$

First, let us evaluate term₁ for a = b = 0 as follows:

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{1}& =\mathrm{Tr}\left[\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\mathcal{A}\,\mathrm{Tr}\left[\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\left({\mathcal{A}}^{2}+2\mathcal{A}\mathcal{B}\right)\mathrm{Tr}\left[\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\left({\mathcal{A}}^{2}+2\mathcal{A}\mathcal{B}\right)\frac{{\mathbf{k}}^{4}}{{m}^{2}}+\frac{{\mathbf{k}}^{8}}{4{m}^{4}}.\hfill \end{align} \tag{ B5 }$$

Similarly, we get:

$$\begin{equation}{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{2}=\left({\mathcal{B}}^{2}+2\mathcal{A}\mathcal{B}\right)\frac{{\mathbf{k}}^{4}}{{m}^{2}}+\frac{{\mathbf{k}}^{8}}{4{m}^{4}}.\end{equation} \tag{ B6 }$$

Evaluating the integrals, we obtain:

$$\begin{equation}{\left\langle {\rho }_{0}{\rho }_{0}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )=0.\end{equation} \tag{ B7 }$$

Let us now evaluate term₁ for a, b > 0 as follows:

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{1}& =\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times {{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\mathcal{A}\,\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times {{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =2\mathcal{A}\mathcal{B}\,\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\right]\hfill \\ \hfill & \quad +{\mathcal{A}}^{2}\,\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =2\mathcal{A}\mathcal{B}{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\mathrm{Tr}\left[{{\Gamma}}_{a}{{\Gamma}}_{b}\right]+\left({\mathcal{A}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\hfill \\ \hfill & \quad \times \mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =8\mathcal{A}\mathcal{B}{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}+4\left({\mathcal{A}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\hfill \\ \hfill & \quad \times \left[2{d}_{a}(\mathbf{k}){d}_{b}(\mathbf{k})-{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}\right]\hfill \\ \hfill & =4\left[2\mathcal{A}\mathcal{B}-\left({\mathcal{A}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\right]{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}\hfill \\ \hfill & \quad +8\left({\mathcal{A}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right){d}_{a}(\mathbf{k}){d}_{b}(\mathbf{k}).\hfill \end{align} \tag{ B8 }$$

Similarly, for a, b > 0, term₂ evaluates to:

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{2}& =\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\mathcal{B}\,\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & =\left({\mathcal{B}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\mathrm{Tr}\left[{{\Gamma}}_{a}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}\right]\hfill \\ \hfill & \quad +2\mathcal{A}\mathcal{B}{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\mathrm{Tr}\left[{{\Gamma}}_{a}{{\Gamma}}_{b}\right]\hfill \\ \hfill & =4\left({\mathcal{B}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\left[2{d}_{a}(\mathbf{k}){d}_{b}(\mathbf{k})-{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}\right]\hfill \\ \hfill & \quad +8\mathcal{A}\mathcal{B}{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}\hfill \\ \hfill & =4\left[2\mathcal{A}\mathcal{B}-\left({\mathcal{B}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right)\right]{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}{\delta }_{ab}\hfill \\ \hfill & \quad +8\left({\mathcal{B}}^{2}+{\left\vert \mathbf{d}(\mathbf{k})\right\vert }^{2}\right){d}_{a}(\mathbf{k}){d}_{b}(\mathbf{k}).\hfill \end{align} \tag{ B9 }$$

Evaluating the integrals, we finally obtain:

$$\begin{equation}{\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{2\;\text{loop}}^{(1)}(\mathrm{i}\omega )=\frac{{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{75{\pi }^{3}c}\left[\frac{1}{\varepsilon }-\frac{1}{2}\,\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)\right]{\delta }_{ab}.\end{equation} \tag{ B10 }$$

The second two-loop-order correction are the vertex-like corrections, which can be expressed as

$$\begin{align}\hfill & {\left\langle {\rho }_{a}{\rho }_{b}\right\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )\hfill \\ \hfill & \qquad =\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{c}\int \frac{\mathrm{d}{k}_{0}\,\mathrm{d}{\ell }_{0}}{{(2\pi )}^{2}}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}\,{\mathrm{d}}^{4-\varepsilon }\boldsymbol{\ell }}{{(2\pi )}^{2(4-\varepsilon )}}\mathrm{Tr}\left[{{\Gamma}}_{a}{G}_{0}(k+q)\frac{1}{{\boldsymbol{\ell }}^{2}}\right.\hfill \\ \hfill & \quad \qquad \left.\times {G}_{0}(k+q+\ell ){{\Gamma}}_{b}{G}_{0}(k+\ell ){G}_{0}(k)\right]\hfill \\ \hfill & \qquad =\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{c}\int \frac{\mathrm{d}{k}_{0}\,\mathrm{d}{\ell }_{0}}{{(2pi)}^{2}}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}\,{\mathrm{d}}^{4-\varepsilon }\boldsymbol{\ell }}{{(2\pi )}^{2(4-\varepsilon )}}\mathrm{Tr}\left[{{\Gamma}}_{a}{G}_{0}({k}_{0}+\omega ,\mathbf{k})\right.\hfill \\ \hfill & \qquad \quad \left.\times \frac{1}{{\boldsymbol{\ell }}^{2}}{G}_{0}({\ell }_{0}+\omega ,\mathbf{k}+\boldsymbol{\ell }){{\Gamma}}_{b}{G}_{0}({\ell }_{0},\mathbf{k}+\boldsymbol{\ell }){G}_{0}({k}_{0},\mathbf{k})\right],\hfill \end{align} \tag{ B11 }$$

with q = (ω, 0) and ℓ = (ℓ₀, ℓ ). After some convenient regrouping of the terms in the integrand, the expression to be evaluated takes the form:

$$\begin{align}\hfill & {\langle {\rho }_{a}{\rho }_{b}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )\hfill \\ \hfill & \qquad =\frac{{e}^{2}{{\Lambda}}^{\varepsilon }}{c}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}\,{\mathrm{d}}^{4-\varepsilon }\boldsymbol{\ell }}{{(2\pi )}^{2(4-\varepsilon )}}\,\mathrm{Tr}\left[\frac{\int \frac{\mathrm{d}{\ell }_{0}}{2\pi }{G}_{0}({\ell }_{0}-\omega ,\boldsymbol{\ell }){{\Gamma}}_{a}{G}_{0}({\ell }_{0},\boldsymbol{\ell })\int \frac{\mathrm{d}{k}_{0}}{2\pi }{G}_{0}({k}_{0}+\omega ,\mathbf{k}){{\Gamma}}_{b}{G}_{0}({k}_{0},\mathbf{k})}{{\left(\mathbf{k}+\boldsymbol{\ell }\right)}^{2}}\right].\hfill \end{align} \tag{ B12 }$$

It is clear from this expression that ${\langle {\rho }_{a}{\rho }_{b}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )\propto {\delta }_{ab}$. Hence, we will consider ${\langle {\rho }_{a}{\rho }_{a}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )$ (no sum over a).

• (a)  
  ${\langle {\rho }_{a}{\rho }_{a}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )=0$.
• (b)  
  ${\langle {\rho }_{a}{\rho }_{0}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )=0$.
• (c)  
  Next, let us consider a, b > 0. Let us evaluate the k₀-integral, such that

  $$\begin{align}\hfill & \int \frac{\mathrm{d}{k}_{0}}{2\pi }{G}_{0}({k}_{0}+\omega ,\mathbf{k}){{\Gamma}}_{b}{G}_{0}({k}_{0},\mathbf{k})\hfill \\ \hfill & \qquad =\int \frac{\mathrm{d}{k}_{0}}{2\pi }\frac{\left\{\mathcal{A}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathcal{B}+\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}}{\left(-{\mathcal{A}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right)\left(-{\mathcal{B}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right)}\hfill \\ \hfill & \qquad =\int \frac{\mathrm{d}{\ell }_{0}}{2\pi }\frac{\mathcal{A}\mathcal{B}{{\Gamma}}_{b}+\mathcal{A}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}+\mathcal{B}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}+\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}}{\left(-{\mathcal{A}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right)\left(-{\mathcal{B}}^{2}+\vert \mathbf{d}(\mathbf{k}){\vert }^{2}\right)}\hfill \\ \hfill & \qquad ={m}^{3}\frac{-\frac{{\mathbf{k}}^{4}}{2{m}^{2}}{{\Gamma}}_{b}-\mathrm{i}\omega {{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}+\mathrm{i}\omega \left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}+2\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}{{\Gamma}}_{b}\left\{\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right\}}{{\mathbf{k}}^{2}\left({\mathbf{k}}^{4}+{m}^{2}{\omega }^{2}\right)}.\hfill \end{align} \tag{ B13 }$$

Plugging this in, we get:

$$\begin{align}\hfill & {\langle {\rho }_{a}{\rho }_{a}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )\hfill \\ \hfill & \qquad =\frac{{e}^{2}{{\Lambda}}^{\varepsilon }{m}^{6}}{c}\int \frac{{\mathrm{d}}^{4-\varepsilon }\mathbf{k}\,{\mathrm{d}}^{4-\varepsilon }\boldsymbol{\ell }}{{(2\pi )}^{2(4-\varepsilon )}}\hfill \\ \hfill & \qquad \quad \times \frac{{\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{3}}{{\left(\mathbf{k}+\boldsymbol{\ell }\right)}^{2}{\boldsymbol{\ell }}^{2}{\mathbf{k}}^{2}\left({\boldsymbol{\ell }}^{4}+{m}^{2}{\omega }^{2}\right)\left({\mathbf{k}}^{4}+{m}^{2}{\omega }^{2}\right)},\hfill \end{align} \tag{ B14 }$$

where

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{3}& =\mathrm{Tr}\left[\left\{-\frac{{\boldsymbol{\ell }}^{4}{{\Gamma}}_{a}}{2{m}^{2}}+\mathrm{i}\omega {{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right)\right.\right.\hfill \\ \hfill & \quad \left.\left.-\mathrm{i}\omega \left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}+2\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right)\right\}\right.\hfill \\ \hfill & \quad \left.\times \left\{-\frac{{\mathbf{k}}^{4}{{\Gamma}}_{a}}{2{m}^{2}}-\mathrm{i}\omega {{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right)\right.\right.\hfill \\ \hfill & \quad \left.\left.+\mathrm{i}\omega \left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}+2\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right)\right\}\right]\hfill \\ \hfill & =\mathrm{Tr}\left[\frac{{\boldsymbol{\ell }}^{4}{\mathbf{k}}^{4}}{4{m}^{4}}-\frac{{\boldsymbol{\ell }}^{4}{{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right)}{{m}^{2}}\right.\hfill \\ \hfill & \quad \left.-\frac{{\mathbf{k}}^{4}{{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right)}{{m}^{2}}\right.\hfill \\ \hfill & \quad \left.+2{\omega }^{2}{{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right)\right.\hfill \\ \hfill & \quad \left.-2{\omega }^{2}\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{d}(\mathbf{k})\right]+\tilde{\mathrm{t}}{\mathrm{e}\mathrm{r}\mathrm{m}}_{3},\hfill \\ \hfill & =4\left[\frac{{\boldsymbol{\ell }}^{4}{\mathbf{k}}^{4}}{4{m}^{4}}-{\boldsymbol{\ell }}^{4}\frac{2{d}_{a}^{2}(\mathbf{k})-\vert \mathbf{d}(\mathbf{k}){\vert }^{2}}{{m}^{2}}\right.\hfill \\ \hfill & \quad \left.-{\mathbf{k}}^{4}\frac{2{d}_{a}^{2}(\boldsymbol{\ell })-\vert \mathbf{d}(\boldsymbol{\ell }){\vert }^{2}}{{m}^{2}}+4{\omega }^{2}\left\{{d}_{a}(\mathbf{k}){d}_{a}(\boldsymbol{\ell })\right.\right.\hfill \\ \hfill & \quad \left.\left.-\mathbf{d}(\mathbf{k})\cdot \mathbf{d}(\boldsymbol{\ell })\right\}\right]+\tilde{\mathrm{t}}{\mathrm{e}\mathrm{r}\mathrm{m}}_{3}.\hfill \end{align} \tag{ B15 }$$

where

$$\begin{align}\hfill \tilde{\mathrm{t}}{\mathrm{e}\mathrm{r}\mathrm{m}}_{3}& =4\,\mathrm{Tr}\left[\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{\Gamma }\right)\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right){{\Gamma}}_{a}\left(\mathbf{d}(\mathbf{k})\cdot \mathbf{\Gamma }\right)\right]\hfill \\ \hfill & =16\left[-\frac{{d}_{a}^{2}(\boldsymbol{\ell }){\mathbf{k}}^{4}+{d}_{a}^{2}(\mathbf{k}){\boldsymbol{\ell }}^{4}}{2{m}^{2}}\right.\hfill \\ \hfill & \quad \left.+\frac{{\mathbf{k}}^{4}{\boldsymbol{\ell }}^{4}}{16{m}^{4}}+4{d}_{a}(\boldsymbol{\ell }){d}_{a}(\mathbf{k})\left\{\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{d}(\mathbf{k})\right\}\right].\hfill \end{align} \tag{ B16 }$$

Hence, using the symmetry of the integral under the exchange k ↔ ℓ , we can rewrite the integrand as

$$\begin{align}\hfill {\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}}_{3}& =\frac{{\boldsymbol{\ell }}^{4}{\mathbf{k}}^{4}}{{m}^{4}}-2{\boldsymbol{\ell }}^{4}\frac{8{d}_{a}^{2}(\mathbf{k})-\frac{{\mathbf{k}}^{4}}{{m}^{2}}}{{m}^{2}}\hfill \\ \hfill & \quad +16{\omega }^{2}\left\{{d}_{a}(\mathbf{k}){d}_{a}(\boldsymbol{\ell })-\mathbf{d}(\mathbf{k})\cdot \mathbf{d}(\boldsymbol{\ell })\right\}\hfill \\ \hfill & \quad -\frac{16{d}_{a}^{2}(\mathbf{k}){\boldsymbol{\ell }}^{4}}{{m}^{2}}+\frac{{\boldsymbol{\ell }}^{4}{\mathbf{k}}^{4}}{{m}^{4}}\hfill \\ \hfill & \quad +64{d}_{a}(\boldsymbol{\ell }){d}_{a}(\mathbf{k})\left\{\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{d}(\mathbf{k})\right\}\hfill \\ \hfill & =\frac{2{\boldsymbol{\ell }}^{4}{\mathbf{k}}^{4}}{{m}^{4}}-2{\boldsymbol{\ell }}^{4}\frac{16{d}_{a}^{2}(\mathbf{k})-\frac{{\mathbf{k}}^{4}}{{m}^{2}}}{{m}^{2}}\hfill \\ \hfill & \quad +16{\omega }^{2}\left\{{d}_{a}(\mathbf{k}){d}_{a}(\boldsymbol{\ell })-\mathbf{d}(\mathbf{k})\cdot \mathbf{d}(\boldsymbol{\ell })\right\}\hfill \\ \hfill & \quad +64{d}_{a}(\boldsymbol{\ell }){d}_{a}(\mathbf{k})\left\{\mathbf{d}(\boldsymbol{\ell })\cdot \mathbf{d}(\mathbf{k})\right\}.\hfill \end{align} \tag{ B17 }$$

Finally, performing the integrals, we get:

$$\begin{equation}{\langle {\rho }_{a}{\rho }_{b}\rangle }_{2\;\text{loop}}^{(2)}(\mathrm{i}\omega )=-\frac{7{e}^{2}{m}^{3-\frac{\varepsilon }{2}}\vert \omega {\vert }^{1-\frac{\varepsilon }{2}}}{600{\pi }^{3}c}\left[\frac{1}{\varepsilon }-\frac{1}{2}\,\mathrm{ln}\left(\frac{m\vert \omega \vert }{{{\Lambda}}^{2}}\right)\right]{\delta }_{ab}.\end{equation} \tag{ B18 }$$