Working with circular GW polarization states is more convenient here than using the plus and cross polarizations as a basis, even though the latter is more common in the GW literature.

Note that at second order in gravitational wave perturbations, the perturbed metric is $\delta g_{ij}=a^2(h_{ij}+\frac{1}{2}{h}_{ik}{h}_{jk})$ [40].

Given that helicity and chirality are equivalent for free fermions, the zeroth-order solution for left-handed Weyl fermions is similar, but with ${\mathbf{E}}_{\mathbf{k}}^{+}$ replaced by ${\mathbf{E}}_{\mathbf{k}}^{-}$.

We have already performed the quantum expectation value over the fermionic fields, and the remaining averaging on GW is purely stochastic.

A long-wavelength GW mode (homogeneous tensor mode) $h_{ij}^{\mathrm{long}}(t,\mathbf{x})\approx C_{ij}(t)$ can be recast as a change in coordinates $h_{ij}^{\mathrm{long}}\mapsto 0$ and $x^i\mapsto x^{i^{\prime }}=e^{h_{ij}^{\mathrm{long}}}{x}^j$ [46, 47]. In other words, it can be considered as a large gauge transformation. As a result, after this coordinate transformation, a short-wavelength fermion mode in the presence of long-wavelength GWs becomes a free fermion, i.e., ${\mathrm{\Psi }}_{\text{short}}(t,\mathbf{x}){|}_{h_{ij}^{\mathrm{long}}}={\mathrm{\Psi }}^{\prime }(t,{\mathbf{x}}^{\prime }){|}_0$, which implies $⟨{\rho }_{\psi ,\text{short}}(x)⟩{|}_{h_{ij}^{\mathrm{long}}}=⟨{\rho }_{\psi ,\text{short}}(x^{\prime })⟩{|}_0=0$.

The energy density of GWs is defined as ${\rho }_{\mathrm{gw}}=(1/32\pi G)⟨{\dot{h}}_{ij}{\dot{h}}^{ij}⟩$, and its spectral energy density is ${\mathrm{\Omega }}_{\mathrm{gw}}=(1/M_{\mathrm{Pl}}^2)[d{\rho }_{\mathrm{gw}}/3H^{2}d\ln q]$.

Note that ${\mathrm{\Omega }}_{\mathrm{peak}}$ may be related to $T_{*}$ and the GW background properties in a model-dependent manner. This work focuses on a broadly applicable framework, leaving a detailed exploration to future numerical studies.