The static hard-loop gluon propagator to all orders in anisotropy

We calculate the (semi-)static hard-loop self-energy and propagator using the Keldysh formalism in a momentum-space anisotropic quark-gluon plasma. The static retarded, advanced, and Feynman (symmetric) self-energies and propagators are calculated to all orders in the momentum-space anisotropy parameter $\xi$. For the retarded and advanced self-energies/propagators, we present a concise derivation and comparison with previously-obtained results and extend the calculation of the self-energies to next-to-leading order in the gluon energy, $\omega$. For the Feynman self-energy/propagator, we present new results which are accurate to all orders in $\xi$. We compare our exact results with prior expressions for the Feynman self-energy/propagator which were obtained using Taylor-expansions around an isotropic state. We show that, unlike the Taylor-expanded results, the all-orders expression for the Feynman propagator is free from infrared singularities. Finally, we discuss the application of our results to the calculation of the imaginary-part of the heavy-quark potential in an anisotropic quark-gluon plasma.


Introduction
When nuclear matter is heated to extremely high temperatures at low net baryon density, it experiences a crossover from ordinary hadronic matter to a new phase in which quark and gluon quasiparticles are the relevant degrees of freedom. The high-temperature phase is called the quark-gluon plasma (QGP) and in order to study it experimentally, physicists use ultrarelativistic heavy-ion collisions (URHICs) performed using, for example, the Relativistic Heavy Ion Collider at Brookhaven National Laboratory and the Large Hadron Collider at the European Organization for Nuclear Research. Due to the short lifetime of the QGP created in URHICs (τ 15 fm/c), in order to study the behavior of the QGP, external probes, such as ultrafast laser pulses or relativistic particle beams, are infeasible. Instead, one must rely on internal probes which are created as a part of system and are subject to the system's full spatiotemporal evolution. If one is interested in probing the full dynamics of the QGP, however, light hadrons which are formed at late times during hadronization do not provide a complete picture. To complement studies of light hadron production, heavy-ion theorists and experimentalists also study particles that are produced at early times and, ideally, do not interact significantly in the final state medium, e.g. dilepton pairs and heavy quarkonia states. Heavy quarkonia production in URHICs, in particular, can provide information about the initial temperature of the QGP, deviations from equilibrium, etc., due to the fact that bound states of heavy quarks are formed very early in the collision and do not melt until the temperature is a few times the phase transition temperature [1,2]. In order to apply this idea to phenomenology, however, one needs to compute the in-medium spectrum of bound states. This step is traditionally done in the context of a finite-temperature equilibrium quantum field theory and the result of this exercise is then folded together with a model for the spatiotemporal evolution of the QGP. However, in practice this approach is insufficient since, at early times, the QGP possesses a significant degree of momentum-space anisotropy.
To address the deviations from isotropic equilibrium, one could attempt to compute the nonequilibrium corrections to the heavy quark potential "perturbatively" as is done in the framework of viscous hydrodynamics, however, at very early times after the collision or near the transverse/longitudinal edges of the QGP this method will become unreliable since such approaches rely on an implicit assumption that the system is near equilibrium. In order to address this problem, one can reorganize the hydrodynamic expansion around a state that is anisotropic in momentum-space at leading-order. This reorganized expansion has been dubbed the "anisotropic hydrodynamics" framework and it is now being widely applied to heavy-ion phenomenology, see e.g. Refs. [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] and references therein. In this approach one assumes a specific form for the momentum-space anisotropic distributions function for the light quarks and gluons which is of Romatschke-Strickland (RS) form [19,20] or generalized RS form [9,10,14,21]. Since momentum-space anisotropies affect e.g. the binding energies and widths of heavy-quarkonia states, there has been a concerted effort to compute their effect on the heavy-quark potential itself [22][23][24][25][26][27] which, when convoluted with a realistic anisotropic hydrodynamics background evolution, has resulted in phenomenological calculations of the nuclear suppression of bottomonia which are in good agreement with available experimental data [28][29][30][31][32].
One of the shortcomings of all previous calculations of the heavy-quark potential in a momentumspace anisotropic QGP is that, although the real part of the potential was obtained to all orders in the plasma momentum-space anisotropy parameter ξ in Ref. [22], calculations of the imaginary part of the heavy-quark potential have relied on a Taylor expansion around ξ = 0 to linear order in ξ [24]. In practice, the in-medium heavy-quark potential is obtained from the Fourier transform of the sum of the static dressed retarded, advanced, and Feynman propagators, i.e.
where g is the strong coupling constant, C F = (N 2 c − 1)/(2N c ) is the quadratic Casimir in the fundamental representation of SU (N c ). The first two terms involving the retarded and advanced propagators are purely real and, consequently, the imaginary part of the potential comes solely from the Fourier transform of the Feynman propagator. As a result, in order to compute the imaginary part of the potential, one must first compute the dressed Feynman propagator to all orders in ξ. In this paper, we perform this computation, laying the groundwork for the all-orders calculation of the imaginary part of the heavy-quark potential in a momentum-space anisotropic QGP. In the process of the derivation, we repeat some previously obtained results for completeness of the presentation, e.g. we revisit the computation of the exact expressions for the retarded and advanced propagators obtained originally in Ref. [22]. In order to carry out the calculation, we make use of the real-time Keldysh formalism and we work in the high-temperature "hard-loop" limit [33]. As a cross-check, our final results for the static limit of the hard-loop Feynman propagator are compared to previously-obtained Taylor-expanded results.

Setup and notation
In this paper, the metric is taken to be "mostly minus", i.e. in Minkowski space η µν ≡ (1, −1). Lower-case letters denote four-vectors and bold lower-case letters denote three vectors. All Greekletter indices stand for the components of four-vectors while Latin indices indicate spatial components of four-vectors. For any two four-vectors x µ ≡ (x 0 , x) and y µ ≡ (y 0 , y) the inner product is defined as x · y ≡ x µ y µ = x 0 y 0 − x · y. The subscripts 'R', 'A', and 'F' for propagators and self-energies stand for retarded, advanced, and Feynman propagators, respectively. According to Feynman slash notation, for any four-vector x, we have / x ≡ γ · x, with γ µ being Dirac matrices. Herein, the anisotropic distribution function is taken to be of Romatschke-Strickland (RS) form [19,20] where n is a three-dimensional unit vector, λ is a temperature-like scale which corresponds to temperature in the equilibrium limit, ξ is the momentum-space anisotropy parameter which, for −1 < ξ < 0, corresponds to a prolate momentum-space anisotropy and, for ξ > 0, corresponds to an oblate momentum-space anisotropy. Herein, we assume that the chemical potential is zero. f iso is an equilibrium distribution function, i.e. f iso (x) ≡ [exp(x) + a] −1 , which depending on the particles statistics can be Boltzmann (a = 0), Fermi-Dirac f F , (a = 1), or Bose-Einstein f B , (a = −1) distribution functions We will perform the calculation initially in QED and then, in the end, we will generalize our results to QCD.

The real-time formalism for non-equilibrium field theories
In this section, we present the basic formalism used to obtain our results in a concise and selfcontained manner. We will use the real-time Keldysh formalism [25,[34][35][36]. The Keldysh formalism is based on contour Green's functions. For a spinor field ψ and a vector field A µ , we can define the fermionic and bosonic QED Green's functions as the following, respectively where {α, β} ∈ {1, 2, 3, 4} are spinor indices, {µ, ν} ∈ {0, 1, 2, 3} are Lorentz indices. The angle brackets · · · denote the quantum expectation value andT is the time-ordering operator defined asT where Θ is the Heaviside step function. Corresponding to different ways of propagation with respect to the contour, one can define four functions based on the contour propagator: Likewise, for bosons, we have In the relation above,T c andT a are time-ordering and anti-time-ordering operators, respectively, which are defined aŝ The upper and lower signs above correspond to bosonic and fermionic cases, respectively. In practice, one can define four different Green's functions with the following meanings: In the Keldysh formulation it is useful to introduce the following 2 × 2 matrix Since, the system under consideration is assumed to be translationally invariant, the two-point function only depends on x − y. Therefore, we can safely set y = 0 and study the two point function function of a single variable (x). The components of the electron propagator S satisfy the following relations The retarded, advanced, and Feynman propagators are defined as Using the above relations one finds some useful identities such as 14) The one-loop photon self-energy is defined as The retarded, advanced, and Feynman self-energies are defined as following with Substituting (3.20) into (3.18), the retarded/advanced self-energy can be written as Using (3.13) and (3.14) we have Thus, Eq. (3.21) gives where in the last line we have used (3.15). Performing the Fourier transform of both sides we have with q ≡ k − p. The Feynman self-energy can be obtained by substituting (3.20) in (3.19) and then using (3.22) After performing the Fourier transform of both sides one has We will calculate the hard loop self-energies and propagators using the Keldysh real-time formalism in the limit of vanishing chemical potential [33]. The "bare" propagators are then 2 × 2 matrices such as for a massless Dirac field and for a massless scalar field. In the above relations, is a small positive number which is sent to zero only at the end of calculation. It should be noted that since Eqs.
with analogous expressions holding for the fermionic propagators. In momentum space, the explicit expressions for the bare propagators as a function of a general momentum k are for fermions and , , for scalar bosons. In the real-time formalism, the following relations hold for the self energies: and Note that, for vector fields, one must add the appropriate Lorentz indices to the propagators and self-energies.

The dressed photon propagator
The dressed photon propagator can be obtained from the Dyson-Schwinger equation where, in the Keldysh formalism, both the propagators and self-energies are 2 × 2 matrices, D is the bare photon propagator andD is the dressed propagator. For the dressed retarded propagator, we haveD The dressed Feynman propagator satisfies which, upon using p = (ω, p) and Note that in the relations introduced so far the distribution functions are general. We will specify the precise forms in the forthcoming sections.

Tensor decomposition in a momentum-space anisotropic plasma
Since the photon propagator is a tensor, we must find a suitable tensor basis and extract the corresponding scalar coefficient functions. For anisotropic systems there are more independent projectors than for the standard equilibrium case due to the fact that there is an additional spacelike vector n which defines the anisotropy direction [19]. Here, we use a four-tensor basis which is appropriate for systems with one anisotropy direction [22]. Specifically, we introduce four tensors Here, m µ is the heat-bath four-velocity, which in the local rest frame is given by is the component of m µ orthogonal to p µ . The direction of anisotropy in momentum space is determined by the vector where n is a three-dimensional unit vector. Likewise,ñ µ is the component of n µ orthogonal to p µ . The self-energies and dressed propagators can be expanded in terms of the tensor basis as We remind the reader that, due to the transversality of the photon self-energy p µ Π µν = 0, not all components of Π µν are independent. One has four equations which can be used, for example, to write the timelike rows/columns of the self-energy tensor in terms of the space-like components Using the symmetry of Π µν = Π νµ , one finds ω 2 Π 00 = p i Π ij p j . These relations show that having Π ii and Π xy , Π xz , Π yz (6 components overall), one can obtain all components of Π µν . Restricting our attention to the spatial block of Π µν , Π ij , one can obtain the expansion tensor coefficients of self-energy using the following projections An alternative method for extracting the coefficient functions, which is based on the four-tensor form of Π µν , is presented in Appendix B.
The dressed retarded/advanced propagators satisfy (3.35), which can be solved to givẽ Using the definition of bare propagator and (4.4) one has where ζ is the gauge fixing parameter. One can obtainD by inverting the above relation [22] D µν (4.14) Comparing to Eq. (4.5), the expansion coefficients are Herein, we take the gauge parameter to be zero, ζ = 0, which is allowed since the static limit of the gauge propagator is gauge invariant.

The retarded and advanced gluon self-energies
In this section, we calculate the gluon self-energy in the hard-loop (HL) limit. Since, in the HLlimit, the photon and gluon self-energies are the same up to definition of Debye mass, we start with Eq. (3.24) for photon self-energy. Note that regarding Eq. (3.24) and (3.30), Π R and Π A are complex conjugates of each other and we only need to find one of them. To proceed, we start with the retarded photon self-energy where, specializing to anisotropic distribution function (2.2) in this section, we made the dependence on the anisotropy parameter, ξ, explicit. Using S R,A, , , and Tr[γ µ γ α γ ν γ β ] = 4(η µα η νβ − η µν η αβ + η µβ η αν ), we have where, in going from the first to the second line, we have used the fact that two terms in the integrand are equal under the transformation k → −k + p (= −q). Using (5.2), one has Note that the first term in ∆ F (k) in Eq. (5.2) corresponds to a divergent vacuum contribution. This contribution is subtracted to obtain the in-medium photon self-energy. Now we can take the HL approximation, that is, taking all internal momenta k to be of order λ (hard) and the external momenta p to be of order eλ (soft) and Taylor-expand the integrand around e = 0. At the leading order of the HL approximation, the quantity in square brackets in (5.4) is Note that the terms containing k 2 are effectively zero due to the delta function which enforces the mass shell condition. Substituting this into the integral (5.4), using integrating over k 0 , and finally setting k → −k in negative-energy contribution, one finds At next-to-leading order we have where by substituting into the integral one obtains the retarded photon self-energy in the HL limit Performing the integral over k 0 and again setting k → −k in the negative energy contribution, one has (5.10) One can show that, for on-shell momentum k µ , so, we after integrating by parts one obtains By specializing to the RS form for the anisotropic distribution function, we can simplify Π µν R a bit more. Using the anisotropic Fermi-Dirac distribution (2.2), one has where v µ ≡ k µ /|k| = (1, k/|k|) and the QED Debye mass is defined as Substituting (5.13) and (5.14) into (5.12), one obtains This precisely corresponds to the result obtained in Ref. [19,35] using relativistic kinetic theory.

The retarded and advanced hard-loop gluon self-energies
As mentioned previously, in the hard-loop limit, the photon and gluon self-energies are the same up to the definition of the Debye mass. The effective QCD distribution function includes contributions from quarks, anti-quarks, and gluons, including the degeneracy factors for the number of quark flavors and gluon color-charge states where f g (k) is gluonic distribution function, f q (k) and fq(k) are quarks and anti-quarks distribution functions. Eq. (5.16) is valid for the gluon self-energy provided that QED Debye mass is replaced with its QCD counterpart, defined as In order to calculate the expansion coefficients (4.4) we only need the spatial block of Π µν R which is In order to obtain the tensor expansion coefficients for the retarded self-energy Π µν R (4.4) we can make use of Eqs. (4.8). For this purpose, we need to choose a frame to be able to define the vectors k, p, and n, and subsequently, the tensor basis matrices. The trivial choice is to take n along z-axis and p in the x-z plane making an angle θ n ≡ arctan(p x /p z ) with z-axis. Based on this coordinate we have n = (0, 0, 1); v = (sin θ cos φ, sin θ sin φ, cos θ);p = p x |p| , 0, p z |p| . (5.20) Up to next-to-leading order in ω, we have

(5.27)
The above expressions apply when n = (0, 0, 1) points along the z-axis and p lies in the x − z plane; in the general case, p z and p x should be replaced by p · n and |p − (p · n)n|, respectively. With this, we can write down an expression for the real part of the potential which is valid to all orders in ξ [22,29] V (r, ξ) = −g 2 C F d 3 p (2π) 3 e ip·r − 1 D 00 R (ω = 0, p, ξ) where we have used the fact that, in this frame, A 00 = C 00 = D 00 = 0 and B 00 = p 2 /ω 2 . Unfortunately, from this point forward one must compute this integral numerically except in some special limiting cases [22,29].

The Feynman gluon self-energy
Starting from the relation (3.26) and using S R,F, one obtains the term inside [...] in (6.1) as Also by definition Using the relations above and performing the integral over k 0 , and setting k → −k in the negativeenergy contribution one has whereq ≡ k + p. Taking the HL limit one obtains the Feynman photon self-energy in the HL limit By specializing to the anisotropic distribution function, we can simplify Π µν F a bit more. Using an anisotropic Fermi-Dirac distribution (2.2) and Eq. (5.14) Using this identity, (6.7) becomes

Feynman gluon self-energy
Once more, the expression for the HL Feynman photon self-energy (6.9) can be used to obtain the HL Feynman gluon self-energy provided that the Debye mass is replaced by (5.18). The spatial block of Feynman gluon self-energy is The integral above can be solved analytically as a function of p, ω, and ξ. For this purpose, it is more convenient to takep along z-axis and n in the x−z plane characterized by (θ n , φ n ) = (arctan(p x /p z ), π), where p x ≡ |p − (p · n)n| and p z ≡ p · n are the perpendicular and parallel components of p respect to n . Using this setup one has ; v = (sin θ cos φ, sin θ sin φ, cos θ);p = (0, 0, 1) , Now by defining and taking the integral over θ in (6.10) we have Note that, for easier evaluation, the Feynman self-energy is calculated in a coordinate system that is different than the one used for retarded and advanced self-energies. Provided that we stay in the same frame until we obtain the scalars (which are coordinate independent), i.e. expansion coefficients of Π µν F , we do not need to do any further transformation to translate between the two choices of coordinate system. Therefore, for example, when needed for calculating Π F coefficients, i.e. Eqs. (4.8), the tensor A should be defined in the frame defined by (6.11). Once we have obtained the coefficients through various tensor contractions, we can use the matrices A, B, C, and D in original frame to construct Π µν F . The resulting HL tensor expansion coefficients for Feynman self-energy in the static limit are where ς ≡ ξp 2 x /p 2 , and K and E are complete elliptic integrals of the first and second kind, respectively, defined by In the isotropic case, the relations above simplify to The first two agree with the isotropic results given by Eqs. (20) and (19) of Ref. [33] for Π T F and Π L F , respectively, upon using β F = (ω 2 /p 2 )Π L F [19] and m 2 D is replaced by its QED definition, i.e. Eq. (5.15).

The hard-loop Feynman propagator
In this section, we obtain the static limit (ω → 0) ofD 00 F using Eq. (3.37) Taking the 00 component using (4.1) and (4.5), the first term becomes where considering (4.1) we have used A 00 = C 00 = D 00 = 0 and B 00 = p 2 /ω 2 . Using (5.21), one can write the tensor basis coefficients of Π R/A as where all new coefficients with subscripts '0' and '1' are independent of ω and can be easily read off from (5.21). Also, specializing to the anisotropic distribution function, Eq. (2.2), one has lim ω→0 (1 + 2f B (p)) sgn(ω) ≈ 2λ Using relations (4.15) for β R/A , one finds (1 + 2f B ) sgn(ω)(D 00 R −D 00 A ) at leading order in the static limit lim ω→0 (1 + 2f B (p))sgn(ω)(D 00 Now, we turn to the next term, (D R Π FDA ) 00 which, in the static limit, becomes The last term can be calculated similarly using the setup above and relations listed in App. A From the relations above (7.5) and (7.7), one can see that the first and the last terms of Eq. (7.1) cancel each other leaving the second term which, using the parametrizations (5.24)-(5.27), gives This is our final result for the static Feynman propagator. Expanding our final result in terms of powers of ξ lim ω→0D 00 which is in agreement with earlier results obtained in the small-ξ limit [24,25,37]. Note that, if one expands (7.8) to higher order in ξ, one finds increasingly negative powers of |p| which result in infrared divergences in the corresponding corrections to the imaginary part of the static heavy quark potential. The full result (7.8) is, however, infrared safe.

Pinch singularity
As mentioned in the introduction, the imaginary part of the heavy-quark potential can be obtained from the Fourier transform of the static limit of the Feynman propagator However, (7.8) contains a pinch singularity which is related to the presence of the (chromo-)Weibel instability in momentum-space anisotropic plasmas [19,35]. This pinch singularity causes the imaginary-part of the potential to be ill-defined. To see that this is the case, we point out that (7.8) can be written more compactly using where with M 2 = m 2 α +m 2 β +m 2 γ [19]. One can show that m 2 + is positive for all ξ and angles of propagation; however, m 2 − can be negative for some propagation angles. This is illustrated in Fig. 2 of Ref. [19] and discussed in the surrounding text. As a result, in unstable regions of phase space, p 2 + m 2 − can go to zero. This occurs already in the integral necessary to obtain the real part of the potential (5.28); however, in this case there is only one power of p 2 + m 2 − in the denominator, which results in a simple pole that can be integrated using a principle part prescription, e.g. p 2 + m 2 − → (|p| + |m − |)(|p| − |m − |). In the case of V I , however, the denominator of the integrand contains (p 2 + m 2 − ) 2 , which results in a double pole in the Fourier transform. To see that this is, in fact, a pinch singularity we note that the prefactor of (7.8) which causes the trouble comes from the product of retarded and advanced propagators, ∆ R G ∆ A G . Keeping track of the i 's, one finds two simple poles shifted by ±i which collapse onto the real axis as → 0, forming a double pole.

Conclusions and outlook
In this paper we presented a calculation of the hard-loop resummed retarded, advanced, and Feynman (symmetric) gluon propagators in a momentum-space anisotropic plasma with a single anisotropy direction, n. We used the real-time formalism throughout and, when available, we compared to previously obtained results. Our main new result is an expression for the Feynman gluon propagator which is accurate to all orders in the anisotropy parameter ξ (7.8). Unlike results obtained using Taylor expansion in ξ, (7.8) is infrared finite, however, it possesses a pinch singularity which formally renders the imaginary part of the heavy-quark potential infinite. The existence of this pinch singularity can be traced back to the existence of unstable modes in a momentum-space anisotropic quark-gluon plasma [19,35,36].
A pinch-singularity emerges because one presumes that the collective modes, which are determined through a linearized treatment, apply at all times. In an equilibrium (stable) situation the field amplitudes are bounded (and small) and such a treatment makes some sense. However, in our case, the system is unstable and some subset of the linearized collective modes grow exponentially for all times, which upon taking the static limit (ω → 0 or t → ∞) results in an infinite effect. As a result, in the presence of unstable modes this scheme is ill-defined and it seems necessary to impose an upper time limit for unstable mode growth. At the most conservative, the upper time limit for unstable mode growth would be set by the lifetime of the QGP, however, in practice one finds that plasma instabilities may saturate on a shorter timescale.
In terms of the calculation presented herein, one could attempt to implement the physics of instability saturation or finite plasma lifetime by imposing an infrared cutoff on the frequency ω 0 ∼ max(τ −1 instability , τ −1 QGP ) where τ instability is the expected timescale for the saturation of unstable field growth and τ QGP ∼ 10 fm/c ∼ 1/(20 MeV) is the typical lifetime of the quark-gluon plasma. Detailed simulations of anisotropic non-abelian plasmas in fixed boxes show that unstable exponential growth terminates when the gauge field amplitude reaches the soft scale and the subsequent gauge field dynamics transform into a much slower turbulent cascade of energy from soft scales to hard scales [38][39][40][41][42]. More recent studies of chromo-Weibel dynamics in an expanding non-Abelian plasmas found that unstable modes saturate on a time scale of 3-4 fm/c at LHC energies [43,44]. Combined with the QGP lifetime estimate, one has ω 0 ∼ 20 − 70 MeV. In practice, an infrared cutoff such as this will lift the poles off the real axis, even in the limit → 0 due to the finite imaginary linear correction in ω to the structure functions α R/A , β R/A , γ R/A , and δ R/A listed in (5.21).
While such a phenomenological prescription may work in practice, it introduces a fundamental problem, since the 00-component gluon propagator which is used to define the potential is not gauge invariant for finite ω [see Eq. (4.12)]. For "reasonable gauges" the dependence may not be large, but nevertheless this is an unsatisfactory resolution of this problem on general grounds. For this reason, one should simultaneously pursue the possibility to measure the potential numerically using classical gauge theory simulations similar to those used to measure the imaginary part of the heavy quark potential in the equilibrium limit [45]. In this method, one determines the imaginary part of the potential by measuring the classical Wilson loop which amounts to a two-point correlation function of two spatial Wilson lines. With this method one would be able to obtain a gaugeinvariant imaginary part of the potential, however, one still would not be able to take the t → ∞ limit due to finite computational resources, break-down of the classical hard loop limit, etc.
Finally, as another path forward, one might consider adding the effect of collisions in the computation of the anisotropic structure functions. Previous studies [46] have shown that at fixed ξ, if the collision rate exceeds a certain threshold, then unstable modes are eliminated from the spectrum. This would provide another way to regulate/eliminate the ill-defined effect of unstable modes in the heavy-quark potential. We leave the investigation of these possibilities for future work(s).

A Tensor relations
In this appendix, we present various tensor identities obeyed by our basis tensors (4.1). Useful identities for the contraction of any two basis tensors are as following A · A · C = A · C · C = A · C · A = C · A · A = C · C · C = C · A · C = C · C · A = C , (A.17) where (A.20) and (A.21) contain the contractions that cannot be expressed in terms of any single basis tensors. In all relations above p ⊥ ≡ |p − (p · n)n| is the component of p perpendicular to n.

B An alternative method for calculating the expansion coefficients
In this section, we present an alternative method to Eq. (4.8) for calculating the expansion coefficients. This method, in comparison to the previous one, is based on the four-vector form of tensors A, B, C, and D. Any arbitrary rank-2 tensor T can be expanded in terms of basis tensors as Let's start with calculating the following terms using the identities presented in App. A As an example of usage of the above formulae, one can set T ≡ Π R,A,F to obtain the expansion coefficients for the retarded, advanced, and Feynman self-energies directly from the their four-tensor expressions prior to restricting to the space-like components.