Nonthermal Gravitino Production after Large Field Inflation

We revisit the nonthermal gravitino production at the (p)reheating stage after inflation. Particular attention is paid to large field inflation models with a $\mathbb{Z}_2$ symmetry, for which the previous perturbative analysis is inapplicable; and inflation models with a stabilizer superfield, which have not been studied non-perturbatively. It is found that in single-superfield inflation models (without the stabilizer field), nonthermal production of the transverse gravitino can be cosmologically problematic while the abundance of the longitudinal gravitino is small enough. In multi-superfield inflation models (with the stabilizer field), production of the transverse and longitudinal gravitinos is significantly suppressed, and they are cosmologically harmless. We also clarify the relation between the background field method used in the preheating context and the standard perturbative decay method to estimate the gravitino abundance.

1 Introduction and summary 1

.1 Introduction
Supersymmetric (SUSY) models are well-motivated as a physics beyond the standard model, since it provides a successful gauge coupling unification, dark matter candidate, a great reduction of the hierarchy problem etc. In supergravity, however, there is a cosmological problem associated with the gravitino, the superpartner of the graviton, called the gravitino problem [1][2][3][4]. If the gravitino is not the lightest SUSY particle (LSP), it can decay into lighter SUSY particles. The lifetime of the gravitino is given by [5] τ 3/2 = 3 8π where M Pl is the reduced Planck scale and m 0 3/2 denotes the present gravitino mass. 1 Here we have assumed that the gravitino decays into only gaugino plus gauge boson pairs. If other decay modes are open, the lifetime becomes slightly shorter. Therefore, if the gravitino is much lighter than 100 TeV, it decays after the beginning of big-bang nucleosynthesis (BBN) and may affect light element abundances [5][6][7][8][9][10]. If it is heavier, the decay itself does not affect BBN but LSPs produced by the gravitino decay can be overabundant compared with the observed dark matter abundance. If the gravitino is LSP and R-parity is conserved, the gravitino itself contributes to the dark matter abundance. In any case, there is a strict upper bound on the gravitino abundance.
There are several processes that produce gravitinos in the early universe. One of the unavoidable production mechanisms is thermal production: in the high-temperature universe, scatterings of high-energy particles produce gravitinos [11][12][13]. 2 The abundance of thermally produced gravitinos is proportional to the reheating temperature after inflation T R , and hence we obtain an upper bound on T R to avoid the gravitino problem.
Gravitinos can also be produced nonthermally. Nonthermal gravitino production by the direct decay of inflaton was extensively studied in a series of works [16][17][18][19][20][21][22][23][24]. It was found that the inflaton generally decays into the gravitino pair with the partial decay rate [18,19,21,22] (1.2) where m φ is the inflaton mass and φ is the vacuum expectation value (VEV) of the inflaton. It gives a stringent constraint on inflation models, although there are some loopholes [25,26]. This expression for the decay rate is valid for small-field inflation models such as new inflation or hybrid inflation. 3 On the other hand, large field inflation models [27] attract lots of attentions in view of recent developments on successful inflation model building in the framework of supergravity [28][29][30][31][32][33][34][35][36]. It is interesting because it can be tested with on-going/future B -mode polarization experiments. In large field inflation models with a 2 symmetry in which the inflaton field oscillates around 1 Throughout this paper, we distinguish the "present gravitino mass" m 0 3/2 and "gravitino mass" m 3/2 , since the notion of gravitino and its mass is time-dependent in a cosmological evolution. The former corresponds to the gravitino mass in the present universe and it is just a constant while the latter is time-dependent. 2 Ref. [14] discussed the gravitino production by the scatterings of energetic inflaton decay products during the process of thermalization [15] and found that it is subdominant compared with the standard thermal production. 3 Chaotic inflation without a 2 symmetry also leads to a similar expression.
the origin φ = 0 after inflation, we cannot use the expression (1.2) as a gravitino production rate. This is simply because the calculations in Refs. [21][22][23][24] assumed the perturbative decay of inflaton around its VEV. In inflation models with the 2 symmetry, however, there is no such decay process due to the 2 symmetry. 4 This does not mean that the inflaton cannot decay into gravitinos as well as other light particles. The inflaton coherent oscillation affects the masses or kinetic terms of coupled particles. The coupled particles, including gravitinos, "feel" the rapid inflaton oscillation and it affects the evolution of their wave functions. It is known that this leads to particle production, often in the context of preheating [37][38][39][40][41]. Therefore, even if the inflaton has the 2 symmetry, its coherent oscillation necessarily transfers its energy to the coupled particles. The question we would like to address is: what amount of gravitinos is produced during the preheating?
Production of gravitinos during the preheating was first discussed in Refs. [42,43] in a singlesuperfield case, in which only one inflaton chiral superfield was introduced. There it was found that in the preheating stage, longitudinal gravitinos are efficiently produced. Later it was recognized that the theory of gravitino preheating is much more involved due to the subtlety of the notion of "gravitino" [44][45][46]. The gravitino becomes massive by "absorbing" the goldstino, but the definition of goldstino is time-dependent in a cosmological background. In the early universe, the inflaton oscillation energy dominantly breaks SUSY and hence the goldstino is almost the inflatino, the fermionic superpartner of the inflaton. At late time, the Polonyi field 5 dominantly breaks SUSY and its fermionic component, Polonyino, becomes the goldstino. Thus the composition of goldstino changes with time. Refs. [45,46] noticed that it is essential to include (at least) two chiral superfields, inflaton and Polonyi, to correctly deal with this problem and concluded that what the preheating efficiently produces eventually becomes the inflatino, which is less harmful than the gravitino.
Still, however, a quantitative/comprehensive analysis of the nonthermal gravitino production rate in such a case is missing. Although Refs. [45,46] revealed that the gravitino production is suppressed than previously thought, it is highly non-trivial how we can extrapolate their numerical results into more realistic setups and parameters both in the inflaton and SUSY breaking sector. Thus we would like to provide general analytic formulae for the nonthermal gravitino abundance that are applicable to any realistic models.

Brief summary
In this paper, we revisit the theory of nonthermal gravitino production in a comprehensive and unified manner. Our purposes and results are summarized below.
• We derive nonthermal gravitino production rates and their resulting abundances quantitatively with useful formulae. We find that in single-superfield inflation models, the production of transverse gravitino is significant and cosmologically problematic, while the production of longitudinal gravitino is less important. This aspect of the nonthermal gravitino production has been overlooked in previous literatures except for a few [47]. 4 Another assumption there was that the SUSY is dominantly broken by the Polonyi field at the end of reheating so that the definition of "gravitino" at that epoch is the same as the present-day gravitino. This assumption is valid as long as we are interested in the gravitino production at the end of reheating. In large field inflation models, however, the gravitino production is often dominated at the epoch just after inflation (preheating) and hence this assumption is not justified, as we will see later. 5 In this paper we call the present-day SUSY breaking field as Polonyi field.
• Recent realistic large field inflation models introduce an additional chiral superfield, called a stabilizer [28][29][30][31][32][33][34][35][36]. We find that in models with the stabilizer field, the production of transverse gravitino is significantly suppressed and it is cosmologically harmless. For the longitudinal component, the production rate is similar to the single-superfield case.
• We show the equivalence between the background field method developed in Refs. [44][45][46] and the perturbative decay method developed in Refs. [21][22][23][24] for evaluating the gravitino abundance in some sense. The former can deal with a broad class of models including 2symmetric large field models, while in models without the 2 symmetry, it gives the same result as the perturbative decay method.
This paper is organized as follows. In Sec. 2, we review the structure of the gravitino Lagrangian to set the stage of discussion in the subsequent Sections. We formulate a general setup to discuss multi-superfield case. Gravitino production in the single-superfield inflation models and in the multi-superfield inflation models are studied in Sec. 3 and Sec. 4, respectively. The analyses include both the transverse and longitudinal components of gravitino. Our conclusion is in Sec. 5, and the gravitino abundance is summarized in Fig. 3. App. A summarizes our notations and conventions. The background field method to evaluate the fermion production rate, often used in the preheating context, is reviewed in App. B. We also briefly cover the gravitino production in small-field inflation models in App. C in order to show the equivalence between our method and the perturbative decay method. In App. D, we review the multi-field scalar dynamics to discuss the induced oscillation of the Polonyi field in the main text. Calculations of mass eigenvalues are given in App. E, which are used in Secs. 3.3 and 4.3.

Master supergravity Lagrangian
We start from the following supergravity Lagrangian: where e is the determinant of the vierbein e a µ , M Pl is the reduced Planck mass, R is the Ricci scalar, φ i and φ * ī are scalar fields and their complex conjugates, χ i L and χj R are left-handed matter fermions and their conjugate right-handed fermions, and g ij = ∂ i ∂j K is the Kähler metric with ∂ i being the derivative with respect to φ i . The scalar potential V is given in terms of the Kähler potential K and the superpotential W as where D i W = ∂ i W + (∂ i K ) W /M 2 Pl , and g ij is the inverse Kähler metric. The hats denote the quantities in the curved space-time. The field strength of the gravitino ψ µ in its kinetic term is given by is the "remnant" of the gauge field of the R-symmetry in the underlying superconformal formulation, and γ a is the Dirac gamma matrix in the flat space. Note that practically the Christoffel symbol does not contribute due to the anti-symmetry of γ µρσ . The gravitino mass m 3/2 and the other fermion mass matrix m i j are respectively defined as where Γ k i j ≡ g kl ∂ i g jl is the Christoffel symbol in the Kähler manifold. The goldstino v is defined as The last term in Eq. (2.1), 4f , denotes the four-fermion interactions originating from the torsion, which we will neglect from now. We consider only gauge singlet components in this paper. For more details on the notation and conventions used in this paper, see App. A. The master Lagrangian (2.1) contains the gauge redundancy as well as unphysical degrees of freedom. From now, we fix the fermionic gauge redundancy by taking the unitary gauge v = 0, (2.9) and integrate out the unphysical degrees of freedom to obtain the physical Lagrangian. In order to do so, we should first find constraint equations, and solve them to express the unphysical degrees of freedom in terms of the physical degrees of freedom, i.e., the transverse and longitudinal modes of the gravitino. Below we give the outline of this procedure. For more details, see Ref. [44]. In the following, we assume that the scalar fields are real for simplicity, and hence m * 3/2 = m 3/2 and A µ = 0.

Constraint equations
In the unitary gauge, the equations of motion for the gravitino are given by From these equations, we can verify that the following equations do not contain the time derivatives with respect to the gravitino, and hence are constraints: The first constraint (2.11) can be solved to give ψ 0 in terms of ψ, χ i L and χj R . Practically, however, we do not need to know the explicit solution for ψ 0 . This is because the Lagrangian (2.1) without 4f depends only linearly on ψ 0 , and hence ψ 0 contributes to the Lagrangian with a combination of ψ 0 Σ 0 . Thus, it automatically drops from the Lagrangian once we impose the second constraint (2.12). For this reason, we concentrate only on the second constraint (2.12), which relates k · ψ and γ · ψ.
So far we took the background metric to be generic. From now on, we take it to be the Friedmann-Lemaître-Robertson-Walker (FLRW) one, with a being the scale factor, since we are interested in the gravitino production in the cosmological background. We also decompose the gravitino as (2.14) where the longitudinal mode is ψ ≡ γ · ψ and the transverse mode satisfies γ · ψ t = k · ψ t = 0, respectively. Here we have moved to the momentum space. Then, we can solve the second constraint (2.12) to obtain k · ψ in terms of ψ as where H ≡ȧ /a is the Hubble parameter. 6 This relation ensures that ψ is actually "longitudinal". By substituting it to the original Lagrangian and sorting things out, we obtain the Lagrangian for the physical gravitino ψ t and ψ . Its explicit form will be shown in the next subsection.

Lagrangian of physical gravitino
The gravitino Lagrangian in terms of the physical degrees of freedom is given by where the first two terms correspond to the kinetic and mass terms of the transverse and longitudinal gravitino, respectively: 7 (2.17) 6 In this paper, the dot denotes the derivative with respect to time t while ∂ 0 represents the derivative with respect to conformal time η. 7 Here we keep terms which are exactly zero due to the Majorana property (e.g.,ψγ 0 ψ = 0), in order for the combinations which appear in the equations of motion to be manifest. and The last term represents the mixing between the longitudinal gravitino and chiral fermions: Note that the transverse mode does not mix with chiral fermions. Here we have defined where the SUSY breaking F -term is defined as F i = −e K /2M 2 Pl g ij Dj W * , for the minimal Kähler potential. 8 The A is an important combination whose phase rotation gives rise to gravitino production as we will see below. The V SB , ρ SB and p SB can be interpreted as the potential energy, energy density and pressure of the matter components which break SUSY while p W measures the time variation of the gravitino mass, which can also be expressed like the "geometric average" of the SUSY breaking kinetic and potential energies. Here and in what follows, we assume the reality of dynamical scalar fields for simplicity (and accordingly, the difference between the subscript and superscript of the field index disappears). Actually, if the parameters in the Kähler and superpotentials are all real and the initial condition of the fields is taken to be real, the subsequent dynamics does not affect the reality of scalar fields. The Friedmann equation reads Here ρ is the energy density of the system. Note that ρ SB is the energy density associated with SUSY breaking, which differs from ρ. It also means that Now let us rewrite these Lagrangians in terms of the canonical field for later convenience. The Lagrangian of the transverse mode in terms of the canonical field ψ t c ≡ a ψ t is given by 9 Note that the mass of the order of the Hubble scale disappears reflecting the conformal invariance. Thus it does not "feel" the Hubble expansion, meaning that the oscillation/non-adiabatic change 8 In this paper the "minimal" Kähler potential means K ij δ ij . Adding a holomorphic function in K does not change the results. Note that the minimal shift-symmetric Kähler potential (K = −(φ − φ † ) 2 /2) also satisfies K ij = δ ij . Also we will introduce a higher order Kähler potential like K ∼ |z | 4 /Λ 2 for a Polonyi field z , but it does not affect the discussion as long as the field value of z is small enough. 9 Note that / Dψ = ( / ∂ + 3 2 a H γ 0 )ψ andψγ 0 ψ = 0 for a Majorana fermion ψ.
of the Hubble parameter does not lead to the production of transverse gravitino. On the other hand, the gravitino mass m 3/2 oscillates rapidly in the inflaton oscillation epoch, which inevitably leads to significant gravitino production. As for the longitudinal mode, let us rewrite the Lagrangian by using the canonical field In the second equality, chiral fermions in the original supergravity Lagrangian χ i are rescaled as χ i ≡ a 3/2 χ i and the prime is dropped in the following discussion so that their kinetic terms without explicit dependence on the scale factor a . In Eq. (2.27), we have defined the generalized gravitino mass term, which is of the order of m 3/2 ∼ (H ) + (m 3/2 ). In the Lagrangian, the coefficient A as well as the gravitino mass m 3/2 oscillates rapidly. The whole structure is slightly complicated, especially in the multi-superfield case. However, in almost all the situations of our interest, the SUSY breaking is dominated by one field and the problem effectively reduces to the single-superfield case, in which the analysis is significantly simplified.

The case of only one chiral superfield
In order to illustrate the structure of the theory, first let us focus on the single-superfield case. In this case, we have a relation Therefore we can take − A ≡ e 2γ 0 θ with θ being a real parameter. Then by defining ψ c ≡ e −γ 0 θ ψ c (and dropping the prime thereafter) the Lagrangian is simplified as The remaining task is to calculate ∂ 0 θ = aθ . 10 Substituting the explicit expression ofθ , we obtain Below we estimate behavior of the mass term in typical one-field dominated cases.

Inflaton-dominated SUSY breaking
Let us consider the case where the inflaton φ dominates the SUSY breaking. As we will see later, if the energy density of the universe is dominated by the inflaton oscillation, and H m 0 3/2 , the SUSY is dominantly broken by the inflaton superfield. The Kähler potential and superpotential are assumed to be where m 3/2 m φ φ 2 /(2M 2 Pl ). Note that we have assumed the reality of φ-dynamics to derive those expressions. Substituting these equations into , we obtain 12 It obtains a mass of the order of the inflatino mass, as expected. 13 Note that in the inflatondominated SUSY breaking case, we roughly have ρ SB ∼ |p SB | ∼ |p W |. Especially p SB and p W are violently oscillating functions with frequency of twice the inflaton oscillation, and m 3/2 itself is also an oscillating function. Thus, both terms proportional to H and m 3/2 contribute to the production of longitudinal gravitino. For the superpotential with a higher power, with n > 2, we obtain (2.42) 11 Useful equations are:ρ = −6H |φ| 2 ,ρ SB =ρ + 3m 3/2 p W ,ṗ SB =ρ − 2V − 3m 3/2 p W . 12 Note thatV = V φφ + V φ * φ * = 2V φφ for real φ. Here V φ should be regarded as ∂ φ V (φ, φ * ). 13 This expression is consistent with Eq. (9.15) of [44].
up to corrections with (φ 2 /M 2 Pl ). The Lagrangian of the longitudinal gravitino becomes (2.43) where m 3/2 λφ n /(n M 2 Pl ). Note that n = 2 reproduces the above result (2.39). In this case with general n , the longitudinal gravitino obtains a mass of ∼ λφ n−2 , which is rapidly oscillating. Other terms proportional to H and m 3/2 are suppressed by φ/M Pl and φ 2 /M 2 Pl , respectively. Therefore, the mass term of ∼ λφ n−2 contributes to the gravitino production dominantly and it is much more efficient than the case of n = 2.
When we start from the minimal shift-symmetric Kähler potential, K ⊃ −(φ−φ † ) 2 /2, the qualitative discussion is maintained. The only differences in the final expressions of the Lagrangian are that 5 and n(n + 3)/2 in the last parenthesis in Eqs. (2.39) and (2.43) respectively are both replaced with −1.

Polonyi-dominated SUSY breaking
Now let us consider the case that the Polonyi field z dominates the SUSY breaking as in the present universe. In a cosmological setup, this approximation is valid at H m 0 3/2 . The Kähler potential and superpotential are assumed to be In the present paper, we always implicitly assume a dynamical SUSY breaking scenario [48][49][50][51][52] in which the Polonyi field z obtains a large SUSY breaking mass, represented by the nonminimal Kähler potential. 14 In this model we obtain where the mass of the Polonyi field is given by m 2 z = 12(m 0 Substituting these equations into , we obtain a simple expression: Here we keep relevant terms in the limit δz → 0, in which we have p SB /ρ SB −1, |p W | ρ SB , etc. One can see that ψ obtains a mass of m 0 3/2 as expected. The last term is responsible for the Polonyi decay into the longitudinal gravitino pair and it is clear that there is an enhancement of ∼ (3m z /m 0 3/2 ) 2 for the Polonyi decay rate into the longitudinal mode compared with that into the transverse mode. We have omitted the Hubble mass term since it is suppressed by the ratio p W /ρ SB . This is because, in the limit δz → 0, the F term of z does not contribute to the Hubble expansion owing to the requirement that the cosmological constant in the present vacuum is (almost) zero.

The case of several chiral superfields 2.4.1 General argument
Let us consider the multi-field case φ i (i = 1, . . . N ). We want to express all N chiral fermions χ i in terms of the goldstino v and those orthogonal to it. In this setup, the goldstino is given by 15 and These are left-handed spinors as indicated by the subscript L. The righthanded counterparts are similarly given by where we have assumed scalar fields are real. It is convenient to define the Majorana spinor where O is an N × N orthogonal matrix and O T is its transpose, whose first row is determined by Eq. (2.53) as from the orthogonality and normalization conditions of v I ⊥ . Also, this matrix fulfills 0 = α i O i J . Using these elements defined above, we can express all terms in the Lagrangian including fermions with v I ⊥ in the unitary gauge v = 0. The fermion mass term becomes The fermion kinetic term becomes The gradient term seems curious. However, by taking into account the mixing of ψ c and v I ⊥ , it can be diagonalized. In particular, the off-diagonal gradient term of ψ c and v ⊥ comes from mix . The longitudinal gravitino-fermion mixing term (2.28) is rewritten in a simple form as The quantity A (2.20) is also expressed as It is easily checked that | A| 2 = 1 in the single-field case. Combined with the gravitino kinetic term (2.27) and mix , we obtain a matrix form of the gradient term: Here and hereafter we suppress the indices of v ⊥ for brevity. This can be written as It is easily seen that | | 2 = 1 and hence can be regarded as a generalization of A in the singlefield case. This leads us to define an N × N matrix θ such that ≡ e −2γ 0 θ . Then we have Therefore we can diagonalize the gradient term by using the rescaled field Using this basis, the kinetic term of ψ c and v ⊥ is completely diagonal. On the other hand, this transformation yields the mass mixing between them. After all, we finally obtain the Lagrangian of the gravitino-fermion system written only by the physical degrees of freedom as For later convenience, we can also estimateθ i aṡ It is also conveniently expressed aṡ where we have decomposed This is a generalization of the single-field case (2.33). Note that there are off-diagonal antisymmetric mass terms proportional to γ 0 . This part can be removed with a transformation by an orthogonal matrix as described in Sec. 3.2 of Ref. [46], and hence does not contribute to the mass eigenvalues.

Two-field example
Now let us consider the two-field case, which is used in Sec. 3. In this case, after removing the goldstino v in the unitary gauge, there remains only one fermion v ⊥ other than the gravitino. The orthogonal matrix O and O are easily found to be The goldstino and the fermionic degree of freedom orthogonal to the goldstino v are, which is inversely solved as In the last equality, we have taken the unitary gauge v = 0. First, the fermion mass term becomes We can substitute this expression into the fermion kinetic term: The matrix of the gradient term is given by where is a 2 × 2 matrix defined by As noted earlier, the matrix may be regarded as a generalization of A. We can express as = e −2γ 0 θ with θ being a real symmetric matrix. Therefore we can diagonalize the gradient term by redefining the fields as in Eq. (2.64). Explicitly, and After the diagonalization of the gradient term, the full fermionic Lagrangian becomes where denotes the mass matrix given by Clearly, in the single-field dominance limit α 1 α 2 or α 2 α 1 , the off-diagonal elements are suppressed by the ratio α 2 /α 1 or α 1 /α 2 . Therefore the mixing between ψ c and v ⊥ is also suppressed by this ratio, and in this limit we effectively recover the single field case studied in Sec. 2.3.
An explicit expression for the mass matrix = 1 + 2 is given by As mentioned above, the terms proportional to γ 0 do not affect the mass eigenvalues.

Gravitino production in single-superfield inflation
Let us consider the two-field case, in which there are two chiral superfields: inflaton φ and Polonyi z . Their fermionic components are denoted by φ and z , respectively. The Kähler and superpotential are assumed to be Although we have imposed an approximate shift symmetry φ → φ +c with c being a real constant, almost all the following arguments do not depend on this specific choice of the Kähler potential as long as K φφ 1. Without loss of generality, we can take all the coupling constants real and positive. In the present vacuum, φ = 0 and 〈z 〉 2 3M Pl (m 0 Pl is the gravitino mass in the present universe. We focus on the case with the 2 symmetry in which there is no linear term in the Kähler potential and the inflaton oscillates around φ = 0. The case without the 2 symmetry will be discussed in Sec. 3.5.
One should note that this theory does not lead to successful chaotic inflation: the simple power-law behavior of the inflaton potential V m 2 φ |φ| 2 is ensured only at the sub-Planckian field range |φ| M Pl . Therefore we need to carefully choose the Kähler and/or superpotential to modify the potential at large field value for successful inflation [56][57][58][59][60][61][62][63][64]. However, we are interested in the behavior after inflation with sub-Planckian field value, and hence such modifications on the potential at large field value do not significantly affect the following discussion.

Dynamics
First let us briefly summarize the scalar dynamics in the present model. The inflaton dynamics is so simple that the inflaton φ just oscillates around the origin φ = 0 along the axis of Re φ after inflation. The oscillation amplitude φ amp decreases as φ amp ∝ a −3/2 ∼ t −1 until it completely decays and the universe is reheated.
The z dynamics is slightly complicated. The oscillation of z is directly induced by the inflaton dynamics through the potential term Although this term is not regarded as a "mixing" between φ and z formally, it inevitably induces a coherent oscillation in the z direction when φ has a finite oscillation amplitude. In that sense, this may be viewed as an effective mixing term between φ and z , with the mixing angle of 16 As shown in App. D, the induced amplitude of z can be estimated as z amp ∼ θ φz φ amp . This nonzero amplitude of z may contribute to the longitudinal gravitino production as will be shown later. Note that this "induced" oscillation of z always exists even for m z H inf . It can be interpreted as a result of tilted axis of the oscillation on (φ, z ) plane due to the effective mixing term (3.3), hence it is just a small mixture of the z into dominantly φ oscillation.
On the other hand, if m z H inf , there is another oscillation mode which dominantly consists of the z direction, which occurs at H ∼ m z . This is roughly the oscillation along the light mass eigenstate, which mostly consists of z . If m z H inf , the adiabatic suppression of the coherent oscillation works [65,66] and the non-induced oscillation is safely neglected. Therefore, taking account of the Hubble expansion, the oscillation amplitude is given by [26] where a (t osc ) denotes the scale factor at H = m z . The resultant oscillation amplitude is the sum of the induced one and the non-induced one: The effect of the non-induced one (3.5) on the gravitino production was considered in Refs. [26,67] and we mainly focus on the effect of the induced one.

Transverse component
First let us evaluate the transverse gravitino production rate. The Lagrangian is given by Eq. (2.25), and the rapid oscillation of the mass m 3/2 contributes to the gravitino production. The oscillating part of the gravitino mass is given by m 3/2 m φ φ 2 /2M 2 Pl . A schematic picture of the time evolution of m 3/2 , H and m φ is shown in the left panel of Fig. 1. Hence, as shown in App. B, we obtain the effective "annihilation" rate of the inflaton as where is a numerical constant of order unity. We take = 1 as a reference [See Eq. (B.29)]. Since m φ m 3/2 always holds after inflation, this process is kinematically accessible. 17 This production rate is comparable to the gravitational particle production rate of a minimal scalar field [68,69], although the gravitino is actually conformal in the massless limit. Note that since the transverse mode does not have mixing with any other fermion, the transverse mode produced in any epoch eventually becomes the present gravitino.
Since the amplitude φ amp is rapidly decreasing with respect to the cosmic expansion, the transverse gravitino production becomes less and less effective as time goes on. This is shown by checking that the gravitino number density produced per Hubble time, , (3.8) decreases faster than a −3 , or Γ(φφ → ψ t ψ t ) decreases faster than H . Thus the dominant contribution comes from those created within a few inflaton oscillation just after inflation. The present gravitino abundance is estimated as where T R is the reheating temperature and H inf denotes the Hubble scale at the end of inflation. This is about three orders of magnitude smaller than the contribution from thermal production.

Longitudinal component
Next let us consider the production of longitudinal gravitino. As already mentioned in Sec. 2, we must be careful on the mixing with fermions to correctly treat the production and evolution of the longitudinal gravitino.
In the present case, there is one physical fermionic degree of freedom v ⊥ in the matter sector, and it mixes with the longitudinal gravitino. The goldstino is given by where Note again that we have assumed that scalar fields have real values: this is justified since we have taken all model parameters and also the initial condition during inflation real. In the unitary gauge, v is set to be zero. The remaining fermionic degree of freedom is that orthogonal to the goldstino, As a rough estimation, we have 14) where denotes the mass matrix defined in Eq. (2.83). Identifying 1 → φ and 2 → z in the notation of Sec. 2.4, is roughly estimated as for H m 0 3/2 and for H m 0 3/2 , where the first 's involve a similarity transformation of the matrix. The analysis in App. E shows that the determinant of is However, the former contribution is negligible for the following reason. From the structure of the mass matrix, the mixing angle between ψ c and v ⊥ is about and we always findθ ψ c −v ⊥ ∼ H θ ψ c −v ⊥ m φ , meaning that the time evolution of the mixing is always adiabatic. Therefore we can neglect the conversion of the heavy mass eigenstate produced at the early epoch (dominantly ψ c ) into the light mass eigenstate at the late epoch (again dominantly ψ c ). 18 Below we analyze the production of the light mass eigenstate and compute the final longitudinal gravitino abundance. Also we will comment on the production of the heavy mass eigenstate. where A very important note here is that m 3/2 appearing in m light does not contain an oscillating part of the order of p SB /ρ SB ∼ (1) despite that m 3/2 contains such a violently oscillating part. Actually a careful calculation of the mass eigenvalues of the full mass matrix in App. E shows that such a violently oscillating term cancels out in the light mass eigenvalue m light (but not for m heavy ). Therefore, the φ dependence of m light is at most m φ φ 2 /M 2 Pl which may be included in m 3/2 , leading to the upper bound on the production rate as It is suppressed by a factor (m 0 3/2 /H ) 4 compared with the transverse gravitino production rate. The origin of this suppression is the fact that the Polonyino z is massless in the global SUSY limit. Although it is absorbed by gravitino to make it massive by the super-Higgs mechanism in the late epoch, goldstino is mainly composed of inflatino in the early epoch and the Polonyino remains light.
Note also that z oscillation also contributes to the production. As long as we are concerned with the induced oscillation explained in Sec. 3.1, its contribution to the light mass eigenstate (∼ v ⊥ ) production rate is at most the same order of the upper bound from φ annihilation mentioned above. 19 Late epoch : The longitudinal gravitino ψ produced at H < m 0 3/2 is essentially the same as the present gravitino. Thus, we can estimate the final longitudinal gravitino abundance by simply evaluating the production rate of ψ at H < m 0 3/2 . Taking the limit α 2 α 1 in the general two-field Lagrangian (2.83), we find the Lagrangian of ψ same as that discussed in Sec. 2.3.2: where Again, it is important to note that m 3/2 in the expression of m light does not contain a violently oscillating part and the φ dependence of m light is at most m φ φ 2 /M 2 Pl . 20 Thus we obtain an upper bound on the production rate of ψ same as that of the transverse gravitino as This is a rapidly time decreasing function and becomes less efficient at the later epoch. We again note that the contribution of the induced z oscillation to the longitudinal gravitino production is at most the same order of the upper bound from φ annihilation mentioned above.

Abundance of longitudinal gravitino
Combining the results of production rate at the early and late epochs, we find that the longitudinal gravitino production is most efficient around H ∼ m 0 3/2 . Thus the abundance of longitudinal gravitino is found to be It is much smaller than the corresponding transverse gravitino abundance. The dotted line in the left panel of Fig. 3 in Sec. 5 summarizes the transverse and longitudinal gravitino abundance.

Abundance of inflatino
Finally, we comment on the abundance of inflatino. One should also be cautious about the notion of inflatino: it is φ if the SUSY breaking is dominated by z , i.e., H m 0 3/2 while it is absorbed by the longitudinal gravitino at H m 0 3/2 . Thus, the heavy mass eigenstate of (ψ c , v ⊥ ) system eventually becomes the inflatino no matter when it is produced.
To estimate the inflatino abundance, we note that the heavy mass eigenvalue is given by where the second term is the oscillating part, see App. E. In contrast to m light , there appears a term of ( m 3/2 ) in m heavy , which can make the production efficient. Just after inflation, | m 3/2 | is of the same order of m φ if φ i ∼ M Pl , hence the production is kinematically accessible around H ∼ H inf and the production becomes less and less efficient after that due both to the kinematical suppression and the decrease of the oscillation amplitude of m 3/2 . Thus, the production of heavy mass eigenstate is dominated at H ∼ H inf and the resultant abundance is of the same order of the transverse gravitino: where we have used m 2 m 3/2 3H inf 3m φ with m being the amplitude of the oscillating mass defined in App. B.
The fate of the inflatino significantly depends on its interactions unspecified here. Since we need to reheat the universe, there must be interactions of the inflaton superfield with the SUSY standard model (SSM) sector, which automatically introduce inflatino-SSM interactions. It may also decay into gravitino if kinematically allowed [70], although such a contribution is smaller than that of the pre-existing transverse gravitino (3.9).

Inflation model with a higher power potential
In this subsection, we briefly discuss the inflation model with higher power. The Kähler and superpotential are assumed to be All parameters are taken to be real and positive without loss of generality. It should be noticed that for n > 2, the presence of constant W 0 eventually leads to the spontaneous breaking of 2 or n so that φ obtains a finite VEV. When discussing a theory of higher n (> 2) hereafter, we implicitly assume that there is also a quadratic term W ∼ m φ φ 2 which is subdominant at the early stage of (p)reheating, and the 2 symmetry is maintained in a practical sense. For consistency with the Planck observation, we need modifications on the inflaton potential in the large field region, and it can also simultaneously affect the potential in the reheating stage. Such effects may be taken into account by including higher powers in the Kähler or superpotential. The theoretical construction is almost parallel to the case of n = 2 discussed so far after replacing m φ → λφ n −2 . One of the significant differences is the background evolution. Let us suppose that the inflaton field φ oscillates around the potential minimum φ = 0 under a potential V λ 2 φ 2(n −1) . The inflaton amplitude and the Hubble parameter decrease as φ amp ∝ a −3/n , H 2 ∝ a −6(n −1)/n . (3.35) This changes the abundance of gravitino, in particular the transverse one. For example, for the quartic inflaton potential n = 3, such as the Higgs-inflation like model, it is much more abundant than the case of n = 2: independently of the reheating temperature. Here we take = 1 as a reference. 21 Therefore in this case the transverse gravitino is problematic even if the reheating temperature is very low. The abundance of longitudinal gravitino does not change much compared with the case of n = 2 after replacing m φ → m eff φ = λφ n −2 amp . This is because the transverse gravitino production is dominated at H = H inf just after inflation while the longitudinal one is dominated at H ∼ m 0 3/2 . Thus the transverse gravitino abundance is much more sensitive to the background evolution than the longitudinal one. See the dashed line in the left panel of Fig. 3 in Sec. 5 for the transverse and longitudinal gravitino abundance for n = 3.
Another comment is that the inflaton "mass" ∼ λφ n −2 itself is a rapidly oscillating function for n > 2. It does not much affect the production of the light mass eigenstate of (ψ , v ⊥ ), but it can significantly change the production of the heavy mass eigenstate, which eventually becomes the inflatino (see Fig. 1). Therefore, the final inflatino abundance is expected to be sensitive to the power n. 22 Although the inflatino seems to be less harmful than the gravitino, precise discussion requires the information of interactions between the inflaton sector and the SSM sector for the reheating.

Comment on the linear term in Kähler potential
Here we comment on the case where the inflaton does not possess the 2 symmetry and there is a linear term in the Kähler potential (3.38) As seen below, the existence of this term drastically changes the picture of the longitudinal gravitino production.
For the moment, we focus on the case of n = 2. First we should notice that there arises a mixing between φ and z if there is a linear term in the Kähler potential as in Eq. (3.38). The mixing comes from The mixing angle is As shown in App. D, the oscillation of z is induced by the φ oscillation, and its amplitude is given by z amp ∼ θ φz φ amp . This induced z oscillation efficiently produces the longitudinal gravitino through the δz term in Eq. (3.27) for H m 3/2 . The longitudinal gravitino production rate is thus given bẏ (3.42) denotes the Polonyi decay width into the longitudinal gravitino pair.
This is consistent with Refs. [24,26]. In contrast to the case of the previous subsection, this process is interpreted as "decay" of φ because the relevant interaction involves a single power of δz and the mixing angle is constant. In other words, ψ "feels" the oscillation of φ linearly. Thus, this process is more and more effective at later time. The longitudinal mode produced during H m 3/2 will eventually become the present gravitino. The resulting abundance is estimated as where Γ inf denotes the total decay width of the inflaton. Cosmological consequences of this kind of scenario is extensively studied in Refs. [24,26]. This abundance is quite huge and it is one of the motivations to consider the 2 -symmetric model. 23 In the 2 -symmetric case c = 0, the decay process does not exist, because there is no mixing between z and φ. 24 A similar result holds for general n(> 2) after replacing m φ with the effective inflaton mass m eff φ ∼ λφ n−2 amp . In this case, however, the partial "annihilation" rate Γ(φφ → ψ ψ ) becomes smaller for late time and the gravitino production is suppressed compared with the case of n = 2.

Gravitino production in multi-superfield inflation
We consider a more realistic three-field case: the inflaton φ, the stabilizer X and the Polonyi field z . The Kähler and superpotential are given by This model has an approximate shift symmetry φ → φ + c with c being a real constant, and the inflaton is identified with 2 Re(φ) [29,30]. Without loss of generality, we can take all the coupling constants real and positive. Simple power-law chaotic inflation models are not favored according to the results from the Planck satellite [71], and hence we need modifications on the models to 23 If the inflaton and Polonyi sectors are effectively sequestered, the enhanced gravitino production can be avoided even when the inflation sector breaks the 2 symmetry. One of the simplest ways is to move the linear term to the superpotential by the Kähler-Weyl transformation before coupling it to the Polonyi sector. 24 In Ref. [20], it was pointed out that the gravitino does not couple to the heavy mass eigenstate of (φ, z ) system in the limit m z m φ . This corresponds to (m z /m φ ) 4 suppression in (3.43) [21]. In the 2 -symmetric case, this statement is trivial because m 3/2 does not explicitly contain φ but only z at the linear level. be consistent with observation [28][29][30][31][32][33][34][35][36]. Again, we are interested in the sub-Planckian regime and these modifications on the potential in the large field region do not much affect the following discussion.
From the viewpoint of gravitino production, the significant difference from the single-field case is that the stabilizer field X is close to zero during inflation. 25 It suppresses the oscillating gravitino mass m 3/2 ∼ λX φ n /M 2 Pl , leading to a suppressed gravitino production rate.

Dynamics
In contrast to the single-field model, m 3/2 m 0 3/2 just after inflation since X is stabilized at |X | |φ|. However, |X | does not remain to be very small because of the mixing between φ and X : Since we are assuming that all parameters in the model are real, the reality of φ and X is also maintained. Thus the potential around the minimum is given by According to App. D, the φ oscillation induces the X oscillation due to the mixing term. The overall magnitude remains to be about m 0 3/2 throughout the whole history of the universe. The oscillating part is the first term. 26 For n = 1, on the other hand, masses of φ and X are almost degenerate. As shown in App. D, for H m 0 3/2 , the amplitude is given by 27 Notably, this is almost constant until H ∼ m 0 3/2 . After that, φ ∼ X have the same order of oscillation amplitudes, meaning the maximal mixing of φ and X . This implies that (4.7) 25 Precisely speaking, we need K ∼ −|X | 4 /M 2 Pl to give X the Hubble mass and stabilize X 0 during inflation, and also it slightly deviates from the origin: X ∼ (m 0 3/2 /m φ )φ. However, it does not affect the following discussion. 26 Even if the term with n > 1 is dominant at the early stage, eventually the term with n = 1 becomes dominant. Otherwise, φ and X would be massless in the present vacuum. The bare mass term n = 1 can change the dynamics described above. 27 The Hubble mass term like ∼ H 2 |X | 2 does not change the discussion. Figure 2: Same as Fig. 1 but for multi-superfield inflation models.
The overall magnitude is about m 0 3/2 throughout the whole history of the universe. The first term is the oscillating part, which is suppressed by m 0 3/2 /H at H > m 0 3/2 compared with the single-field case.
The dynamics of Polonyi field z is similar to the single-superfield inflation case studied in Sec. 3.1, so we do not repeat the analysis here.
Based on this picture, below we estimate the gravitino abundance for n = 1 in which W = m φ X φ. The case of higher n is discussed in Sec. 4.4.

Transverse component
Let us estimate the transverse gravitino production rate for n = 1. As already shown, it is solely determined by the field-dependence of the gravitino mass m 3/2 . According to Eq. (4.7), it is suppressed by the small factor m 0 3/2 /H compared with the single-field case studied in Sec. 3.2. Thus we obtain the effective "annihilation" rate of the inflaton into the transverse gravitino pair as In contrast to the single-superfield inflation case, the production rate is time-independent at the early epoch, hence the transverse gravitino production is maximized at H ∼ m 0 3/2 . Thus the abundance of transverse gravitino is found to be This is similar to the abundance of longitudinal gravitino in the single-superfield inflation case (3.30) and much smaller than the gravitino abundance from thermal production. Thus it does not give significant effects on cosmology.

Longitudinal component
Next, let us discuss the longitudinal gravitino production. The discussion is parallel to the previous single-superfield inflation case, but one care is needed since it is |F X | ∼ |φ| ∼ H M Pl that dominantly breaks SUSY at H > m 0 3/2 . It may be useful to rewrite the Kähler and superpotential as where In this basis, we can follow the same method as that in Sec. 2.4 to derive the Lagrangian of the gravitino-fermion system. The goldstino is given by where We can define the 3 × 3 orthogonal matrix O that transforms the original basis into (v, v ⊥ ) basis: e.g., where s 2 1 + c 2 1 = s 2 2 + c 2 2 = 1 and O is the 3 × 2 matrix. It is easily checked that we have s 2 0 and v s 1 e γ 0 θ + Φ + + c 1 e γ 0 θ − Φ − and v (2) ⊥ ∼ e γ 0 θ z z for H m 0 3/2 , while s 2 1 and v e γ 0 θ z z and v ⊥ ∼ −s 1 e γ 0 θ + Φ + − c 1 e γ 0 θ − Φ − in the opposite limit H m 0 3/2 . Following the same procedure as in Sec. 2.4, we obtain the Lagrangian of the longitudinal gravitino and fermion system as where denotes the mass matrix defined in Eq. (2.66), whose determinant is Although a full expression of is lengthy, we can deduce the mass eigenvalues as Therefore the expression for m light is similar to the single-superfield inflation case. The lightest mass eigenstate, which should be mostly composed of z , is v (2) ⊥ for H m 0 3/2 and ψ c for H m 0 3/2 . A more detailed explanation for the mass eigenvalues is given in App. E. The right panel of Fig. 2 shows the schematic picture of the time evolution of the absolute values of the mass eigenvalues.

Abundance of longitudinal gravitino
The lightest mass eigenstate eventually becomes the present-day longitudinal gravitino. Hence, we here would like to estimate the production rate of the lightest mass eigenstate in the universe. Recalling that m 3/2 in the inflation model with a stabilizer field is given by Eq. (4.7), and repeating the same analysis for the production rate as that in Sec. 3.3, we find that an upper bound on the longitudinal gravitino production rate as for H > m 0 3/2 , and for H < m 0 3/2 . We again find that the resulting abundance is dominated by those created around H ∼ m 0 3/2 , with an amount of It is the same order of the abundance of the transverse gravitino and also comparable to the longitudinal gravitino abundance in the case of single-superfield inflation model (3.30). It is too small to give significant phenomenological effects. As noted in the single-superfield infaltion case, the contribution from the induced z oscillation to the longitudinal gravitino abundance is at most comparable to this upper bound. After all, in inflation models with a stabilizer field X , it is safe to neglect the nonthermal gravitino production after inflation. The dotted line in the right panel of Fig. 3 in Sec. 5 summarizes the transverse and longitudinal gravitino abundance in models with a stabilizer field.

Abundance of inflatino and stabilizino
There are two heavy mass eigenstates in the present model. They are roughly (ψ , v ⊥ ) at the late epoch (H < m 0 3/2 ). They are finally regarded as inflatino/stabilizino fields, since both v (1) ⊥ and v (2) ⊥ are composed of φ and X at late epoch. The heavy mass eigenvalues are given by m ± heavy ±m φ + 2α 2 ± m ± 3/2 , (4.24) and m ± 3/2 contains an oscillating term of (H ) (see App. E). Thus the production rate of the heavy mass eigenstates at H m 0 3/2 is similar to the case without the stabilizer field. This process is accessible only just after inflation when H ∼ m φ . The production of heavy states is expected to be dominated at H ∼ H inf , and the resulting inflatino/stabilizino abundance at late time is estimated to be 10 GeV , (4.25) which is similar to the inflatino abundance in the single-superfield inflation case. Again we note that the fate of the inflatino and stabilizino depends on the inflaton-SSM interactions which is not specified here. Generally the inflatino decay into the gravitino plus inflaton may be kinematically allowed depending on the soft SUSY breaking mass of the inflaton, but the branching fraction is expected to be small. However, depending on the inflatino/stabilizino branching fractions into gravitino, this channel can be the dominant source of nonthermal gravitino.

Inflation model with a higher power potential
So far we have focused on the case of quadratic inflaton potential n = 1. Now let us briefly discuss the case with a higher power n > 1. As described in Sec. 3.4, there are mainly two differences from the quadratic case. One is the change of the background evolution: see Eq. (3.35). The other is that the inflaton mass itself becomes a rapidly oscillating function. Although the full analysis is complicated, we can qualitatively discuss these effects. First, the transverse gravitino abundance is dominated at H ∼ H inf , but it is suppressed by a factor (m 0 3/2 /m φ ) 2 compared to the single-superfield case because of the suppression of the gravitino mass (4.5). The production of the longitudinal component is dominated at H ∼ m 0 3/2 . For example, for the quartic inflaton potential n = 2, we obtain where Γ φ denotes the total decay width of the inflaton. See the dashed line in the right panel of Fig. 3 in Sec. 5 for the transverse and longitudinal gravitino abundance for n = 2. The oscillation of inflaton mass itself also does not much affect the final gravitino abundance, since the mass eigenvalue of the light state is determined by m 3/2 , not the inflaton mass. On the other hand, the inflatino abundance can be enhanced since its production is dominated at the early epoch H ∼ H inf , and the oscillation of the effective inflaton mass itself directly contributes to the heavy mass eigenstates.

Conclusions
We have studied the nonthermal gravitino production during (p)reheating paying particular attention to the case of 2 symmetric large field inflation models. The result crucially depends on inflation models. In single-superfield inflation without a stabilizer field, production of the transverse gravitino is efficient and it can cause cosmological problems depending on the power law index of the inflaton potential. The longitudinal gravitino production is safely neglected. On the other hand, in multi-superfield inflation models with a stabilizer field, the transverse gravitino production is significantly suppressed and nonthermal gravitino production plays no important role in cosmology. Fig. 3 shows the gravitino abundance as a function of the reheating temperature T R for the single-superfield inflation models (left) and multi-superfield inflation models with a stabilizer field (right). The solid line shows a contribution from thermal production [11][12][13], while dashed (dotted) lines show nonthermally produced ones for the inflaton potential V ∝ φ p with p = 2 (p = 4). Here we have taken H inf = 10 13 GeV (left) and m 0 3/2 = 10 6 GeV (right). If the inflaton potential changes from φ 4 to φ 2 at some point, the prediction lies between these two lines. It is clearly seen that inflation models with a stabilizer predicts negligibly small nonthermal gravitino abundance. Therefore inflation models with a stabilizer field is motivated not only from the viewpoint of model building, but also from the requirement to avoid the nonthermal gravitino overproduction. Note that in this plot we have ignored some other nonthermal gravitino production processes such as those from Polonyi/inflatino decay since they are rather model-dependent [26,70]. However, inclusion of them does not much affect this conclusion.
Some comments are in order. In the most part we assumed the (nearly) minimal Kähler potential for the inflaton superfield for simplicity. The production rate is significantly modified if there is a 2 -symmetric Kähler potential of the form for the single-superfield inflaton and multi-superfield inflaton case, respectively. Since these terms directly give the large oscillating Polonyino mass like mz ∼ m φ φ 2 /M 2 Pl , the longitudinal gravitino production rate is expected to be significantly enhanced to the same level as the transverse production rate, if the coefficients of these terms are (1). Moreover, some inflation models, especially those constructed from the Jordan frame action, have a nonminimal Kähler potential of the inflaton sector itself which can potentially induce violent phenomena [72]. Also we assumed that the inflaton is a gauge singlet: for a gauge non-singlet inflaton, the structure becomes more complicated. We will come back to these issues in future.

A Notations and conventions
In this appendix, we summarize the notations and conventions used in this paper. We follow those of the textbook [73].

Gamma matrix
Here we summarize the notations and conventions related to the gamma matrices. We take the metric to be "almost plus", i.e., η µν = diag(−1, +1, +1, +1) for the flat space, and similar way for the FLRW one. Then, the Clifford algebra is defined as where [A, B ] + = A B + B A, and the hats denote the quantities in the curved space-time as it is also noted in the main text. We express the gamma matrices in the flat space without the hats. Through the vierbein e a µ , they are related: Here the vierbein is defined as The sign convention of Eq. (A.1) as well as that of the metric determine the (anti-)hermiticity of the gamma matrices. In our convention, γ 0 is anti-hermitian, while γ is hermitian: We define the short hand notations for the product of the gamma matrices as Note that they are anti-symmetric with respect to the indices. We also define the hermitian matrix γ * as where we have used the gamma matrices in the flat space, not in the curved space. In terms of it, the projection operators are defined as

Dirac/Majorana conjugation
The Dirac conjugation of a fermion ψ is defined as In order to define the Majorana conjugation, we need the charge conjugation matrix C . It is a unitary matrix that satisfies With this matrix, the Majorana conjugation is defined as For the Majorana fermion λ that satisfies the Dirac and Majorana conjugations are equivalent. 28 The explicit form of C in terms of γ depends on the representation of the gamma matrices.

Curvature
First we define the spin connection as where we have neglected the torsion since it contributes only to the four Fermi interaction. Here [...] denotes the anti-symmetrization with respect to those indices. The normalization of the antisymmetrization of n indices includes the factor of 1/n!. In terms of it, the Riemann tensor is defined as The Ricci tensor is defined as Finally we define the Ricci scalar as 16)

B Fermion production in the background field method
In this appendix, we discuss production of a Majorana fermion due to a time-dependent mass term: Though the time dependence of mass term typically arises from a coherently oscillating scalar field in the cosmological context, here we do not specify the origin of m (t ) and just assume that it oscillates with frequency of Ω. We follow the discussion given in Refs. [74][75][76]. See also Appendix A in Ref. [77].

B.1 Quantization
Let us start with the mode expansion of the Majorana field: The mode function obeys 29 puts the following constraints on the mode function: To quantize the Majorana field, we introduce the creation/annihilation operator: where v k ,s (t ) = −C −1 u T − k ,s (t ). One can see that the Majorana condition is fulfilled. We take the following normalization of the spinor: The quantization condition comes from the equal time anti-commutators: Together with the normalization of spinors, we obtain the following algebras for creation/annihilation operators: 29 Here we match the definition of the charge conjugation matrix C with that of Ref. [73].
where we have used the equation of motion for ψ in the first equality, and the bra-ket stands for the expectation value with respect to vacuum. One can check that the energy density only contains vacuum fluctuations for the initial condition given by Eqs. (B.15) and (B.16): where the phase space density is given by We have omitted the helicity subscript, h, in the wave function, u ± , because f ψ does not depend on helicity. The factor two in Eq. (B.19) represents the degrees of freedom for the Majorana fermion.
We assume the following form of the solution to Eq. (B.14): Initial values of A k ,h and B k ,h are obtained as Let us estimate the growth rate of f ψ at the very beginning, f ψ 1. At that time, we expect A k ω k + m and B k 0 at the leading order. We take into account the growth of B k perturbatively by using Eq. (B.22). Under this assumption, one can easily obtain In the second similarity, we have performed integration by parts and assumed k 2 m 2 . For given time t , the integration cancels out due to oscillations of the phase except for Ω−∆Ω 2ω k Ω+∆Ω with ∆Ω ∼ 1/t . Within this frequency range, the phase of m (t ) and e −2i ω k t cancels and B k grows linearly with time: Here and hereafter m stands for the amplitude of oscillating m (t ). Similar arguments lead to Plugging those approximated solutions into Eq. (B.20), we get This expression is valid as long as f ψ 1, namely q Ωt 1 with the resonance parameter being q ≡ m 2 /Ω 2 1. Finally, performing the phase space integral in Eq. (B.19), we obtain the master equation for Majorana fermion production due to its oscillating mass term m (t ): 16π Here we have introduced an order one factor , which depends on how m (t ) oscillates. For instance, in the case of m (t ) ∝ cos (Ωt ), we have = 1. Now let us assume that the oscillation of m (t ) is caused by oscillating (canonical) scalar field φ φ amp cos(m φ t ). The result (B.28) may be interpreted as the decay or annihilation of φ into ψ, depending on how m (t ) oscillates with time. Suppose that m (t ) ∝ φ n (t ), which involves Ω = j m φ with j = n, n − 2, n − 4, . . . . One can see that allowed processes depend on the parity of n. To avoid unnecessary complications, let us consider the case of n = 1 and n = 2 as an illustration. In this case, the relevant frequency is Ω = m φ (2m φ ) for n = 1 (2). For Ω = m φ (2m φ ), the decay (annihilation) rate is estimated as Here we have substituted n φ ∼ m φ φ 2 amp . Actually if m (t ) ∝ φ ∝ cos (m φ t ), this coincides with the perturbative decay rate of φ into two ψ particles calculated in the standard method in quantum field theory. In this case, the particle production can be reasonably interpreted as the decay of φ. On the other hand, if m (t ) ∝ φ 2 , this may be rather regarded as the annihilation of φ into ψ particles, with an annihilation rate of Γ(φφ → ψψ) ∼ n φ σv with σv ∼ ( m /φ 2 amp ) 2 . In the main text, we use the formula (B.29) for the gravitino production rate rather than (B.28), since reinterpreting the gravitino production as if it is caused by the inflaton decay may help readers understand the underlying physics.
Here are some comments. The calculations here are assuming that the background field (inflaton) is spatially homogeneous. This assumption is not always valid particularly if it is steeper than the quadratic one [78,79]. 31 For the quartic inflaton potential V ∝ φ 4 , for example, inflaton fluctuations with momenta of its effective mass scale develop due to the parametric resonance of the inflaton fluctuation itself and the initially homogeneous configuration may mostly become semi-relativistic waves. However, it does not much affect the estimate of resulting fermion abundance, since the equation of state of the universe remains the same no matter how the parametric resonance is efficient [79] and the fermion production rate is also roughly the same even if the production is caused by the decay/annihilation of the inflaton quanta, as explicitly shown above. Note also that in the most relevant case, i.e., the transverse gravitino production in the single-superfield inflaton case, the abundance is dominated by just the first few oscillations after inflation during which the coherence of the inflaton is maintained. Therefore, the gravitino abundance shown in Fig. 3 is not much affected by this subtlety. For a higher power V ∝ φ n with n > 4, we need much more care on the particle production rate, since in such a case the inflaton fluctuation develops and the equation of state may approach to the radiation-dominated one [79], which significantly modifies the naive estimate obtained by the assumption that the background is dominated by the homogeneous inflaton condensation. It also means that inflaton quanta may be highly relativistic so that the production rate may be suppressed by the Lorentz factor. A careful investigation of particle production rate in this situation is beyond the scope of this paper.

C Gravitino production in small-field inflation models
In this section we briefly comment on the gravitino production in small-field inflation models, such as new-inflation or hybrid inflation models. The known results of Refs. [21,24] are reproduced in our framework.

C.1 Single-superfield model
Let us consider a single-field new inflation model [93] as an example of small-field inflation models. The Kähler and superpotential are assumed to be K = |φ| 2 + |z | 2 − |z | 4 Λ 2 , (C.1) In this model, at the potential minimum φ = (M 2 /λ) 1/n , the superpotential takes a finite value: After a few Hubble time after inflation ends, the oscillation amplitude of the inflaton becomes much smaller than φ . Expanding φ = φ + δφ, Kähler and superpotentials are expressed as where the inflaton mass around the potential minimum is given by m φ = n M 2 / φ . This is the same as the single-field chaotic inflation model with linear term in the Kähler potential studied in Sec. 3.5 after the identification c → φ , except for the small linear δφ term in the superpotential. 32 Although the linear δφ term in the superpotential can induce the inflaton decay, the dominant contribution comes from the inflaton-induced z oscillation as studied in Sec. 3.5. The inflaton decay rate into the longitudinal gravitino pair and the resultant gravitino abundance is consistent with [21,24].

C.2 Multi-superfield model
Next, we consider multi-field new inflation model [94,95]: The potential minimum is 〈X 〉 0 and φ (M 2 /λ) 1/n . Expanding the field as φ = φ + δφ, the Kähler and superpotentials can be written as where the inflaton mass around the potential minimum is given by m φ = n M 2 / φ . Thus the gravitino production in this theory is also the same as that in the chaotic inflation model with linear term in the Kähler potential after the identification c → φ in Sec. 3.5. Although the model of Sec. 3.5 does not have a stabilizer X , the discussion is almost parallel, considering that there is a mixing term ∼ X * z in the scalar potential and also φ and X are maximally mixed with each other at least for H m 3/2 . As a result, the inflaton decay rate into the longitudinal gravitino pair and its abundance is consistent with [21,24].

(D.4)
In the mass eigenstate basis, the potential looks like Hence after the time t ∼ 2m 1 /m 2 12 , φ 1 and φ 2 may be regarded as maximally mixed with each other.

E Calculation of mass eigenvalues
In this appendix, we estimate the mass eigenvalues of the system of the longitudinal gravitino and matter fermions discussed in Secs. 3.3 and 4.3.