# Chiral gravitational waves from thermalized neutrinos in the early Universe

Philipp Gubler,  Naoki Yamamoto  and Di-Lun Yang

###### Abstract

We investigate polarized gravitational waves generated by chiral fermions in
the early Universe. In particular, we focus on the contribution from left-
handed neutrinos in thermal equilibrium with finite temperature and chemical
potential in the radiation dominated era. We compute the correlation functions
of gravitational fields pertinent to the Stokes parameter $V$ characterizing
the circular polarization of gravitational waves in the Minkowski and
expanding spacetime backgrounds. In the expanding universe, we find that the
thermalized neutrinos induce a non-vanishing $V$ linear to the neutrino
degeneracy parameter and wavenumber of gravitational waves in the long
wavelength region. While the magnitude of the gravitational waves generated by
thermal neutrinos is too small to be detectable by current and planned third
generation gravitational wave detectors, their observations by future
generation detectors for ultra-high frequency regimes could provide
information on the neutrino degeneracy parameter in the early Universe.

##  1 Introduction

With the experimental observation [1] of gravitational waves in the ground
based detectors LIGO [2] and Virgo [3], while KAGRA [4] is starting its
operation, we are now in the era of gravitational wave astronomy. The
gravitational waves observed so far, provide information about violent
astronomical events such as black hole mergers or binary inspiraling neutron
stars, but have also been used for constraining the energy density of the
stochastic gravitational wave background [5]. With increasing precision of the
existing facilities and future ground and space based detectors, it will
furthermore become possible to measure the parity-odd components of
gravitational waves [6, 7, 8, 9] (see also ref. [10] for a recent analysis of
the already available data). In the context of the early Universe, mechanisms
that generate such parity violating gravitational waves have been discussed
for instance in relation to pseudoscalar (or axion) inflation models with the
gravitational topological coupling $\varphi
R_{\mu\nu\alpha\beta}\tilde{R}^{\mu\nu\alpha\beta}$ [11, 12, 13, 14, 15] or
gauge topological coupling $\varphi F_{\mu\nu}\tilde{F}^{\mu\nu}$ [16, 17, 18]
and its non-Abelian version [19, 20, 21], the existence of chiral fermions
[23, 24, 25], the helical magnetic field [26] due to the chiral plasma
instability [27, 28], and helical turbulence [29].

We will in this paper focus on (nearly) massless left-handed neutrinos in the
early Universe that were thermalized after the inflation period.111Note that
thermalized matter can emit gravitational waves due to the thermal
fluctuations of energy-momentum tensors; see, e.g., ref. [30]. This is
analogous to the generation of thermal photons from a neutral plasma [31, 32].
Because these neutrinos have fixed left-handed chirality and violate parity
within the Standard Model, they generate chiral gravitational waves with an
imbalance between the two possible circular polarizations, which can be
characterized by the Stokes parameter $V$. We especially consider the role of
thermal effects, and demonstrate that in the regime where both the neutrino
chemical potential $\mu$ and the momentum $k$ of the gravitational wave are
much smaller than the temperature $T$ of the thermalized neutrino gas
($\mu/T$, $k/T\ll 1$), $V$ depends linearly on both $\mu/T$ and $k/T$; see eq.
(4.44) for our main result. It is also straightforward to extend the results
of the present paper to other types of chiral fermions or Dirac fermions with
finite chiral chemical potential in thermal equilibrium. A deeper
understanding of these chiral gravitational waves would also provide a
baseline to study the possible parity-violating effects beyond the Standard
Model.

This paper is organized as follows. In section 2, we briefly review the Stokes
parameters for electromagnetic and gravitational waves. In section 3, we
introduce a general framework for the production of gravitational waves from
fermionic sources and obtain the correlations of gravitational fluctuations
from chiral fermions in thermal equilibrium. In section 4, we adopt the
developed formalism to investigate the birefringence of gravitational waves
coming from thermalized left-handed neutrinos in both the Minkowski and
expanding universes, where the Stokes parameters are evaluated and the
relation to experimental measurements is discussed. Finally, a brief summary
and outlook are presented in section 5. The details of calculations are shown
thereafter in multiple appendices.

Throughout the paper, we use the natural units $\hbar=c=1$ and mostly minus
metric $\eta_{\mu\nu}={\rm diag}(1,-1,-1,-1)$. We also employ the shorthand
notations $A_{(\mu}B_{\nu)}\equiv(A_{\mu}B_{\nu}+A_{\nu}B_{\mu})/2$ and
$\int_{\bm{p}}\equiv\int\frac{{\rm d}^{3}p}{(2\pi)^{3}}$.

##  2 The Stokes parameters for electromagnetic and gravitational waves

We here briefly review the Stokes parameters, which were initially proposed
for characterizing the polarization of electromagnetic waves [33, 34], but can
also be used to describe the polarization of gravitational waves [35].

First, let us for illustration consider an electromagnetic wave moving in the
direction of the z-axis, partially following ref. [36]. It suffices to
consider only the properties of the electric field $\bm{E}$, as the magnetic
field will be simply pointing to the direction perpendicular to the direction
of the wave and $\bm{E}$. The electric field can be generally parametrized as

$\displaystyle\bm{E}=\begin{pmatrix}E_{x}\\\ E_{y}\\\
0\end{pmatrix}=\begin{pmatrix}A_{x}\cos(z-t)\\\ A_{y}\cos(z-t+\beta)\\\
0\end{pmatrix}$ (2.1)

Eliminating the time variable, $E_{x}$ and $E_{y}$ can be shown to satisfy the
equation

$\displaystyle\frac{E_{x}^{2}}{A_{x}^{2}}+\frac{E_{y}^{2}}{A_{y}^{2}}-2\frac{E_{x}E_{y}}{A_{x}A_{y}}\cos\beta=\sin^{2}\beta,$
(2.2)

which corresponds to a rotated ellipse in the $E_{x}$-$E_{y}$ plane. Defining
$A_{x}=A\cos\theta$ and $A_{y}=A\sin\theta$, the rotation angle $\alpha$ of
the ellipse is determined from

$\displaystyle\tan 2\alpha=\tan 2\theta\cos\beta,$ (2.3)

and the semi-major and semi-minor axes $a$ and $b$ are given as

$\displaystyle a^{2}$
$\displaystyle=\frac{A^{2}}{2}\Bigl{(}1+\sqrt{1-\sin^{2}2\theta\sin^{2}\beta}\Bigr{)},$
(2.4) $\displaystyle b^{2}$
$\displaystyle=\frac{A^{2}}{2}\Bigl{(}1-\sqrt{1-\sin^{2}2\theta\sin^{2}\beta}\Bigr{)}.$
(2.5)

For $\beta=0$, the wave is linearly polarized, while it is circularly
polarized for $\beta=\pi/2$ and $\theta=\pi/4$.

The nature of the polarization can be expressed using the Stokes parameters
$I$, $Q$, $U$ and $V$, which are defined as

$\displaystyle I$ $\displaystyle\equiv a^{2}+b^{2}=A^{2},$ (2.6)
$\displaystyle Q$ $\displaystyle\equiv(a^{2}-b^{2})\cos 2\alpha=A^{2}\cos
2\theta,$ (2.7) $\displaystyle U$ $\displaystyle\equiv(a^{2}-b^{2})\sin
2\alpha=A^{2}\sin 2\theta\cos\beta,$ (2.8) $\displaystyle V$
$\displaystyle\equiv 2abh=A^{2}h\sin 2\theta\sin\beta,$ (2.9)

where $h=\pm 1$ parametrizes the direction of rotation for the circularly
polarized component of the wave. It is clear from the above definitions that
$I$ gives the overall magnitude of the wave, $V$ the degree and direction of
its circular polarization, and $Q$ and $U$ the degree and direction of the
linear polarization of the wave. The four Stokes parameters are not
independent, but are related to each other via

$\displaystyle Q^{2}+U^{2}+V^{2}=I^{2}.$ (2.10)

Using the complex representation of an electric field,

$\displaystyle\bm{E}=\begin{pmatrix}E_{x}\\\ E_{y}\\\ 0\end{pmatrix}{\rm
e}^{{\rm i}(z-t)},$ (2.11)

where $E_{x}$ and $E_{y}$ are complex and contain relative phases between the
$x$ and $y$ component of the field, the Stokes parameters can be expressed as

$\displaystyle I$ $\displaystyle\equiv|E_{x}|^{2}+|E_{y}|^{2},$ (2.12)
$\displaystyle Q$ $\displaystyle\equiv|E_{x}|^{2}-|E_{y}|^{2},$ (2.13)
$\displaystyle U$ $\displaystyle\equiv 2\mathrm{Re}(E_{x}E_{y}^{\ast}),$
(2.14) $\displaystyle V$
$\displaystyle\equiv-2\mathrm{Im}(E_{x}E_{y}^{\ast}),$ (2.15)

which can readily be confirmed to satisfy eq. (2.10) and to be consistent with
the definitions of eqs. (2.6)–(2.9) if one considers $\beta$ to be the
relative phase of $E_{x}$ and $E_{y}$ and $\theta$ an angle determining the
relative size of $|E_{x}|$ and $|E_{y}|$. Note that for a non-purely polarized
wave, the above definitions should be regarded as expectation values over time
of different components of $E_{x}$, $E_{y}$ and their complex conjugates. This
can be cast in a concise form using the matrix

$\displaystyle\langle
E_{i}E_{j}^{\ast}\rangle=\frac{1}{2}\begin{pmatrix}I+Q&U-{\rm i}V\\\ U+{\rm
i}V&I-Q\end{pmatrix},$ (2.16)

with $i$, $j=x$, $y$.

In analogy to eq. (2.16), the Stokes parameters of gravitational waves can be
defined as [6]

$\displaystyle\langle
h_{i}(\bm{k})h_{j}^{\ast}(\bm{k}^{\prime})\rangle=\frac{1}{2}\delta(\bm{k}-\bm{k}^{\prime})\begin{pmatrix}I(\bm{k})+Q(\bm{k})&U(\bm{k})-{\rm
i}V(\bm{k})\\\ U(\bm{k})+{\rm i}V(\bm{k})&I(\bm{k})-Q(\bm{k})\end{pmatrix},$
(2.17)

where the subscripts $i$ and $j$ stand for the $+$ and $\times$ polarization
modes of the gravitational wave in the transverse-traceless gauge. The $+$ and
$\times$ components can be rewritten using right- and left-handed modes as

$\displaystyle h_{\rm R}$ $\displaystyle=\frac{1}{\sqrt{2}}(h_{+}-{\rm
i}h_{\times}),$ (2.18) $\displaystyle h_{\rm L}$
$\displaystyle=\frac{1}{\sqrt{2}}(h_{+}+{\rm i}h_{\times}),$ (2.19)

which modifies the above matrix to

$\displaystyle\langle
h_{\lambda}(\bm{k})h_{\lambda^{\prime}}^{\ast}(\bm{k}^{\prime})\rangle=\frac{1}{2}\delta(\bm{k}-\bm{k}^{\prime})\begin{pmatrix}I(\bm{k})+V(\bm{k})&Q(\bm{k})-{\rm
i}U(\bm{k})\\\ Q(\bm{k})+{\rm i}U(\bm{k})&I(\bm{k})-V(\bm{k})\end{pmatrix},$
(2.20)

with $\lambda$, $\lambda^{\prime}={\rm R}$, ${\rm L}$. In particular, the $V$
parameter is explicitly given as the circular polarization of the
gravitational wave:

$\displaystyle\langle h_{\rm R}(\bm{k},t)h_{\rm
R}^{\ast}(\bm{k}^{\prime},t)\rangle-\langle h_{\rm L}(\bm{k},t)h_{\rm
L}^{\ast}(\bm{k}^{\prime},t)\rangle$ $\displaystyle=$
$\displaystyle\delta(\bm{k}-\bm{k}^{\prime})V(\bm{k}).$ (2.21)

Starting from the next section, we will for convenience of notation employ
$\lambda$, $\lambda^{\prime}=\pm 1$ instead of $\lambda$,
$\lambda^{\prime}={\rm R}$, ${\rm L}$. Especially, $\lambda$,
$\lambda^{\prime}=+1$ should not be confused with the $+$ polarization
mentioned above. As we are mostly interested in parity violation components of
the gravitational wave spectrum, our focus will be laid on the Stokes
parameter $V$, which parametrizes the difference between the magnitudes of the
right- and left-handed spectrum.

##  3 Production of gravitational waves from thermal chiral fermions

In this section, we will demonstrate how to derive the two-point functions of
gravitational fields induced by chiral fermions in the helicity basis. We will
first review the general framework in subsection 3.1 and the application to
the vacuum polarization in subsection 3.2. Then, we generalize the results to
compute thermal corrections in subsection 3.3.

###  3.1 General framework

In setting up the general framework to study gravitational waves of differing
helicities emitted from left-handed fermionic sources, we partly follow ref.
[23], but generalize it to take into account effects of finite temperature and
chemical potential. The key objects that we wish to calculate are helical mode
correlators of gravitational perturbations generated by left-handed neutrinos.

We start with the equation of motion (EOM) for gravitational perturbations
$h_{\mu\nu}$ in the Friedmann-Lemaître-Robertson-Walker (FLRW) background in
the conformal coordinates, ${\rm d}s^{2}=a(\tau)^{2}({\rm d}\tau^{2}-{\rm
d}\bm{x}^{2})$. Here, $\tau$ is the conformal time, which is related to the
cosmic time $t$ via ${\rm d}t=a{\rm d}\tau$, and $a(\tau)$ is the scale
factor. (Later, we will also consider the Minkowski universe, which
corresponds to setting $\tau=t$ and $a(\tau)=1$.) We adopt the transverse-
traceless (TT) gauge such that
$h_{00}=h_{0i}=h^{i}_{i}=\partial_{j}h^{j}_{i}=0$, where $i,j$ denote the
spatial coordinates. The EOM now reads222Here we neglect the damping term due
to free-streaming neutrinos [22], which may reduce the correlation of
gravitational waves by several $10\%$ while this suppression does not change
the order-of-magnitude estimate in our study.

$\displaystyle\big{(}\partial_{\tau}^{2}+2\mathcal{H}\partial_{\tau}+k^{2}\big{)}h_{\mu\nu}({\bm{k}},\tau)=\frac{2}{m_{\rm
P}^{2}}T_{\mu\nu}({\bm{k}},\tau),$ (3.1)

where ${\cal H}\equiv({\rm d}a/{\rm d}\tau)/a$ is the conformal Hubble
parameter, $T_{\mu\nu}$ is the physical energy-momentum tensor (whose explicit
form will be specified below), $m_{\rm P}=1/\sqrt{8\pi G}$ is the reduced
Planck mass, and $k=|{\bm{k}}|$, $\bm{k}$ being the three momentum.

In the TT gauge, $h_{\mu\nu}$ can be decomposed into circular helicity modes
$h_{\pm}$ as

$\displaystyle
h^{\mu\nu}({\bm{k}},\tau)=\sum_{\lambda=\pm}\epsilon_{\lambda}^{\mu}({\bm{k}})\epsilon_{\lambda}^{\nu}({\bm{k}})h_{\lambda}({\bm{k}},\tau),$
(3.2)

where $\epsilon^{\mu}_{\pm}$ denote the polarization four vectors with $\pm$
helicity satisfying

$\displaystyle\epsilon_{\lambda}^{\mu~{}*}(\bm{k})$ $\displaystyle=$
$\displaystyle\epsilon_{-\lambda}^{\mu}(\bm{k})\,,\quad\epsilon_{\lambda}^{\mu}(-\bm{k})=-\epsilon_{-\lambda}^{\mu}(\bm{k})\,,\quad
k\cdot\epsilon_{\lambda}(\bm{k})=0\,,\quad$
$\displaystyle\epsilon_{\lambda}(\bm{k})\cdot\epsilon_{\lambda^{\prime}}^{*}(\bm{k})$
$\displaystyle=$
$\displaystyle-\delta_{\lambda,\lambda^{\prime}},\quad\epsilon_{\lambda}(\bm{k})\cdot\epsilon_{\lambda^{\prime}}(\bm{k})=-\delta_{\lambda,-\lambda^{\prime}}.$
(3.3)

Given the decomposition in eq. (3.2), the EOM for $h_{\lambda}$ yields

$\displaystyle\big{(}\partial_{\tau}^{2}+2\mathcal{H}\partial_{\tau}+k^{2}\big{)}h_{\lambda}({\bm{k}},\tau)=\frac{2}{m_{\rm
P}^{2}}\epsilon^{\alpha\beta}_{-\lambda}(\bm{k})\Pi^{\mu\nu}_{\quad,\alpha\beta}(\bm{k})T_{\mu\nu}(\bm{k},\tau),$
(3.4)

where

$\displaystyle\epsilon^{\alpha\beta}_{\lambda}(\bm{k})$ $\displaystyle\equiv$
$\displaystyle\epsilon_{\lambda}^{\alpha}(\bm{k})\epsilon_{\lambda}^{\beta}(\bm{k}),\quad\Pi_{\mu\nu,\alpha\beta}(\bm{k})\equiv{\cal
P}_{\mu\alpha}(\bm{k}){\cal P}_{\nu\beta}(\bm{k})-\frac{1}{2}{\cal
P}_{\mu\nu}(\bm{k}){\cal P}_{\alpha\beta}(\bm{k}),\quad$ $\displaystyle{\cal
P}_{\mu\nu}(\bm{k})$ $\displaystyle\equiv$
$\displaystyle\delta_{\mu\nu}-\frac{k_{\mu}k_{\nu}}{k^{2}}.$ (3.5)

The solution for $h_{\lambda}(\bm{k},\tau)$ from the EOM of eq. (3.4) can be
written as

$\displaystyle
h_{\lambda}(\bm{k},\tau)=h_{\lambda}^{\text{hom}}(\bm{k},\tau)-\frac{2}{m_{\rm
P}^{2}}\epsilon^{\mu\nu}_{-\lambda}(\bm{k})\int{\rm
d}\tau^{\prime}G_{k}(\tau,\tau^{\prime})T_{\mu\nu}(\bm{k},\tau^{\prime}),$
(3.6)

where
$\epsilon^{\alpha\beta}_{\lambda}(\bm{k})\Pi^{\mu\nu}_{\quad,\alpha\beta}(\bm{k})=\epsilon^{\mu\nu}_{\lambda}(\bm{k})$
is used, $h_{\lambda}^{\text{hom}}(\bm{k},t)$ is the homogeneous solution and
$G_{k}(\tau,\tau^{\prime})$ is the retarded Green’s function satisfying

$\displaystyle\big{(}\partial_{\tau}^{2}+2\mathcal{H}\partial_{\tau}+k^{2}\big{)}G_{k}(\tau,\tau^{\prime})=\delta(\tau-\tau^{\prime}).$
(3.7)

As $h_{\lambda}^{\text{hom}}(\bm{k},\tau)$ is independent of $\lambda$ [23],
we can focus on just the inhomogeneous solution. Accordingly, the correlation
function of the helical modes of the gravitational waves reads

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle=\frac{4}{m_{\rm
P}^{4}}\int{\rm d}\tau^{\prime}{\rm
d}\tau^{\prime\prime}G_{k}(\tau,\tau^{\prime})G_{k^{\prime}}(\tau,\tau^{\prime\prime})M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime}),$
(3.8)

where

$\displaystyle
M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})=\epsilon^{\mu\nu}_{-\lambda}(\bm{k})\epsilon^{\alpha\beta}_{-\lambda^{\prime}}(\bm{k}^{\prime})\langle
T_{\mu\nu}(\bm{k},\tau^{\prime})T_{\alpha\beta}(\bm{k}^{\prime},\tau^{\prime\prime})\rangle.$
(3.9)

For the symmetric energy-momentum tensor, we consider the contribution of
left-handed fermions (neutrinos),

$\displaystyle T_{\mu\nu}(x)=\frac{1}{a^{2}}\left(\frac{\rm
i}{2}\big{[}\psi^{\dagger}(x)\bar{\sigma}_{(\mu}\partial_{\nu)}\psi(x)-\partial_{(\mu}\psi^{\dagger}(x)\bar{\sigma}_{\nu)}\psi(x)\big{]}-{\rm
i}\eta_{\mu\nu}\psi^{\dagger}(x)\bar{\sigma}^{\lambda}\partial_{\lambda}\psi(x)\right)\,,$
(3.10)

where $\sigma^{\mu}=(1,{\bm{\sigma}})$,
$\bar{\sigma}^{\mu}=(1,-{\bm{\sigma}})$, $\sigma^{i}$ being Pauli matrices,
$\psi$ the rescaled fermionic field, and $x^{\mu}=(\tau,\bm{x})$. Note that
the term proportional to $\eta_{\mu\nu}$ will be irrelevant when we consider
the transverse part of $T_{\mu\nu}$ in eq. (3.4) in momentum space. The mode
expansion of the $\psi$ field takes the form

$\displaystyle\psi(x)=\int_{\bm{p}}\left[\xi_{-}(\bm{p})u(\tau,\bm{p})a_{\bm{p}}+\xi_{-}(-\bm{p})v(\tau,-\bm{p})b^{\dagger}_{-\bm{p}}\right]{\rm
e}^{{\rm i}\bm{p}\cdot\bm{x}},$ (3.11)

where $\xi_{-}(\bm{p})$ is the negative eigenstate of the helicity operator
such that

$\displaystyle{\bm{\sigma}}\cdot{\bm{p}}\xi_{-}(\bm{p})$ $\displaystyle=$
$\displaystyle-p\xi_{-}(\bm{p}),\quad\xi_{-}^{\dagger}(\bm{p})\xi_{-}(\bm{p})={\bm{1}},$
$\displaystyle\xi_{-}(\bm{p})\xi_{-}^{\dagger}(\bm{p})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sigma\cdot n(\bm{p}),\quad n^{\mu}=(1,{\bm{p}}/p),$
(3.12)

and $a_{\bm{p}}$ and $b^{\dagger}_{\bm{p}}$ denote the annihilation operator
for fermions and creation operator for anti-fermions, respectively. They
follow the fermionic commutation relations

$\displaystyle\\{a_{\bm{q}},a^{\dagger}_{\bm{p}}\\}=(2\pi)^{3}\delta^{3}({\bm{q}}-{\bm{p}}),\quad\\{b_{\bm{q}},b^{\dagger}_{\bm{p}}\\}=(2\pi)^{3}\delta^{3}({\bm{q}}-{\bm{p}}),\quad\\{a_{\bm{q}},a_{\bm{p}}\\}=\\{b_{\bm{q}},b_{\bm{p}}\\}=0.$
(3.13)

Using the mode expansion in eq. (3.11), we find

$\displaystyle T_{\mu\nu}(\bm{k},\tau)=\frac{1}{2a(\tau)^{2}}\int{\rm
d}^{3}x\int_{\bm{p}}\int_{\bm{p}^{\prime}}L^{\dagger}(\bm{p}^{\prime},\tau)\bar{\sigma}_{(\mu}\big{(}p_{\nu)}+p^{\prime}_{\nu)}\big{)}L(\bm{p},\tau){\rm
e}^{{\rm i}({\bm{p}-\bm{p}^{\prime}-\bm{k}})\cdot\bm{x}},$ (3.14)

where

$\displaystyle
L(\bm{p},\tau)\equiv\xi_{-}(\bm{p})u(\tau,\bm{p})a_{\bm{p}}+\xi_{-}(-\bm{p})v(\tau,-\bm{p})b^{\dagger}_{-\bm{p}}.$
(3.15)

As shown in ref. [23], by inserting eq. (3.14) into eq. (3.8), one obtains

$\displaystyle
M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})$
$\displaystyle=$
$\displaystyle\frac{1}{4a(\tau^{\prime})^{2}a(\tau^{\prime\prime})^{2}}\int_{\bm{p}}\int_{\bm{p}^{\prime}}\Big{\langle}L^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})\epsilon^{\mu\nu}_{-\lambda}(\bm{k})$
(3.16)
$\displaystyle\times\big{[}\bar{\sigma}_{\mu}p_{\nu}+\bar{\sigma}_{\mu}(p-k)_{\nu}\big{]}L(\bm{p},\tau^{\prime})L^{\dagger}({\bm{p}^{\prime}-\bm{k}^{\prime}},\tau^{\prime\prime})\epsilon^{\rho\sigma}_{-\lambda^{\prime}}(\bm{k}^{\prime})$
$\displaystyle\times\big{[}\bar{\sigma}_{\rho}p^{\prime}_{\sigma}+\bar{\sigma}_{\rho}(p^{\prime}-k^{\prime})_{\sigma}\big{]}L(\bm{p}^{\prime},\tau^{\prime\prime})\Big{\rangle}.$

###  3.2 Polarization in vacuum

When considering just the vacuum expectation value, the only non-vanishing
combination of the creation and annihilation operators (besides the trivial
$\bm{k}=\bm{k}^{\prime}=\bm{0}$ case, which we ignore here) is

$\displaystyle\langle
0|b_{\bm{k}-\bm{p}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}|0\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}({\bm{p}-\bm{p}^{\prime}+\bm{k}^{\prime}})\delta^{3}({\bm{p}^{\prime}-\bm{p}+\bm{k}})$
(3.17) $\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{p}^{\prime}-\bm{p}+\bm{k}}).$

We hence have

$\displaystyle
M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})$
$\displaystyle=(2\pi)^{3}\frac{\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{a(\tau^{\prime})^{2}a(\tau^{\prime\prime})^{2}}\int_{\bm{p}}v^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})u(\bm{p},\tau^{\prime})u^{\dagger}({\bm{p}^{\prime}-\bm{k}^{\prime}},\tau^{\prime\prime})v(\bm{p}^{\prime},\tau^{\prime\prime})\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})\,,$
(3.18)

where

$\displaystyle\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}\epsilon^{\mu}_{-\lambda}(\bm{k})\epsilon^{\rho}_{\lambda^{\prime}}(\bm{k})\xi^{\dagger}_{-}({\bm{k}-\bm{p}})\bar{\sigma}_{\mu}\xi_{-}(\bm{p})\xi^{\dagger}_{-}({\bm{p}})\bar{\sigma}_{\rho}\xi_{-}({\bm{k}-\bm{p}}).$
(3.19)

Making use of eq. (3.1), $\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
can be further rewritten as

$\displaystyle\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{1}{4}\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}\epsilon^{\mu}_{-\lambda}(\bm{k})\epsilon^{\rho}_{\lambda^{\prime}}(\bm{k})n^{\alpha}({\bm{k}-\bm{p}})n^{\beta}(\bm{p})\mathrm{tr}\,\big{(}\sigma_{\alpha}\bar{\sigma}_{\mu}\sigma_{\beta}\bar{\sigma}_{\rho}\big{)}$
(3.20) $\displaystyle=$
$\displaystyle-\frac{1}{2}\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}\Bigg{[}\frac{2\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}}{p|{\bm{k}-\bm{p}}|}-\delta_{\lambda,\lambda^{\prime}}\left(1+\frac{p^{2}-{\bm{k}\cdot\bm{p}}}{p|{\bm{k}-\bm{p}}|}\right)$
$\displaystyle+{\rm
i}\epsilon_{\mu\beta\rho\alpha}\epsilon^{\mu}_{-\lambda}(\bm{k})\epsilon^{\rho}_{\lambda^{\prime}}(\bm{k})n^{\alpha}({\bm{k}-\bm{p}})n^{\beta}(\bm{p})\Bigg{]},$

where we utilized

$\displaystyle\mathrm{tr}\,\big{(}\sigma_{\alpha}\bar{\sigma}_{\mu}\sigma_{\beta}\bar{\sigma}_{\rho}\big{)}=2\big{(}\eta_{\mu\beta}\eta_{\rho\alpha}-\eta_{\mu\rho}\eta_{\beta\alpha}+\eta_{\mu\alpha}\eta_{\beta\rho}-{\rm
i}\epsilon_{\mu\beta\rho\alpha}\big{)}$ (3.21)

with $\epsilon_{0123}=-1$. The last (imaginary) term in eq. (3.20) is critical
to break parity.

To integrate over $\bm{p}$, one has to explicitly fix the coordinate system.
Following ref. [23], we take $\bm{k}$ along $z$ (or equivalently $x^{3}$)
direction and explicitly define the polarization vectors as

$\displaystyle\epsilon^{\mu}_{\lambda}(\bm{k})=\frac{1}{\sqrt{2}}(0,\mathop{\mathrm{sgn}}(k^{z}),{\rm
i}\lambda,0),$ (3.22)

where $\mathop{\mathrm{sgn}}(k^{z})$ is introduced to satisfy the relation
$\epsilon_{\lambda}^{\mu}(-\bm{k})=-\epsilon_{-\lambda}^{\mu}(\bm{k})$ in eq.
(3.1). We will focus on $k^{z}>0$ such that $k^{\mu}=(k,0,0,k)$ without loss
of generality. $\bm{p}$ is parametrized as
$\bm{p}=p(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\phi)$. Now,
$\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ in spherical coordinates
can be given as

$\displaystyle\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})=$
$\displaystyle\frac{p^{2}}{4}\sin^{2}\theta\big{(}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{)}\Bigg{[}\delta_{\lambda,\lambda^{\prime}}-\frac{\lambda+\lambda^{\prime}}{2}\cos\theta$
(3.23)
$\displaystyle+\frac{1}{|{\bm{k}-\bm{p}}|}\bigg{(}\delta_{\lambda,\lambda^{\prime}}\big{(}p-k\cos\theta\big{)}+\frac{\lambda+\lambda^{\prime}}{2}\big{(}k-p\cos\theta\big{)}$
$\displaystyle-p\sin^{2}\theta\big{(}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{)}\bigg{)}\Bigg{]},$

which yields333Here, $\beta_{\lambda\lambda}(\bm{k},\bm{p})$ is different from
the result found in ref. [23] with opposite signs in front of $\lambda$.

$\displaystyle\beta_{\lambda\lambda}({\bm{k},\bm{p}})$
$\displaystyle=\frac{p^{2}\sin^{2}\theta}{4}\Bigg{[}1-\lambda\cos\theta+(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}\,,$
(3.24) $\displaystyle\beta_{\lambda,-\lambda}({\bm{k},\bm{p}})$
$\displaystyle=-\frac{p^{3}\sin^{4}\theta\big{(}\cos 4\phi-{\rm i}\lambda\sin
4\phi\big{)}}{4|{\bm{k}-\bm{p}}|}\,,$ (3.25)

where we used $\lambda^{2}=1$. Note that
$\beta_{\lambda,-\lambda}({\bm{k},\bm{p}})$, which is the only term depending
on $\phi$ in the integrand of eq. (LABEL:M_in_beta), does not contribute to
$M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})$
by symmetry. Now, inserting eq. (3.24) into eq. (LABEL:M_in_beta), the
correlation function in eq. (3.8) becomes

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle$
$\displaystyle=\delta_{\lambda,\lambda^{\prime}}\frac{\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{m_{\rm
P}^{4}}\int{\rm d}\tau^{\prime}{\rm
d}\tau^{\prime\prime}\frac{G_{k}(\tau,\tau^{\prime})G_{k^{\prime}}(\tau,\tau^{\prime\prime})}{a(\tau^{\prime})^{2}a(\tau^{\prime\prime})^{2}}\int{\rm
d}^{3}pv^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})u(\bm{p},\tau^{\prime})$
$\displaystyle\quad\times
u^{\dagger}({\bm{p}^{\prime}-\bm{k}^{\prime}},\tau^{\prime\prime})v(\bm{p}^{\prime},\tau^{\prime\prime})p^{2}\sin^{2}\theta\Bigg{[}1-\lambda\cos\theta+(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}.$
(3.26)

###  3.3 Thermal corrections

So far, we have closely followed ref. [23] in our derivation of the
correlation function. The purpose of this work, however, is to study not the
effect of vacuum fluctuations of left-handed neutrinos, but to evaluate their
thermal contributions to the generation of gravitational waves in the early
Universe before decoupling. Despite our focus on neutrinos, the following
approach can be applied to chiral fermions in general.

We first introduce the thermal average, which is defined for a generic
operator $\hat{O}$ as

$\displaystyle\langle\hat{O}\rangle$ $\displaystyle\equiv$
$\displaystyle\mathrm{Tr}\,\big{(}\hat{O}\hat{\rho}\big{)},\quad\hat{\rho}=\frac{1}{Z}{\rm
e}^{-\beta\big{[}(E_{q}-\mu)\hat{n}_{\bm{q}}+(E_{q}+\mu)(\hat{\tilde{n}}_{\bm{q}}-1)\big{]}}\,,$
$\displaystyle Z$ $\displaystyle\equiv$ $\displaystyle\mathrm{Tr}\,\left({\rm
e}^{-\beta\left[(E_{q}-\mu)\hat{n}_{\bm{q}}+(E_{q}+\mu)(\hat{\tilde{n}}_{\bm{q}}-1)\right]}\right),$
(3.27)

where $\hat{n}_{\bm{q}}\equiv(a^{\dagger}_{\bm{q}}a_{\bm{q}})/{\cal V}$ and
$\hat{\tilde{n}}_{\bm{q}}\equiv(b^{\dagger}_{\bm{q}}b_{\bm{q}})/{\cal V}$ with
${\cal V}\equiv(2\pi)^{3}\delta^{3}(\bm{0})$ denote the normalized density
operators for neutrinos and anti-neutrinos and $\beta=1/T$. For the case of
the expanding universe that we will consider later, all dimensionful
quantities, such as the physical temperature $T$ and chemical potential $\mu$
(as well as momenta) will be rescaled by the scale factor $a$ as $\bar{T}=aT$
and $\bar{\mu}=a\mu$, when taking the thermal average. As $T(\tau)\propto
1/a(\tau)$ and $\mu(\tau)\propto 1/a(\tau)$, $\bar{T}$ and $\bar{\mu}$ are
constants as a function of $\tau$. Nonetheless, for notational simplicity, we
will use ${\bm{k}}$ etc. to denote the rescaled momentum. The physical
momentum depending on $\tau$ is given by ${\bm{k}}_{\rm
phys}(\tau)={\bm{k}}/a(\tau)$. The traces in eq. (3.27) are taken for vacuum
and one-particle (and anti-particle) states, while states with more than one
particle are forbidden from the Pauli exclusion principle. Specifically, we
define

$\displaystyle\mathrm{Tr}\,(\hat{O})=\langle
0|\hat{O}|0\rangle+\int_{\bm{q}}\langle
0|a_{\bm{q}}\hat{O}a_{\bm{q}}^{\dagger}|0\rangle+\int_{\bm{q}}\langle
0|b_{\bm{q}}\hat{O}b_{\bm{q}}^{\dagger}|0\rangle+\frac{1}{\cal
V}\int_{\bm{q}}\langle
0|a_{\bm{q}}b_{\bm{q}}\hat{O}b_{\bm{q}}^{\dagger}a_{\bm{q}}^{\dagger}|0\rangle$
(3.28)

and $\int_{\bm{q}}{\cal V}\equiv 1$. For this, we find

$\displaystyle Z=1+{\rm e}^{-\beta(E_{q}-\mu)}+{\rm e}^{\beta(E_{q}+\mu)}+{\rm
e}^{-\beta(E_{q}-\mu)+\beta(E_{q}+\mu)}=\frac{1}{(1-f_{E_{q}-\mu})f_{E_{q}+\mu}}\,,$
(3.29)

where $f_{q}\equiv 1/({\rm e}^{\beta q}+1)$. Furthermore, one can easily
confirm that the expectation values of density operators give rise to the
Fermi-Dirac distribution functions,

$\displaystyle\langle\hat{n}_{\bm{p}}\rangle$ $\displaystyle=$
$\displaystyle\int_{\bm{q}}(2\pi)^{3}\delta^{3}({\bm{q}-\bm{p}})\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}\bigl{[}{\rm
e}^{\beta(E_{q}+\mu)}+1\bigr{]}=\frac{1}{{\rm e}^{\beta(E_{p}-\mu)}+1}\,,$
$\displaystyle\langle\hat{\tilde{n}}_{\bm{p}}\rangle$ $\displaystyle=$
$\displaystyle\int_{\bm{q}}(2\pi)^{3}\delta^{3}({\bm{q}-\bm{p}})\frac{1}{Z}\bigl{[}{\rm
e}^{-\beta(E_{q}-\mu)}+1\bigr{]}=\frac{1}{{\rm e}^{\beta(E_{p}+\mu)}+1}\,.$
(3.30)

We now replace the vacuum expectation value of eq. (3.16) with the trace of
eq. (3.27), to incorporate the thermal excitations and the effects of the
chemical potential, and look for combinations of operators, which have non-
trivial expectation values. This will include the vacuum part that we have
already computed in the previous subsection. As the details of this
calculation are rather involved, we relegate them to appendix A, and here show
only the final result, which can be divided into the vacuum and thermal parts,

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle\equiv\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\text{vac}}+\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\text{th}},$
(3.31)

where $\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\text{vac}}$
is given in eq. (3.2), which generally suffers from an ultraviolet divergence
and may have to be renormalized by subtracting proper counter terms or by
introducing a cutoff [23]. For the thermal part, after the rescaling of
dimensionful quantities, we have

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\text{th}}$
$\displaystyle=$
$\displaystyle\delta_{\lambda,\lambda^{\prime}}\frac{\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{m_{\rm
P}^{4}}\int{\rm d}\tau^{\prime}{\rm
d}\tau^{\prime\prime}\frac{G_{k}(\tau,\tau^{\prime})G_{k^{\prime}}(\tau,\tau^{\prime\prime})}{a(\tau^{\prime})^{2}a(\tau^{\prime\prime})^{2}}\int{\rm
d}^{3}p\ p^{2}$ (3.32)
$\displaystyle\times\Big{[}\Gamma_{+}(\bm{p},\theta;\bm{k},\lambda;\tau^{\prime},\tau^{\prime\prime})f_{|{\bm{k}-\bm{p}}|-{\bar{\mu}}}+\Gamma_{-}(\bm{p},\theta;\bm{k},\lambda;\tau^{\prime},\tau^{\prime\prime})f_{p+\bar{\mu}}$
$\displaystyle-\Gamma_{0}(\bm{p},\theta;\bm{k},\lambda;\tau^{\prime},\tau^{\prime\prime})\big{(}f_{|{\bm{k}-\bm{p}}|+\bar{\mu}}+f_{p-\bar{\mu}}\big{)}\Big{]},$

where

$\displaystyle\Gamma_{\pm}(\bm{p},\theta;\bm{k},\lambda;\tau^{\prime},\tau^{\prime\prime})$
$\displaystyle={\rm e}^{{\rm i}(|{\bm{k}-\bm{p}}|-p)\tau^{\prime}}{\rm
e}^{{\rm
i}(p-|{\bm{k}-\bm{p}}|)\tau^{\prime\prime}}\tilde{\Gamma}_{\pm}(\bm{p},\theta;\bm{k},\lambda),$
(3.33)
$\displaystyle\Gamma_{0}(\bm{p},\theta;\bm{k},\lambda;\tau^{\prime},\tau^{\prime\prime})$
$\displaystyle={\rm e}^{-{\rm i}(|{\bm{k}-\bm{p}}|+p)\tau^{\prime}}{\rm
e}^{{\rm
i}(p+|{\bm{k}-\bm{p}}|)\tau^{\prime\prime}}\tilde{\Gamma}_{0}(\bm{p},\theta;\bm{k},\lambda),$
(3.34)

with

$\displaystyle\tilde{\Gamma}_{\pm}(\bm{p},\theta;\bm{k},\lambda)$
$\displaystyle=\sin^{2}\theta\Bigg{[}1\mp\lambda\cos\theta-(\cos\theta\mp\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]},$
(3.35) $\displaystyle\tilde{\Gamma}_{0}(\bm{p},\theta;\bm{k},\lambda)$
$\displaystyle=\sin^{2}\theta\Bigg{[}1-\lambda\cos\theta+(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}.$
(3.36)

##  4 Gravitational birefringence from primordial neutrinos

Utilizing the formalism introduced in the previous section, we next explicitly
evaluate the Stokes parameter due to primordial neutrinos in thermal
equilibrium. As a heuristic study, we first investigate the case of the
Minkowski universe in subsection 4.1. Subsequently, we analyze the more
practical case for an expanding universe in the radiation dominated period in
subsection 4.2.

###  4.1 Minkowski universe

In the Minkowski spacetime, where $\tau=t$, $a(\tau)=1$, and
$\mathcal{H}(\tau)=0$, eq. (3.7) simply results in

$\displaystyle G_{k}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle-\int\frac{{\rm d}k_{0}}{2\pi}{\rm e}^{-{\rm
i}k_{0}(t-t^{\prime})}\frac{1}{(k_{0}+i\epsilon)^{2}-k^{2}}$ (4.1)
$\displaystyle=$ $\displaystyle\theta(t-t^{\prime})\frac{1}{k}\sin
k(t-t^{\prime})\,.$

Moreover, the fermionic wave functions reduce to the plane waves,

$\displaystyle u(t,p)={\rm e}^{-{\rm i}p\cdot x},\quad v(t,p)={\rm e}^{{\rm
i}p\cdot x}.$ (4.2)

One can then easily evaluate the integrations with respect to $t^{\prime}$ and
$t^{\prime\prime}$ in eq. (3.2) and derive the correlation function in vacuum,

$\displaystyle\langle
h_{\lambda}(\bm{k},t)h_{\lambda^{\prime}}(\bm{k}^{\prime},t)\rangle$
$\displaystyle=$
$\displaystyle\delta_{\lambda,\lambda^{\prime}}\frac{\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{4m_{\rm
P}^{4}}\int{\rm
d}^{3}p\frac{\sin^{2}\theta}{\big{(}p-k\cos\theta+|{\bm{k}-\bm{p}}|\big{)}^{2}}$
(4.3)
$\displaystyle\times\Bigg{[}1-\lambda\cos\theta+(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]},$

where $|{\bm{k}-\bm{p}}|=\sqrt{p^{2}-2pk\cos\theta+k^{2}}$. The correlation
function is in fact independent of $t$ since now we work in the static
background spacetime.

For the thermal part, by carrying out similar computations, eq. (3.32) becomes

$\displaystyle\langle
h_{\lambda}(\bm{k},t)h_{\lambda^{\prime}}(\bm{k}^{\prime},t)\rangle_{\text{th}}$
$\displaystyle=$
$\displaystyle\delta_{\lambda,\lambda^{\prime}}\frac{\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{4m_{\rm
P}^{4}}\int{\rm
d}^{3}p\Big{[}\Pi_{+}(\bm{p},\theta;\bm{k},\lambda)f_{|{\bm{k}-\bm{p}}|-\mu}+\Pi_{-}(\bm{p},\theta;\bm{k},\lambda)f_{p+\mu}$
(4.4)
$\displaystyle-\Pi_{0}(\bm{p},\theta;\bm{k},\lambda)\big{(}f_{|{\bm{k}-\bm{p}}|+\mu}+f_{p-\mu}\big{)}\Big{]}\,,$

where

$\displaystyle\Pi_{\pm}(\bm{p},\theta;\bm{k},\lambda)$
$\displaystyle=\frac{\sin^{2}\theta}{\big{(}p-k\cos\theta-|{\bm{k}-\bm{p}}|\big{)}^{2}}\left[1\mp\lambda\cos\theta-(\cos\theta\mp\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\right]\,,$
(4.5) $\displaystyle\Pi_{0}(\bm{p},\theta;\bm{k},\lambda)$
$\displaystyle=\frac{\sin^{2}\theta}{\big{(}p-k\cos\theta+|{\bm{k}-\bm{p}}|\big{)}^{2}}\left[1-\lambda\cos\theta+(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\right]\,.$
(4.6)

When considering the long-wavelength (small $k$) limit, the $p$ integral can
be computed analytically. In this limit, one finds
$\Pi_{\pm}(\bm{p},\theta;\bm{k},\lambda)$ to be the dominant terms. In
particular, we find

$\displaystyle\Pi_{+}(\bm{p},\theta;\bm{k},\pm
1)\approx\left[\frac{k^{4}}{4p^{2}}\pm\frac{\cos^{2}(\theta/2)k^{5}}{2p^{3}}+\mathcal{O}\left(\frac{k^{6}}{p^{4}}\right)\right]^{-1}\approx\frac{4p^{2}}{k^{4}}\mp\frac{8p\cos^{2}(\theta/2)}{k^{3}}\,,$
(4.7)

and can approximate $f_{|{\bm{k}-\bm{p}}|-\mu}\rightarrow f_{p-\mu}$.

Interestingly, when taking $\mu\gg T$ despite being an unrealistic condition
in the early Universe, the Fermi-Dirac distribution function reduces to a
unit-step function, which leads to $f_{p-\mu}\rightarrow\Theta(\mu-p)$ and
$f_{p+\mu}\rightarrow\Theta(\mu+p)=0$ for $\mu>0$. In such a case, only
$\Pi_{+}(\bm{p},\theta;\bm{k},\pm 1)$ contributes for $k\ll\mu$ and eq. (4.4)
approximately becomes

$\displaystyle\langle
h_{\pm}(\bm{k},t)h_{\pm}(\bm{k}^{\prime},t)\rangle_{\text{th}}$
$\displaystyle\approx$
$\displaystyle\frac{2\pi\delta^{3}({\bm{k}+\bm{k}^{\prime}})}{4m_{\rm
P}^{4}}\int^{\infty}_{0}{\rm d}p\int^{\pi}_{0}{\rm d}\theta
p^{2}\sin\theta\left[\frac{4p^{2}}{k^{4}}\mp\frac{8p\cos^{2}(\theta/2)}{k^{3}}\right]\Theta(\mu-p)$
(4.8) $\displaystyle=$
$\displaystyle\frac{\pi\delta^{3}({\bm{k}+\bm{k}^{\prime}})\mu^{4}}{5m_{\rm
P}^{4}k^{4}}(4\mu\mp 5k)\,.$

Note that the primary contributions for both integrals come from $p=\mu$. The
results above are accordingly insensitive to the infrared cutoff of $p$ and
consistent with the assumption $p\gg k$ even though we take
$p_{\text{min}}=0$. The above result furthermore has an IR singularity,
behaving as $\sim 1/k^{4}$ from small $k$, which can be traced back to the
infinite ranges of the $t^{\prime}$ and $t^{\prime\prime}$ integrals in eq.
(3.2). It can be explicitly verified that this IR singularity vanishes after
introducing a finite cutoff to the before-mentioned time integrals. As will be
seen in the following section, such a singularity is indeed absent in the
expanding universe case, where all time integral ranges are finite.

In the general case with finite $T$ and $\mu$, we find

$\displaystyle\Pi_{+}(\bm{p},\theta;\bm{k},\lambda)f_{|{\bm{k}-\bm{p}}|-\mu}$
$\displaystyle\approx$
$\displaystyle\frac{4p^{2}}{k^{4}}f_{p-\mu}-\frac{4p}{k^{3}}\left(\lambda+\cos\theta\left[1-(1-f_{p-\mu})\frac{p}{T}\right]\right)f_{p-\mu}\,,$
$\displaystyle\Pi_{-}(\bm{p},\theta;\bm{k},\lambda)f_{p+\mu}$
$\displaystyle\approx$
$\displaystyle\frac{4p^{2}}{k^{4}}f_{p+\mu}+\frac{4p}{k^{3}}(\lambda-\cos\theta)f_{p+\mu}\,,$
(4.9)

up to $\mathcal{O}(k^{-3})$ when $k\ll T,\mu$. From the definition of eq.
(2.21), we can obtain the Stokes parameter $V(\bm{k})$ as

$\displaystyle V(\bm{k})$ $\displaystyle\approx$
$\displaystyle\frac{\pi}{2m_{\rm P}^{4}}\int^{\infty}_{0}{\rm
d}p\int^{\pi}_{0}{\rm d}\theta
p^{2}\sin\theta\frac{8p}{k^{3}}(f_{p+\mu}-f_{p-\mu})$ (4.10) $\displaystyle=$
$\displaystyle-\frac{48\pi T^{4}}{m_{\rm
P}^{4}k^{3}}\Big{[}\rm{Li}_{4}\bigl{(}-{\rm
e}^{-\mu/T}\bigr{)}-\rm{Li}_{4}\bigl{(}-{\rm e}^{\mu/T}\bigr{)}\Big{]},$

where we used

$\displaystyle\int^{\infty}_{0}{\rm
d}pp^{3}f_{p-\mu}=-6T^{4}\rm{Li}_{4}\bigl{(}-{\rm e}^{\mu/T}\bigr{)}$ (4.11)

with $\mathrm{Li}_{s}(z)=\sum_{k=1}^{\infty}z^{k}/k^{s}$ being the
polylogarithm function. When $\mu\ll T$, eq. (4.10) reduces to

$\displaystyle V(\bm{k})\approx-\frac{72\pi T^{4}}{m_{\rm
P}^{4}k^{3}}\zeta(3)\frac{\mu}{T}+\mathcal{O}\left(\left(\frac{\mu}{T}\right)^{2}\right),$
(4.12)

where $\zeta(x)$ is the Riemann zeta function with $\zeta(3)\simeq 1.202$.
Accordingly, one finds that $V(\bm{k})$ vanishes for the
$\Pi_{\pm}(\bm{p},\theta;\bm{k},\lambda)$ terms when $\mu=0$. Even when
including the higher-order corrections up to $\mathcal{O}(k^{0})$, this
conclusion remains unchanged.

We, however, find a non-vanishing $V(\bm{k})$ for $\mu=0$ by including terms
up to $\mathcal{O}(k)$. The nonzero contribution here comes only from the
$\Pi_{0}(\bm{p},\theta;\bm{k},\lambda)$ term, as

$\displaystyle V(\bm{k})\approx-\frac{\pi k}{15m_{\rm P}^{4}}\,,$ (4.13)

which shows that $V(\bm{k})$ can be non-vanishing even for $\mu=0$. A detailed
analysis going beyond the long wavelength limit is, however, out of the scope
of the present study.

###  4.2 Expanding universe

In this section, we consider the expanding universe in the radiation dominated
period. The scale factor $a(\tau)$ in this period behaves as

$\displaystyle a(\tau)=a_{\rm i}\frac{\tau}{\tau_{\rm i}}\,,$ (4.14)

where $\tau_{\rm i}$ stands for the initial time at which thermalized left-
handed neutrinos are realized and $a_{\rm i}$ is the corresponding scale
factor. Furthermore, the solution of eq. (3.7) is known to be444Note that the
solution is distinct from the one in ref. [23] due to the different periods
considered in the early Universe, which give rise to the different forms of
$\mathcal{H}$ in eq. (3.7). Here we have $\mathcal{H}=1/\tau$ for the
radiation dominated period as opposed to $\mathcal{H}=-1/\tau$ for the end of
inflation considered in ref. [23].

$\displaystyle
G(\tau,\tau^{\prime})=\Theta(\tau-\tau^{\prime})\tau^{\prime}\frac{\sin
k(\tau-\tau^{\prime})}{k\tau}\,.$ (4.15)

With these ingredients, we next need to evaluate the $\tau^{\prime}$,
$\tau^{\prime\prime}$ and $\bm{p}$ integrals in eq. (3.32). Unlike the case in
the Minkowski spacetime, it is not possible to compute them analytically for
arbitrary $\tau$, $\bm{k}$, $\bar{T}=aT$ and $\bar{\mu}=a\mu$ (where these
dimensionful quantities are rescaled by the scale factor $a$ as we mention in
subsection 3.3), while analytical expressions can be obtained for certain
physically relevant limits. Specifically, we will consider the regime
$\bar{T}\gg k,\bar{\mu}$. It will, however, become clear below that both $k$
and $\bar{\mu}$ need to be non-zero to obtain a non-zero $V$ parameter.
Furthermore, $\bar{T}$ is considered large enough such that $\tau$, $\tau_{\rm
i}\gg 1/\bar{T}$ (see appendix E). As will be demonstrated later, these are
reasonable assumptions for the temperature ranges of the radiation dominated
period considered in this work.

Let us start with the $\tau^{\prime}$ and $\tau^{\prime\prime}$ integrals. We
define

$\displaystyle\bar{\mathcal{I}}(\tau,\tau_{\rm i},b)=\int_{\tau_{\rm
i}}^{\infty}{\rm
d}\tau^{\prime}\frac{G(\tau,\tau^{\prime})}{a^{2}(\tau^{\prime})}{\rm e}^{{\rm
i}b\tau^{\prime}},$ (4.16)

where

$\displaystyle b=|\bm{k}-\bm{p}|-p\equiv b_{-}$ (4.17)

for the $\Gamma_{\pm}$ terms and

$\displaystyle b=|\bm{k}-\bm{p}|+p\equiv b_{+}$ (4.18)

for the $\Gamma_{0}$ term in eq. (3.32). We will treat each case separately in
the small $k$ limit and keep terms linear in $k$. Note that the
$\tau^{\prime\prime}$ integral in eq. (3.32) corresponds to
$\bar{\mathcal{I}}(\tau,\tau_{\rm
i},-b)=\bar{\mathcal{I}}^{\ast}(\tau,\tau_{\rm i},b)$. Hence, we have

$\displaystyle\int_{\tau_{\rm i}}^{\infty}{\rm d}\tau^{\prime}{\rm
d}\tau^{\prime\prime}\frac{G(\tau,\tau^{\prime})G(\tau,\tau^{\prime\prime})}{a^{2}(\tau^{\prime})a^{2}(\tau^{\prime\prime})}{\rm
e}^{{\rm i}b\tau^{\prime}}{\rm e}^{-{\rm
i}b\tau^{\prime\prime}}=|\bar{\mathcal{I}}(\tau,\tau_{\rm
i},b)|^{2}=|\bar{\mathcal{I}}(\tau,\tau_{\rm i},-b)|^{2},$ (4.19)

and it is therefore sufficient to consider either the $\tau^{\prime}$ or
$\tau^{\prime\prime}$ integral. The correlation function is furthermore
guaranteed to be real. Leaving the details of the calculation to appendix B,
we here directly show the final results for the $\Gamma_{\pm}$ and
$\Gamma_{0}$ terms. For $b=b_{-}$ ($\Gamma_{\pm}$ terms) we have

$\displaystyle|\bar{\mathcal{I}}(\tau,\tau_{\rm
i},b_{-})|^{2}\simeq\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}+\mathcal{O}(k^{2})\,.$
(4.20)

For $b=b_{+}$ ($\Gamma_{0}$ term), the result reads

$\displaystyle|\bar{\mathcal{I}}(\tau,\tau_{\rm i},b_{+})|^{2}\simeq$
$\displaystyle\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\Biggl{\\{}\Bigl{(}\mathrm{ci}(2\tau p)-\mathrm{ci}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\cos[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}$ (4.21) $\displaystyle+\Bigl{(}\mathrm{Si}(2\tau
p)-\mathrm{Si}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\sin[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle+2\frac{k}{p}\cos\theta\sin[(\tau-\tau_{\rm
i})p]\Bigr{(}\sin[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)\bigr{]}$ $\displaystyle-\cos[(\tau+\tau_{\rm
i})p]\bigl{[}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)\bigr{]}\Bigr{)}\Biggr{\\}}+\mathcal{O}(k^{2})\,,$

where $\mathrm{ci}(x)$ and $\mathrm{Si}(x)$ are the cosine and sine integral
functions defined as

$\displaystyle\mathrm{ci}(x)$ $\displaystyle\equiv-\int_{x}^{\infty}{\rm
d}t\frac{\cos t}{t}=\gamma_{\,\mathrm{E}}+\log x+\int_{0}^{x}{\rm
d}t\frac{\cos t-1}{t}\,,$ (4.22) $\displaystyle\mathrm{Si}(x)$
$\displaystyle\equiv\int_{0}^{x}{\rm d}t\frac{\sin t}{t}\,,$ (4.23)

respectively, with $\gamma_{\,\mathrm{E}}$ standing for the Euler-Mascheroni
constant.

Next, we consider the angular integrals of $\bm{p}$, $\phi$ and $\theta$.
While the $\phi$ integral is trivial, the $\theta$ integral requires some
work, but can nevertheless be done analytically in the small $k$ limit.
Relegating the details of this calculation to appendix C, we here give the
final results individually for the contributions of both $\Gamma_{\pm}$ and
$\Gamma_{0}$ terms to eq. (3.32). For the $\Gamma_{\pm}$ terms, we have

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{\pm}}$
$\displaystyle\simeq\frac{32\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}\int_{0}^{\infty}{\rm
d}p\
p^{4}\left(1\mp\lambda\frac{k}{p}\right)f_{p\mp\bar{\mu}}+\mathcal{O}(k^{2}).$
(4.24)

The expression of the $\Gamma_{0}$ term turns out to be much lengthier. We
therefore give the leading-order (LO) and next-to-leading-order (NLO) terms in
terms of $k$ separately. They read

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{LO}}$
$\displaystyle\simeq\frac{16\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\int_{0}^{\infty}{\rm d}p\
p^{4}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})$
$\displaystyle\quad\times\Biggl{\\{}\Bigl{(}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\cos[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle\quad+\Bigl{(}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\sin[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}\Biggr{\\}}\,,$ (4.25)

and

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{NLO}}$
$\displaystyle\simeq-\frac{16\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\lambda k\int_{0}^{\infty}{\rm
d}p\ p^{4}$
$\displaystyle\quad\times\Biggl{\\{}\frac{2}{p}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})\sin[(\tau-\tau_{\rm
i})p]\Bigr{(}\sin[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)\bigr{]}$
$\displaystyle\quad-\cos[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{Si}(2\tau
p)-\mathrm{Si}(2\tau_{\rm i}p)\bigr{]}\Bigr{)}$
$\displaystyle\quad+\left[-\frac{2}{p}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})+\frac{1}{T}f_{p+\bar{\mu}}(1-f_{p+\bar{\mu}})\right]$
$\displaystyle\quad\times\biggl{[}\Bigl{(}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\cos[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle\quad+\Bigl{(}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\sin[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}\biggr{]}\Biggr{\\}}\,.$ (4.26)

As the last step, we compute the remaining integral of $p$. Let us first
consider the $\Gamma_{\pm}$ terms which are simpler than $\Gamma_{0}$. We
leave the calculational details to appendix D and give here the result for the
sum of the $\Gamma_{+}$ and $\Gamma_{-}$ terms, which reads

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{+}+\Gamma_{-}}$
$\displaystyle\simeq-\frac{256\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}\bar{T}^{5}$
$\displaystyle\quad\times\Biggl{\\{}\mathrm{Li}_{5}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}+\mathrm{Li}_{5}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}-\frac{\lambda}{4}\frac{k}{\bar{T}}\Bigl{[}\mathrm{Li}_{4}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}-\mathrm{Li}_{4}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}\Bigr{]}\Biggr{\\}}$
$\displaystyle\quad+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\Bigr{)}.$
(4.27)

Expanding the above expression in $\bar{\mu}/\bar{T}$, we have

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{+}+\Gamma_{-}}$
$\displaystyle\simeq\frac{96\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm
i}}{\tau}\right]^{2}\bar{T}^{5}\left[5\zeta(5)-\zeta(3)\lambda\frac{k}{\bar{T}}\frac{\bar{\mu}}{\bar{T}}\right]$
$\displaystyle\quad+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\mathrm{,}\,\Bigl{(}\frac{\bar{\mu}}{\bar{T}}\Bigr{)}^{2}\Bigr{)},$
(4.28)

with $\zeta(5)\simeq 1.037$.

To compute the $\Gamma_{0}$ term, one needs to integrate eqs. (4.25) and
(4.26), which does not appear to be possible analytically in general. However,
considering the time and temperature regions that are relevant for the
radiation dominant period, some of the functions $\mathrm{Si}(x)$ and
$\mathrm{ci}(x)$ appearing in these equations can be simplified as follows.
After the change of integration variable from $p$ to $y=p/\bar{T}$, the
arguments of these functions become either $x=2\tau{\bar{T}}y$ or $2\tau_{\rm
i}{\bar{T}}y$, while the integral can be expected to be dominated by
contributions from the region $y\sim 1$. Noting that $\tau{\bar{T}}\gg 1$ and
$\tau_{\rm i}{\bar{T}}\gg 1$ in the radiation dominant period (see appendix E
for details), the functions $\mathrm{Si}(2\tau{\bar{T}}y)$,
$\mathrm{Si}(2\tau_{\rm i}{\bar{T}}y)$, $\mathrm{ci}(2\tau{\bar{T}}y)$ and
$\mathrm{ci}(2\tau_{\rm i}{\bar{T}}y)$ can be approximated using their
asymptotic forms

$\displaystyle\mathrm{Si}(x)$ $\displaystyle=$
$\displaystyle\frac{\pi}{2}-\frac{\cos
x}{x}+\mathcal{O}\left(\frac{1}{x^{2}}\right),$ (4.29)
$\displaystyle\mathrm{ci}(x)$ $\displaystyle=$ $\displaystyle\frac{\sin
x}{x}+\mathcal{O}\left(\frac{1}{x^{2}}\right).$ (4.30)

Within these approximations, the $p$ integral in eqs. (4.25) and (4.26) can be
computed as (see appendix D for details)

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0}}$
$\displaystyle\simeq-\frac{8\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}\bar{T}^{3}\Biggl{\\{}3\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}\Bigl{[}\mathrm{Li}_{3}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}+\mathrm{Li}_{3}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}\Bigr{]}$
$\displaystyle\quad+\lambda\frac{k}{\bar{T}}\Bigl{[}\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}\mathrm{Li}_{2}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}-2\bar{T}{\rm
e}^{\bar{\mu}/\bar{T}}\mathrm{Im}\,\Phi^{\ast}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}},3,1-2{\rm i}(\tau-\tau_{\rm i})\bar{T}\bigr{)}$
$\displaystyle\quad-2\bar{T}{\rm
e}^{-\bar{\mu}/\bar{T}}\mathrm{Im}\,\Phi^{\ast}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}},3,1-2{\rm i}(\tau-\tau_{\rm
i})\bar{T}\bigr{)}\Bigr{]}\Biggr{\\}}+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\Bigr{)}\,,$
(4.31)

where $\Phi^{\ast}(z,s,c)$ is the Lerch transcendent function, defined as

$\displaystyle\Phi^{\ast}(z,s,c)\equiv\sum_{k=0}^{\infty}\frac{z^{k}}{\bigl{[}(c+k)^{2}\bigr{]}^{s/2}}\,.$
(4.32)

Expanding up to first order in $\bar{\mu}/\bar{T}$, we have

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0}}$
$\displaystyle\simeq\frac{4\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}\bar{T}^{3}\Biggl{\\{}9\zeta(3)\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}$
$\displaystyle\quad+\lambda\frac{k}{\bar{T}}\left(\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm
i}}\left[\frac{\pi^{2}}{6}+2\log(2)\frac{\bar{\mu}}{\bar{T}}\right]-\frac{1}{2(\tau-\tau_{\rm
i})^{3}\bar{T}^{2}}+\pi^{3}\bar{T}\frac{3+\cosh\bigl{[}4\pi(\tau-\tau_{\rm
i})\bar{T}\bigr{]}}{\sinh^{3}\bigl{[}2\pi(\tau-\tau_{\rm
i})\bar{T}\bigr{]}}\right)\Biggr{\\}}$
$\displaystyle\quad+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\mathrm{,}\,\Bigl{(}\frac{\bar{\mu}}{\bar{T}}\Bigr{)}^{2}\Bigr{)}\,.$
(4.33)

It is interesting to observe that the difference between the two helicities
here does not vanish even for $\mu=0$, which is different from the situation
in the $\Gamma_{\pm}$ terms [see eq. (4.28)].

Let us now add all terms to get the full correlator, which in the small
$k/\bar{T}$ and small $\bar{\mu}/\bar{T}$ expansion reads

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\tau_{\rm
i}{\bar{T}}\gg 1}$
$\displaystyle\simeq\frac{8\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\bar{T}^{5}\Biggl{\\{}180\zeta(5)\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm
i}}{\tau}\right]^{2}+\frac{9}{2}\zeta(3)\frac{(\tau-\tau_{\rm
i})^{2}}{\tau^{2}\tau_{\rm i}^{2}\bar{T}^{2}}$
$\displaystyle\quad+\lambda\frac{k}{\bar{T}}\Biggl{[}\frac{\pi^{2}}{12}\frac{(\tau-\tau_{\rm
i})^{2}}{\tau^{2}\tau_{\rm i}^{2}\bar{T}^{2}}-\frac{1}{4(\tau-\tau_{\rm
i})^{2}\tau\tau_{\rm i}\bar{T}^{4}}+\frac{\pi^{3}}{2}\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}\bar{T}}\frac{3+\cosh\bigl{[}4\pi(\tau-\tau_{\rm
i})\bar{T}\bigr{]}}{\sinh^{3}\bigl{[}2\pi(\tau-\tau_{\rm i})\bar{T}\bigr{]}}$
$\displaystyle\quad+\frac{\bar{\mu}}{\bar{T}}\biggl{(}-36\zeta(3)\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm
i}}{\tau}\right]^{2}+\log(2)\frac{(\tau-\tau_{\rm i})^{2}}{\tau^{2}\tau_{\rm
i}^{2}\bar{T}^{2}}\biggr{)}\Biggr{]}\Biggr{\\}}$
$\displaystyle\quad+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\mathrm{,}\,\Bigl{(}\frac{\bar{\mu}}{\bar{T}}\Bigr{)}^{2}\Bigr{)}\,.$
(4.34)

One can further clarify the importance of the different terms by considering
the $\tau\gg\tau_{\rm i}$ limit. Furthermore, as we concluded in eq. (E.7),
factors such as $1/(\tau\bar{T})$ or $1/\sinh^{3}\bigl{[}2\pi(\tau-\tau_{\rm
i})\bar{T}\bigr{]}$ are strongly suppressed and $1/(\tau_{\rm i}\bar{T})$ can
in the limit considered here also be regarded as small. The correlator thus
becomes

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\tau_{\rm
i}{\bar{T}}\gg 1}$
$\displaystyle\simeq\frac{8\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\bar{T}^{5}\Biggl{\\{}180\zeta(5)\log^{2}\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}+\lambda\frac{k}{\bar{T}}\biggl{[}\frac{\pi^{2}}{12}\frac{1}{\tau_{\rm
i}^{2}\bar{T}^{2}}-36\zeta(3)\frac{\bar{\mu}}{\bar{T}}\log^{2}\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}\biggr{]}\Biggr{\\}}$
$\displaystyle\quad+\mathcal{O}\Bigl{(}\Bigl{(}\frac{k}{\bar{T}}\Bigr{)}^{2}\mathrm{,}\,\Bigl{(}\frac{\bar{\mu}}{\bar{T}}\Bigr{)}^{2}\Bigr{)}\,.$
(4.35)

In terms of the Stokes parameters, the above results can be summarized as

$\displaystyle I(\bm{k})$ $\displaystyle\simeq 192\pi\zeta(5)\frac{1}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\bar{T}^{5}\log^{2}\Bigl{(}\frac{\tau}{\tau_{\rm i}}\Bigr{)}\,,$ (4.36)
$\displaystyle V(\bm{k})$ $\displaystyle\simeq\frac{192\pi}{5}\frac{1}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\bar{T}^{4}k\bigg{[}\frac{\pi^{2}}{432}\frac{1}{\tau_{\rm
i}^{2}\bar{T}^{2}}-\zeta(3)\frac{\bar{\mu}}{\bar{T}}\log^{2}\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}\bigg{]}\,.$ (4.37)

For $V(\bm{k})$ we have kept the leading $1/(\tau_{\rm i}\bar{T})^{2}$ term,
having the situation of vanishing chemical potential in mind, in which case it
becomes the leading contribution to the difference between the two helicities.
We will, however, for simplicity ignore this term in the following
discussions, as for currently available estimates of $\bar{\mu}/\bar{T}$ the
second term will be dominant.

###  4.3 Comparisons and discussions

In general, the leading-order correlation functions of gravitational fields in
the Minkowski and expanding universes in the long wavelength limit are
different. In the Minkowski spacetime, the correlation functions diverge at
$k\rightarrow 0$. From the contributions of neutrinos both in vacuum and in
thermal equilibrium, the parity-even and parity-odd components of correlation
functions become $\mathcal{O}(k^{-4})$ and $\mathcal{O}(k^{-3})$,
respectively. On the contrary, the correlation functions in the FLRW
background are finite at $k\rightarrow 0$. From neutrinos in thermal
equilibrium, the parity-even and parity-odd components of correlation
functions emerge at $\mathcal{O}(k^{0})$ and $\mathcal{O}(k)$, respectively.
Compared with eq. (4.10) in the Minkowski spacetime, the $V$ parameter in eq.
(4.37) in the expanding universe is thus much suppressed in the long
wavelength regime. Note that both the parity-even and parity-odd components
are at $\mathcal{O}(k^{0})$ in ref. [23] in the long wavelength limit, while a
different mechanism for the helicity imbalance generated in the inflation
period is considered therein.

When focusing on the primary contribution from a non-zero chemical potential
with $\tau_{\rm i}\bar{T}\gg 1$ and recast the temperature into physical one
at $\tau_{\rm i}$, $T_{\rm i}\equiv T(\tau_{\rm i})=\bar{T}/a_{\rm i}$, eq.
(4.37) becomes

$\displaystyle V(\bm{k})\simeq-\frac{192\pi\zeta(3)}{5}\frac{\tau_{\rm
i}^{4}}{m_{\rm P}^{4}}T_{\rm
i}^{4}k\mathcal{N}_{\nu}\xi\log^{2}\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}\,,$ (4.38)

where $\mathcal{N}_{\nu}=3$ is the number of neutrinos and we define the
flavor averaged neutrino degeneracy parameter $\xi$ as

$\displaystyle\xi\equiv\frac{1}{\mathcal{N}_{\nu}}\sum_{j}\xi_{j},$ (4.39)

where $\xi_{j}\equiv\mu_{j}/T=\bar{\mu}_{j}/\bar{T}$, in which the subscript
“$j$” represents the neutrino flavor. Note that $V(\bm{k})$ itself depends on
the scale factor through $k=a(\tau)k_{\rm phys}(\tau)$. In order to make a
direct comparison with observations, we may introduce an observable
independent of the scale factor,

$\displaystyle\Delta\chi({\bm{k}})\equiv\frac{{V}(\bm{k})}{{I}(\bm{k})}=\frac{\langle
h_{+}(\bm{k},t)h_{+}(\bm{k}^{\prime},t)\rangle_{\text{th}}-\langle
h_{-}(\bm{k},t)h_{-}(\bm{k}^{\prime},t)\rangle_{\text{th}}}{\langle
h_{+}(\bm{k},t)h_{+}(\bm{k}^{\prime},t)\rangle_{\text{th}}+\langle
h_{-}(\bm{k},t)h_{-}(\bm{k}^{\prime},t)\rangle_{\text{th}}}\,.$ (4.40)

From eq. (4.2), we find

$\displaystyle\Delta\chi({\bm{k}})=-\frac{\zeta(3)}{5\zeta(5)}\xi\frac{k}{\bar{T}}=-\frac{\zeta(3)}{5\zeta(5)}\xi\frac{k_{\rm
g}}{T_{\rm g}}\,,$ (4.41)

where we take $k_{\rm g}=k_{\rm phys}(\tau_{\rm g})$ and $T_{\rm
g}=T(\tau_{\rm g})$ with $k_{\rm phys}(\tau)=k/a(\tau)$ and $\tau_{\rm g}$ the
conformal time when gravitational waves are generated from the thermal plasma.
Note that $\Delta\chi({\bm{k}})$ is not suppressed by $m_{\rm P}$ and depends
on three parameters, $k_{\rm g}$, $T_{\rm g}$ and $\xi$; the parameter $k_{\rm
g}$ can be related to the momentum (wavenumber) of gravitational waves
measured at present time via $k_{\rm g}=a(\tau_{0})k_{\rm
phys}(\tau_{0})/a(\tau_{\rm g})$, where $\tau_{0}$ denotes the present
conformal time and $k_{\rm phy}(\tau_{0})$ corresponds to the measured
momentum of gravitational waves; $T_{\rm g}$ may be chosen as $T_{\rm
g}\gtrsim 1\,{\rm MeV}$ when neutrinos are still in equilibrium in a
primordial plasma during the radiation dominant era; $\xi$ is not well
constrained, but at least for electron neutrinos known to be small, the recent
ref. [37] for instance obtaining a value of $\xi_{\rm e}=0.05\pm 0.03$ (see
also the recent ref. [38], which reports a similar result and refs. [27, 42,
41, 39] for earlier estimates). In ref. [40], a possible theoretical model
that could generate a sufficiently large lepton asymmetry consistent with the
estimate in ref. [37] is also proposed. Yet the exact value of $\xi$ is still
not fully determined from present observations for all neutrino types. In this
work, we assume that the chemical potentials of all neutrino flavors are
equilibrated and hence $\xi_{\rm e}=\xi_{\mu}=\xi_{\tau}=\xi$ [39] and employ
the numerical value of ref. [37] quoted above.

The computation for $\Delta\chi({\bm{k}})$ above only incorporates the
gravitational waves produced from thermalized neutrinos. In practice, we
should further include other sources that contribute to $I(\bm{k})$ and
$V(\bm{k})$. Considering only the gravitational waves produced by matter in
thermal equilibrium from a primordial plasma, we may simply multiply
$I(\bm{k})$ in eq. (4.36) by the effective degrees of freedom. Recall that
$I(\bm{k})$ is dominated by the $\Gamma_{\pm}$ terms associated with the
integral $I_{1}(k,\pm\bar{\mu},\bar{T})$ in eq. (D.1) for $\mu\ll T$. We hence
simply replace the fermionic distribution function in
$I_{1}(k,\pm\bar{\mu},\bar{T})$ by its bosonic counterpart to roughly estimate
the similar bosonic contributions. In addition, we neglect particle masses,
assuming that they are much smaller than the considered temperature scales.
Consequently, the effective degrees of freedom for $I(\bm{k})$ here can be
characterized by

$\displaystyle\mathcal{N}_{\rm d.o.f.}=\mathcal{N}_{\rm
f}+\frac{16}{15}\mathcal{N}_{\rm b}\,,$ (4.42)

due to the relation

$\displaystyle\int_{0}^{\infty}{\rm d}x\frac{x^{4}}{{\rm
e}^{x}-1}=\frac{16}{15}\int_{0}^{\infty}{\rm d}x\frac{x^{4}}{{\rm
e}^{x}+1}\,,$ (4.43)

where $\mathcal{N}_{\rm f}$ and $\mathcal{N}_{\rm b}$ denote the total numbers
of fermions and bosons in thermal equilibrium, respectively. For $V(\bm{k})$,
the LO contribution arises due to the non-vanishing neutrino degeneracy
parameter $\xi$. We eventually obtain a generic expression

$\displaystyle\Delta\chi({\bm{k}})=-\frac{\zeta(3)r(T_{\rm
g})\mathcal{N}_{\nu}}{5\zeta(5)\mathcal{N}_{\rm d.o.f.}}\frac{k_{\rm
phys}(\tau_{0})}{T_{\rm g}}\xi\,.$ (4.44)

Here, $r(T_{\rm g})\equiv a(\tau_{0})/a(\tau_{\rm g})$, for which we may
utilize [43]

$\displaystyle r(T_{\rm g})\simeq 1.25\times 10^{13}\left(\frac{g_{\ast
s}(T_{\rm g})}{100}\right)^{1/3}\left(\frac{T_{\rm g}}{\rm GeV}\right)\,,$
(4.45)

where $g_{\ast s}(T_{\rm g})$ represents effective entropic degrees of freedom
at $T_{\rm g}$. As an example, we choose $T_{\rm g}=200\,{\rm MeV}$ just above
the QCD phase transition. Accordingly, we take $\mathcal{N}_{\rm f}=50$,
$\mathcal{N}_{\rm b}=18$, and $g_{\ast s}(T_{\rm g})\approx 61.75$ [44], which
yield $\mathcal{N}_{\rm d.o.f.}\approx 69.2$ and $r(T_{\rm g})\approx
2.13\times 10^{12}$. In such a case, we obtain

$\displaystyle\Delta\chi({\bm{k}})\approx-5.35\times 10^{9}\left(\frac{k_{\rm
phys}(\tau_{\rm 0})}{\rm GeV}\right)=-2.21\times 10^{-14}\left(\frac{f}{\rm
Hz}\right)\,$ (4.46)

by taking $\xi=0.05$, where we have defined $f\equiv 2\pi k_{\rm
phys}(\tau_{\rm 0})$.

Let us next estimate the gravitational wave spectrum emitted from thermalized
primordial neutrinos that can be measured at the present time. For this
purpose (and for making a direct comparison to the existing literature
simpler), we define the dimensionless power spectra $h^{2}_{I}(\bm{k},\eta)$
and $h^{2}_{V}(\bm{k},\eta)$ as

$\displaystyle\langle
h_{\lambda}(\bm{k},\eta)h_{\lambda^{\prime}}^{\ast}(\bm{k}^{\prime},\eta)\rangle$
$\displaystyle=$
$\displaystyle\frac{8\pi^{5}}{k^{3}}\delta(\bm{k}-\bm{k}^{\prime})\delta_{\lambda\lambda^{\prime}}[h^{2}_{I}(\bm{k},\eta)+\lambda
h^{2}_{V}(\bm{k},\eta)]\,,$ (4.47)

where we have

$\displaystyle h^{2}_{I}(\bm{k},\tau)$ $\displaystyle\simeq$
$\displaystyle\frac{12\zeta(5)}{\pi^{4}}\frac{1}{m_{\rm P}^{4}}\frac{\tau_{\rm
i}^{4}}{a_{\rm i}^{4}}k^{3}\bar{T}^{5}\log^{2}\left(\frac{\tau}{\tau_{\rm
i}}\right)\,,$ (4.48) $\displaystyle h^{2}_{V}(\bm{k},\tau)$
$\displaystyle\simeq$
$\displaystyle-\frac{12\zeta(3)}{5\pi^{4}}\frac{1}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}k^{4}\bar{T}^{4}\xi\log^{2}\left(\frac{\tau}{\tau_{\rm i}}\right)\,.$
(4.49)

To express the gravitational wave spectrum as a function of frequency $f$ at
the present time, we use [43]

$\displaystyle h^{2}\Omega^{I/V}_{\rm{GW}}(f)$ $\displaystyle=$
$\displaystyle\frac{2\pi^{2}}{3}\frac{h^{2}}{H_{0}^{2}}f^{2}h_{I/V}^{2}(f)$
(4.50) $\displaystyle=$
$\displaystyle\frac{2\pi^{2}}{3}10^{-4}\left(\frac{\rm{Mpc}}{\rm{km}}\right)^{2}\left(\frac{f}{\rm{Hz}}\right)^{2}h_{I/V}^{2}(f)\,,$

where $H_{0}$ is the Hubble constant today, which is related to the parameter
$h$ as $H_{0}=100\,h\,\rm{km}\,\rm{s}^{-1}\,\rm{Mpc}^{-1}$. Furthermore,
$h_{I/V}^{2}(f)$ and $f$ are obtained from $h^{2}_{I/V}(\bm{k},\tau)$ at time
$\tau$ and the momentum $\bm{k}$ as [45, 46, 18]

$\displaystyle h^{2}_{I/V}(f)$ $\displaystyle=$
$\displaystyle\left(\frac{\Omega_{m}}{\Omega_{\Lambda}}\right)^{2}\left(\frac{g_{\ast}(T)}{g_{\ast}(T_{0})}\right)\left(\frac{g_{\ast
s}(T_{0})}{g_{\ast
s}(T)}\right)^{4/3}\overline{\left(\frac{3j_{1}(k\tau_{0})}{k\tau_{0}}\right)^{2}}T_{1}^{2}\left(\frac{k_{\mathrm{phys}}}{k_{\mathrm{eq}}}\right)h^{2}_{I/V}(\bm{k},\tau)\,,$
(4.51) $\displaystyle f$ $\displaystyle=$
$\displaystyle\frac{a(\tau)}{a(\tau_{0})}\frac{k_{\mathrm{phys}}(\tau)}{2\pi}\,.$
(4.52)

Here, the factor $\Omega_{m}/\Omega_{\Lambda}$ expresses the ratio of matter
and dark energy in the current universe ($\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$). $g_{\ast}(T)$ denotes, similar to the entropic
degrees of freedom $g_{\ast s}(T)$, the effective number of relativistic
degrees of freedom contributing to the energy density at temperature $T$.
$j_{1}(x)=(\sin x/x-\cos x)/x$ is the spherical Bessel function and the bar
above it represents the averaging over multiple periods. Considering its
argument $x=k\tau_{0}$, we use $\tau_{0}=1.4\times 10^{3}\,{\rm Mpc}$ and
$k=k_{\mathrm{phys}}(\tau_{0})/a(\tau_{0})$ with $a(\tau_{0})=1$ to obtain

$\displaystyle k\tau_{0}$ $\displaystyle=$ $\displaystyle 9.1\times
10^{17}\left(\frac{f}{\mathrm{Hz}}\right)\gg 1\,,$ (4.53)

for all frequencies $f$ relevant for realistic measurements. Therefore, we can
approximate $j_{1}(x)\simeq-\cos x/x$ and the barred factor becomes

$\displaystyle\overline{\left(\frac{3j_{1}(k\tau_{0})}{k\tau_{0}}\right)^{2}}$
$\displaystyle\simeq$ $\displaystyle\frac{9}{2}\frac{1}{(k\tau_{0})^{4}}\,.$
(4.54)

$T_{1}(x)$ is the so-called transfer function which describes the transfer
effect of the universe from the radiation dominant to matter dominant eras and
is obtained by numerically solving the equation governing
$h_{\lambda}(\bm{k},\tau)$ in the expanding universe. It can be given as

$\displaystyle T^{2}_{1}(x)=1+1.57x+3.42x^{2}.$ (4.55)

The momentum $k_{\mathrm{eq}}$ appearing in eq. (4.51) is determined by the
conformal time $\tau_{\mathrm{eq}}$ corresponding to the matter-radiation
equality, $k_{\mathrm{eq}}=1/\tau_{\mathrm{eq}}=7.1\times
10^{-2}\Omega_{m}h^{2}\,\mathrm{Mpc}^{-1}$ with $h=0.7$. Moreover, the ratio
$a(\tau)/a(\tau_{0})$ appearing in eq. (4.52) can, similar to eq. (4.45), be
related to the physical temperatures at times $\tau$ and $\tau_{0}$ from the
principle of conservation of entropy, giving

$\displaystyle\frac{a(\tau)}{a(\tau_{0})}=\frac{T(\tau_{0})}{T(\tau)}\left(\frac{g_{\ast
s}(T_{0})}{g_{\ast s}(T)}\right)^{1/3}\,.$ (4.56)

It should be noted that eq. (4.51) is only applicable for gravitational waves
generated during the radiation dominant period that do not cross the horizon.
This leads to a constraint of the $\tau_{\rm i}$ values (or their
corresponding temperatures) that can be used for a specific frequency $f$. The
respective condition reads

$\displaystyle k\tau_{\rm i}>1.$ (4.57)

To make the physical meaning of this condition more transparent, let us
rewrite it in terms of the temperature $T_{\rm i}$ and the present day
frequency $f$. $k\tau_{\rm i}$ can be obtained by multiplying the factor
$\tau_{\rm i}/\tau_{0}$ to eq. (4.53). We have

$\displaystyle\frac{\tau_{\rm i}}{\tau_{0}}$ $\displaystyle=$
$\displaystyle\frac{\tau_{\rm
i}}{\tau_{\mathrm{eq}}}\frac{\tau_{\mathrm{eq}}}{\tau_{0}}=\frac{a(\tau_{\rm
i})}{a(\tau_{\mathrm{eq}})}\left(\frac{a(\tau_{\mathrm{eq}})}{a(\tau_{0})}\right)^{1/2}\,,$
(4.58) $\displaystyle=$ $\displaystyle\left(\frac{g_{\ast
s}(T_{\mathrm{eq}})g_{\ast s}(T_{0})}{g_{\ast s}^{2}(T_{\rm
i})}\right)^{1/6}\left(\frac{T_{\mathrm{eq}}T_{0}}{T^{2}_{\rm
i}}\right)^{1/2}\,,$ $\displaystyle=$ $\displaystyle 4.6\times
10^{-12}\left(\frac{T_{\rm i}}{\mathrm{GeV}}\right)^{-1}\,,$

where $T_{\mathrm{eq}}=0.8\,{\rm eV}$ is the temperature at the matter-
radiation equality with $g_{\ast s}(T_{\mathrm{eq}})=3.91$. Here we have used
eq. (4.14) in the radiation dominated period and $a(\tau)\propto\tau^{2}$ in
the matter dominated period to derive the second equality of the first line in
eq. (4.58). Furthermore, $T_{0}=2.35\times 10^{-13}\,{\rm GeV}$, $g_{\ast
s}(T_{0})=3.91$ and $g_{\ast s}(T_{\rm i})=106.75$, the latter implying that
$T_{\rm i}$ is large enough to thermalize all particles of the standard model.
In combination with eq. (4.53), eq. (4.57) can hence be recast as

$\displaystyle 4.2\times 10^{6}\left(\frac{T_{\rm
i}}{\mathrm{GeV}}\right)^{-1}\left(\frac{f}{\mathrm{Hz}}\right)>1,$ (4.59)

which shows that for, say, $f\sim 1\,{\rm Hz}$, $T_{\rm i}$ should not be much
larger than $10^{6}\,{\rm GeV}$.

Rewriting $h^{2}_{I/V}(\bm{k},\tau)$ in terms of physical quantities,
combining eqs. (4.50)–(4.56), we arrive at555Here we again use eq. (4.14) in
the radiation dominated period and the relation (4.56) based on the entropy
conservation to evaluate $\ln^{2}(\tau/\tau_{\rm i})$. Such an estimation
could be different from evaluating $\tau/\tau_{\rm i}$ by using the estimation
in eqs. (E.7) and (E.8) in light of the energy conservation, while their
difference in the numerical results is negligible.

$\displaystyle h^{2}\Omega_{\rm{GW}}^{I}(f)$ $\displaystyle=$
$\displaystyle\frac{23328\zeta(5)}{\pi^{3}}\mathcal{N}_{\nu}\left(\frac{\rm{Mpc}}{\rm{km}}\right)^{2}\left(\frac{\rm{Hz}}{\rm{GeV}}\right)^{3}\frac{1}{g_{\ast}(T)g_{\ast}(T_{0})}\left(\frac{g_{\ast
s}(T_{0})}{g_{\ast
s}(T)}\right)^{1/3}\left(\frac{T_{0}}{\rm{GeV}}\right)^{-3}$ (4.60)
$\displaystyle\times\left(\frac{\Omega_{m}}{\Omega_{\Lambda}}\right)^{2}(k\tau_{0})^{-4}\,T_{1}^{2}\left(\frac{k_{\mathrm{phys}}}{k_{\mathrm{eq}}}\right)\log^{2}\left[\frac{T_{\rm
i}}{T}\left(\frac{g_{\ast s}(T_{\rm i})}{g_{\ast
s}(T)}\right)^{1/3}\right]\left(\frac{f}{\rm{Hz}}\right)^{5}\,,$
$\displaystyle h^{2}\Omega_{\rm{GW}}^{V}(f)$ $\displaystyle=$
$\displaystyle\frac{46656\zeta(3)}{5\pi^{2}}\mathcal{N}_{\nu}\left(\frac{\rm{Mpc}}{\rm{km}}\right)^{2}\left(\frac{\rm{Hz}}{\rm{GeV}}\right)^{4}\frac{1}{g_{\ast}(T)g_{\ast}(T_{0})}\left(\frac{T_{0}}{\rm{GeV}}\right)^{-4}$
(4.61)
$\displaystyle\times\left(\frac{\Omega_{m}}{\Omega_{\Lambda}}\right)^{2}(k\tau_{0})^{-4}\,T_{1}^{2}\left(\frac{k_{\mathrm{phys}}}{k_{\mathrm{eq}}}\right)\xi\log^{2}\left[\frac{T_{\rm
i}}{T}\left(\frac{g_{\ast s}(T_{\rm i})}{g_{\ast
s}(T)}\right)^{1/3}\right]\left(\frac{f}{\rm{Hz}}\right)^{6}\,.$

For making a numerical estimate, we need to fix $\tau$, the respective
temperature $T$ and the corresponding active degrees of freedom. While keeping
$\tau_{\rm i}$ and $T_{\rm i}$ as free parameters (with $g_{\ast}(T_{\rm
i})=g_{\ast s}(T_{\rm i})=106.75$ as mentioned above), $\tau$ in our approach
marks the end of the period during which neutrinos emit gravitational waves,
for which we take $T(\tau)=1\,{\rm MeV}$ at the neutrino decoupling
temperature and $g_{\ast}(T)=g_{\ast s}(T)=14.25$. We note that neutrinos will
in principle continue to generate gravitational waves even after $\tau$. For
times between the neutrino decoupling and the end of the radiation dominant
period, the thermal distribution of neutrinos is frozen, with momenta that are
rescaled due to the expanding universe. This rescaling allows one to define a
time-dependent temperature analogous to the situation before $\tau$. The
rescaled momenta $\bm{k}$ and $\bm{p}$ in eq. (3.32) and $\bar{T}$, can
however not both be considered as constants and our analytical results derived
above can therefore no longer be applied. As the universe enters the matter
dominant era, $a(\tau)$ can furthermore no longer be described by eq. (4.14).

Substituting the above numerical values, $\xi=0.05$, and $g_{\ast}(T_{0})=2$,
we have

$\displaystyle h^{2}\Omega_{\rm{GW}}^{I}(f)$ $\displaystyle=$ $\displaystyle
1.4\times 10^{-32}\times\left[1+\frac{1}{49.9}\log\left(\frac{T_{\rm
i}}{m_{\mathrm{P}}}\right)\right]^{2}\left(\frac{f}{\rm{Hz}}\right)^{3}\,,$
(4.62) $\displaystyle h^{2}\Omega_{\rm{GW}}^{V}(f)$ $\displaystyle=$
$\displaystyle 2.9\times
10^{-44}\times\left[1+\frac{1}{49.9}\log\left(\frac{T_{\rm
i}}{m_{\mathrm{P}}}\right)\right]^{2}\left(\frac{f}{\rm{Hz}}\right)^{4}\,,$
(4.63)

which shows that the parity-odd component is (at $f=1\,\mathrm{Hz}$)
suppressed by more than 11 orders of magnitude compared to its parity-even
counterpart. This is caused primarily by one additional factor of $k/\bar{T}$
in its power spectrum and secondarily by the further suppression due to the
smallness of $\xi$. The parity-even gravitational spectrum due to hydrodynamic
fluctuations of the primordial plasma has been previously obtained from the
emission rate of gravitons based on a perturbative Standard Model calculation
[31, 32, 47].

Another frequently used measure to quantify the strength of a gravitational
wave signal is the power spectral density $S_{h}(f)$ [48], which has the
dimension of $\rm{Hz}^{-1}$ and is related to $h^{2}\Omega_{\rm{GW}}(f)$ as

$\displaystyle\frac{S_{h}(f)}{\rm{Hz}^{-1}}=\frac{3}{2\pi^{2}}10^{4}\left(\frac{\rm{Mpc}}{\rm{km}}\right)^{-2}\left(\frac{f}{\rm{Hz}}\right)^{-3}h^{2}\Omega_{\rm{GW}}(f).$
(4.64)

In figure 1, we show $\sqrt{S_{h}(f)/\rm{Hz}^{-1}}$ and
$h^{2}\Omega_{\rm{GW}}(f)$ as functions of the frequency $f/\rm{Hz}$ together
with the sensitivity curves of a few GW detectors [48, 49]. $T_{\rm
i}=10^{2}\,{\rm GeV}$ was used in this figure, corresponding to the
temperature just below the electroweak (EW) phase transition [50]. For
temperatures above the EW phase transition, there may be further chiral
asymmetries related to different conserved charges (see, for instance, ref.
[51]), which will inevitably be model dependent. We however note that the
results of this paper can be easily generalized to any model before the EW
phase transition by summing up all chiral asymmetries related to its available
chemical potentials.

Figure 1: The power spectral densities of the gravitational waves generated
from primordial thermal neutrinos. The even (odd) parity component is shown by
a thick red (blue) line. For comparison, the sensitivity curves of various GW
detectors, extracted from ref. [49], are shown as thin black lines.

As seen in the plots of figure 1, the GW signals of thermal neutrinos are
likely too weak to be measured in the foreseeable future. Let us, however,
mention here the possibility of measuring high frequency GWs (see refs. [52,
53, 54, 55, 56] for recent developments in this direction), which might be
sensitive to the signals discussed in this work. To make more accurate
statements about this possibility, one would however need to go beyond the
small $k/\bar{T}$ expansion used here to study the behavior of
$h^{2}\Omega_{\rm{GW}}(f)$ in the high frequency limit.

As a last point, we check the validity of the small $k/\bar{T}$ expansion,
that we have used to derive our analytic results in view of the experimentally
measurable frequency ranges illustrated in figure 1. As shown in this plot,
typical operating or planned GW detectors can measure frequency ranges up to
about $10^{4}\,{\rm Hz}$. Let us translate this into the corresponding range
of $k/\bar{T}$. Using eq. (4.52), we have

$\displaystyle\frac{k}{\bar{T}}$ $\displaystyle=$ $\displaystyle
2\pi\frac{a(\tau_{0})}{a(\tau)}\frac{f}{T}$ (4.65) $\displaystyle=$
$\displaystyle
2\pi\left(\frac{\rm{Hz}}{\rm{GeV}}\right)\left(\frac{T(\tau_{0})}{\rm{GeV}}\right)^{-1}\left(\frac{g_{\ast
s}(T)}{g_{\ast s}(T_{0})}\right)^{1/3}\left(\frac{f}{\rm{Hz}}\right)$
$\displaystyle=$ $\displaystyle 1.8\times 10^{-11}\left(\frac{g_{\ast
s}(T)}{g_{\ast s}(T_{0})}\right)^{1/3}\left(\frac{f}{\rm{Hz}}\right),$

where we have in the second line used eq. (4.56). Due to the factor involving
the degrees of freedom $g_{\ast s}(T)$, the coefficient of $f/\rm{Hz}$ weakly
depends on $\tau$. This factor takes values ranging from about 3 above
$200\,{\rm GeV}$, down to 1 at eV scales. For whatever time or temperature is
considered, the ratio $k/\bar{T}$ is small for frequencies $f\lesssim
10^{10}\,{\rm Hz}$. Hence, for frequencies shown in figure 1 this requirement
is satisfied, which shows that the small $k/\bar{T}$ expansion is valid. The
results of this work do not allow us to make predictions about what happens
for frequencies of the order of $f\sim 10^{11}\,{\rm Hz}$ or larger due to the
adopted approximation. It is, however, possible to make the educated guess
that because of the thermal distribution factors appearing in eq. (3.32), the
GW spectrum will be exponentially suppressed roughly as ${\rm e}^{-k/\bar{T}}$
in the high frequency limit. This exponential suppression can be expected to
become dominant for $f\gtrsim 10^{11}\,{\rm Hz}$. We note that a similar
behavior was obtained in the explicit calculations of refs. [31, 47].

##  5 Summary and outlook

We have studied the birefringence of gravitational waves triggered by
primordial left-handed neutrinos in thermal equilibrium. A generic expression
for correlation functions of gravitational fluctuations with different
helicities due to chiral fermions at finite temperature and density is derived
and further applied to analyze the Stokes parameters of gravitational waves
induced by thermalized neutrinos in the early Universe. We investigated both
the cases of Minkowski and expanding spacetime backgrounds and obtained a non-
vanishing $V$ parameter characterizing gravitational birefringence for non-
zero neutrino chemical potentials. In the expanding universe, the general
expression for the parity-violating parameter $\Delta\chi$ is given by eq.
(4.44).

While we have in this work focused on the thermalized neutrinos in the early
Universe, our argument and derivation may be extended to the thermalized
neutrino as chiral matter in core-collapse supernovae [57], which is expected
to produce polarized gravitational waves as well. More generically, chiral
gravitational waves from neutrinos or other chiral fermions can also be
generated in out of equilibrium situations both in the early Universe and
supernovae. To study such situations, one needs to use the radiation transport
theory taking into account the effects of chirality of fermions, termed the
chiral radiation transport theory [58], in which the energy-momentum tensor
with chiral corrections is expressed in terms of the distribution function of
chiral fermions. These questions are deferred to future work.

## Acknowledgments

We thank Kohei Kamada, Jun’ya Kume, Kin-Wang Ng, and Yusuke Yamada for useful
discussions and comments. P.G. is supported by the Grant-in Aid for Scientific
Research (C) (JSPS KAKENHI Grant Number JP20K03940) and the Leading Initiative
for Excellent Young Researchers (LEADER) of the Japan Society for the
Promotion of Science (JSPS). N.Y. is supported by the Keio Institute of Pure
and Applied Sciences (KiPAS) project at Keio University and JSPS KAKENHI Grant
Number JP19K03852. D.-L.Y. is supported by the National Science and Technology
Council (Taiwan) under Grant Number MOST 110-2112-M-001-070-MY3.

##  Appendix A Derivation of eqs. (3.32)–(3.36)

In this appendix, we demonstrate all necessary steps needed to obtain eqs.
(3.32)–(3.36). For simplicity of notation, we here perform computations
without rescaling dimensionful quantities with respect to the scaling factor
$a$ (which can be done later). We start from eq. (3.16) and replace the vacuum
expectation value with its thermal counterpart of eq. (3.27). First, let us
define

$\displaystyle\hat{M}_{\lambda\lambda^{\prime}}$ $\displaystyle\equiv$
$\displaystyle
L^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})\epsilon^{\mu\nu}_{-\lambda}(\bm{k})\big{[}\bar{\sigma}_{\mu}p_{\nu}+\bar{\sigma}_{\mu}(p-k)_{\nu}\big{]}L(\bm{p},\tau^{\prime})$
(A.1) $\displaystyle\times
L^{\dagger}({\bm{p}^{\prime}-\bm{k}^{\prime}},\tau^{\prime\prime})\epsilon^{\rho\sigma}_{-\lambda^{\prime}}(\bm{k}^{\prime})\big{[}\bar{\sigma}_{\rho}p^{\prime}_{\sigma}+\bar{\sigma}_{\rho}(p^{\prime}-k^{\prime})_{\sigma}\big{]}L(\bm{p}^{\prime},\tau^{\prime\prime}),$

and compute $\mathrm{Tr}\,(\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho})$
according to eqs. (3.27) and (3.28). Throughout the following discussion,
trivial terms that lead to $\bm{k}=\bm{k}^{\prime}=\bm{0}$ will be neglected.

Considering first the pure vacuum term in the trace of eq. (3.28), $\langle
0|\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho}|0\rangle$, the non-vanishing
operator combination reads

$\displaystyle\Omega_{0}\equiv\langle
0|b_{\bm{k}-\bm{p}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}|0\rangle=(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\frac{1}{Z}{\rm
e}^{\beta(E_{q}+\mu)}\,.$ (A.2)

Next, for $\langle
0|a_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho}a_{\bm{q}}^{\dagger}|0\rangle$,
the non-trivial terms are given by

$\displaystyle\Omega^{\text{th}(1)}_{a}$ $\displaystyle\equiv$
$\displaystyle\langle
0|a_{\bm{q}}a^{\dagger}_{\bm{p}-\bm{k}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}a_{\bm{p}^{\prime}}\hat{\rho}a^{\dagger}_{\bm{q}}|0\rangle$
(A.3) $\displaystyle=$
$\displaystyle(2\pi)^{9}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{p}^{\prime}-\bm{p}+\bm{k}})\delta^{3}({\bm{q}-\bm{p}^{\prime}})\frac{1}{Z}{\rm
e}^{\beta\left[(E_{q}+\mu)-(E_{q}-\mu)\right]}$

and

$\displaystyle\Omega^{(2)}_{a}$ $\displaystyle\equiv$ $\displaystyle\langle
0|a_{\bm{q}}b_{\bm{k}-\bm{p}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}a^{\dagger}_{\bm{q}}|0\rangle$
(A.4) $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta({\bm{k}-\bm{p}+\bm{p}^{\prime}})\langle
0|a_{\bm{q}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}a_{\bm{q}}^{\dagger}|0\rangle\frac{1}{Z}{\rm
e}^{\beta\left[(E_{q}+\mu)-(E_{q}-\mu)\right]}$ $\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\Big{[}\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}}){\cal
V}-(2\pi)^{3}\delta^{3}({\bm{q}-\bm{p}^{\prime}+\bm{k}^{\prime}})\delta^{3}({\bm{p}-\bm{q}})\Big{]}$
$\displaystyle\times\frac{1}{Z}{\rm
e}^{\beta\left[(E_{q}+\mu)-(E_{q}-\mu)\right]}\,.$

Note that part of eq. (A.4) is similar to eq. (A.2), which does not correlate
${\bm{q}}$ and other momenta. Such terms actually contribute to the vacuum
expectation value in the end. We hence separate the vacuum and thermal parts
as $\Omega^{(2)}_{a}=\Omega^{\text{vac}}_{a}+\Omega^{\text{th}(2)}_{a}$, where

$\displaystyle\Omega^{\text{vac}}_{a}$
$\displaystyle=(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}}){\cal
V}\frac{1}{Z}{\rm e}^{\beta\left[(E_{q}+\mu)-(E_{q}-\mu)\right]}\,,$ (A.5)
$\displaystyle\Omega^{\text{th}(2)}_{a}$
$\displaystyle=-(2\pi)^{9}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{p}^{\prime}-\bm{p}+\bm{k}})\delta^{3}({\bm{q}-\bm{p}})\frac{1}{Z}{\rm
e}^{\beta\left[(E_{q}+\mu)-(E_{q}-\mu)\right]}\,.$ (A.6)

For $\langle
0|b_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho}b_{\bm{q}}^{\dagger}|0\rangle$,
we find

$\displaystyle\Omega^{\text{th}(1)}_{b}$ $\displaystyle\equiv$
$\displaystyle\langle
0|b_{\bm{q}}b_{\bm{k}-\bm{p}}b^{\dagger}_{-\bm{p}}b_{\bm{k}^{\prime}-\bm{p}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}|0\rangle$
(A.7) $\displaystyle=$ $\displaystyle\frac{1}{Z}\langle
0|b_{\bm{k}-\bm{p}}\big{[}(2\pi)^{3}\delta^{3}({\bm{q}+\bm{p}})-b^{\dagger}_{-\bm{p}}b_{\bm{q}}\big{]}\big{[}(2\pi)^{3}\delta^{3}({\bm{k}^{\prime}-\bm{p}^{\prime}-\bm{q}})-b^{\dagger}_{\bm{q}}b_{\bm{k}^{\prime}-\bm{p}^{\prime}}\big{]}b^{\dagger}_{-\bm{p}^{\prime}}|0\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{Z}(2\pi)^{9}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\delta^{3}({\bm{q}+\bm{p}})\,.$

Furthermore,

$\displaystyle\Omega^{(2)}_{b}\equiv\langle
0|b_{\bm{q}}b_{\bm{k}-\bm{p}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}|0\rangle=\Omega^{\text{vac}}_{b}+\Omega^{\text{th}(2)}_{b},$
(A.8)

where

$\displaystyle\Omega^{\text{vac}}_{b}$
$\displaystyle=\frac{1}{Z}(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}}){\cal
V}\,,$ (A.9) $\displaystyle\Omega^{\text{th}(2)}_{b}$
$\displaystyle=-\frac{1}{Z}(2\pi)^{9}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\delta^{3}({\bm{q}+\bm{p}^{\prime}})\,.$
(A.10)

Finally, for $\langle
0|a_{\bm{q}}b_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho}b_{\bm{q}}^{\dagger}a_{\bm{q}}^{\dagger}|0\rangle$,
the first two non-vanishing and non-trivial terms are

$\displaystyle\Omega^{\text{th}(1)}_{ab}$ $\displaystyle\equiv$
$\displaystyle\langle
0|a_{\bm{q}}b_{\bm{q}}a^{\dagger}_{\bm{p}-\bm{k}}b^{\dagger}_{-\bm{p}}b_{\bm{k}^{\prime}-\bm{p}^{\prime}}a_{\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}a^{\dagger}_{\bm{q}}|0\rangle$
(A.11) $\displaystyle=$
$\displaystyle(2\pi)^{12}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-2\bm{p}})\delta^{3}({\bm{q}-\bm{p}^{\prime}})\delta^{3}({\bm{p}+\bm{p}^{\prime}})\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}$

and

$\displaystyle\Omega^{(2)}_{ab}$ $\displaystyle\equiv$ $\displaystyle\langle
0|a_{\bm{q}}b_{\bm{q}}b_{\bm{k}-\bm{p}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}a^{\dagger}_{\bm{q}}|0\rangle$
(A.12) $\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\Big{[}{\cal
V}^{2}-{\cal V}(2\pi)^{3}\delta^{3}({\bm{p}-\bm{q}})-{\cal
V}(2\pi)^{3}\delta({\bm{q}+\bm{p}^{\prime}})$
$\displaystyle+(2\pi)^{6}\delta^{3}({\bm{q}-\bm{p}})\delta^{3}({\bm{q}-\bm{p}^{\prime}-\bm{k}^{\prime}})\Big{]}\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}\,,$

which lead to
$\Omega^{(2)}_{ab}=\Omega^{\text{vac}}_{ab}+\Omega^{\text{th}(2)}_{ab}$, where

$\displaystyle\Omega^{\text{vac}}_{ab}$
$\displaystyle=(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}}){\cal
V}^{2}\frac{1}{Z}{\rm e}^{-\beta(E_{q}-\mu)}\,,$ (A.13)
$\displaystyle\Omega^{\text{th}(2)}_{ab}$
$\displaystyle=(2\pi)^{12}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}+2\bm{q}})\delta^{3}({\bm{k}-\bm{p}-\bm{q}})\delta^{3}({\bm{q}-\bm{p}})\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}$ $\displaystyle\quad-(2\pi)^{9}{\cal
V}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\Big{[}\delta^{3}({\bm{p}-\bm{q}})+\delta^{3}({\bm{p}^{\prime}+\bm{q}})\Big{]}\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}\,.$ (A.14)

Further non-vanishing terms are found as

$\displaystyle\Omega^{\text{th}(3)}_{ab}$ $\displaystyle\equiv$
$\displaystyle\langle
0|a_{\bm{q}}b_{\bm{q}}a^{\dagger}_{\bm{p}-\bm{k}}a_{\bm{p}}a^{\dagger}_{\bm{p}^{\prime}-\bm{k}^{\prime}}a_{\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}a^{\dagger}_{\bm{q}}|0\rangle$
(A.15) $\displaystyle=$ $\displaystyle(2\pi)^{9}{\cal
V}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\delta^{3}({\bm{p}^{\prime}-\bm{q}})\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}\,,$

and

$\displaystyle\Omega^{\text{th}(4)}_{ab}$ $\displaystyle\equiv$
$\displaystyle\langle
0|a_{\bm{q}}b_{\bm{q}}b_{\bm{k}-\bm{p}}b^{\dagger}_{-\bm{p}}b_{\bm{k}^{\prime}-\bm{p}^{\prime}}b^{\dagger}_{-\bm{p}^{\prime}}\hat{\rho}b^{\dagger}_{\bm{q}}a^{\dagger}_{\bm{q}}|0\rangle$
(A.16) $\displaystyle=$ $\displaystyle(2\pi)^{9}{\cal
V}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}})\delta^{3}({\bm{p}+\bm{q}})\frac{1}{Z}{\rm
e}^{-\beta(E_{q}-\mu)}\,.$

It is clear that the combination of eqs. (A.2), (A.5), (A.9) and (A.13) yields
the vacuum expectation value in eq. (3.17),

$\displaystyle\Omega_{0}+\int_{\bm{q}}\Omega^{\text{vac}}_{a}+\int_{\bm{q}}\Omega^{\text{vac}}_{b}+\frac{1}{\cal
V}\int_{\bm{q}}\Omega^{\text{vac}}_{ab}=(2\pi)^{6}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\delta^{3}({\bm{k}-\bm{p}+\bm{p}^{\prime}}).$
(A.17)

The wave-function parts for all these terms take the same form and thus their
full contribution just matches the result in vacuum. We can accordingly
decompose

$\displaystyle
M_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})=M^{\text{vac}}_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})+M^{\text{th}}_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime}),$
(A.18)

where
$M^{\text{vac}}_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})$
takes the same form as eq. (LABEL:M_in_beta). In the end, only
$\Omega^{\text{th}(1)}_{a/b}$, $\Omega^{\text{th}(2)}_{a/b}$,
$\Omega^{\text{th}(3,4)}_{ab}$, and the term linear to ${\cal V}$ in
$\Omega^{\text{th}(2)}_{ab}$ contribute to the thermal expectation part,
$M^{\text{th}}_{\lambda\lambda^{\prime}}(\bm{k},\bm{k}^{\prime};\tau^{\prime},\tau^{\prime\prime})$,
while the contributions from $\Omega^{\text{th}(1)}_{ab}$ and the remaining
part of $\Omega^{\text{th}(2)}_{ab}$ vanish as will be shown later. We first
compute the thermal contributions from $\Omega^{\text{th}(1)}_{a}$ and
$\Omega^{\text{th}(2)}_{a}$,

$\displaystyle\int_{\bm{q}}\langle
0|a_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}\hat{\rho}a^{\dagger}_{\bm{q}}|0\rangle_{\text{th}}$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\int_{\bm{p}}\Big{[}u^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})u(\bm{p},\tau^{\prime})u^{\dagger}(\bm{p},\tau^{\prime\prime})u({\bm{p}-\bm{k}},\tau^{\prime\prime})\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
(A.19)
$\displaystyle+v^{\dagger}({\bm{k}-\bm{p}},\tau^{\prime})u(\bm{p},\tau^{\prime})u^{\dagger}(\bm{p},\tau^{\prime\prime})v({\bm{k}-\bm{p}},\tau^{\prime\prime})\chi^{(a2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})\Big{]}.$

The specific expressions for
$\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ and
$\chi^{(a2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ can be derived in
analogy to the computation for
$\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ in vacuum. By definition,
$\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ takes the form

$\displaystyle\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{1}{4}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\xi^{\dagger}_{-}({\bm{p}-\bm{k}})\epsilon^{\mu\nu}_{-\lambda}(\bm{k})\big{[}\bar{\sigma}_{\mu}p_{\nu}+\bar{\sigma}_{\mu}(p-k)_{\nu}\big{]}\xi_{-}(\bm{p})$
(A.20)
$\displaystyle\times\xi_{-}^{\dagger}(\bm{p})\epsilon^{\rho\sigma}_{-\lambda^{\prime}}(-\bm{k})\big{[}\bar{\sigma}_{\rho}(p-k)_{\sigma}+\bar{\sigma}_{\rho}p_{\sigma}\big{]}\xi_{-}({\bm{p}-\bm{k}}).$

Using the properties of $\xi_{-}$ and $\epsilon_{\lambda}^{\mu}$, we have

$\displaystyle\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}\Bigg{[}\frac{2\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}}{p|{\bm{k}-\bm{p}}|}$ (A.21)
$\displaystyle+\delta_{\lambda,\lambda^{\prime}}\bigg{(}1-\frac{p^{2}-{\bm{k}\cdot\bm{p}}}{p|{\bm{k}-\bm{p}}|}\bigg{)}-{\rm
i}\epsilon_{\mu\beta\rho\alpha}\epsilon^{\mu}_{-\lambda}(\bm{k})\epsilon^{\rho}_{\lambda^{\prime}}(\bm{k})n^{\alpha}({\bm{p}-\bm{k}})n^{\beta}(\bm{p})\Bigg{]}.$

Using the same fixed coordinate system as in the vacuum calculations, we
eventually derive

$\displaystyle\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{p^{2}}{4}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\sin^{2}\theta\big{[}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{]}$ (A.22)
$\displaystyle\times\Bigg{\\{}\delta_{\lambda,\lambda^{\prime}}-\frac{\lambda+\lambda^{\prime}}{2}\cos\theta+\frac{1}{|{\bm{k}-\bm{p}}|}\bigg{(}\delta_{\lambda,\lambda^{\prime}}\big{(}k\cos\theta-p\big{)}+\frac{\lambda+\lambda^{\prime}}{2}\big{(}p\cos\theta-k\big{)}$
$\displaystyle+p\sin^{2}\theta\big{[}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{]}\bigg{)}\Bigg{\\}}\,,$

which yields

$\displaystyle\chi^{(a1)}_{\lambda\lambda}({\bm{k},\bm{p}})$
$\displaystyle=\frac{p^{2}}{4}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\sin^{2}\theta\Bigg{[}1-\lambda\cos\theta-(\cos\theta-\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}\,,$
(A.23) $\displaystyle\chi^{(a1)}_{\lambda,-\lambda}({\bm{k},\bm{p}})$
$\displaystyle=\frac{p^{3}}{4}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\frac{\sin^{4}\theta\big{(}\cos
4\phi-{\rm i}\lambda\sin 4\phi\big{)}}{|{\bm{k}-\bm{p}}|}\,.$ (A.24)

Since $|{\bm{k}-\bm{p}}|$ is independent of $\phi$, only
$\chi^{(a1)}_{\lambda\lambda}({\bm{k},\bm{p}})$ contributes to the correlation
function, in analogy to the vacuum case. We hence have

$\displaystyle\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\delta_{\lambda\lambda^{\prime}}\frac{p^{2}}{4}(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}\sin^{2}\theta$
(A.25)
$\displaystyle\times\Bigg{[}1-\lambda\cos\theta-(\cos\theta-\lambda)\frac{p\cos\theta-k}{\sqrt{p^{2}-2pk\cos\theta+k^{2}}}\Bigg{]}.$

$\chi^{(a2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ is obtained as

$\displaystyle\chi^{(a2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle-(1-f_{p+\mu})f_{p-\mu}\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}}).$
(A.26)

Next, we evaluate the thermal contribution from $\Omega^{\text{th}}_{b}$,
which gives

$\displaystyle\int_{\bm{q}}\langle
0|b_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}b^{\dagger}_{\bm{q}}|0\rangle_{\text{th}}$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\int_{\bm{p}}\Big{[}v^{\dagger}({\bm{k}-\bm{p}},\tau^{\prime})v(-\bm{p},\tau^{\prime})v^{\dagger}(-\bm{p},\tau^{\prime\prime})v({\bm{k}-\bm{p}},\tau^{\prime\prime})\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
(A.27)
$\displaystyle+v({\bm{k}-\bm{p}},\tau^{\prime})u(\bm{p},\tau^{\prime})u^{\dagger}(\bm{p},\tau^{\prime\prime})v^{\dagger}({\bm{k}-\bm{p}},\tau^{\prime\prime})\chi^{(b2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})\Big{]}.$

$\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ here is defined as

$\displaystyle\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{1}{4}(1-f_{p-\mu})f_{p+\mu}\xi^{\dagger}_{-}({\bm{k}-\bm{p}})\epsilon^{\mu\nu}_{-\lambda}(\bm{k})\big{[}\bar{\sigma}_{\mu}p_{\nu}+\bar{\sigma}_{\mu}(p-k)_{\nu}\big{]}\xi_{-}(-\bm{p})$
(A.28)
$\displaystyle\times\xi_{-}^{\dagger}(-\bm{p})\epsilon^{\rho\sigma}_{-\lambda^{\prime}}(-\bm{k})\big{[}\bar{\sigma}_{\rho}(p-k)_{\sigma}+\bar{\sigma}_{\rho}p_{\sigma}\big{]}\xi_{-}({\bm{k}-\bm{p}}).$

Following a computation similar to that for
$\chi_{\lambda\lambda^{\prime}}^{(a1)}({\bm{k},\bm{p}})$, we find

$\displaystyle\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{1}{2}(1-f_{p-\mu})f_{p+\mu}\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}\Bigg{[}\frac{2\big{(}\epsilon_{-\lambda}(\bm{k})\cdot
p\big{)}\big{(}\epsilon_{\lambda^{\prime}}(\bm{k})\cdot
p\big{)}}{p|{\bm{k}-\bm{p}}|}$ (A.29)
$\displaystyle+\delta_{\lambda,\lambda^{\prime}}\bigg{(}1-\frac{p^{2}-{\bm{k}\cdot\bm{p}}}{p|{\bm{k}-\bm{p}}|}\bigg{)}-{\rm
i}\epsilon_{\mu\beta\rho\alpha}\epsilon^{\mu}_{-\lambda}(\bm{k})\epsilon^{\rho}_{\lambda^{\prime}}(\bm{k})n^{\alpha}({\bm{k}-\bm{p}})n^{\beta}(-\bm{p})\Bigg{]}\,.$

Employing again the same fixed coordinate system as in the vacuum calculation,
we eventually derive

$\displaystyle\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\frac{p^{2}}{4}(1-f_{p-\mu})f_{p+\mu}\sin^{2}\theta\big{[}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{]}$ (A.30)
$\displaystyle\times\Bigg{[}\delta_{\lambda,\lambda^{\prime}}+\frac{\lambda+\lambda^{\prime}}{2}\cos\theta+\frac{1}{|{\bm{k}-\bm{p}}|}\bigg{(}\delta_{\lambda,\lambda^{\prime}}\big{(}k\cos\theta-p\big{)}+\frac{\lambda+\lambda^{\prime}}{2}\big{(}k-p\cos\theta\big{)}$
$\displaystyle+p\sin^{2}\theta\big{[}\cos^{2}\phi+\lambda\lambda^{\prime}\sin^{2}\phi+{\rm
i}(\lambda^{\prime}-\lambda)\cos\phi\sin\phi\big{]}\bigg{)}\Bigg{]},$

which yields

$\displaystyle\chi^{(b1)}_{\lambda\lambda}({\bm{k},\bm{p}})$
$\displaystyle=\frac{p^{2}}{4}(1-f_{p-\mu})f_{p+\mu}\sin^{2}\theta\Bigg{[}1+\lambda\cos\theta-(\cos\theta+\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}\,,$
(A.31) $\displaystyle\chi^{(b1)}_{\lambda,-\lambda}({\bm{k},\bm{p}})$
$\displaystyle=\frac{p^{3}}{4}(1-f_{p-\mu})f_{p+\mu}\frac{\sin^{4}\theta\big{(}\cos
4\phi-{\rm i}\lambda\sin 4\phi\big{)}}{|{\bm{k}-\bm{p}}|}\,.$ (A.32)

Again, only $\chi^{(b1)}_{\lambda\lambda}({\bm{k},\bm{p}})$ contributes to the
correlation function, leading to

$\displaystyle\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle\delta_{\lambda\lambda^{\prime}}\frac{p^{2}}{4}(1-f_{p-\mu})f_{p+\mu}\sin^{2}\theta\Bigg{[}1+\lambda\cos\theta-(\cos\theta+\lambda)\frac{p\cos\theta-k}{|{\bm{k}-\bm{p}}|}\Bigg{]}\,.$
(A.33)

Furthermore, $\chi_{\lambda\lambda^{\prime}}^{(b2)}({\bm{k},\bm{p}})$ is
obtained as

$\displaystyle\chi^{(b2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=$
$\displaystyle-(1-f_{|{\bm{k}-\bm{p}}|-\mu})f_{|{\bm{k}-\bm{p}}|+\mu}\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}}).$
(A.34)

Finally, the contribution from $\Omega^{\text{th}(1-4)}_{ab}$ reads

$\displaystyle\frac{1}{\cal V}\int_{\bm{q}}\langle
0|a_{\bm{q}}b_{\bm{q}}\hat{M}_{\lambda\lambda^{\prime}}b_{\bm{q}}^{\dagger}a^{\dagger}_{\bm{q}}|0\rangle_{\text{th}}$
$\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}({\bm{k}+\bm{k}^{\prime}})\int_{\bm{p}}\Big{[}u^{\dagger}({\bm{p}-\bm{k}},\tau^{\prime})v(-\bm{p},\tau^{\prime})v^{\dagger}(-\bm{p},\tau^{\prime\prime})u({\bm{p}-\bm{k}},\tau^{\prime\prime})\chi^{(ab1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
(A.35)
$\displaystyle+v^{\dagger}({\bm{k}-\bm{p}},\tau^{\prime})u(\bm{p},\tau^{\prime})u^{\dagger}(\bm{p},\tau^{\prime\prime})v({\bm{k}-\bm{p}},\tau^{\prime\prime})\chi^{(ab2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})\Big{]},$

where

$\displaystyle\chi^{(ab1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=\frac{f_{|{\bm{k}-\bm{p}}|-\mu}f_{|{\bm{k}-\bm{p}}|+\mu}\chi^{(a1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})}{(1-f_{|{\bm{k}-\bm{p}}|+\mu})f_{|{\bm{k}-\bm{p}}|-\mu}}+\frac{f_{p-\mu}f_{p+\mu}\chi^{(b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})}{(1-f_{p-\mu})f_{p+\mu}}\,,$
(A.36) $\displaystyle\chi^{(ab2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$
$\displaystyle=-\big{(}f_{|{\bm{k}-\bm{p}}|-\mu}f_{|{\bm{k}-\bm{p}}|+\mu}+f_{p-\mu}f_{p+\mu}\big{)}\beta_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})\,.$
(A.37)

Here $\chi^{(ab1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ comes from
$\Omega^{\text{th}(3,4)}_{ab}$, while
$\chi^{(ab2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ originates from the
terms linear in ${\cal V}$ in $\Omega^{\text{th}(2)}_{ab}$.

In the conformal coordinates, the wave functions reduce to plane waves such
that $u(\bm{p},\tau)={\rm e}^{-{\rm i}p\tau}$ and $v(\bm{p},\tau)={\rm
e}^{{\rm i}p\tau}$. It is then easy to recognize that
$\chi^{(ab1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ and
$\chi^{(ab2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ take similar forms
as $\chi^{(a,b1)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$ and
$\chi^{(a,b2)}_{\lambda\lambda^{\prime}}({\bm{k},\bm{p}})$, respectively.
Furthermore, one finds that the contributions from
$\Omega^{\text{th}(1)}_{ab}$ and the first term of
$\Omega^{\text{th}(2)}_{ab}$ vanish due to their on-shell conditions,
$\bm{p}=\bm{k}/2$ and $\bm{p}=-\bm{k}/2$, respectively.

Adding all non-vanishing terms and replacing the wave functions with the
abovementioned plane wave solutions, we obtain eqs. (3.32)–(3.36).

##  Appendix B Derivation of eqs. (4.20) and (4.21)

We here present the details of the calculations related to the $\tau^{\prime}$
integrals. Substituting eqs. (4.14) and (4.15) into eq. (4.16) and rearanging
some terms, we obtain

$\displaystyle\bar{\mathcal{I}}(\tau,\tau_{\rm i},b)=\frac{1}{2{\rm
i}}\frac{\tau_{\rm i}^{2}}{a_{\rm i}^{2}k\tau}\left[{\rm e}^{{\rm i}\tau
k}\int_{\tau_{\rm i}}^{\tau}\frac{{\rm d}\tau^{\prime}}{\tau^{\prime}}{\rm
e}^{{\rm i}\tau^{\prime}(b-k)}-{\rm e}^{-{\rm i}\tau k}\int_{\tau_{\rm
i}}^{\tau}\frac{{\rm d}\tau^{\prime}}{\tau^{\prime}}{\rm e}^{{\rm
i}\tau^{\prime}(b+k)}\right],$ (B.1)

which will be the starting point for the following considerations.

###  B.1 $b=b_{-}$ case

To derive eq. (4.20), we first expand $|\bm{k}-\bm{p}|$, which gives

$\displaystyle|\bm{k}-\bm{p}|\simeq
p\left(1-\frac{k}{p}\cos\theta+\frac{1}{2}\frac{k^{2}}{p^{2}}\sin^{2}\theta\right),$
(B.2)

where $\theta$ is the angle between $\bm{k}$ and $\bm{p}$. Thus

$\displaystyle{\rm e}^{{\rm i}\tau^{\prime}(b_{-}\mp k)}\simeq 1-{\rm
i}\tau^{\prime}k(\cos\theta\pm 1)+\frac{\rm
i}{2}\tau^{\prime}\frac{k^{2}}{p}\sin^{2}\theta-\frac{1}{2}\tau^{\prime
2}k^{2}(\cos\theta\pm 1)^{2}\,.$ (B.3)

Hence, we have

$\displaystyle\int_{\tau_{\rm i}}^{\tau}\frac{{\rm
d}\tau^{\prime}}{\tau^{\prime}}{\rm e}^{{\rm i}\tau^{\prime}(b_{-}\mp k)}$
$\displaystyle\simeq$ $\displaystyle\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-{\rm i}k(\tau-\tau_{\rm i})\Bigl{(}\cos\theta\pm
1-\frac{1}{2}\frac{k}{p}\sin^{2}\theta\Bigr{)}$ (B.4)
$\displaystyle-\frac{1}{4}(\tau^{2}-\tau^{2}_{0})k^{2}(\cos\theta\pm
1)^{2}\,.$

Note that because the LO term cancels in eq. (B.1), terms up to second order
in $k$ have to be kept up to this point of the calculation. Substituting the
above results into eq. (B.1) and expanding the factors ${\rm e}^{{\rm i}\tau
k}$ and ${\rm e}^{-{\rm i}\tau k}$ in $k$, we obtain

$\displaystyle\bar{\mathcal{I}}(\tau,\tau_{\rm i},b_{-})\simeq\frac{\tau_{\rm
i}^{2}}{a_{\rm i}^{2}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}-\frac{\rm
i}{2}k\frac{(\tau-\tau_{\rm i})^{2}}{\tau}\right]+\mathcal{O}(k^{2})\,.$ (B.5)

As the zeroth-order term in $k$ is real while the first-order term is
imaginary, the latter becomes a second-order contribution to
$|\bar{\mathcal{I}}(\tau,\tau_{\rm i},b_{-})|^{2}$ and can hence be neglected.
Therefore, we have

$\displaystyle|\bar{\mathcal{I}}(\tau,\tau_{\rm
i},b_{-})|^{2}\simeq\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}+\mathcal{O}(k^{2})\,.$
(B.6)

###  B.2 $b=b_{+}$ case

Using eq. (B.2), we have

$\displaystyle{\rm e}^{{\rm i}\tau^{\prime}(b_{+}\mp k)}\simeq{\rm e}^{2{\rm
i}\tau^{\prime}p}\left[1-{\rm i}\tau^{\prime}k(\cos\theta\pm 1)+\frac{\rm
i}{2}\tau^{\prime}\frac{k^{2}}{p}\sin^{2}\theta-\frac{1}{2}\tau^{\prime
2}k^{2}(\cos\theta\pm 1)^{2}\right],$ (B.7)

which leads to

$\displaystyle\int_{\tau_{\rm i}}^{\tau}\frac{{\rm
d}\tau^{\prime}}{\tau^{\prime}}{\rm e}^{{\rm i}\tau^{\prime}(b_{+}\mp k)}$
$\displaystyle\simeq$ $\displaystyle\mathrm{Ei}(2{\rm i}\tau
p)-\mathrm{Ei}(2{\rm i}\tau_{\rm i}p)+\frac{k}{p}{\rm e}^{{\rm
i}(\tau+\tau_{\rm i})p}\Biggl{\\{}-{\rm i}\sin[(\tau-\tau_{\rm
i})p](\cos\theta\pm 1)$ (B.8) $\displaystyle+\frac{\rm
i}{2}\frac{k}{p}\sin[(\tau-\tau_{\rm
i})p]\sin^{2}\theta-\frac{1}{4}\frac{k}{p}\Bigl{(}({\rm i}+(\tau+\tau_{\rm
i})p)\sin[(\tau-\tau_{\rm i})p]$ $\displaystyle-{\rm i}(\tau-\tau_{\rm
i})p\cos[(\tau-\tau_{\rm i})p]\Bigr{)}(\cos\theta\pm 1)^{2}\Biggr{\\}}\,,$

where $Ei(x)$ is the exponential integral function. Substituting the above
results into eq. (B.1) and expanding the factors ${\rm e}^{{\rm i}\tau k}$ and
${\rm e}^{-{\rm i}\tau k}$ in $k$, we obtain

$\displaystyle\bar{\mathcal{I}}(\tau,\tau_{\rm i},b_{+})$
$\displaystyle\simeq$ $\displaystyle\frac{\tau_{\rm i}^{2}}{a_{\rm
i}^{2}}\left[\mathrm{Ei}(2{\rm i}\tau p)-\mathrm{Ei}(2{\rm i}\tau_{\rm
i}p)-\frac{1}{\tau p}{\rm e}^{{\rm i}(\tau+\tau_{\rm i})p}\sin[(\tau-\tau_{\rm
i})p](1+{\rm i}\tau k\cos\theta)\right]$ (B.9)
$\displaystyle+\mathcal{O}(k^{2}).$

For the next step, we make use of

$\displaystyle\mathrm{Ei}({\rm i}x)=\mathrm{ci}(x)+{\rm i}\mathrm{Si}(x),$
(B.10)

for $x$ real and positive, where $\mathrm{ci}(x)$ and $\mathrm{Si}(x)$ are the
cosine and sine integral functions defined in eqs. (4.22) and (4.23),
respectively. Via eq. (B.10), we can now divide
$\bar{\mathcal{I}}(\tau,\tau_{\rm i},b)$ into its real and imaginary parts,
and hence obtain

$\displaystyle|\bar{\mathcal{I}}(\tau,\tau_{\rm i},b_{+})|^{2}$
$\displaystyle\simeq$ $\displaystyle\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\Biggl{\\{}\Bigl{(}\mathrm{ci}(2\tau p)-\mathrm{ci}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\cos[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}$ (B.11) $\displaystyle+\Bigl{(}\mathrm{Si}(2\tau
p)-\mathrm{Si}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\sin[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle+2\frac{k}{p}\cos\theta\sin[(\tau-\tau_{\rm
i})p]\Bigr{(}\sin[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)\bigr{]}$ $\displaystyle-\cos[(\tau+\tau_{\rm
i})p]\bigl{[}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)\bigr{]}\Bigr{)}\Biggr{\\}}+\mathcal{O}(k^{2})\,.$

##  Appendix C Derivation of eqs. (4.24)–(4.26)

We in this appendix provide the details of the angular integrals of $\phi$ and
$\theta$ of $\bm{p}$ in eq. (3.32). As the $\phi$ integral is trivial, we only
need to discuss the $\theta$ integral explicitly.

###  C.1 $\Gamma_{\pm}$ terms

Before considering the $\theta$ integrals of these terms, we first need to
expand eq. (3.35) in small $k$. Making use of eq. (B.2), the result reads

$\displaystyle\tilde{\Gamma}_{\pm}(\bm{p},\theta;\bm{k},\lambda)\simeq\sin^{4}\theta\left[1+\frac{k}{p}(\cos\theta\mp\lambda)\right]+\mathcal{O}(k^{2}).$
(C.1)

Furthermore, for the $\Gamma_{+}$ term, we also need to expand
$f_{|\bm{k}-\bm{p}|-\bar{\mu}}$, which gives

$\displaystyle f_{|\bm{k}-\bm{p}|-\bar{\mu}}\simeq
f_{p-\bar{\mu}}\left[1+\frac{k}{\bar{T}}\cos\theta(1-f_{p-\bar{\mu}})\right]+\mathcal{O}(k^{2}).$
(C.2)

Next, substituting the results of eqs. (4.20), (C.1) and (C.2) into eq.
(3.32), we obtain for the $\Gamma_{+}$ term

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{+}}$
$\displaystyle\simeq$ $\displaystyle
2\pi\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}\int_{0}^{\infty}{\rm
d}p\ p^{4}f_{p-\bar{\mu}}$ (C.3) $\displaystyle\times\int_{0}^{\pi}{\rm
d}\theta\sin^{5}\theta\left[1+\frac{k}{p}(\cos\theta-\lambda)\right]\left[1+\frac{k}{\bar{T}}\cos\theta(1-f_{p-\bar{\mu}})\right]+\mathcal{O}(k^{2})\,,$

Using

$\displaystyle\int_{0}^{\pi}{\rm d}\theta\sin^{5}\theta$ $\displaystyle=$
$\displaystyle\frac{16}{15}\,,$ (C.4) $\displaystyle\int_{0}^{\pi}{\rm
d}\theta\sin^{5}\theta\cos\theta$ $\displaystyle=$ $\displaystyle 0\,,$ (C.5)

we have

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{+}}$
$\displaystyle\simeq$
$\displaystyle\frac{32\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}\int_{0}^{\infty}{\rm
d}p\ p^{4}\left(1-\lambda\frac{k}{p}\right)f_{p-\bar{\mu}}$ (C.6)
$\displaystyle+\mathcal{O}(k^{2}).$

For the $\Gamma_{-}$ term, a similar calculation leads to

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{-}}$
$\displaystyle\simeq$
$\displaystyle\frac{32\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm
i}^{4}}\left[\log\Bigl{(}\frac{\tau}{\tau_{\rm
i}}\Bigr{)}-\frac{\tau-\tau_{\rm i}}{\tau}\right]^{2}\int_{0}^{\infty}{\rm
d}p\ p^{4}\left(1+\lambda\frac{k}{p}\right)f_{p+\bar{\mu}}$ (C.7)
$\displaystyle+\mathcal{O}(k^{2}).$

###  C.2 $\Gamma_{0}$ term

For handling the $\theta$ integral of this term, we first need to expand
$\tilde{\Gamma}_{0}$ in eq. (3.36) and $f_{|\bm{k}-\bm{p}|+\bar{\mu}}$ in
small $k$,

$\displaystyle\tilde{\Gamma}_{0}(\bm{p},\theta;\bm{k},\lambda)\simeq\sin^{2}\theta\Bigl{[}(1-\lambda\cos\theta)^{2}-\frac{k}{p}\sin^{2}\theta(\cos\theta-\lambda)\Bigr{]}+\mathcal{O}(k^{2}),$
(C.8) $\displaystyle f_{|\bm{k}-\bm{p}|+\bar{\mu}}\simeq
f_{p+\bar{\mu}}\Bigl{[}1+\frac{k}{\bar{T}}\cos\theta(1-f_{p+\bar{\mu}})\Bigr{]}+\mathcal{O}(k^{2}).$
(C.9)

As the analytic structure of this term is very involved, we will consider the
LO ($k^{0}$) and NLO ($k^{1}$) terms separately.

####  C.2.1 Leading-order contribution

Substituting the LO parts of eqs. (4.21), (C.8) and (C.9) into eq. (3.32),
using

$\displaystyle\int_{0}^{\pi}{\rm
d}\theta\sin^{3}\theta(1-\lambda\cos\theta)^{2}=\frac{8}{5}\,,$ (C.10)

and noting that the $\lambda$-dependent term vanishes, we obtain

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{LO}}$
$\displaystyle\simeq$
$\displaystyle\frac{16\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\int_{0}^{\infty}{\rm d}p\
p^{4}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})$ (C.11)
$\displaystyle\times\Biggl{\\{}\Bigl{(}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\cos[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle+\Bigl{(}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\sin[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}\Biggr{\\}}\,.$

####  C.2.2 Next-to-leading-order contribution

Collecting all terms linear in $k$ from eqs. (4.21), (C.8) and (C.9), and
performing the $\theta$ integrals

$\displaystyle\int_{0}^{\pi}{\rm
d}\theta\sin^{3}\theta\cos\theta(1-\lambda\cos\theta)^{2}$ $\displaystyle=$
$\displaystyle-\frac{8}{15}\lambda\,,$ (C.12) $\displaystyle\int_{0}^{\pi}{\rm
d}\theta\sin^{5}\theta(\cos\theta-\lambda)$ $\displaystyle=$
$\displaystyle-\frac{16}{15}\lambda\,,$ (C.13)

one gets

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{NLO}}$
$\displaystyle\simeq$
$\displaystyle-\frac{16\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\lambda k\int_{0}^{\infty}{\rm
d}p\ p^{4}$ (C.14)
$\displaystyle\times\Biggl{\\{}\frac{2}{p}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})\sin[(\tau-\tau_{\rm
i})p]\Bigr{(}\sin[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)\bigr{]}$ $\displaystyle-\cos[(\tau+\tau_{\rm
i})p]\bigl{[}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm i}p)\bigr{]}\Bigr{)}$
$\displaystyle+\left[-\frac{2}{p}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})+\frac{1}{T}f_{p+\bar{\mu}}(1-f_{p+\bar{\mu}})\right]$
$\displaystyle\times\biggl{[}\Bigl{(}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\cos[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$
$\displaystyle+\Bigl{(}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\sin[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}\biggr{]}\Biggr{\\}}\,.$

##  Appendix D Derivation of eqs. (4.27) and (4.31)

The remaining $p$ integrals in eqs. (4.24), (4.25) and (4.26) are discussed in
this appendix.

###  D.1 $\Gamma_{\pm}$ terms

Starting from eq. (4.24), we need to evaluate the following integral

$\displaystyle I_{1}(k,\pm\bar{\mu},\bar{T})=\int_{0}^{\infty}{\rm d}p\
p^{4}\left(1\pm\lambda\frac{k}{p}\right)f_{p\pm\bar{\mu}}\,.$ (D.1)

Changing the integration variable from $p$ to $y=p/\bar{T}$ and making use of
the following general integration formula

$\displaystyle\int_{0}^{\infty}{\rm d}yy^{n}\frac{1}{{\rm
e}^{y\pm\bar{\mu}/\bar{T}}+1}=-n!\mathrm{Li}_{n+1}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}}\bigr{)},$ (D.2)

where $\mathrm{Li}_{s}(z)=\sum_{k=1}^{\infty}z^{k}/k^{s}$ is the polylogarithm
function, we obtain

$\displaystyle
I_{1}(k,\pm\bar{\mu},\bar{T})=-24\bar{T}^{5}\Bigl{[}\mathrm{Li}_{5}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}}\bigr{)}\pm\frac{\lambda}{4}\frac{k}{\bar{T}}\mathrm{Li}_{4}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}}\bigr{)}\Bigr{]},$ (D.3)

which leads to eq. (4.27).

###  D.2 $\Gamma_{0}$ term

As we discussed in the main text, it is possible to recast the functions
$\mathrm{Si}(x)$ and $\mathrm{ci}(x)$ into simpler expressions for the
temperature and time regions considered in this work. By making use of eqs.
(4.29) and (4.30), the parts of the integrands involving $\mathrm{Si}(x)$ and
$\mathrm{ci}(x)$ can be simplified as

$\displaystyle\Bigl{(}\mathrm{ci}(2\tau p)-\mathrm{ci}(2\tau_{\rm
i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm i})p]\cos[(\tau+\tau_{\rm
i})p]\Bigr{)}^{2}$ $\displaystyle+\Bigl{(}\mathrm{Si}(2\tau
p)-\mathrm{Si}(2\tau_{\rm i}p)-\frac{1}{\tau p}\sin[(\tau-\tau_{\rm
i})p]\sin[(\tau+\tau_{\rm i})p]\Bigr{)}^{2}$ $\displaystyle\simeq$
$\displaystyle\frac{(\tau-\tau_{\rm i})^{2}}{4\tau_{\rm
i}^{2}\tau^{2}p^{2}}+\mathcal{O}\Bigl{(}\frac{1}{(\tau p)^{3}}\Bigr{)}\,,$
(D.4)

and

$\displaystyle\sin[(\tau+\tau_{\rm i})p]\bigl{[}\mathrm{ci}(2\tau
p)-\mathrm{ci}(2\tau_{\rm i}p)\bigr{]}-\cos[(\tau+\tau_{\rm
i})p]\bigl{[}\mathrm{Si}(2\tau p)-\mathrm{Si}(2\tau_{\rm i}p)\bigr{]}$
$\displaystyle=$ $\displaystyle-\frac{\tau-\tau_{\rm i}}{2\tau_{\rm i}\tau
p}\cos[(\tau-\tau_{\rm i})p]+\mathcal{O}\Bigl{(}\frac{1}{(\tau
p)^{2}}\Bigr{)}\,.$ (D.5)

Using eq. (D.4), the LO $\Gamma_{0}$ term of eq. (4.25) becomes

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{LO}}$
$\displaystyle\simeq$
$\displaystyle\frac{4\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\frac{(\tau-\tau_{\rm
i})^{2}}{\tau^{2}\tau_{\rm i}^{2}}\int_{0}^{\infty}{\rm d}p\
p^{2}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})$ (D.6) $\displaystyle=$
$\displaystyle\,\,-\frac{8\pi}{5}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\frac{(\tau-\tau_{\rm
i})^{2}}{\tau^{2}\tau_{\rm i}^{2}}\bar{T}^{3}$
$\displaystyle\times\Bigl{[}\mathrm{Li}_{3}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}+\mathrm{Li}_{3}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}\Bigr{]}\,,$

which corresponds to the LO term in eq. (4.31). Moreover, employing both eqs.
(D.4) and (D.5), the NLO term can be simplified as

$\displaystyle\langle
h_{\lambda}(\bm{k},\tau)h_{\lambda^{\prime}}(\bm{k}^{\prime},\tau)\rangle_{\mathrm{th},\,\Gamma_{0},\,\mathrm{NLO}}$
$\displaystyle\simeq$
$\displaystyle\frac{4\pi}{15}\delta_{\lambda,\lambda^{\prime}}\frac{\delta(\bm{k}+\bm{k}^{\prime})}{m_{\rm
P}^{4}}\frac{\tau_{\rm i}^{4}}{a_{\rm i}^{4}}\frac{\tau-\tau_{\rm
i}}{\tau\tau_{\rm i}}\lambda k\int_{0}^{\infty}{\rm d}p\ p^{2}$ (D.7)
$\displaystyle\times\Biggl{\\{}2(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})\sin[2(\tau-\tau_{\rm
i})p]+\frac{\tau-\tau_{\rm i}}{\tau_{\rm i}\tau}$
$\displaystyle\times\left[\frac{2}{p}(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})-\frac{1}{\bar{T}}f_{p+\bar{\mu}}(1-f_{p+\bar{\mu}})\right]\Biggr{\\}}\,.$

The integrals appearing in eq. (D.7) can be computed as

$\displaystyle\int_{0}^{\infty}{\rm d}p\
p^{2}f_{p\pm\bar{\mu}}\sin[2(\tau-\tau_{\rm i})p]$ $\displaystyle=$
$\displaystyle\bar{T}^{3}\int_{0}^{\infty}{\rm d}x\
x^{2}\frac{\sin[2(\tau-\tau_{\rm i})\bar{T}x]}{{\rm
e}^{x\pm\bar{\mu}/\bar{T}}+1}$ $\displaystyle=$ $\displaystyle-{\rm
i}\bar{T}^{3}{\rm e}^{\mp\bar{\mu}/\bar{T}}\Bigl{[}\Phi^{\ast}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}},3,1-2{\rm i}(\tau-\tau_{\rm
i})\bar{T}\bigr{)}-\Phi^{\ast}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}},3,1+2{\rm i}(\tau-\tau_{\rm
i})\bar{T}\bigr{)}\Bigr{]}$ $\displaystyle=$ $\displaystyle 2\bar{T}^{3}{\rm
e}^{\mp\bar{\mu}/\bar{T}}\mathrm{Im}\,\Phi^{\ast}\bigl{(}-{\rm
e}^{\mp\bar{\mu}/\bar{T}},3,1-2{\rm i}(\tau-\tau_{\rm i})\bar{T}\bigr{)},$
(D.8)

where $\Phi^{\ast}(z,s,a)$ is the Lerch transcendent function, defined in eq.
(4.32), and

$\displaystyle\int_{0}^{\infty}{\rm d}p\
p(f_{p-\bar{\mu}}+f_{p+\bar{\mu}})=-\bar{T}^{2}\Bigl{[}\mathrm{Li}_{2}\bigl{(}-{\rm
e}^{\bar{\mu}/\bar{T}}\bigr{)}+\mathrm{Li}_{2}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}\Bigr{]}\,,$ (D.9)
$\displaystyle\int_{0}^{\infty}{\rm d}p\
p^{2}f_{p+\bar{\mu}}(1+f_{p+\bar{\mu}})=-2\bar{T}^{3}\mathrm{Li}_{2}\bigl{(}-{\rm
e}^{-\bar{\mu}/\bar{T}}\bigr{)}\,,$ (D.10)

leading to the term linear in $k$ in eq. (4.31).

##  Appendix E Estimate of $\tau\bar{T}$ and $\tau_{\rm i}\bar{T}$

We here derive an order-of-magnitude estimate of $\tau\bar{T}$ and $\tau_{\rm
i}\bar{T}$ in the radiation dominated period, in which the scale factor
evolves as (see, e.g., ref. [59])

$\displaystyle a(t)=At^{1/2},$ (E.1)

where

$\displaystyle A=\left[\frac{32\pi G}{3}(\epsilon a^{4})\right]^{1/4},$ (E.2)

with $\epsilon a^{4}=$ constant. Here, $\epsilon$ is the energy density, which
is expressed at leading order of the $\mu/T$ expansion as

$\displaystyle\epsilon=g_{\ast}(T)\frac{\pi^{2}}{30}T^{4}\,,$ (E.3)

where $g_{\ast}(T)$ is the number of relativistic degrees of freedom at
temperature $T$.

From the equations above, $t$ can be related to $T$ as

$\displaystyle
t=\left(\frac{45}{16\pi^{3}Gg_{\ast}(T)}\right)^{1/2}\frac{1}{T^{2}}\,.$ (E.4)

Furthermore, we can obtain $\tau$ from ${\rm d}t=a{\rm d}\tau$ and eq. (E.1)
as

$\displaystyle\tau-\tau_{\rm i}=\int^{t}_{t_{\rm i}}\frac{{\rm
d}t}{a(t)}=\frac{2}{A}(t^{1/2}-t_{\rm i}^{1/2})\,,$ (E.5)

and hence,

$\displaystyle\tau\bar{T}-\tau_{\rm i}\bar{T}$ $\displaystyle=$ $\displaystyle
2\left(\frac{2025}{256\pi^{6}G^{2}g^{2}_{\ast}(T)}\right)^{1/4}\frac{1}{T}\left(1-\frac{\sqrt{g_{\ast}(T)}T}{\sqrt{g_{*}(T_{\rm
i})}T_{\rm i}}\right).$ (E.6)

By extracting the $\tau$-dependent component above, we find

$\displaystyle\tau\bar{T}=2\left(\frac{2025}{256\pi^{6}G^{2}g^{2}_{\ast}(T)}\right)^{1/4}\frac{1}{T}\sim
10^{20}\left(\frac{T}{\mathrm{GeV}}\right)^{-1}.$ (E.7)

We can therefore conclude that $\tau\bar{T}\gg 1$, for relevant temperatures
in the radiation dominant period.

Considering next $\tau_{\rm i}\bar{T}$, based on the $\tau$-independent part
in eq. (E.6), we can derive

$\displaystyle\tau_{\rm i}{\bar{T}}$ $\displaystyle=$ $\displaystyle
2\left[\frac{2025}{256\pi^{6}G^{2}g^{2}_{\ast}(\tau_{\rm
i})}\right]^{1/4}\frac{1}{T(\tau_{\rm i})},$ (E.8)

which can further be simplified as

$\displaystyle\tau_{\rm i}{\bar{T}}$ $\displaystyle=$
$\displaystyle\frac{3\sqrt{10}}{\pi g^{1/2}_{\ast}(\tau_{\rm
i})}\frac{m_{\mathrm{P}}}{T(\tau_{\rm i})}.$ (E.9)

The numerical factor in front of $m_{\mathrm{P}}/T(\tau_{\rm i})$ is
approximately 1 and we hence have

$\displaystyle\tau_{\rm i}{\bar{T}}$ $\displaystyle\simeq$
$\displaystyle\frac{m_{\mathrm{P}}}{T(\tau_{\rm i})},$ (E.10)

which shows that $\tau_{\rm i}{\bar{T}}\gg 1$ as long as $T(\tau_{\rm i})\ll
m_{\rm P}$.

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, no.6, 061102 (2016) [arXiv:1602.03837 [gr-qc]].
  * [2] J. Aasi et al. [LIGO Scientific], “Advanced LIGO,” Class. Quant. Grav. 32, 074001 (2015) [arXiv:1411.4547 [gr-qc]].
  * [3] F. Acernese et al. [VIRGO], “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32, no.2, 024001 (2015) [arXiv:1408.3978 [gr-qc]].
  * [4] T. Akutsu et al. [KAGRA], “Overview of KAGRA : KAGRA science,” Prog. Theor. Exp. Phys. 2021, 05A101 (2021). [arXiv:2008.02921 [gr-qc]].
  * [5] B. P. Abbott et al. [LIGO Scientific and Virgo], “Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run,” Phys. Rev. D 100, no.6, 061101 (2019) [arXiv:1903.02886 [gr-qc]].
  * [6] N. Seto and A. Taruya, “Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode,” Phys. Rev. D 77, 103001 (2008) [arXiv:0801.4185 [astro-ph]].
  * [7] V. Gluscevic and M. Kamionkowski, “Testing Parity-Violating Mechanisms with Cosmic Microwave Background Experiments,” Phys. Rev. D 81, 123529 (2010) [arXiv:1002.1308 [astro-ph.CO]].
  * [8] K. Hayama, T. Kuroda, K. Nakamura and S. Yamada, “Circular Polarizations of Gravitational Waves from Core-Collapse Supernovae: A Clear Indication of Rapid Rotation,” Phys. Rev. Lett. 116, no.15, 151102 (2016) [arXiv:1606.01520 [astro-ph.HE]].
  * [9] V. Domcke, J. Garcia-Bellido, M. Peloso, M. Pieroni, A. Ricciardone, L. Sorbo and G. Tasinato, “Measuring the net circular polarization of the stochastic gravitational wave background with interferometers,” JCAP 05, 028 (2020) [arXiv:1910.08052 [astro-ph.CO]].
  * [10] K. Martinovic, C. Badger, M. Sakellariadou and V. Mandic, “Searching for parity violation with the LIGO-Virgo-KAGRA network,” Phys. Rev. D 104, no.8, L081101 (2021) [arXiv:2103.06718 [gr-qc]].
  * [11] A. Lue, L. M. Wang and M. Kamionkowski, “Cosmological signature of new parity violating interactions,” Phys. Rev. Lett. 83, 1506-1509 (1999) [arXiv:astro-ph/9812088 [astro-ph]].
  * [12] K. Choi, J. c. Hwang and K. W. Hwang, “String theoretic axion coupling and the evolution of cosmic structures,” Phys. Rev. D 61, 084026 (2000) [arXiv:hep-ph/9907244 [hep-ph]].
  * [13] S. Alexander and J. Martin, “Birefringent gravitational waves and the consistency check of inflation,” Phys. Rev. D 71, 063526 (2005) [arXiv:hep-th/0410230 [hep-th]].
  * [14] D. H. Lyth, C. Quimbay and Y. Rodriguez, “Leptogenesis and tensor polarisation from a gravitational Chern-Simons term,” JHEP 03, 016 (2005) [arXiv:hep-th/0501153 [hep-th]].
  * [15] M. Satoh, S. Kanno and J. Soda, “Circular Polarization of Primordial Gravitational Waves in String-inspired Inflationary Cosmology,” Phys. Rev. D 77, 023526 (2008) [arXiv:0706.3585 [astro-ph]].
  * [16] L. Sorbo, “Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton,” JCAP 06, 003 (2011) [arXiv:1101.1525 [astro-ph.CO]].
  * [17] J. L. Cook and L. Sorbo, “Particle production during inflation and gravitational waves detectable by ground-based interferometers,” Phys. Rev. D 85, 023534 (2012) [erratum: Phys. Rev. D 86, 069901 (2012)] [arXiv:1109.0022 [astro-ph.CO]].
  * [18] S. Okano and T. Fujita, JCAP 03, 026 (2021) [arXiv:2005.13833 [astro-ph.CO]].
  * [19] P. Adshead, E. Martinec and M. Wyman, “Gauge fields and inflation: Chiral gravitational waves, fluctuations, and the Lyth bound,” Phys. Rev. D 88, no.2, 021302 (2013) [arXiv:1301.2598 [hep-th]].
  * [20] A. Maleknejad, “Chiral Gravity Waves and Leptogenesis in Inflationary Models with non-Abelian Gauge Fields,” Phys. Rev. D 90, no.2, 023542 (2014) [arXiv:1401.7628 [hep-th]].
  * [21] I. Obata and J. Soda, “Chiral primordial Chiral primordial gravitational waves from dilaton induced delayed chromonatural inflation,” Phys. Rev. D 93, no.12, 123502 (2016) [arXiv:1602.06024 [hep-th]].
  * [22] S. Weinberg, “Damping of tensor modes in cosmology,” Phys. Rev. D 69 (2004), 023503 [arXiv:astro-ph/0306304 [astro-ph]].
  * [23] M. M. Anber and E. Sabancilar, “Chiral Gravitational Waves from Chiral Fermions,” Phys. Rev. D 96, no.2, 023501 (2017) [arXiv:1607.03916 [hep-th]].
  * [24] N. D. Barrie and A. Kobakhidze, “Gravitational Instabilities of the Cosmic Neutrino Background with Non-zero Lepton Number,” Phys. Lett. B 772, 459-463 (2017) [arXiv:1701.00603 [hep-ph]].
  * [25] K. Kamada, J. Kume and Y. Yamada, “Chiral gravitational effect in time-dependent backgrounds,” JHEP 05, 292 (2021) [arXiv:2104.00583 [hep-ph]].
  * [26] S. Anand, J. R. Bhatt and A. K. Pandey, “Chiral plasma instability and primordial gravitational waves,” Eur. Phys. J. C 79, no.2, 119 (2019) [arXiv:1801.00650 [astro-ph.CO]].
  * [27] M. Joyce and M. E. Shaposhnikov, “Primordial magnetic fields, right-handed electrons, and the Abelian anomaly,” Phys. Rev. Lett. 79, 1193-1196 (1997) [arXiv:astro-ph/9703005 [astro-ph]].
  * [28] Y. Akamatsu and N. Yamamoto, “Chiral Plasma Instabilities,” Phys. Rev. Lett. 111, 052002 (2013) [arXiv:1302.2125 [nucl-th]].
  * [29] T. Kahniashvili, A. Brandenburg, G. Gogoberidze, S. Mandal and A. Roper Pol, “Circular polarization of gravitational waves from early-Universe helical turbulence,” Phys. Rev. Res. 3, no.1, 013193 (2021) [arXiv:2011.05556 [astro-ph.CO]].
  * [30] S. Weinberg, “Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity,” (John Wiley & Sons, New York, 1971).
  * [31] J. Ghiglieri and M. Laine, “Gravitational wave background from Standard Model physics: Qualitative features,” JCAP 07 (2015), 022 [arXiv:1504.02569 [hep-ph]].
  * [32] J. Ghiglieri, G. Jackson, M. Laine and Y. Zhu, “Gravitational wave background from Standard Model physics: Complete leading order,” JHEP 07 (2020), 092 [arXiv:2004.11392 [hep-ph]].
  * [33] G. G. Stokes, “On the Composition and Resolution of Streams of Polarized Light from Different Sources,” Trans. Cambridge Phil. Soc., 9, 399 (1852).
  * [34] S. Chandrasekhar, “Radiative Transfer,” Clarendon Press, Oxford (1950).
  * [35] R. A. Breuer, J. Tiomno and C. V. Vishveshwara, “Polarization of gravitational synchrotron radiation,” Lett. Nuovo Cim. 4S2, 857-860 (1972).
  * [36] O. F. Piattella, “Lecture Notes in Cosmology,” [arXiv:1803.00070 [astro-ph.CO]].
  * [37] A. Matsumoto et al. “EMPRESS. VIII. A New Determination of Primordial He Abundance with Extremely Metal-Poor Galaxies: A Suggestion of the Lepton Asymmetry and Implications for the Hubble Tension,” [arXiv:2203.09617 [astro-ph.CO]].
  * [38] O. Seto and Y. Toda, “Hubble tension in lepton asymmetric cosmology with an extra radiation,” [arXiv:2104.04381 [astro-ph.CO]].
  * [39] G. Steigman, “Neutrinos And Big Bang Nucleosynthesis,” Adv. High Energy Phys. 2012 (2012), 268321 [arXiv:1208.0032 [hep-ph]].
  * [40] M. Kawasaki and K. Murai, “Lepton Asymmetric Universe,” [arXiv:2203.09713 [hep-ph]].
  * [41] S. Pastor, T. Pinto and G. G. Raffelt, “Relic density of neutrinos with primordial asymmetries,” Phys. Rev. Lett. 102, 241302 (2009) [arXiv:0808.3137 [astro-ph]].
  * [42] P. D. Serpico and G. G. Raffelt, “Lepton asymmetry and primordial nucleosynthesis in the era of precision cosmology,” Phys. Rev. D 71, 127301 (2005) [arXiv:astro-ph/0506162 [astro-ph]].
  * [43] C. Caprini and D. G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav. 35 (2018) no.16, 163001 [arXiv:1801.04268 [astro-ph.CO]].
  * [44] L. Husdal, “On Effective Degrees of Freedom in the Early Universe,” Galaxies 4 (2016) no.4, 78 [arXiv:1609.04979 [astro-ph.CO]].