Gravitational waves from fermion production during axion inflation

TL;DR

This paper analytically investigates gravitational waves generated during axion inflation due to fermion production, finding that the signal is suppressed and does not surpass standard vacuum fluctuation mechanisms, remaining within observational bounds.

Contribution

It provides the first analytic calculation of the gravitational wave spectrum from fermion production during axion inflation, highlighting suppression of polarization and compatibility with current bounds.

Findings

Polarized gravitational wave component is parametrically suppressed.

Gravitational wave amplitude cannot exceed standard vacuum fluctuation signals.

Produced gravitational waves are consistent with current observational constraints.

Abstract

We present analytic results for the gravitational wave power spectrum induced in models where the inflaton is coupled to a fermionic pseudocurrent. We show that although such a coupling creates helically polarized fermions, the polarized component of the resulting gravitational waves is parametrically suppressed with respect to the non-polarized one. We also show that the amplitude of the gravitational wave signal associated to this production cannot exceed that generated by the standard mechanism of amplification of vacuum fluctuations. We previously found that this model allows for a regime in which the backreaction of the produced fermions allows for slow-roll inflation even for a steep inflaton potential, and still leads to Gaussian primordial scalar perturbations. The present analysis shows that this regime also results in a gravitational wave signal compatible with the current bounds.

Figures

Equations

$$\Delta \mathcal{L}=\frac{\phi }{f}F\tilde{F}+\frac{{\partial }_{\mu }\phi }{f}\bar{X}{\gamma }^{\mu }{\gamma }_{5}X,$$

$$S=\int d^{4}x\sqrt{-g}\left[\frac{M_{\mathrm{Pl}}^2}{2}R+\frac{g^{\mu \nu }}{2}{\partial }_{\mu }\phi {\partial }_{\nu }\phi -V(\phi )+\bar{X}\left(i{\gamma }^{\mu }{D}_{\mu }X-m-\frac{1}{f}{\partial }_{\mu }\phi {\gamma }^{\mu }{\gamma }_5\right)X\right].$$

$$X\to \Psi =e^{i{\gamma }_5\frac{\phi }{f}}X,$$

$$S_{\Psi }=\int d^{4}x\sqrt{-g}\bar{\Psi }\left[i{\gamma }^C{e}_C{}^{\mu }\left({\partial }_{\mu }+\frac{1}{2}{\omega }_{\mu AB}{\Sigma }^{AB}\right)\Psi -m\cos \left(\frac{2\phi }{f}\right)+im\sin \left(\frac{2\phi }{f}\right){\gamma }_5\right]\Psi ,$$

$$d{s}^2=N^{2}d{\tau }^2-h_{ij}(d{x}^i+N^{i}d\tau )(d{x}^j+N^{j}d\tau ),$$

$$S_{\mathrm{B}}=$$

$$K_{ij}\equiv$$

$$S_F=\int d^{4}x{\mathcal{L}}_F,$$

$$\displaystyle\mathcal{L}_{F}=a^{3}\Bigg{\{}$$

$$+i\bar{\Psi }\left({\gamma }^{a}N{e}_a{}^k-{\gamma }^0{N}^k\right)\left[{\partial }_k-K_{ik}{e}_b{}^i{\Sigma }^{0b}+\frac{1}{2}\left(e^c{}_i{\partial }_k{e}_b{}^i{+}^{(3)}{\Gamma }^m{}_{ik}{e}^c{}_m{e}_b{}^i\right){\eta }_{ac}{\Sigma }^{ab}\right]\Psi$$

$$\displaystyle-Nm\bar{\Psi}\left[\cos\left(\frac{2\varphi}{f}\right)+i\sin\left(\frac{2\varphi}{f}\right)\gamma_{5}\right]\Psi\Bigg{\}}\,.$$

$$0=\frac{\delta S}{\delta N}=$$

$$+\frac{i}{2}{e}_a{}^i(\bar{\Psi }{\gamma }^a{\partial }_i\Psi -{\partial }_i\bar{\Psi }{\gamma }^a\Psi )-\frac{1}{4}{e}_a{}^i{e}_{bj}{}^{(3)}{\nabla }_i{e}_c{}^j{ϵ}^{abc}\bar{\Psi }{\gamma }^0{\gamma }^5\Psi$$

$$\displaystyle-m\,\bar{\Psi}\Big{[}\mathrm{cos}\left(\frac{2\varphi}{f}\right)+i\,\mathrm{sin}\left(\frac{2\varphi}{f}\right)\gamma^{5}\Big{]}\Psi\,,$$

$$0=\frac{\delta S}{\delta N_i}-{}^{(3)}{\nabla }_j\frac{\delta S}{\delta ({}^{(3)}{\nabla }_j{N}_i)}=$$

$$\displaystyle+M_{\rm Pl}^{2}{}^{(3)}\nabla^{j}\left(K_{ij}-h_{ij}\,K\right)+\frac{1}{4}{}^{(3)}\nabla^{j}\Big{(}e^{a}{}_{i}\,e^{b}{}_{j}\,\epsilon_{abc}\,\bar{\Psi}\gamma^{c}\gamma^{5}\Psi\Big{)}\;.$$

$$h_{ij}=a^2{\left({\mathrm{e}}^{\gamma }\right)}_{ij}=a^2\left[{\delta }_{ij}+{\gamma }_{ij}+\frac{1}{2}{\gamma }_{im}{\gamma }_{mj}+\dots \right],$$

$$e^a{}_i=a{\delta }^{ak}{e}^{\frac{1}{2}{\gamma }_{ki}}=a{\delta }^{ak}\left[1+\frac{1}{2}{\gamma }_{ki}+\frac{1}{8}{\gamma }_{kj}{\gamma }_{ji}\dots \right].$$

$$N$$

$$N_i$$

$${\mathcal{H}}^2=\frac{1}{3M_{\mathrm{Pl}}^2}\left(\frac{{\phi }_0^{\prime }{}^2}{2}+a^{2}V(\phi )\right),\mathcal{H}\equiv \frac{a^{\prime }}{a},$$

$${\beta }^{(1)}=0,{\alpha }^{(1)}=\frac{{\phi }_0^{\prime }}{2\mathcal{H}{M}_P^2}\delta \phi ,\Delta {\theta }^{(1)}=-\frac{{\phi }_0^{\prime 2}}{2M_P^2{\mathcal{H}}^2}{\left(\frac{\mathcal{H}\delta \phi }{{\phi }_0^{\prime }}\right)}^{\prime }.$$

$${\alpha }^{(2)}$$

$${\beta }_j^{(2)}$$

$$+\frac{ia}{M_{\mathrm{Pl}}^2}{\partial }_j{\Delta }^{-1}\left[\bar{\Psi }{\gamma }^0\Delta \Psi -(\Delta \bar{\Psi }){\gamma }^0\Psi \right]$$

$$-\frac{a}{M_{\mathrm{Pl}}^2}\left[i\left(\bar{\Psi }{\gamma }^0{\partial }_j\Psi -({\partial }_j\bar{\Psi }){\gamma }^0\Psi \right)-\frac{1}{2}{ϵ}_{ijk}{\partial }_i(\bar{\Psi }{\gamma }^k{\gamma }^5\Psi )\right]\},$$

$${\theta }^{(2)}$$

$$-\frac{a^2}{M_{\mathrm{Pl}}^2\mathcal{H}}V({\phi }_0){\Delta }^{-1}\left\{\frac{1}{8\mathcal{H}}{\partial }_j\left[({\partial }_j{\gamma }_{\ell i}){\gamma }_{i\ell }^{\prime }\right]+\frac{ia}{4M_{\mathrm{Pl}}^2\mathcal{H}}\left[\bar{\Psi }{\gamma }^0\Delta \Psi -(\Delta \bar{\Psi }){\gamma }^0\Psi \right]\right\}\},$$

$$S_{\gamma }^{(2)}=\frac{M_{\mathrm{Pl}}^2}{8}\int d^{4}x{a}^2\left[{\gamma }_{ij}^{\prime }{\gamma }_{ij}^{\prime }-{\partial }_k{\gamma }_{ij}{\partial }_k{\gamma }_{ij}\right],$$

$$S_F^{(2)}=\int d^{4}x\left[i\bar{\psi }\left({\gamma }^0{\partial }_0+{\gamma }^a{\partial }_a\right)\psi -ma\cos \left(\frac{2\phi }{f}\right)\bar{\psi }\psi +ima\sin \left(\frac{2\phi }{f}\right)\bar{\psi }{\gamma }_5\psi \right],$$

$$S_F^{(3)}=-\frac{i}{2}\int d^{4}x{\gamma }_{ij}\bar{\psi }{\gamma }^i{\partial }_j\psi \equiv \int d^{4}x{\mathcal{L}}^{(3)},$$

$$S_F^{(4)}=$$

$$+\frac{i}{4}\left(1-\frac{V}{4H^2{M}_{\mathrm{Pl}}^2}\right){\Delta }^{-2}{\partial }_m\left({\partial }_m{\gamma }_{kn}{\gamma }_{kn}^{\prime }\right)\left(\bar{\psi }{\gamma }^0\Delta \psi -(\Delta \bar{\psi }){\gamma }^0\psi \right)$$

$$-\frac{{\Delta }^{-1}}{8}({\gamma }_{jk}^{\prime }{\partial }_j{\gamma }_{ik}-{\gamma }_{jk}{\partial }_j{\gamma }_{ik}^{\prime }-{\gamma }_{kj}^{\prime }{\partial }_i{\gamma }_{kj})\left[{ϵ}_{aic}{\partial }_a(\bar{\psi }{\gamma }^c{\gamma }^5\psi )+2i\left(\bar{\psi }{\gamma }^0{\partial }_i\psi -{\partial }_i\bar{\psi }{\gamma }^0\psi \right)\right]$$

$$-\frac{i}{32aH}\left({\partial }_i{\gamma }_{jk}{\partial }_i{\gamma }_{jk}+{\gamma }_{ij}^{\prime }{\gamma }_{ij}^{\prime }\right){\Delta }^{-1}\left(\bar{\psi }{\gamma }^0\Delta \psi -\Delta \bar{\psi }{\gamma }^0\psi \right)$$

$$\displaystyle-\frac{i}{16\,aH}\left(\bar{\psi}\gamma^{0}\psi^{\prime}-\bar{\psi}^{\prime}\,\gamma^{0}\psi\right)\,\Delta^{-1}\partial_{i}\left(\partial_{i}\gamma_{kj}\,\gamma^{\prime}_{kj}\right)\Bigg{\}}\equiv\int d^{4}x\,{\cal L}^{(4)}\equiv\sum_{i=1}^{7}\,\int d^{4}x\,{\cal L}_{i}^{(4)}\,,$$

$$H_{\mathrm{int}}^{(3)}\left(\tau \right)$$

$$H_{\mathrm{int},i}^{(4)}\left(\tau \right)$$

Full text

1 Introduction

There is strong observational evidence supporting the hypothesis that the early universe went through a rapid period of accelerated expansion dubbed inflation [1, 2]. As well as resolving the horizon, flatness, curvature, and gravitino/monopole problems, inflation provides a simple explanation for the observed red-tilted, approximately Gaussian and adiabatic density fluctuations [3, 4]. The inflationary scenario also generically predicts primordial gravitational waves, which can be measured or constrained through the $B$-mode polarization of the cosmic microwave background. Current measurements restrict the tensor-to-scalar ratio to $r\lesssim 0.09$ [5, 6], which tightens to $r\lesssim 0.06$ when the consistency relation $n_{t}=-r/8$, as appropriate for vacuum fluctuations from slow-roll inflation, is imposed. Future experiments are expected to reach $\sigma_{r}\sim 0.001$ [7].

Slow-roll inflation requires a flat potential; in axion (or natural) inflation this flatness is a natural consequence of an approximate shift symmetry [8]. Axion, or pseudoscalar, models can be motivated within string theory, in which case monodromy [9, 10, 11, 12, 13] or alignment [14, 15, 16, 17] effects can give rise to vacuum primordial gravitational waves within the reach of future experiments.

Due to the approximate shift symmetry, the axion inflaton $\varphi$ must couple derivatively to matter fields. At mass-dimension five, the possible couplings are

[TABLE]

to gauge fields and fermions $X$, respectively. Here $F$ is the usual gauge-field field-strength tensor, $\tilde{F}$ is its dual, and $f$ is a scale known as the axion decay constant. This coupling of the axion to gauge fields leads to exponentially large gauge field amplification, with several possible phenomenological consequences (see [18] for a review). These include steep inflation [19], thermalized inflation [20, 21], magnetic field production [22, 23, 24, 25, 26, 27], large non-Gaussianity [28, 29, 30], chiral gravitational wave production [31, 32, 33, 34, 35, 36], instantaneous preheating [37, 38], and the generation of primordial black holes [39, 40, 41].

The fermionic coupling has not attracted as much attention, first being studied by one of us in reference [42], and, more recently, in references [43, 44, 45]. Due to Pauli blocking, fermions cannot undergo the same exponential amplification as the gauge fields. Furthermore, on the one hand, massless fermions are conformal and therefore cannot be created gravitationally through the expansion of the Universe [46]. On the other hand, very heavy fermions decouple, so that, in the absence of the coupling (1.1), only fermions with mass $m\approx H$ are produced in sizable quantities. With only one scale in the problem, the Hubble scale, the energy density is $\sim H^{4}$, which is too small to produce observable effects (with the possible exception of super-heavy dark matter [47, 48]). During axion inflation, however, this conclusion does not follow due to the presence of the additional scale $\dot{\varphi}/f$.

We have previously studied the regime $\dot{\varphi}/f\gg H$ in reference [43]. Because fermion modes can be populated up to $k\sim\dot{\varphi}/f$, the energy density can be parametrically larger than $H^{4}$, as first noticed in reference [42]. In reference [43] we also identified a regime in which the sourced contribution to the power spectrum dominates the vacuum contribution, yet the non-Gaussianity is beneath current observational bounds. This behavior is in striking contrast to the analogous effect in systems with strong bosonic particle production. This difference can be understood by noting that phenomenologically interesting results require one to populate a large number of fermion modes and, since the occupation of these modes is restricted by Pauli exclusion, their sum is uncorrelated and becomes increasingly Gaussian by the central limit theorem. Bosonic systems, conversely, allow for large occupation numbers per mode, which add coherently leading to sourced $n$-point functions that are generically related to the two-point function by $\langle\delta\varphi^{n}\rangle\sim\left(\langle\delta\varphi^{2}\rangle\right)^{n/2}$ [49]. As the fermion-axion coupling is increased, eventually one enters a regime of strong backreaction where the evolution of the inflaton zero-mode is controlled by particle production; this is the fermionic analogue of steep inflation studied in reference [19]. We expect that the above argument that non-Gaussianity remains small still holds in the regime of strong backreaction.

Given the rich phenomenology of the scalar perturbations found in reference [43], it is important to characterize how fermions source gravitational waves in the regime where $\dot{\varphi}/f\gg H$. As was noticed in reference [42] and confirmed in reference [43], the produced fermions have a helicity asymmetry, which was used for leptogenesis in reference [50]. This helicity asymmetry raises the possibility that the spectrum of sourced gravitational waves has a chiral component. Sourced production of gravitational waves in this context was previously studied111Gravitational wave production by non-chiral fermions has also been studied in reference [51]. in reference [52]. However, the fermion basis used in that work was leading to pathologies as $m\to 0$, whereas in this work we use the basis introduced in reference [43], in which perturbation theory remains valid as the fermion mass $m$ becomes small. We review this basis while introducing our model in section 2.

This paper is organized as follows. In section 2, we introduce our theory, and working in the Arnowitt-Deser-Misner (ADM) formalism, we solve the gravitational constraint equations to second order. We then use these solutions to obtain the interaction Lagrangian to fourth order in fluctuations. From this interaction Lagrangian, we obtain one ${\cal O}(\gamma\Psi^{2})$ vertex and seven ${\cal O}(\gamma^{2}\Psi^{2})$ vertices. In section 3, we use these interactions to compute eight loop diagrams in the in-in formalism. The ${\cal O}(\gamma^{2}\Psi^{2})$ interactions lead to seven one-loop one-vertex diagrams, which we evaluate in section 3.2.1, while the ${\cal O}(\gamma\Psi^{2})$ generates a two-vertex loop, which we evaluate in section 3.2.2. We discuss our results in section 3.3; we show that the chirally asymmetric contribution is subdominant, and that the total sourced contribution to the tensor-to-scalar ratio is beneath the vacuum component. Details of our calculations can be found in the various appendices. We work in natural units where $\hbar=c=1$, and $M_{\rm Pl}=1/\sqrt{8\pi G}$ is the reduced Planck mass.

2 Fermion-graviton interactions during axion inflation

The aim of this work is to compute the production of gravitational waves by the fermions. At leading order in perturbation theory, fermions source gravitational wave power at one-loop. At one-loop, diagrams of two topologies are possible. The first topology—the cubic loop—is a two-vertex diagram that is generated by two cubic order vertices consisting of one gravitational wave and a fermion bilinear. The second topology—the quartic loop—is a one-vertex diagram generated by a quartic-order vertex consisting of two gravitational waves and a fermion bilinear (see figure 1 below). To find the required interactions, we therefore need to expand the full action to quartic order in fluctuations.

2.1 Starting action

We consider a theory containing a pseudoscalar inflaton $\varphi$ with a shift-symmetric coupling to a fermion $X$ and minimally coupled to gravity, so that our action, in mostly minus convention, reads

[TABLE]

As discussed in reference [43] (in particular, see appendix B of that work for details), the use of these fields makes apparent the shift-symmetric nature of the inflaton-fermion interaction, but it obscures the fact that such interaction vanishes as $m\to 0$. It is therefore convenient to redefine the fermion according to

[TABLE]

which puts the fermion action in the form

[TABLE]

where greek letters are spacetime indices $\mu,\nu\in\{0,1,2,3\}$, capital roman letters are 4D Lorentz indices $A,B,C\in\{0,1,2,3\}$, lower case roman letters from the start of the alphabet are spatial Lorentz indices $a,b,c\in\{1,2,3\}$, and finally roman letters from the middle of the alphabet are spatial spacetime indices, $i,j,k\in\{1,2,3\}$. The generator of local Lorentz transformations is $\Sigma^{AB}=\frac{1}{4}\left[\gamma^{A},\gamma^{B}\right]$, and the spin connection is $\omega_{\mu}{}^{AB}=e^{A}{}_{\nu}\nabla_{\mu}e^{B\nu}$, where $e^{A}{}_{\nu}$ is the vierbein.

2.2 The action in ADM form

Because certain components of the metric are constrained degrees of freedom whose values depend on the fermion bilinears (as well as the other dynamical degrees of freedom), gravitationally-mediated fermion-graviton couplings are generated when these constraints are eliminated from the action. In order to perform this analysis, it is convenient to decompose the metric using the ADM formalism. The key advantage of this formulation is that the constrained degrees of freedom enter the theory algebraically; their equations of motion are algebraic constraints. The metric in ADM form reads

[TABLE]

where $N$ is the lapse and $N^{i}$ is the shift. For the background metric we choose $N=a$, so that $\tau$ denotes conformal time. Derivatives with respect to $\tau$ are represented with primes. The spatial indices, $i,\,j,\,k,...$ are raised and lowered with $h_{ij}$, so that $N^{i}\equiv h^{ij}\,N_{j}$, $h^{ij}\,h_{jk}=\delta^{i}_{k}$. Finally, ${\rm det}\left[g\right]=-N^{2}\,{\rm det}\left[h\right]$.

In these coordinates, the action for the purely bosonic sector of the theory (involving gravity and the inflation) becomes

[TABLE]

where $\pi_{\varphi}\equiv{\varphi}^{\prime}-N^{j}\partial_{j}\varphi$ and

[TABLE]

The fermionic action, $S_{F}$, in these coordinates reads

[TABLE]

where (see appendix A)

[TABLE]

In these expressions, ${}^{(3)}\nabla_{i}$ denotes the three dimensional covariant derivative, $\eta_{ab}$ is the spatial part of the Minkowski metric, $\eta_{AB}={\rm diag}[1,-1,-1,-1]$, and the spatial vielbeins satisfy $\delta_{ab}e^{a}\phantom{}_{i}e^{b}\phantom{}_{j}=h_{ij}$. The total action is the sum $S=S_{B}+S_{F}$.

2.3 Constraints

When written in terms of the ADM decomposition, one can see that the lapse $N$ and the shift $N_{i}$ enter in the action, eqs. (2.5) and (2.8), without time derivatives (in the case of the lapse, spatial derivatives are also missing from the action). This implies that the corresponding Euler-Lagrange equations are constraints. The equation of motion for the lapse is the Hamiltonian constraint

[TABLE]

while the equation of motion for the shift is the momentum constraint

[TABLE]

In these expressions, ${\epsilon}^{abc}$ is the “flat” three-dimensional Levi-Civita tensor, with the convention ${\epsilon}^{123}=+1$.

We work in the spatially flat gauge where $\det[h_{ij}]=a^{6}$ and the dynamical scalar fluctuation degrees of freedom are in the fluctuations of the inflaton. We parametrize the tensor perturbations of the metric as222Repeated lower roman indices are summed with the Kronecker delta: $x_{i}x_{i}\equiv\sum_{i,j=1}^{3}\delta^{ij}x_{i}x_{j}$

[TABLE]

where $a(\tau)$ is the scale factor. The transverse-traceless nature of the tensor modes, $\delta^{ij}\gamma_{ij}=\gamma_{ij,j}=0$, implies that $\det\left[e^{\gamma}\right]=1$. Similarly, the spatial vielbeins are expanded in terms of the tensor perturbations as

[TABLE]

Our goal is to determine the effective cubic $\gamma\bar{\Psi}\Psi$ and quartic $\gamma\gamma\bar{\Psi}\Psi$ component of the Lagrangian, where $\gamma$ schematically denotes the graviton. As is well known, in order to determine the action to $n$-th order in the fluctuations, the solutions to the constraint equations are required at order $n-2$; terms of order $n$ and $n-1$ simply multiply lower-order constraint equations [53]. Thus, to obtain the action up to fourth order in the fluctuations, we require solutions for the constraints (the lapse and shift) up to quadratic order. Note that Lorentz invariance means that the fermion fields only begin to contribute to the constraint equations at quadratic order. We therefore solve the above constraints in eqs. (2.3) and (2.3) perturbatively, to second order in the fluctuations, before plugging them back into the original action.

To facilitate a perturbative solution, we expand the lapse and shift functions as

[TABLE]

where the superscript denotes the order of the expansion, and where $\beta_{i}^{(1,\,2)}$ are transverse, $\partial_{i}^{(1,\,2)}\beta_{i}=0$. We also expand the inflaton as $\varphi(\tau,\,\vec{x})=\varphi_{0}(\tau)+\delta\varphi(\tau,\,\vec{x})$, and we treat the fermions as first order quanties, so that fermion bilinears are of second order.

To zeroth order, the Hamiltonian constraint reduces to the Friedmann equation

[TABLE]

while the momentum constraint is automatically satisfied. At first order, we obtain [54]

[TABLE]

Since we are not interested in the perturbations sourced by fluctuations of the inflaton, we drop these from now on. Ignoring inflaton fluctuations, the second order constraints read

[TABLE]

where $\Delta=\partial_{i}\partial_{i}$ is the spatial Laplacian, $\Delta^{-1}$ is its inverse, and we note that $\epsilon_{123}=-1$. We have used the linear equation of motion for the fermion to simplify the solution for $\theta^{(2)}$; the details of this calculation are given in appendix A.

2.4 Explicit form of the fermion action, and fermion-GW interactions

We insert the solutions to the constraint equations for $N$ and $N^{i}$ to second order (eqs. (2.14), (2.15), and (2.16)) into the action, eq. (2.5) + (2.7), and then expand order by order in the fluctuations. This gives the quadratic action $S^{(2)}=S^{(2)}_{\gamma}+S_{F}^{(2)}$ for the free gravitons and fermions, the cubic action $S_{F}^{(3)}$ describing the ${\cal O}\left(\Psi^{2}\gamma\right)$ interactions, and the quartic action $S_{F}^{(4)}$ describing the ${\cal O}\left(\Psi^{2}\gamma^{2}\right)$ interactions.

The quadratic action for the gravitons reads

[TABLE]

while the quadratic action for the fermions is

[TABLE]

where we have rescaled the fermion field according to ${\psi}\equiv a^{3/2}\Psi$.

At cubic order we find

[TABLE]

and some straightforward, but lengthy, algebra leads to the quartic order action

[TABLE]

where $H=\dot{a}/a$ is the Hubble parameter and $\epsilon_{123}=-1$. ${\cal L}_{1}^{(4)},\,{\cal L}_{2}^{(4)},\,{\cal L}_{3}^{(4)}$ refer to the three terms in the first line of eq. (2.20), while the remaining ${\cal L}_{i}^{(4)}$ refer to the other four lines (one term per line). Note that the interactions ${\cal L}_{4}^{(4)},\,{\cal L}_{5}^{(4)},\,{\cal L}_{6}^{(4)},\,{\cal L}_{7}^{(4)}$ arise from integrating out the non-dynamical constraints (the second order parts of the lapse and shift). From the cubic and quartic Lagrangian densities we find the interaction Hamiltonian densities

[TABLE]

To proceed, we expand the tensors in Fourier space as

[TABLE]

where the field ${\rm t}^{\lambda}_{\bf k}$ is canonically normalized. The sum is over the right-handed $(\lambda=+1)$ and left-handed $(\lambda=-1)$ tensor polarizations, with the polarization tensors satisfying

[TABLE]

We also Fourier transform the fermions according to

[TABLE]

In terms of the fields ${\rm t}^{\lambda}_{\bf k}$ and $\psi_{\bf k}$, the quadratic action is

[TABLE]

and we note that the kinetic terms are canonically normalized.

Inserting eqs. (2.22) and (2.24) in the interaction Hamiltonians eq. (2.21), we obtain the Fourier space Hamiltonian densities we report in appendix B.

3 Fermion contributions to the tensor power spectrum

In this section, making use of the interaction Hamiltonians derived in the previous section, we compute the fermion contribution to the gravitational wave two-point correlation function. After quantizing the free theory, we introduce the in-in formalism and compute the cubic and quartic loops generated by the interactions derived in section 2. Finally, we end this section by showing how simple scaling arguments concur with our results.

3.1 Quantization

We canonically quantize the theory by expanding the fields into modes

[TABLE]

where the creation-annihilation operators for the tensor modes satisfy the commutation relations

[TABLE]

and the fermionic operators satisfy anti-commutation relations

[TABLE]

The mode functions ${\rm t}_{k}^{\lambda}\left(\tau\right)$, and the spinors $U^{r}_{\bf k}\left(\tau\right)$ and $V^{r}_{-{\bf k}}\left(\tau\right)$ are solutions of the Euler-Lagrange equations of motion that follow from the action in eq. (2.25). We further decompose the 4-component fermionic spinors into helicity states

[TABLE]

where $C=i\gamma^{0}\gamma^{2}$ is the charge-conjugation operator, and the spinors $\chi_{r}\left({{\bf k}}\right)$ are explicitly given by

[TABLE]

where $k_{z}$ is the $z-$component of ${\bf k}$, and $\sigma_{i}$ are the Pauli matrices. Note that $\chi_{r}\left({{\bf k}}\right)$ are helicity eigenspinors which satisfy ${\bf k}\cdot\bm{\sigma}\chi_{r}\left({{\bf k}}\right)=rk\chi_{r}\left({{\bf k}}\right)$. We use the Dirac representation for the $\gamma$ matrices,333In these expressions, $\mathds{1}$ denotes the $2\times 2$ identity matrix.

[TABLE]

To obtain solutions to the classical mode equations, we approximate the background inflationary spacetime as de Sitter space and take the evolution of the inflaton to be rolling at a constant speed in cosmic time. This implies $\varphi_{0}(\tau)/f=\varphi_{0}^{\rm{in}}/f-2\xi\log\left(x/x_{\rm in}\right)$, with $x\equiv-k\tau$ and $x_{\rm in}\equiv-k\tau_{\rm in}$, where $\tau_{\rm in}$ is some reference time. With these approximations, the mode functions for the fermion field are given by

[TABLE]

which satisfy the normalization condition $|u^{r}|^{2}+|v^{r}|^{2}=2$, with [42, 43]

[TABLE]

where $W_{\mu,\,\lambda}(z)$ denotes the Whittaker W-function and

[TABLE]

In the same approximation, the tensor mode functions read

[TABLE]

In both cases, the integration constants have been chosen so that the solutions match onto the appropriate Bunch-Davies vacuum solution at early times, $x\to\infty$.

3.2 Fermion loop-corrections to the gravitational wave power spectrum

The interaction Hamiltonians derived above allow us to compute the leading order contributions from the produced fermions to the two-point function of the graviton. These are computed using the in-in formalism, where the correlation function of an operator ${\hat{O}}_{1}\dots\hat{O}_{n}\left(\tau\right)$ at time $\tau$ is given by

[TABLE]

The interactions in section 2 result in two classes of diagrams: there are seven quartic loop diagrams, illustrated in the left panel of figure 1, one for each of the seven vertices generated by the quartic action (2.20), and one cubic loop diagram with two vertices generated by the cubic action (2.19), illustrated in the right panel of figure 1. We discuss these diagrams in the next two subsections.

3.2.1 Quartic loops

We begin with the left diagram of figure 1. The seven terms in the quartic action lead to seven quartic contributions to the graviton spectrum of the form

[TABLE]

where the interaction Hamiltonians $H^{(4)}_{{\rm int},i}$ are given in eq. (B). Remarkably, all these diagrams can be computed exactly. The details of the calculation, as well as the exact results, are presented in appendix C. Here we summarize the main issues one encounters when performing this calculation. At the end of this section we present the expression of the sum of the quartic loops in the limit $\mu\ll 1\ll\xi$.

First, several of the terms in the interaction Hamiltonian contain the nonlocal operator $\Delta^{-1}$, the inverse of the Laplacian. When evaluating eq. (3.16) one often encounters the expectation value of quantities evaluated at vanishing momentum that, when acted upon by $\Delta^{-1}$, lead to a undetermined “$0/0$” that needs to be regularized. To deal with this limit we follow the prescription given in [55]: these undetermined quantities are schematically given by

[TABLE]

where ${\bf p}_{1}$ is an external momentum. We regularize eq. (3.17) by setting ${\bf q}_{1}={\bf q}_{2}+\bm{\epsilon}$, where we eventually send $\bm{\epsilon}\to 0$. Since $f$ is a scalar, it depends only on $\bm{\epsilon}\cdot{\bf p}_{1}$, $p_{1}^{2}$, and ${\epsilon}^{2}$. We then impose that $\bm{\epsilon}$ approaches zero in a direction that is orthogonal to ${\bf p}_{1}$, so that $\bm{\epsilon}\cdot{\bf p}_{1}={\cal O}({\epsilon}^{2})$. With this convention, all the operators containing $\Delta^{-1}$ give finite and unambiguous results.

Secondly, many integrals contributing to the graviton two-point function are divergent in the ultraviolet. We deal with these divergences as we did in [43], by introducing a ultraviolet cutoff $\Lambda$ and by subtracting all the terms that are divergent as $\Lambda\to\infty$. As we have discussed in [43], we expect the result of this procedure to be equivalent to that obtained by adiabatic subtraction in the limit $\xi\gg 1$.

After long calculations, which we outline in appendix C, we obtain the leading contribution from the quartic diagrams

[TABLE]

in the limit $\mu\ll 1\ll\xi$ and for superhorizon modes $-k\tau\ll 1$. We note that this contribution is parity-even. Parity-odd terms are associated to the operators $H_{\mathrm{int},2}^{(4)}$, $H_{\mathrm{int},3}^{(4)}$ and $H_{\mathrm{int},5}^{(4)}$, which contain the Levi-Civita symbol. However, the contributions from $H_{\mathrm{int},2}^{(4)}$ and from the parity-odd part of $H_{\mathrm{int},5}^{(4)}$ vanish identically after angular integrations, so that the only parity-odd contribution to the tensor power spectrum is given by $H_{\mathrm{int},3}^{(4)}$ and yields

[TABLE]

which is sub-leading, by a factor $1/\xi$, with respect to the parity-even component.

3.2.2 Cubic loop

Next we consider the cubic loop, shown on the right side of figure 1. There is a single contribution to this diagram, given by

[TABLE]

where $H^{(3)}_{\rm int}(\tau)$ is given by eq. (B.1). Unlike those appearing in the quartic loops, the integrals in the cubic loop are prohibitively difficult to evaluate exactly. The expressions appearing in this diagram, however, are very similar to those which appeared in the cubic loop contribution to the spectrum of scalar perturbations considered in reference [43]. Therefore, we apply the same sequence of approximations developed in that work to the present calculation. Here we outline these approximations; the details of the calculation are presented in appendix D.

We start by setting the external momenta to zero; as discussed in [43] we expect this approximation to generate at most a ${\cal O}(1)$ error. Next, since the functions appearing in the integrals are rapidly oscillating, we perform a Wick rotation on the time integration variables, so that the Whittaker functions appearing in the fermion mode functions now have real argument and are exponentially increasing or decreasing. Next, we approximate those Whittaker functions as linear combinations of monomials times exponentials, with special attention given to the branch cuts. The explicit form of those approximations are given in eqs. (D) and (D.6), and we have verified their validity in the regime $\xi\gg 1$ we are interested in. These approximate expressions contain a part that behaves like positive frequency (we schematically denote the coefficient of this part by ${\cal A}$) and a part that behaves like negative frequency (whose coefficient is denoted schematically by ${\cal B}$). Explicit expressions for $\cal{A}$ and $\cal{B}$ can be found in appendix D.

Once the above approximations are in place, the integrals can be computed analytically. We find a divergence in the limit $\tau_{2}\to\tau_{1}$, although it is only present in the ${\cal A}^{2}$ term. Since this term corresponds to positive frequency, “vacuum only” modes, we subtract them. Once this component is subtracted, we are left with the final result

[TABLE]

which has the same parametric dependence as the contribution from the quartic loop.

3.3 Scaling of our result

The gravitational wave power spectrum is related to the two point function as

[TABLE]

The contribution from the produced fermions, using the results eqs. (3.18) and (3.21) we derived in the previous subsections, is

[TABLE]

We now compare this sourced gravitational wave signal to the vacuum contribution, ${\cal P}_{t}^{\rm{vacuum}}\sim 0.1\,H^{2}/M_{\rm Pl}^{2}$. Their ratio can be written as

[TABLE]

As discussed in reference [43], the quantity $\mu^{2}\,\xi^{3}\,H^{2}/M_{\rm Pl}^{2}$ corresponds to the ratio between the energy density in fermions and the total energy density in the Universe, which must be much smaller than unity. We thus conclude that an axion-like inflaton coupled to fermions cannot produce tensor modes that dominate over the spectrum of vacuum fluctuations, in contrast to the scenario in which an axionic inflaton is coupled to gauge fields.

Finally, in reference [43] we determined the overall scaling, as a function of the parameter $\xi$ and $\mu$, of the diagrams that are relevant for the scalar (bi)spectrum. Here we apply analogous arguments to the diagrams that led to the result in eq. (3.23) for the tensor spectrum.

Our first observation is that the cubic interaction (2.19) gives a contribution $\sim\gamma_{ij}\,\bar{\psi}\,\Gamma\,p\,\psi$, whereas each quartic interaction in eq. (2.20) gives a contribution $\sim\gamma_{ij}^{2}\,\bar{\psi}\,\Gamma\,p\,\psi$, where $p$ schematically denotes a quantity that scales as fermion momentum and $\Gamma$ denotes some combination of the Dirac $\gamma$-matrices (exceptions are ${\cal L}_{2}^{(4)}$ and ${\cal L}_{3}^{(4)}$ which contain no dependence on the fermion momentum, but only on the graviton momentum). The contribution from $\langle\bar{\psi}\,\gamma^{i}\,\hat{p}_{i}\,\gamma^{5}\,\psi\rangle$ is ultimately zero, due to the asymmetry of the Levi-Civita symbol. A numerical evaluation shows that the combinations $\langle\bar{\psi}\,\gamma^{i}\,\hat{p}_{i}\,\psi\rangle$ and $\langle\bar{\psi}\,\gamma^{5}\,\gamma^{0}\,\psi\rangle$ (which appear in both the cubic and quartic gravitational vertices) oscillate with amplitude $\mu^{2}/\xi$, for momenta up to $-k\tau\simeq\xi$. Moreover, the fermion part of the operators in ${\cal L}^{(4)}_{4,\,5,\,6,\,7}$ can always be brought to a form $\langle{\Im}\left(\partial_{k}\bar{\psi}\,\gamma^{0}\psi\right)\rangle$ (see eqs. (C) and (C) – here ${\Im}$ denotes the imaginary part), which can also be seen to scale as $\mu^{2}/\xi$ for momenta up to $-k\tau\simeq\xi$. (We note for comparison that the bilinears which appear in the diagrams involving fluctuations of the inflaton, $\langle\bar{\psi}\,\psi\rangle$ and $\langle\bar{\psi}\,\gamma^{5}\psi\rangle$, oscillate instead with amplitude $\mu/\xi$. These were considered in reference [43].) Since all fermion bilinears in the $\gamma\bar{\psi}\psi$ and the $\gamma\gamma\bar{\psi}\psi$ sector scale as $\mu^{2}/\xi$, each fermionic line in the diagrams of figure 1 contributes $\mu^{2}/\xi$.

We recall that each interaction Hamiltonian (with the exception of ${\cal L}_{2}^{(4)}$ and ${\cal L}_{3}^{(4)}$, as noted above) also carries a power of $p$, and therefore each vertex gives an additional power of $\xi$. Furthermore, every fermionic loop integral, which goes as $d^{3}k$, gives a contribution $\sim\xi^{3}$.

Once we apply these scalings to the quartic diagrams, we have a scaling $(\mu^{2}/\xi)$ (one fermion line) times $\xi$ (one vertex) times $\xi^{3}$, giving an overall scaling $(\mu^{2}/\xi)\times(\xi)\times(\xi^{3})\sim\mu^{2}\,\xi^{3}$. The quartic diagrams involving ${\cal L}_{2}^{(4)}$ and ${\cal L}_{3}^{(4)}$ are respectively vanishing (as a consequence of the symmetries of the operator) and scaling as $(\mu^{2}/\xi)\times(1)\times(\xi^{3})\sim\mu^{2}\,\xi^{2}$—the factor of $(1)$ instead of $(\xi)$ originates from the fact that this vertex does not contain a power of the fermion momentum. These results are in agreement with the direct calculations presented in appendix C.

For the cubic gravitational diagram, a naive implementation of these scalings would read $(\mu^{2}/\xi)^{2}$ (two fermion lines) times $\xi^{2}$ (two vertices) times $\xi^{3}$, giving an overall scaling $\sim\mu^{4}\,\xi^{3}$. This would disagree by a factor of $\mu^{2}$ with the result obtained by the direct (albeit approximate) calculation of appendix D. The different scaling with $\mu$ can, however, be understood as follows. Each fermion line in fact scales as $\mathcal{O}(1)+\mathcal{O}(\mu^{2})$, although in the quartic diagrams the $\mathcal{O}(1)$ contribution is always divergent and therefore is removed by regularization. The cubic diagram, with two fermion lines, has terms of order $\mathcal{O}(1)$, $\mathcal{O}(\mu^{2})$, and $\mathcal{O}(\mu^{4})$. The first is again removed by renormalization, leaving the $\mathcal{O}(\mu^{2})$ term which arises from interference.

Given the parity violating nature of the system, one can expect a parity violating tensor spectrum $\delta{\cal P}_{t}^{+1}\neq\delta{\cal P}_{t}^{-1}$, which we did indeed find. As mentioned above, however, this is subdominant by a factor $1/\xi$ with respect to the parity-even part.

4 Conclusion

Axion, or natural inflation is a class of models for slow-roll inflation where the required flatness of the potential is protected from radiative corrections by an approximate shift symmetry. This shift symmetry means that any axion-matter couplings must be via derivatives, and the lowest dimension couplings of an axion inflaton to gauge fields and fermions are given by eq. (1.1). These couplings are typically employed for reheating in these models. The recent literature (see [18] for a review) has shown that the coupling of the axion to gauge fields can lead to a rich phenomenology during inflation. Analogous studies for fermions are more scarce [42, 52, 43, 44]. The present work is a direct continuation of our previous paper [43], where it was pointed out that it is more convenient to perform computations in the redefined basis ${\rm e}^{i\gamma_{5}\frac{\varphi}{f}}\,X$ (where $\varphi$ and $X$ are, respectively, the inflaton and the fermionic fields appearing in eq. (1.1)), as perturbation theory in the $X$ basis breaks down as the fermion mass $m$ approaches zero. Our previous work considered the non-Gaussianity of inflaton perturbations sourced by the coupled fields. In the fermionic case, Pauli blocking limits the occupation number in each individual mode, and therefore phenomenologically interesting results occur only when the Fermi sphere of produced modes becomes large. The many modes produced are uncorrelated sources of inflaton perturbations, and therefore summing over the modes produces Gaussian inflaton perturbations by the central limit theorem. This is in contrast to the scenario in which the axion inflaton is coupled to gauge fields. In that case, the inflaton perturbations of any given scale $k$ are mostly produced by a few (highly amplified) gauge field modes of comparable scale, thus resulting in a strongly non-Gaussian sourced signal [28].

In this paper we complemented our previous study with the computation of the sourced gravitational waves. We worked in the ADM basis, in which the nondynamical metric perturbations $N$ and $N^{i}$ are integrated out via the corresponding energy and momentum constraints. We have solved these constraints perturbatively in the fermionic field $\psi$ and in the gravitational waves $\gamma$. We thus obtained the ${\cal O}\left(\gamma\bar{\Psi}\Psi\right)$ and ${\cal O}\left(\gamma\gamma\bar{\Psi}\Psi\right)$ interactions, which we used to compute the contribution to the gravitational wave spectrum from the diagrams shown in figure 1. The left diagram is technically simpler; after regularizing it we were able to evaluate it exactly, as described in appendix C. The right diagram is much more involved, necessitating the approximations discussed in appendix D, analogous to those made for the cubic diagram in reference [43]. Although the computations were very involved, in section 3.3 we presented some simple scaling arguments that correctly capture the scaling of the result with the parameters, $\mu$ and $\xi$, of the model.

Our main conclusion is that, in contrast to the scenario in which the axion inflaton is coupled to vector fields, the gravitational waves sourced by the fermions cannot be greater than the vacuum gravitational waves. This conclusion holds also in the regime, studied in reference [43], where the strong backreaction of the fermion degrees of freedom controls the dynamic of the zero mode of the inflaton, and we plan to further explore the implications of our studies for inflationary model building.

Acknowledgments

The work of P.A. and L.P. was supported by the US Department of Energy through grant DE-SC0015655. L.P. gratefully acknowledges support from a Fortner Fellowship at the University of Illinois at Urbana-Champaign. L.P. thanks the Aspen Center for Physics for hospitality and support through National Science Foundation grant PHY-1607611. The work of M.R. and L.S. is partially supported by the US-NSF grant PHY-1820675.

Appendix A Computation of the fermion-gravitational wave interactions

In this appendix we discuss the steps that lead from eq. (2.3) to eq. (2.8). We begin by writing down the vielbeins for the line element in eq. (2.4), $g_{\mu\nu}=\eta_{AB}e^{A}{}_{\mu}e^{B}{}_{\nu}$. These are given by

[TABLE]

where the spatial components $e^{a}\phantom{}_{i}$ satisfy $\delta_{ab}e^{a}\phantom{}_{i}e^{b}\phantom{}_{j}=h_{ij}$. Starting from (A.1), we can also write

[TABLE]

and one can indeed verify that the product $g^{\mu\nu}=e_{A}^{\mu}e_{B}^{\nu}\eta^{AB}$ is the inverse of the metric in eq. (2.4).

In terms of the vielbein, the spin connections are obtained from

[TABLE]

where $\Gamma^{\nu}{}_{\sigma\mu}$ are the usual Christoffel symbols associated with the metric $g_{\mu\nu}$. Using eq. (A.1) and (A.2), we find the components of the spin connection in ADM coordinates

[TABLE]

Inserting these in the action (2.3) we obtain an expanded form of the fermion action in ADM coordinates

[TABLE]

where

[TABLE]

and $\Sigma^{AB}=[\gamma^{A},\gamma^{B}]/4$. The full action is the sum of the bosonic part in eq. (2.5) and the fermionic part in eq. (A.5).

We are interested in the interactions between the fermions and the tensor modes of the metric, and so we expand the spatial part of the vielbein as

[TABLE]

which leads to the following components of the spin connection (expanded here to quadratic order in tensors)

[TABLE]

These relations, along with the expansion of the lapse and shift in eqs. (2.13), are inserted into the constraint equations (2.3) and (2.3). The part of the Hamiltonian constraint equation quadratic in field fluctuations is

[TABLE]

while the quadratic part of the momentum constraint reads

[TABLE]

Eqs. (A) and (A) can be solved to find the quadratic order lapse and shift

[TABLE]

where $\Delta=\partial_{i}\partial_{i}$ is the spatial Laplacian, and $\Delta^{-1}$ is its inverse. In deriving these solutions, we have disregarded fluctuations of the inflaton field because we are only interested in the interactions between gravitational waves and fermions. The inclusion of inflaton fluctuations introduces terms that are quadratic in the inflaton fluctuations, as well as terms quadratic in the first order perturbation to the lapse and shift.444Our solutions for the part of the second order constraints that is quadratic in tensors differs from the results presented in reference [56]; however, this does not affect the conclusions of that work. We have also made use of the linear order equation of motion for the fermion. This induces corrections to the action that begin at fifth order in fluctuations and are thus irrelevant here.

We are now ready to evaluate the action, eq. (A.6), on the constraint surface and eliminate the non-dynamical lapse and shift. Inserting the solutions to the constraints (A.15) into the full action, eq. (2.5) + (2.8), we expand the result to quartic order. This results in an action for the dynamical fields, $\psi$ and $\gamma$, consisting of a quadratic (free) part, $S^{(2)}$ and cubic and quartic parts, $S_{F}^{(3)}$ and $S_{F}^{(4)}$, which describe the interactions of a fermion bilinear with one and two gravitational waves, respectively. The action that results from the procedure just described is very long and complicated. However, in practice there are a number of facts that considerably simplify the result. Because the constraint equations are derived from the variation of the action with respect to the lapse and shift, one can use their equations of motion before substituting in their solutions. This results in the cancellation of a large number of terms, and leaves the result we report above in eqs. (2.19) and (2.20).

Appendix B Interaction Hamiltonian

In this appendix we write the explicit forms of the interaction Hamiltonian terms eq. (2.21), obtained by inserting the decompositions (2.22) and (2.24) into the actions in eqs. (2.19) and (2.20). For the cubic term, we find

[TABLE]

The various terms contributing at quartic order are

[TABLE]

[TABLE]

Appendix C Details of the quartic loop computation

Having identified the seven quartic vertices in eq. (B), we now proceed to evaluate the left loop diagram of figure 1, using the mode functions presented in section 3.1. Evaluating the fermion expectation values produces terms of the form $\delta({\bf k}_{3}-{\bf k}_{4})$, which sets the denominator in several of the vertices to zero. To ensure mathematically sensible results, we follow an approach similar to that of [55], in which we set ${\bf k}_{4}={\bf k}_{3}+\bm{\epsilon}$ (and therefore ${\bf p}_{2}=-{\bf p}_{1}+\bm{\epsilon}$) and take the limit $\bm{\epsilon}\rightarrow 0$, as discussed below.

It is convenient to note that, to quadratic order in $\bm{\epsilon}$,

[TABLE]

which can be found using the explicit expressions for the spinors given in [57]; the lack of a linear term allows for several derivations to be simplified. Similarly, using explicit expressions for spinors one can show that

[TABLE]

We calculate the two point correlation function with a single loop correction; the subscript indicates which of the seven vertices of eq. (B) was used. Four of the results are explicitly independent of ${\bm{\epsilon}}$ at lowest order,

[TABLE]

However, the other three are not explicitly independent of the direction of ${\bm{\epsilon}}$ at lowest order,

[TABLE]

where

[TABLE]

Following reference [55], we impose that the limit is approached from an orthogonal direction, ${\bf p}_{1}\cdot{\bm{\epsilon}}\sim{\epsilon}^{2}$, which makes the problematic terms subdominant.

To separate the integrals, we introduce the new variables $x_{1}=-p_{1}\tau_{1}$ and $y_{1}=-k_{3}\tau_{1}$; we also take $p_{1}\tau\rightarrow 0$, to compute superhorizon quantities. We arrive at

[TABLE]

where $m_{c}\equiv m\cos(\varphi_{0}/f)$ and $m_{s}\equiv m\sin(\varphi_{0}/f)$. Several of the $x_{1}$ integrals are logarithmically divergent and are regulated by the finite amount of e-foldings between when the mode leaves the horizon and the end of inflation. This logarithmic divergence is a consequence of the fact that the fermions have a nonzero average density that continues to source graviton fluctuations even when outside the horizon; analogous behavior was observed in the sourced inflaton perturbations of reference [43].

The remaining integrals are UV-divergent and need to be regularized. We thus introduce a ultraviolet cutoff $\Lambda$. (We discuss renormalization further below.) After some algebra, the expectation values can be expressed in terms of the three integrals

[TABLE]

where the functions $s^{r}(y)$ and $d^{r}(y)$ are defined in eq. (3.1). We have

[TABLE]

In reference [43], integrals $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ were evaluated analytically and exactly. Furthermore, it was found that in the $\mu\ll 1\ll\xi$ range, adiabatic subtraction agreed with the result obtained by simply dropping off the divergent pieces. Outside of this regime, the regularizating term generated by adiabatic regularization dominates the calculated contribution even outside of the UV limit, and furthermore, the value of the term depends on the order of adiabatic regularization. Therefore, as in reference [43], we regularize by simply dropping the divergent pieces, noting that in the regime in which adiabatic regularization is well-behaved, these approaches agree.

We reproduce here the exact analytic results for $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ from reference [43] along with their finite piece in the regime $\mu\ll 1\ll\xi$. The third integral, $\mathcal{I}_{3}$, can be evaluated with the same techniques. The first integral is

[TABLE]

The second integral reads

[TABLE]

Finally, the third integral evaluates to

[TABLE]

In these equations $H_{n}$ denotes the $n$-th harmonic number. From this we find final results for each loop,

[TABLE]

from which we conclude that the dominant contribution in the $\mu\ll 1\ll\xi$ limit reads

[TABLE]

which we also write as eq. (3.18) of the main text. We conclude this appendix by noting that the parity violating vertex that originates from the portion of interaction Hamiltonian denoted by $H_{\rm{int},3}^{(4)}$. The parity violation is evident because this term gives a contribution to the two-point function of the graviton which is proportional to the sign of the graviton helicity; however, the contribution is sub-dominant, scaling as ${\cal O}(\xi^{2})$.

Appendix D Details of the cubic loop computation

In this appendix, we calculate the contribution from the one loop diagram with two cubic vertices shown in the right panel of figure 1, using the interaction term (B.1) and the mode functions calculated in section 3.1. We follow the approach, including the approximations, used in reference [43]. To start with, we drop the external momentum within the loop integral (${\bf k}_{1}+{\bf p}_{1}\rightarrow{\bf k}_{1}$); then after taking the $\tau\rightarrow 0$ we have

[TABLE]

The evaluation of these integrals is complicated by the fact that the integrand is oscillating rapidly. Therefore, we perform a Wick rotation, noting that $\tau_{1}>\tau_{2}$ and, because of the nested integrals, both must be Wick rotated in the same direction. This gives

[TABLE]

where $x_{1}=-p_{1}\tau_{1}$, $x_{2}=-p_{1}\tau_{2}$, and $y_{1}=-k_{1}\tau_{1}$, and functions such as $s^{r}(-iy_{1})^{*}$ should be interpreted as first complex conjugating $s^{r}(x)$ and then substituting $x=-iy_{1}$.

Next we introduce polar coordinates in the $x_{1}-x_{2}$ plane,

[TABLE]

The true region of integration is shown on the left side of figure 2; however, we approximate the integral by using the region of integration shown on the right side. We show below that at fixed $y_{1}$, the integral is exponentially suppressed as $\alpha\rightarrow\pi/2$, and therefore we expect this to contribute at most a $\mathcal{O}(1)$ factor.

After further substituting $\beta=\tan(\alpha)$, we can perform the integral in $d\rho$ to find

[TABLE]

in the $\rho_{*}\rightarrow 0$ limit. This is dominated by the $\log(\rho_{*})$ piece, which is regulated by the finite number of e-foldings between when the mode leaves the horizon and the end of inflation; the same behavior was also observed in the quartic loops in appendix C above. Physically, this logarithmic behavior is a consequence of the fact that the nonzero fermion energy density continues to source graviton fluctuations even once it is outside the horizon.

Finally, we approximate the Whittaker functions as in reference [43]. Along the positive axis we use (defining $M\equiv\sqrt{\mu^{2}+4\,\xi^{2}}$)

[TABLE]

and we note that the Whittaker functions are even in their second index. Along the negative axis, we use

[TABLE]

where

[TABLE]

With these approximations, we find

[TABLE]

where we have defined

[TABLE]

Eq. (D) shows that, as claimed, at fixed $y_{1}$, this contribution to the graviton two-point function is exponentially suppressed as $\beta\rightarrow\infty$ ($\alpha\rightarrow\pi/2$). The $dy_{1}$ integral may now be performed analytically, leaving the $d\beta$ integral. We separate this into parts involving $\mathcal{A}$ only, $\mathcal{B}$ only, and mixed. From the form (D.6), we recognize that the $\mathcal{A}$ terms correspond to the vacuum part of the modes, which are nonzero and positive frequency as $p\rightarrow\infty$. We subtract this part by hand; we also note that it is only the purely $\mathcal{A}$ part which has a divergent piece. We regularize this by writing the lower limit of the $d\beta$ integral as $1+\epsilon$ and taking the limit as $\epsilon\rightarrow 0$.

Taking the large $\xi$ limit, the purely $\mathcal{A}$ piece which we renormalize away is

[TABLE]

where we have dropped a term that, in the large $\xi$ limit, is proportional to $\xi^{-2}\epsilon^{-4}$, and where $H_{n}$ is the $n$-th harmonic number. The mixed term is

[TABLE]

while the purely $\mathcal{B}$ contribution is

[TABLE]

where $B_{n}(x,\,y)$ denotes the incomplete beta function. This term has a piece which scales as $\mu^{2}/\xi$ and a piece which scales approximately as $\mu^{4}\xi$. In the regime $\mu\lesssim 1$ and $\xi\gg 1$ either can dominate, depending on how small $\mu$ is. Regardless, the mixed term (D.11) dominates, and the leading contribution from the cubic diagram is

[TABLE]

where we have noted an overall uncertainty due to the approximations used. We also note that $\rho_{*}\sim-p_{1}\tau$ from the change to polar coordinates. This expression is reported as eq. (3.21) of the main text. This is the same order of magnitude as the quartic result in eq. (C.13), although with the opposite sign.