can be performed straightforwardly. The integration over time yields $\delta$ functions and the spatial integrations are of Gaussian type. As the resulting expressions are not very elucidating we do not reproduce them here.

For the following discussion it is convenient to switch to spherical momentum coordinates $\vec{k}={\rm k}\hat{\vec{k}}$, where ${\rm k}=\sqrt{k_{\rm x}^2+k_{\rm y}^2+k_{\rm z}^2}$ and $\hat{\vec{k}}=(\cos\varphi\sin\vartheta,\sin\varphi\sin\vartheta,\cos\vartheta)$, with $\varphi\in[0\ldots2\pi)$ and $\vartheta\in[0\ldots\pi]$. The orthogonal vectors to $\vec{k}$ can then be parameterized by a single angle $\beta\in[0\ldots2\pi)$, $$\vec{e}_{\perp,\beta}=
\left(\begin{array}{c}
  \cos\varphi\cos\vartheta\cos\beta-\sin\varphi\sin\beta \\
  \sin\varphi\cos\vartheta\cos\beta+\cos\varphi\sin\beta \\
  -\sin\vartheta\cos\beta
 \end{array}\right) . \label{eq:e_perpbeta}$$ These vectors live in the tangent space of the unit sphere. Correspondingly, the two transverse polarization modes of photons with wave vector $\vec{k}$ can be spanned by two orthonormalized four-vectors $\epsilon^\mu_{(p)}(\vec{k})$, with $p\in\{1,2\}$, $$\epsilon^\mu_{(1)}(\vec{k})=(0,\vec{e}_{\perp,\beta}) \quad\textrm{and}\quad \epsilon^\mu_{(2)}(\vec{k})=\epsilon^\mu_{(1)}(\vec{k})\big|_{\beta\to\beta+\frac{\pi}{2}}\,, \label{eq:epsilons}$$ representing linear polarization states in the specific basis characterized by a particular choice of $\beta$. In this work we exclusively focus on linear polarization modes. Polarizations other than linear can be obtained through linear combinations of the vectors . Resorting to these definitions, the $10$ and $31$ entries of $\hat f^{\mu\nu}_{(p)}(k)$ entering [Eq. ]{} read $$\begin{gathered}
 \hat f^{10}_{(1)}(k)
  = {\rm k} \bigl(\sin\varphi\sin\beta-\cos\varphi \,u \cos\beta\bigr),
  \quad \hat f^{10}_{(2)}(k)=\hat f^{10}_{(1)}(k)\big|_{\beta\to\beta+\frac{\pi}{2}} , \\
 \hat f^{31}_{(1)}(k)
  = {\rm k}\bigl( \cos\varphi\cos\beta-\sin\varphi \,u \sin\beta \bigr),
  