Fermion production during and after axion inflation

We study derivatively coupled fermions in axion-driven inflation, specifically $m_\phi^2\phi^2$ and monodromy inflation, and calculate particle production during the inflationary epoch and the post-inflationary axion oscillations. During inflation, the rolling axion acts as an effective chemical potential for helicity which biases the gravitational production of one fermion helicity over the other. This mechanism allows for efficient gravitational production of heavy fermion states that would otherwise be highly suppressed. Following inflation, the axion oscillates and fermions with both helicities are produced as the effective frequency of the fermion field changes non-adiabatically. For certain values of the fermion mass and axion-fermion coupling strength, the two helicity states are produced asymmetrically, resulting in unequal number-densities of left- and right-helicity fermions.


Introduction
In this work we study the behavior of fermionic degrees of freedom during and immediately following inflation driven by an axionic field (for a recent review of axion inflation see [1] and references therein). Specifically, we consider Majorana fermions coupled derivatively to the inflaton field via a dimension-five operator. When the fermions have degenerate masses and are paired to make Dirac fermions, this is the coupling to the axial-vector current, a conserved current in the massless limit (in the absence of anomalous gauge couplings).
In local thermal equilibrium, this coupling between the fermion fields and the rolling axion acts as a chemical potential for helicity, allowing the system to reduce its free energy by converting some fermions of one helicity into the other. The universe is not in thermal equilibrium during inflation, however, in this work we demonstrate that this coupling has a similar effect on the gravitational production of fermions. During inflation, the exponential expansion of spacetime produces particles from the vacuum as their wavelength becomes comparable to the Hubble length. In the absence of the axion coupling, all spin, or helicity states are populated with a probability that depends only on the fermion mass, m ψ , and the Hubble rate, H. This gravitational pair production is most efficient for small masses, with the production rate for larger masses being exponentially suppressed by the ratio m ψ /H. In the limit that the fermions are massless, no particle production occurs because the theory is conformally equivalent to a Minkowski space theory (see e.g. [2]). We show that the -1 -

JCAP11(2015)021
coupling to the rolling axion leads to the biasing of this gravitational fermion-production, which enhances the production rate of one of the helicities while suppressing the other. Furthermore, this enhancement can partially counter the exponential suppression of heavy masses allowing for efficient production of fermion states that would otherwise be heavily suppressed. On the other hand, in the limit that the fermion mass vanishes, the axionfermion coupling reduces to a boundary term which has no effect on the equations of motion and the effect vanishes.
The significant feature arising from considering a pseudo-scalar axion field as the inflaton, that we will discuss here, is the nature of the derivative coupling to fermions. The specific form of the axion potential appears to play a role in the details, but does not change the mechanism and the qualitative results of the paper. It is shown that even the final results only differ by an O(1) factor for a significant parameter range in two different forms of the inflaton potential we considered.
Derivative couplings of vector gauge-field degrees of freedom to rolling pseudo-scalar fields [3][4][5][6][7][8][9] lead to strong growth of one of the polarizations of the gauge fields during de Sitter expansion, but the effect of such couplings on fermion fields has been little studied. We demonstrate that in the fermionic case there is a similar effect, however, Pauli blocking prevents the occupation numbers from exceeding unity and consequently the growth of the fermion field is constrained.
Following the inflationary epoch, the axion oscillates about the minima of its potential which causes the effective frequency of the fermion fields to vary non-adiabatically. These oscillations cause both helicities to be produced in turn as the effective momentum of the leftand right-helicity fermions vanishes. However, the axion oscillations damp as the universe expands which reduces the efficiency of successive particle-production events. In combination with the inflationary production, this damping can lead to asymmetric number densities of fermions with left and right helicities. Moreover, since neutrino helicity is equal to lepton number in the standard model, this effect may be relevant for the generation of the matterantimatter asymmetry in the universe -we pursue this question in a companion paper [10].
It is well established that fermions can be produced via Yukawa couplings to oscillating scalar fields after inflation [11][12][13][14][15]. Preheating in these models is, in some instances, known to produce fermions with masses much larger than the inflaton mass and may be important for leptogenesis scenarios where the baryon asymmetry is produced by the decay of right-handed neutrinos with masses near the grand unified energy (GUT) scale [16,17]. Further, ref. [18], studied neutrino production during the relaxation and oscillations of the Higgs condensate following inflation as a mechanism for the generation of a matter-antimatter asymmetry via leptogenesis. The cosmological consequences of the gravitational production of heavy fermionic states during inflation is well appreciated [19], and gravitational production has been considered as a mechanism for producing super-heavy fermions during inflation [20,21]. The behavior of fermions coupled to axions in the post-inflationary universe was first considered in the context of spontaneous baryogenesis [22][23][24][25][26] and has more recently been revived for leptogenesis scenarios [27,28]. Our work on fermion production differs from these existing works because we consider both the production during the inflationary period due to the axion coupling, which has significantly different phenomenology, as well production following inflation during the reheating phase.
We work in natural units where c = = 1 and denote the reduced Planck mass by -2 -

JCAP11(2015)021
2 Background and conventions Our notation and conventions for the axion, gravitational, and fermion sectors of the theory are as follows. We consider an axion inflationary sector minimally coupled to Einstein gravity where V (φ) is a potential capable of supporting an extended period of slow-roll inflation. For concreteness we consider the simplest chaotic inflation model with a quadratic potential [29,30] as well as the potential from the simplest type of axion-monodromy inflation [31,32] In addition to the inflaton sector, we consider a set of fermions 1χ i , described by the action where m ij is hermitian, C i is real, and f is a mass scale associated with the axion. The covariant derivative of the spinor fields is where the spin connection is We work in cosmic time and ignore metric fluctuations as well as inflaton fluctuations. While fluctuations of the axion and gravitational fields are likely to be interesting, our primary interest in this work is the behavior of the fermion fields in the classical, homogeneous axion and gravitational backgrounds. We work in 'mostly plus' convention for the metric with where a(t) is the scale factor and e a µ are vierbeins. We rescale our fermion fields by a −3/2 to absorb the factor of √ −g in the action and turn the covariant derivative into a partial derivative. We write χ i = a 3/2χ i so that the fermion action reads We work with two-component spinors and consider only fields that transform under the left-handed representation of the Lorentz group. Right-handed fields can be obtained simply by complex conjugation, that is, if χi is a left-handed field, χ † i is right-handed. See appendix A for some key notations used in this work. Refer to [33] for an extensive review.

JCAP11(2015)021
For simplicity, we consider only a pair of fermions, i ∈ {1, 2}. Further, we treat separately the case where the mass matrix m ij has degenerate eigenvalues, in which case there is an O(2) flavor symmetry and the fields can be paired to make a Dirac fermion.
Above and in what follows, Greek letters from the middle of the alphabet represent four-dimensional space-time indices, while Greek letters from the beginning of the alphabet denote spinor indices. Roman letters from the start of the alphabet denote four-dimensional Lorentz indices, while Roman indices from the middle of the alphabet denote spatial indices as well as flavor indices. Einstein summation convention is used for spacetime, Lorentz and repeated flavor indices, unless otherwise noted. Paired ascending dotted indices and descending undotted spinor indices will also be summed.

Majorana and Dirac fermions
There are subtle differences between the cases where the fermions can be combined into a single Dirac fermion and can carry charged currents compared with the case where they are Majorana fermions. We treat these cases separately in what follows and then consider their vacuum states and quasi-particle numbers in expanding Friedmann-Robertson-Walker spacetimes.

Non-degenerate fermions: Majorana fermions
First we consider the case where the fermions have non-degenerate masses and there are no conserved currents. Without loss of generality, we work in a basis where m ij is diagonal, and the fermions have masses m i . Variation of the action at eq. (2.8) with respect to χ † iα yields the equations of motion for these Majorana fermions It is convenient to expand each field into a Fourier basis, where we have introduced creation and annihilation operators, a λ ik and a †λ ik , which satisfy the anti-commutation relations with all other anti-commutators vanishing, as usual. Here and throughout we use λ = ±1 to denote the spin or helicity states of the fermions. It follows that We quantize the fields by imposing the anti-commutation relations where the canonical momenta of the fermion fields are found in the usual way where here and throughout an overdot on a field denotes a derivative with respect to cosmic time, i.e.χ iβ = ∂ t χ iβ . After inserting the mode expansions, the fields are canonically quantized provided that the fermion wavefunctions are normalized according to (no sum on i) Substituting the mode expansions, eq. (3.2) and (3.4), into the equation of motion, eq. (3.1), we find the coupled equations for x β and y †β We work in a basis of helicity eigenspinors, which satisfy where ι λ k is a phase that satisfies Writing the spinors as (3.12) canonical quantization requires that (no sum on i) In terms of these wavefunctions, the equations of motion (3.8) and (3.9) become These coupled equations can be decoupled to give two second order ordinary differential equations for each wavefunction where We are interested in particle creation, and so we compute the Hamiltonian density After Fourier transforming, inserting the mode expansions, and some straightforward algebra, the Hamiltonian is given by where . (3.20) Note that the presence of the phase ι λ k prevents the off-diagonal terms in the Hamiltonian in eq. (3.19) from vanishing. The diagonal and off-diagonal parts of the Hamiltonian are related by

Degenerate fermions: Dirac fermions
In the case where the mass matrix has degenerate eigenvalues there is a global internal O(2) flavor symmetry χ i → O i j χ j , where O is an orthogonal matrix, provided we assume that C i = C, that is, that the axion couples universally to the fermions. There is a conserved hermitian Noether current associated with this symmetry However, the flavor basis we have been using is off-diagonal with respect to this current and thus the flavor states are not eigenstates of the charge operator Q = d 3 xJ 0 . It is more useful to work in a basis of states of definite charge. We can achieve this by diagonalizing the Noether current. Redefining our fields diagonalizes the current in eq. (3.22). In this basis, the action reads

JCAP11(2015)021
and the SO(2) flavor symmetry (the part continuously connected to the identity) is realized in this basis as the U(1) symmetry ϕ → e iθ ϕ and η → e −iθ η, where θ is the angle that defines the SO(2) rotation matrix [33]. The spinors ϕ and η satisfy the free-field Dirac equations Together ϕ and η † compose a Dirac fermion, which we can expand into plane waves as where a † , b † , a, and b are creation and annihilation operators that satisfy the anticommutation relations with all others vanishing. It follows that Substituting these mode expansions into the equations of motion, we find the same equations of motion as in the non-degenerate Majorana case We proceed in the same fashion as the non-degenerate mass case, and work in a basis of helicity eigenspinors, as above Then, the equations of motion for the wavefunctions of the Dirac fields are identical to eq. (3.14) derived earlier for each flavor in the non-degenerate system,

JCAP11(2015)021
However, note that for two fermions there is only a single pair of equations, unlike the Majorana case where there is a pair of equations for each flavor. This is ultimately because Dirac particles are related to their antiparticles by charge conjugation; the U(1) symmetry relates the fields to each other. We can construct the Hamiltonian for this theory analogously to before where we have defined , (3.38) which is simply the sum over the flavors of the Hamiltonian at eq. (3.19) where m i = m. The relations between these functions for Majorana and Dirac Hamiltonians are (3.40) and the corresponding consistency relation in the Dirac case becomes

Vacuum states
The Hamiltonian in both eq. (3.19) and eq. (3.37) can be diagonalized at any given time t 0 by choosing the wavefunctions X λ k and Y λ k so that the function F λ k vanishes. A straightforward calculation shows that F λ k = 0 whenever where φ is an arbitrary phase. Here, we are interested in modes that start in the Bunch-Davies vacuum state during inflation. In order to determine the initial condition, we look for Wentzel-Kramers-Brillouin (WKB) solutions to the equations of motion. Assuming that |ω/ω 2 | 1, it is straightforward to derive the lowest order WKB solution to eqs. (3.14)-(3.16) where as above

JCAP11(2015)021
Substituting these solutions into eq. (3.20) we obtain for the Majorana case, while substituting them into eq. (3.38) gives for the Dirac case.
On the other hand, canonical quantization requires that the coefficients satisfy Thus, in order for the fields to be canonically normalized and in the lowest energy state (the state with the lowest vacuum energy density), we are required to take Notice that in the absence of the coupling to the axion (whenk λ = kλ/a), the first term in the vacuum energy in eq. (3.45)-(3.46) vanishes when summed over spins. When the coupling to the axion is turned on, this term no longer vanishes, and leads to a splitting of the vacuum levels for the λ = +1 and λ = −1 helicity states. The choice of initial conditions at eq. (3.48) diagonalizes the Hamiltonian at some initial time, but in general the equations of motion drive F λ i (k) away from zero at times t > t 0 . At these later times, the Hamiltonian can be re-diagonalized after performing a (time-dependent) Bogoliubov transformation on the creation and annihilation operators.
We begin by constructing the Bogoliubov transformation for the Majorana case, by first writing the Hamiltonian in the following form This Hamiltonian can be diagonalized by choosing where α k (t) and β k (t) are complex functions of time that depend on the magnitude of the wavenumber. Under this Bogoliubov transformation, the creation-annihilation operator pairs becomeã †

JCAP11(2015)021
The off-diagonal part proportional to a k a −k is where we supressed helicity and flavor indices. Diagonalizing the Hamiltonian amounts to setting the off-diagonal terms equal to zero. Solving for α is now straightforward, setting eq. (3.52) to zero yields where we used ι −k = −ι k , and took E 2 + 4|F | 2 = ω 2 /4. Using the relation |α| 2 + |β| 2 = 1 leads to where we used the fact that |A iλ | 2 + |B iλ | 2 = 1/4 when not summed over spins. We can proceed in a similar fashion for the Dirac case, where the Hamiltonian can be written as where Note that the Dirac Hamiltonian is comprised of two identical copies of the Majorana Hamiltonian, as expected. The form of the Bogoliubov transformation can be read off as as, for example, in ref. [13]. The complex phases ι k are absorbed into the parameter β in this case, since -in contrast to the Majorana case -we do not need to invert the wavenumber of the operators as defined in the Bogoliubov transformation, so the antisymmetric nature of ι k for the relabelling k → −k is irrelevant. Performing the Bogoliubov transformations on each part of the Dirac Hamiltonian, we obtain (suppressing flavor and helicity indices for clarity) The off-diagonal terms (proportional to b −k a k ) in this expression are

JCAP11(2015)021
Diagonalizing the Hamiltonian requires that these terms vanish and gives Using the relation |α| 2 + |β| 2 = 1 we get as in the Majorana case. We use these time-dependent operators (theã † andb † ) to define (time-dependent) Fock spaces built from the zero (quasi)particle states As noted above, a state that starts off at time t 0 as the vacuum state, i.e. a state containing no particles, evolves to a state which at a generic later time contains particles. This is because the particle number at a given moment is provided by operators which instantaneously diagonalize the Hamiltonian. This instantaneous particle number can be computed by projecting the full mode-functions onto an instantaneous WKB basis (eq. (3.43)), and extracting the negative frequency parts. For both Majorana and Dirac Fermions, the particle number can be evaluated as We express this quasi-particle number in a more useful form, by noting that the particle number can be extracted by comparing a given solution of the Dirac equation to an instantaneous WKB solution of the form eq. (3.43). The square of the coefficient of the negative frequency mode is the particle number and can be extracted via Note, however, that this particle number is only well defined where the WKB approximation is valid. That is, during times where |ω/ω 2 λ | 1.

Inflationary solutions
During the inflationary epoch, the accelerated expansion of spacetime leads to particle production of all fields that are not conformally invariant [2]. In this section we derive analytic solutions for the quasi-particle number by directly solving the Dirac equation for the fermion wavefunctions.
During slow-roll inflation we can treat the inflationary spacetime as de Sitter space to a good approximation. In this section, we work with conformal time, which is taken to be a negative increasing quantity during inflation We assume that the axion rolls at an approximately constant rate in cosmic time during inflationφ Hφ, and we introduce the parameter

JCAP11(2015)021
In conformal time, the equation of motion, eq. (3.16), for the spinor fields becomes where we have treated H ≈ constant andφ ≈ constant, dropping time derivatives of these quantities, and we have rescaled Redefining the time coordinate to u = 2ikτ puts the equation into standard form of the Whittaker equation We make use of the fact that in the limit z → ∞, the Whittaker functions have limiting forms [34] where θ 1 and θ 2 are irrelevant phases. Thus, choosing the Bunch-Davies vacuum implies we should set B λ ik = 0 in eq. (4.6). Next, we find Y λ * ik in terms of X λ ik and thus Φ λ ik by using the Dirac equation, eq. (3.14).
To solve this equation, we make use of the fact that the Whittaker functions satisfy the identities [34] z d dz z where (a) n is Pochhammer's symbol, which satisfies (4.14) For n = 1, eqs. (4.12) and (4.13) are 16) and note that, for λ = −1, The full solutions to the Dirac equation that match onto the vacuum in the UV are then given by It remains to fix the overall normalization of these modes according to (no sum on i) Making use of eqs. (4.12) and (4.13), the fact that the Whittaker functions satisfy the connection formula 20) and have Wronskian, the normalization condition at eq. (4.19) reduces to the condition Therefore, our canonically normalized fermion wavefunction solutions are where θ and θ are arbitrary phases.  Figure 1. We display the dependence of the quasi-particle number (for modes with k < aH) on the mass of the fermion (in units of the hubble rate) and the coupling strength to the rolling axion. In the left panel we show the particle number for fixed values of ϑ = 10, 0.05, 10 −3 (curves from from right to left) as the fermion mass is varied. In the right plot we show the particle number for fixed m/H = 0.01, 0.1, 0.5 (curves from top to bottom) as the coupling to the rolling axion is varied. In both cases, we have chosen ϑ > 0 so that n + k (red curves) is enhanced while n − k (blue curves) is suppressed.
In the absence of the axion coupling, there are well-known solutions to the Dirac equation in de Sitter space, see for example refs. [19,[35][36][37]. These results can be obtained after some algebra from eq. (4.23) by making use of the relation between the Whittaker and Hankel functions [34] (4.24) and the recursion relations for the Whittaker functions Using the analytic expressions in eq. (4.23), we derive analytic expressions for the quasiparticle number during inflation. We are interested in modes late in inflation that have left the horizon. For modes which satisfy k/aH 1, and assuming m i = 0, we find 2 Eqs. (4.23) and (4.27) are the main result of this section.

JCAP11(2015)021
In the absence of the coupling to the axion, production of both fermion helicities is symmetric, as expected, and highly suppressed for fermions with masses larger than the Hubble rate. For small masses, the occupation number approaches its maximum value of 1/2 as m i /H → 0. However, note that for m i = 0, the theory is conformally equivalent to Minkowski space, and no particle production occurs. 3 When the coupling to the axion is switched on, the particle production here is asymmetric between the helicity states. For ϑ > 0 (ϑ < 0), particle production of the λ = + (λ = −) helicity state is enhanced while particle production of the λ = − (λ = +) mode is suppressed. In figure 1 we display the dependence of the occupation probability on the quantity ϑ, defined in eq. (4.2), which encodes the effective axion-fermion coupling, and the mass of the fermion in units of the Hubble rate.
In the case where the fermions are degenerate in mass and combined to make a single Dirac fermion, we note that while the axion coupling leads to the asymmetric production of helicity states, this coupling keeps the number of particles equal to the number of antiparticles, so that the overall Noether charge is conserved. In this case, if left-handed particles are produced, so are equal numbers of right-handed anti-particles.

Post-inflationary evolution and reheating in quadratic chaotic inflation
So far, we have only solved for long-wavelength modes of the fermion field in a basis of quasiparticles states. These states do not coincide with an observable basis until the Universe is in a matter-or radiation-dominated phase. 4 In the absence of instantaneous reheating, at the end of inflation the inflaton oscillates about the minima of its potential. As the axion oscillates, the effective frequency of the fermion modes (see eq. (3.17)) changes non-adiabatically when the effective wavenumber,k λ (t), vanishes Around these times we expect particle production to occur. Note that for one of the helicity states, this non-adiabatic evolution occurs for the first time during inflation and corresponds to the particle production considered earlier in section 4. While exact analytic calculations are difficult in this regime, we can gain some insight by making use of the WKB approximation. This is the calculation we turn to next. The equations of motion for the Majorana and Dirac cases are identical, therefore we treat them uniformly, unless otherwise noted. In order to simplify notation for the rest of this section, the comoving wavenumber k and the fermion mass m i are measured in units of the axion mass m φ . Similarly, time is measured in units of m −1 φ and f and φ are measured in units of M Pl . We also drop the subscript i on the fermion mass, because we focus on only one flavor of fermion.

Static universe approximation
In order to gain intuition about the post-inflationary processes in this model, we first work in the limit where the universe does not expand. In this static universe limit, the oscillations

JCAP11(2015)021
of the axion about the minima of its potential are harmonic, φ(t) = φ 0 sin(t), and we can use the methods outlined by Peloso and Sorbo [13] to study the evolution of the fermion particle number. Peloso and Sorbo used a Yukawa coupling of the inflaton to the fermion sector, however, the equations for the axial coupling considered here can be mapped directly onto the Yukawa case by identifying Proceeding analogously to ref. [13] we look at the points of non-adiabatic behavior, wherẽ k(t * ) = 0. These points occur whenever kλ = −(C/f )φ for the wavenumber of interest. For the rest of the static universe analysis we choose λ = 1 without loss of generality. We thus consider modes for which k < (C/f )φ 0 , wherek(t * ) = 0 is possible, and work in the limit C f φ 0 1 in order for the WKB approximation to be valid. Defining the time t * byk(t * ) = 0, we expand the effective wavenumber near these points as where A = C f φ 0 . These zero-crossings happen multiple times and differ based on helicity and wavenumber, so we should be writing t * = t λ i (k). However, we keep t * , to not overly complicate our notation. Introducing the parameters the equation of motion for the fermion wavefunctions is approximated near t * as If the fermion field started in the vacuum state initially (t < t * ), then the particle spectrum after the point wherek = 0 for the first time (t > t * ) is given by [13] n (1) where the superscript (1) denotes the firstk zero-crossing, henceforth also called a production event. We call the combination πp 2 i the fermion production exponent at the i'th production event. This is independent of i in the static universe case, but not once we take into account the effects of the expansion. In the static universe approximation, the particle number after the first production event becomes

JCAP11(2015)021
Fermion production is exponentially sensitive to the square of the fermion mass and the inverse of the axion-fermion coupling constant. Production is therefore significant when πm 2 /A 1 and effectively shuts off when πm 2 /A 3. A important characteristic of the present model is that each wavenumber undergoes particle production at a different time, defined byk = 0, which for the static universe case is simply t * = arccos (−k/A). In the Yukawa-coupled model studied in ref. [13], particle production happens uniformily for all wavenumbers k when the effective fermion mass crosses zero.
We have thus far only considered the first instance of fermion production where many modes are uniformly populated up to a maximum wavenumber set by the axion coupling. There are many oscillations of the background axion field, and we are primarily interested in the occupation numbers after a very long time. For bosons, occupation numbers grow exponentially due to Bose-enhancement and parametric resonance. However, fermions obey Fermi-Dirac statistics which limits their occupation number to be less than unity. Parametric resonance 5 and Pauli blocking in the fermion case leads to chaotic filling and emptying of fermion states. In the static universe approach, the time averaged particle number count can be evaluated semi-analytically. Following the analysis of Green and Kofman [11,12], we define a locally averaged particle number count where T is the period of φ(t). The amplitude F k and frequency ν k are (5.10) The function Z k t=0 = 0. The WKB method used in ref. [13] provides a different construction for the particle number after one full inflaton oscillation. The two different analytical approximations to the fermion number produce a final result of the same form n k ≈ F k sin 2 (ν k t). However, these analytic approximations are originally constructed for different regimes. Peloso and Sorbo worked in the broad resonance regime 6 and used the WKB approximation, which produces results only for k < (C/f )φ 0 , which is not a restriction for the analysis of Green and Kofman. However, the method of Peloso and Sorbo extends to the expanding universe case, as we see in later sections. The agreement between the two formulas is excellent for (C/f )φ 0 > 1 (the "broad resonance" regime) and there is an increasing difference as the coupling decreases (the "narrow resonance" regime). We tested these approximations against results obtained by solving the Dirac equation numerically on the homogeneous axion background for the broad resonance regime, where both formulae give identical results. Figure 2 shows that the agreement between the approximate and numerical results is exact at the points t = 2πn (where n is an integer). That is, the approximation exactly reproduces the particle number after an integer number of inflaton oscillations. Figure 3 shows the evolution of the fermion JCAP11(2015)021  Left panel: evolution of the particle number for m = 1, C f φ 0 = 1 and k = 1 (blue), k = 0.7 (red). The black lines show the particle number envelope given by eq. (5.9), (5.10). The vertical lines correspond to points where t is an integer multiple of 2π. Right panel: envelope function F k (black) for m = 1, C f φ 0 = 1 along with the particle spectra for t = 2π (blue), t = 4π (red) and t = 20π (green). In the n k (π) plot the green curve is multiplied by 3, the black curve is multiplied by 10 and the brown curve is multiplied by 300. In the n k (7π) and n k (8π) plots the brown curves are multiplied by 10. All curves are shifted vertically by 1, 2, 3, and 4 for visual clarity. The vertical line at k = A = 50 shows the maximum produced wavenumber according to the WKB method.
particle number for a system that is set at the threshold between broad and narrow resonance (C/f )φ 0 = 1. In figure 4 we show the results of numerically solving for the evolution of the particle spectra for different fermion masses with A = C/f φ 0 = 50, chosen to be well within the region of validity of the WKB results. From comparison of the blue and red curves in the upper right panel of figure 4 with the corresponding curves in the lower panels, it appears that the mass dependence of the final particle spectrum is weaker than that of the initial particle spectrum (the particle spectrum after the first zero crossing of the effective wavenumber, k = 0). Further, for a range of masses, in our example m 5, the particle spectra after an even number of zero-crossings ofk are more similar than their initial spectra would suggest. Lighter fermions are more easily produced in the first oscillations, yet the final occupation numbers appear to be similar, with the occupation number of several wavelengths reaching unity. This effect is due to the nature of fermion production. The occupation number of each mode n k cannot exceed unity, therefore systems with initial particle spectra of n k ∼ 1 can be either reduced or stay almost constant. This behavior causes the particle spectrum to develop fine band structure as parametric resonance develops, which can be averaged out when calculating the total particle number.
For light fermions with m 1, n k ≈ 1 for all excited modes after the first production event. In these cases, the second zero-crossing can only de-populate modes, since modes that are already at n k ≈ 1 cannot be further populated. This behavior results in a decrease in  The blue curve corresponds to the first zero-crossing, while the red and green curves correspond to n k (t = 17π) and n k (t = 18π) respectively. the averaged particle occupation number after the second production event; modes that were n k ≈ 1 after the first event end up with occupation numbers n k 1 after the second event. This effect can be seen in figure 4, especially for the blue dotted curve.
Systems with heavier fermions have an initial particle production rate that is much less efficient, resulting in smaller initial particle spectra, n k 1. However, for these modes, each particle production event is more likely to create more particles than destroy them, resulting in an occupation number that grows with time. Given enough time these models can produce heavy fermions more effectively (at least for certain wavenumbers) than their initial production rate would suggest. In these cases, the final result is a particle spectrum that has broadly similar features to its lighter counterparts, as can be seen from the green curve in figure 4. The important difference between the final particle spectra for fermions of different masses is in the range of wavenumbers that reach occupation numbers of n k ≈ 1 rather then in the amplitude of the occupation number (which is n k ≤ 1 by definition for fermions).
This picture changes dramatically once the expansion of the universe is taken into account. We show in the next section that the expansion introduces new effects which provide the means for generating a significant left-right helicity asymmetry in heavy fermions. Before proceeding to the expanding universe case, we calculate the average occupation number per unit occupied volume in momentum space where k max is the maximum wavenumber that is excited according to the WKB approximation. For large values of the coupling, to a good approximation no particle production occurs for k > k max , and eq. (5.11) provides an estimate of the (weighted) average occupation number of each bin in k-space. We plot this quantity in figure 5 for A = 50. However, we caution that this quantity is dependent on k max itself. For small masses, where particle production is more effective, we see a large variation in the particle number after integer and half-integer numbers of axion oscillations. However, this is not true for larger masses where, as discussed above, after each axion-crossing the particle number is predominantly rising.

JCAP11(2015)021
In the expanding universe case, this oscillating behavior of the average particle number is suppressed, since particle production after each production event is less efficient than the previous one due to the fact that the amplitude, and thus velocity of the axion oscillations, is damped by Hubble friction. We end this section on the static universe approximation with a comment about the dependence of the maximum produced wavenumber on the coupling strength. According to the WKB approximation, particle production occurs only for k/A ≤ 1. However, from figures 2, 3, and 4 note that this limit is increasingly violated as the coupling gets smaller. For couplings approaching unity there is significant particle production for k/A > 1, which is absent for larger couplings A 1. Parametric resonance allows for the production of these high-k modes.

Expanding universe
The static universe approximation is a useful simplification that allows us to build intuition about the relevant physical processes without the complications of expanding space. However, away from the individual particle production events, it is not a good description of the processes occurring during and after inflation. Following inflation, the oscillation frequency of the background inflaton field is typically not much faster than the expansion rate, and consequently it is important to take expansion into account. As we demonstrate in this section, there are two main differences between the static-and expanding-universe cases. Both the range of wavenumbers that get excited and the particle production strength are progressively reduced by the expansion, in contrast to the static universe, where they remain constant.
The main goal of this section is to estimate the final occupation number of the fermion states. Additionally, we calculate the initial and final fermion production asymmetry between the two helicity states. Since one helicity is being produced for the first time during inflation, while the other is produced only during preheating, we will demonstrate that unequal numbers of each helicity are produced in the overall evolution. We show that, because the behavior of both helicity states during the preheating phase is comparable and highly stochastic, the early behavior determines the ultimate fate of the asymmetry. Provided this initial asymmetry is large enough, we demonstrate that, at least in the model at hand, it can be preserved throughout the preheating period. This is quantified in section 5.2.3.
Analogously to the WKB-based static universe approach, particle production in expanding spacetime occurs in distinct incidents at the times whenk(t * ) = 0, corresponding to Fermion production occurs in a very narrow region aroundk = 0, therefore the expansion of the universe can be neglected during the production event and the static universe result of eq. (5.7) can be immediately applied in this case. Different wavenumbers are excited at different times, according to eq. (5.12) and as shown schematically in figure 6. The dependence of the particle production exponent on helicity and wavelength is hidden in the term where the time t * depends on k through eq. (5.12).

JCAP11(2015)021
Despite describing particle production starting from a vacuum state, the form of n k given in eq. (5.7) is useful for calculating statistical properties of the particle spectrum, even after many production events and the emergence of a fine band structure, as shown in figure 4 for the static universe case. We generalize this form in the following sections to include the case of a(t) = const and multiple productions. We consider n k as given in eq. (5.7) to be an indicator of the strength of any single particle production event.

During inflation
As discussed earlier and in section 4, the axion-fermion coupling is active during the inflationary epoch. For one of the helicities (in our case, withφ < 0 or equivalently ϑ > 0, the λ = +1 helicity)k λ may vanish during inflation resulting in particle production during the inflationary epoch itself. In this case, the time of the production event, wherek(t * ) = 0, is found to be where we used the approximation of exact de-Sitter space (H = const.) together with the slow-roll approximationφ = 0. We see that forφ < 0 only the positive helicity state (λ = 1) can be excited during inflation. The resulting fermion number with momentum k can then be evaluated from eq. (5.7) yielding We compare this result to the exact solution derived for particle production during inflation given in eq. (4.27). Working in the limit ϑ 1 (actually we do not need it to be much bigger than one, simply a few times larger), as well as m = O(H), eq. (4.27) becomes eq. (5.15) for the growing (λ = +1) state. This WKB analysis agrees with the exact solution obtained for the modes during inflation, in the appropriate limit of large coupling (or ϑ 1), where the WKB method is applicable, leading to a unified treatment of the inflationary and postinflationary fermion production for the parameter space of interest. It is worth noting that the WKB method employed here is valid only in the large-coupling regime, and therefore it cannot capture the particle production occurring solely due to the expanding space-time.
The mode that is excited at the end of inflation can be read off from eq. (5.12) to be k = −(C/f )λφ(t end )a(t end ). In the case of chaotic inflation with a quadratic potential −φ(t end )a(t end ) ≈ 0.7.

After inflation
Following the inflationary epoch, the axion oscillates about the minima of its potential. The particle production due to each zero-crossing ofk (neglecting the emerging band structure due to stochasticity) is governed by eq. (5.7) up to a maximum wavenumber given by We consider first the λ = −1 helicity state which, for chaotic inflation with our conventions, is excited only after inflation has ended and the axion has crossed zero. The maximum particle number is produced for low wavenumbers, since |dk/dt| t * is a decreasing function of k. For low wavenumbers, k C/f , the conditionk = 0 becomesφ(t * ) ≈ 0 and particle production occurs predominantly around points where the axion velocity vanishes.  The red-dashed and brown-dot-dashed lines correspond to modes for λ = 1 that get excited for the first time during and after inflation respectively. Particle production occurs wheneverk crosses zero. Note that these particle production events occur in pairs which are closely spaced for high wavenumbers.
inflation with a quadratic potential,φ(t * ) ≈ 1/3 during the first instance ofφ(t * ) = 0 and a simple application of eq. (5.13) and eq. (5.7) leads to In the region of C/f 1 the particle spectrum n k approaches a step function in wavenumberspace with a sharp cut-off at k = k − max , hence the above formula, which suppresses the k−dependence, provides a useful estimate. Despite being a crude approximation, it provides an increasingly good estimate of the particle number for decreasing values of the combination m 2 /(C/f ) 1 and it can provide an upper limit on the particle production efficiency. We now consider the positive helicity state which becomes excited both during and after inflation. The range of modes excited during inflation is k 0.7C/f for chaotic inflation, as shown in the previous section, while the total range of modes excited is k k + max ≈ 0.9C/f , obtained from eq. (5.16). The difference of k + max and k − max plays a key role in our discussion of asymmetry in section 5.2.3.
An important characteristic of fermion preheating is the evolution of the Fermi sphere, defined in this case as the range of comoving wavenumbers that are (partially) filled with fermions as the universe expands during preheating. We again consider only the large coupling regime C/f 1, in order for the WKB approximation to be valid. For a matter dominated universe, arising from an axion potential that is predominantly quadratic during preheating, the inflaton amplitude decays as φ ∼ a −3/2 and we may approximatė

JCAP11(2015)021
where δθ is an arbitrary constant phase. We then find a simple expression for the maximum excited wavenumber during a given production event Since the scale factor increases, the highest-momentum state that can be excited decreases as time progresses. For a locally quadratic minimum, the derivative of the effective wavenumber is given by (5.20) and thus the particle number can be written Note that the particle production rate decreases as a function of time. When taken together with eq. (5.19), eq. (5.21) implies the first few production events are the most important for determining the characteristics of the particle spectrum produced in this model. These events produce the most particles out to the largest wavenumber with the greatest efficiency. Subsequent production events can only affect smaller and smaller wavenumbers and with a decreasing strength. These subsequent events simply introduce a band structure that can be averaged out for all practical purposes, analogously to the static universe case. Finally, we note the difference between the case considered here and the case of Yukawacoupled light fermions (m m φ ). In the limit of large Yukawa coupling, the maximum comoving wavenumber that is excited grows with the expansion as k max ∼ a 1/4 where a(t) is the scale-factor, normalized as a = 1 at the end of inflation [12]. This result can be derived by noting that the particle production probability is given by and thus the time dependence of the maximum wavenumber excited can be read off as k max ∼ a 1/4 . In the case of heavy fermions (m > m φ ) the result is more complicated, giving the final form of n k as a shifted Gaussian in k (when studied numerically, this Gaussian is found to contain even more structure), rather than the smooth fermi sphere of the light fermion case. As time progresses, the successive production events at each point wherek = 0 change the particle number with decreasing effect. This is illustrated in the upper panels of figure 8 where we show the average particle occupation number of eq. (5.11) as a function of time. We have so far only considered particle production from the first time wherek = 0. In order to accurately estimate the final number densities for each helicity state, we extend our analysis to successive events. The problem of several particle production events can be mapped to the well-known quantum-mechanical problem of a particle undergoing multiple scatterings. In the static universe the constructive or destructive interference creates a well-defined band structure. However, in the expanding universe case, successive phases can be thought of as random.  The calculation is analogous to that in [13], and we simply quote the result here, referring the reader to the original works for more details. After N successive production events the smoothed particle occupation number is given by This result shows that eq. (5.7) can be used as to characterize particle production at all times, despite being strictly true only for the first particle production event. We test its accuracy in the next section. Figure 8 shows the average particle number as a function of the successivek zerocrossings for a fixed value of the coupling C/f = 50 and varying fermion mass. The first result is the fact that particle production is suppressed for larger values of the mass, and effectively shuts off for m 2 /(C/f ) 0.5, and m 2 /(C/f ) 0.3, equivalently m 5 and m 4, for the positive and negative helicity mode respectively. A similar behavior was seen in figure 5 for the static universe analysis.
The second important feature of figure 8 is the decline of the final particle number for small values of the mass. This presents a major departure from the results of the static universe analysis. In the static universe approximation, the particle number oscillates between two extrema for small masses, as shown in figure 5. From figure 6, note that particle production events always come in pairs of similar strength, especially for wavenumbers near k max . This is shown in figure 6 by the two pairs of brown and green dots for the λ = +1 and λ = −1 states respectively, which represent successive production events. For lighter fermions, particle production is very efficient, and the particles that are created by the first event tend to be destroyed by the second event. In the static universe the production efficiency is constant, but in the expanding universe subsequent events are much less efficient, as seen from eq. (5.21). Furthermore, as show in section 5.2.2, particles are produced over a smaller range of wavenumbers as the universe expands. Taken together, somewhat surprisingly, the net result is a low overall particle number for light fermions. We also demonstrate this result analytically. Applying eq. (5.23) and taking the first two particle production events to be -25 -

JCAP11(2015)021
very efficient, that is e −πp i = 1 − i ≈ 1, we get n (2) meaning that (on average) modes that start almost fully occupied become almost empty after the second production event.
An important qualitative feature of the expanding universe case is the similarity between the particle spectra after two production events, n (2) k , and the result after a very long time, n (∞) k . Successive production events become increasingly less efficient and affect a decreasing range of wavenumbers, therefore these first two production events give the dominant contribution to the final particle spectra. We have already seen in the case of light fermions that these two events lead to the creation of a large number of particles followed by the immediate destruction of many of them. The upper panels of figure 8 show that particle number is effectively constant after the first two oscillations. This is important because it means that any helicity asymmetry will be generated and finalized almost immediately after the end of inflation.
Production events can be always described in pairs, therefore simple intuition tells us that, on average, modes with n k < 0.5 initially tend to grow their population towards unity after the secondk = 0 point, while on average, modes with n k > 0.5 initially tend to reduce their population (away from unity). At fixed coupling C/f , production of lower mass fermions is more efficient, and they are mostly destroyed by the double production event. In the large mass region on the other hand, the efficiency of each event is lower, but successive scatterings are more likely to be constructive. Taken together, this suggests that the first two production events will produce a maximum of n k ≈ 0.5. This is exactly what the blue and red lines in the lower panels of figure 8 show.

Asymmetric helicity production
From our discussion of the range of excited wavenumbers and the evolution of the Fermi sphere, we already see a first hint of helicity asymmetry, simply from the fact that the positive helicity mode is produced up to a larger wavenumber k + max > k − max . With the average particle number in eq. (5.11) and the total particle number density n = d 3 k n k , (5.25) the particle asymmetry between the two helicity states becomes The two terms can therefore be comparable, depending on the structure and values of n ± k . Having given estimates of the range of excited wavenumbers, we can now estimate the particle number and study the cases where it is possible to generate large asymmetries in the different fermion helicities. . Upper panels: average particle number after successive production events for the positive (left) and negative (right) helicity states for C/f = 50 and m fermion = 0.5, 1, 2, 3, 4 (blue, green, red, black and brown respectively). The 0'th production event occurs during inflation. Lower panels: average particle number after the first (blue) and second (red) production event and after multiplek zero-crossings (green).
Returning to figure 8, we note that positive and negative helicity states have comparable average particle occupation numbers. However, when converting this to a total particle number, one has to take into account the fact that each helicity state occupies a different volume in phase space. For m 2 φ φ 2 inflation, for example, k + max = 0.9C/f for the positive helicity state and k − max = 0.7C/f for the negative helicity state, which means that the particle number of the positive helicity state is more than twice that of the negative helicity state for n + = n − . Thus the positive helicity mode is dominant, and a helicity asymmetry is present. We will quantify this asymmetry further using numerical simulations.

Numerical results
To test our analytic and semi-analytic results, we numerically evolve the equations of motion for the fermions in the homogeneous axion and gravitational field. We evaluate the evolution of the fermion number for a range of couplings 2 ≤ M Pl C/f ≤ 500, and fermion masses, such that 0 < m 2 /(C/f ) < 1. We neglect the back-reaction of the produced fermions on the axion evolution; back-reaction can be safely ignored for values of the coupling that satisfy (C/f )M Pl < 10 3 , at least during inflation and for the first stages of preheating, where all interesting effects are expected to occur.
Later axion oscillations (as shown in figure 8) become increasingly irrelevant, and so we focus on the first couple of fermion production events, which determine the final result to a good approximation. For the first two production events eq. (5.23) becomes for the two -27 -JCAP11(2015)021 helicity states which provide an increasingly good approximation for the numerically calculated result for C/f 1, both after the second production event as well as after multiple ones. The important feature of the WKB-based results is their dependence solely on the combination m 2 /(C/f ), rather than on the fermion mass and coupling separately. However, as shown in figure 9, this is not strictly true for the particle number obtained from a numerical evolution of the system. There is a significant departure from the WKB result that occurs for wavenumbers that are near the maximum excited wavenumber. In particular, smaller couplings tend to have a larger "tail" in the n k distribution. This is the same behavior manifest in the static universe case. Because there is much more phase space at high wave number, the volume factor d 3 k weighs these modes more, and this discrepancy widens when the total particle number density or average particle occupation number n k is evaluated. This phenomenon occurs for the positive and negative helicity states alike, leading to an enhanced n k at smaller values of the couplings C/f , for the same ratio m 2 /(C/f ).
For order unity values of the coupling, the departure from the semi-analytic result is even larger -this is perhaps not surprising as the WKB approximation breaks down in this region. After several production events, these smaller couplings exhibit a slow but steady parametric excitation of larger wavenumbers, reaching k max ∼ 2C/f . Despite the occupation number being small n k 1 in this range, the particle number is again enhanced due to the d 3 k phase-space factor. Because we calculate the total particle density by integrating n k up to a sufficiently large wavenumber in order to enclose all particle production, and normalize it by k 3 max (given by the WKB approximation), the average particle number density n ± can exceed unity for late times and small couplings. Furthermore, because of the steady excitation of these short-wavelength modes, the particle number keeps rising as time progresses, and there is no well-defined asymptotic particle number for late times, contrary to the case of larger couplings. A further effect of the behavior of these modes is the restoration of the symmetry between negative and positive helicity states. This leads to ∆n = n + − n − taking both positive and negative values for C/f = 2 in our simulations. However, the range of C/f = O(1) is not of particular interest. This is due to the fact that the total number of particles (and the corresponding helicity asymmetry) is proportional to the cube of the maximum wave number k 3 max ∝ (C/f ) 3 , which implies that for a considerable effect to be achieved, we need C/f 1. This makes model-building more robust, since ∆n/(C/f ) 3 is largely insensitive to the value of the coupling in this region. Figure 10 shows the total particle number density asymmetry in helicity states ∆n. This quantity scales with the cube of the coupling, so we have divided out by a factor of (C/f ) 3 to put all curves on the same axes. The main results of this calculation are the dependence of the normalized helicity asymmetry ∆n/(C/f ) 3 on the number of axion oscillations, as well as on the coupling C/f and the fermion mass m, or equivalently on one of these two parameters and the ratio m 2 /(C/f ). Figure 10 illustrates that, for large values of the coupling C/f 1, the normalized helicity asymmetry is largely insensitive to the value of the coupling and is only determined by the ratio m 2 /(C/f ). Furthermore, there is no significant evolution of the normalized helicity asymmetry after the second production event. We omit the case C/f = 2 JCAP11(2015)021 from the right panel of figure 10, due to the absence of a well defined late-time asymmetry, as explained above.
From a model-building perspective, for a fixed fermion mass m, the resulting asymmetry can be easily read from figure 10 for different values of the coupling. The maximum value of the asymmetry occurs for m 2 /(C/f ) ≈ 0.07 and the full width at half maximum is given by which serves as a condition for the existence of a strong helicity-asymmetry. In order to provide a rough, but intuitive estimate of the maximum asymmetry strength, we start with the observation that the average occupation number for each helicity state, when considered as a function of m 2 /(C/f ), peaks at n ± k max ≈ 0.5. The total particle number is where we used the peak value n ≈ 0.5. By assuming that the two distributions for the n + k and n − k peak at the same value of m 2 /(C/f ) the asymmetry is simply Thus, for the peak value of the asymmetry, we see that the first term in eq. (5.26) vanishes, while the particle number density average in the second term can be taken as approximately 0.5. For chaotic inflation with a quadratic potential, the values of k ± max lead to This formula provides a slight underestimation of the peak helicity asymmetry by about 10%, as seen in figure 10. However, both the accuracy, as well as the simplicity of this formula make it very useful for exploring the fermion production capability of different inflationary models, like the axion monodromy potential which will be discussed in the next section.

Axion monodromy inflation
In this section, we extend the results derived above in section 5 for quadratic chaotic inflation to the simplest model of axion monodromy inflation with a potential given by eq. (2.3). All the machinery we have developed in section 5 can be transferred directly to this case. The only differences are in the specific numerical values of |dk/dt| at the various production events and the range of excited wavenumbers. We begin by rescaling our parameters to rewrite the background inflaton equation of motion in a dimensionless form. To generate a spectrum of curvature fluctuations with an amplitude that matches the observed microwave background, we choose log 10 (µ/M Pl ) = −3.22 [38]. Measuring the inflaton field in units of the Planck mass, and time in units of 1/ M Pl /µ 3 , the background equation of motion for the inflaton becomes In what follows, wavenumbers will be measured in terms of µ 3 /M Pl .  Most of the qualitative features of the m 2 φ φ 2 case can be directly applied to monodromy inflation. For monodomy inflation, gravitational particle production during inflation is biased by the rolling axion so that only a single helicity is produced. After inflation both helicity states are excited at the points wherek(t) = 0, which differ for each helicity state and wavenumber. The range of excited wavenumbers is smaller for each subsequent production event, due to the decay of the inflaton oscillations, making the first two production events the dominant ones. For axion monodromy inflation, there is an additional parameter, φ c , which controls the shape of the potential near the end of inflation and quantitatively changes the results.
For φ c = 0 the first few oscillations dominate and define the range of the Fermi sphere, as we found in the m 2 φ φ 2 case. The corresponding maximum wavenumbers excited during the first production event for each of the two helicity states are shown in figure 11, along with JCAP11(2015)021 which tracks the difference in phase space volume that is populated in each case. Note that for m 2 φ 2 inflation (∆k max ) 3 ≈ 0.4, which is similar to the values for axion monodromy inflation, in the range of 0.1 < φ c < 1.
For large values of φ c ≥, the minimum of the potential behaves more and more like a quadratic model during the oscillation phase, and we find similar quantitative results to m 2 φ φ 2 inflation. Specifically for φ c = 1 the form of a(t)φ(t) is very similar for quadratic and axion monodromy inflation, with a difference in amplitude of about 25%. This leads to a difference in the number density of produced fermions of about 1.25 3 ≈ 2. The maximum wavenumber excited in each axion oscillation, k ± max , grows with decreasing φ c , as shown in figure 11. However, the combination (k + max ) 3 − (k − max ) 3 decreases for φ c 0.3. For small values of φ c , both helicities are produced in increasingly equal amounts, resulting in smaller helicity asymmetries. The reason for this behavior can be seen from the (unphysical) limit φ c → 0. In this limit,φ(t)a(t) oscillates with a constant amplitude generating a Fermi sphere with a constant radius in momentum space. Each axion oscillation therefore produces particles out to the same maximum wavenumber with an efficiency that remains essentially constant. Hence for φ c < 0.1 we have an increased particle production (due to the increased range of produced wavenumbers) and a decreased level of helicity asymmetry.
The average particle number after the second production event for each helicity state is shown in figure 12. This is a very good approximation for the final particle number, especially near the peak of the curve n ± as a function of m 2 /(C/f ), where both the maximum particle production and the maximum asymmetry occur. We present results for the range 0.2 < φ c < 0.6, where the asymmetry is largest, as expected from figure 11, but simulations in the whole range φ c ∈ [0.1, 1] show similar results. As expected, the maximum average occupation number is n ± ≈ 0.5. Making use of eq. (5.30) and figure 11, for axion monodromy inflation, the total helicity asymmetry in the number density of left-versus right-helicity fermions is given approximately by where the subscript refers to the maximum asymmetry. The value of m 2 /(C/f ) for which the maximum asymmetry occurs is a model-dependent quantity, but for the model of axion

Conclusion
In this work we have studied the behavior of derivatively-coupled, massive fermionic degrees of freedom during, and immediately following axion-driven inflation. When the masses of the fermions are degenerate, these fermions pair up to make Dirac fermions and the coupling is to the axial-vector current. During inflation, the motion of the inflaton is monotonic and the coupling of the fermions to the rolling axion acts as an effective chemical-potential for helicity. This chemical potential biases the gravitational production of helicity states of the fermions. In the absence of the coupling, massive fermions are produced gravitationally by the expansion of spacetime equally in all helicity states. This production is most efficient for light fermions; the production rate is exponentially suppressed by the ratio of the fermion mass to the Hubble rate. When the axion coupling is turned on, we have shown that one helicity state is produced more efficiently while the other is strongly suppressed. This mechanism allows for the gravitational production of heavy fermion states that would otherwise be heavily suppressed.
Following inflation, as the axion oscillates, the effective frequency of the fermion field varies non-adiabatically resulting in particle production. These oscillations also mean the axion velocity changes sign, resulting in the production of both fermion helicities. The damping of the inflaton amplitude after the end of inflation, due to Hubble friction, results in biased fermion production in favor of the helicity state that is excited during inflation. This is due to the declining efficiency of the production rate as the axion damps combined with the fact that production of the second helicity does not begin until the axion velocity changes sign. For certain combination of the fermion mass and axion-fermion coupling, the -32 -

JCAP11(2015)021
produced fermions can have a significant helicity asymmetry. This helicity asymmetry is robust to the details of the reheating history, at least in the limited models we considered here, because it requires the amplitude of the axion oscillations to decay faster than 1/a(t), which does not need any significant tuning of the model.
We studied the post-inflationary evolution of the particle number in both quadratic chaotic inflation -m 2 φ φ 2 inflation -and in the simplest model of axion monodromy inflation. While there are minor numerical differences between the models, broadly we find very similar behavior. In both cases the average particle number is solely defined by the combination m 2 /(C/f ) in the regime of large coupling, while there is increasing deviation from this behavior as we decrease the coupling. Further, in both models the range of excited wavenumber shrinks with time as well as the production efficiency of any single event, leading to the first production events giving the dominant contribution to the final particle spectrum. The average particle number remains almost constant after the second production event for C/f 1. Finally, the maximum achievable helicity asymmetry in the number density of produced fermions, ∆n, scales as (C/f ) 3 , with a proportionality factor of O(1), which depends on the specific inflationary model.
Neutrino helicity is equal to lepton number in the standard model of particle physics, therefore the asymmetric production of helicity may be important for the generation of the matter-antimatter asymmetry in the Universe. We apply the results of this work in a companion paper [10], where we propose that an axion-inflaton coupled to left-handed standard-model neutrinos as a mechanism for the observed baryon asymmetry in the Universe via leptogenesis.

JCAP11(2015)021
and we formally define αβ = ( αβ ) * and αβ = ( αβ ) * . Note that these conventions imply For convenience, we work only with left-handed fields, that is, fields that transform under the ( 1 2 , 0) representation of the Lorentz group. Since these representations are related by hermitian conjugation, right-handed spinors are then simply conjugates of left-handed spinors.