Heavy Fermion Production and Primordial N-Spectra

TL;DR

This paper investigates how heavy fermions produced during inflation, despite their high mass, can leave detectable imprints on primordial N-spectra, highlighting differences from their bosonic partners due to quantum effects.

Contribution

It provides a detailed analysis of heavy fermion production during inflation and compares their signatures in primordial spectra with those of bosonic partners, including quantum interference effects.

Findings

Fermions can produce detectable signatures in CMB despite high mass.

Classical contributions to N-spectra are similar for fermions and bosons.

Quantum interference causes distinguishable effects between fermions and bosons.

Abstract

We compute the non-adiabatic production of heavy fermion during inflation due to its coupling with inflaton. The coupling, partly inspired by axion monodromy, comes from the modulation of the fermion mass by the inflaton field. Even though the fermion mass is always much higher than the Hubble scale and the density of the produced fermions is low, they can still have detectable signatures in the cosmic microwave background. Their contributions to the primordial $N$-spectra are then analyzed in detail and compared with those from the fermion's bosonic super-partner. At the classical level, where the produced particles are treated as classical sources, the effect on the $N$-spectra is proportional to the density of the produced particles and the fermion and boson cases have the same contribution. Quantum interference, however, leads to distinction between the two cases. Implications of this similarity and distinction are discussed before making general remarks about the limitations of our calculation and possible ways of overcoming them.

Equations

$$|M{|}^2={\mu }_a^2+{\hat{\mu }}_a^2[a(\varphi )-2\pi n{]}^2\simeq {\mu }_a^2+g_a^2(\varphi -2\pi nf),$$

$$|M{|}^2={\mu }^2+2g^2{f}^2\cos \frac{\varphi }{f}.$$

$$gf<\frac{\mu }{\sqrt{2}}.$$

$$|M{|}^2\simeq {\mu }^2-2g^2{f}^2+g^2(\varphi -{\varphi }_n{)}^2={\mu }_b^2+g^2{\dot{\varphi }}^2(t-t_n{)}^2,$$

$$\frac{{\mu }_b}{\sqrt{2}g\dot{\varphi }}<\frac{\pi }{\dot{\varphi }/f}\Rightarrow gf>\frac{\sqrt{2}{\mu }_b}{\pi },$$

$$S=\int d\eta d^{3}x{a}^4\left[i\bar{\psi }{e}_a^{\mu }{\gamma }^a{\nabla }_{\mu }\psi -M(\varphi )\bar{\psi }\psi \right],$$

$$S=\int d\eta d^{3}x{a}^4\bar{\psi }\left[\frac{i}{a}{\gamma }^{\mu }{\partial }_{\mu }+i\frac{3}{2}H{\gamma }^0-M\right]\psi ,$$

$$S=\int d\eta d^{3}x\bar{\chi }\left[i{\gamma }^{\mu }{\partial }_{\mu }-aM\right]\chi .$$

$${\mathcal{L}}_0$$

$${\mathcal{L}}_I$$

$$(i{\gamma }^{\mu }{\partial }_{\mu }-aM)\chi =0.$$

$$\chi ={\int }_{\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{x}}\sum _r\left[u_r(\eta ,\mathbf{k})a_r(\mathbf{k})+v_r(\eta ,\mathbf{k})b_r^{†}(-\mathbf{k})\right],$$

$$u_r=\frac{1}{\sqrt{2}}\left(\begin{array}{c}{u}_{+}(\eta ,\mathbf{k}){\psi }_r(\mathbf{k})\\ u_{-}(\eta ,\mathbf{k}){\psi }_r(\mathbf{k})\end{array}\right),v_r=\frac{1}{\sqrt{2}}\left(\begin{array}{c}{v}_{+}(\eta ,\mathbf{k}){\psi }_r(\mathbf{k})\\ v_{-}(\eta ,\mathbf{k}){\psi }_r(\mathbf{k})\end{array}\right)$$

$${\gamma }^0=\left(\begin{array}{cc}I& 0\\ 0& -I\end{array}\right),{\gamma }^i=\left(\begin{array}{cc}0& {\sigma }^i\\ -{\sigma }^i& 0\end{array}\right),{\gamma }^5=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right).$$

$${\gamma }^0=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right),{\gamma }^5=\left(\begin{array}{cc}-I& 0\\ 0& I\end{array}\right).$$

$$v_{+}=-u_{-}^{\ast },v_{-}=u_{+}^{\ast },$$

$$i{u}_{+}^{\prime }+rk{u}_{-}-aM{u}_{+}=0,i{u}_{-}^{\prime }+rk{u}_{+}+aM{u}_{-}=0,$$

$$u_{\pm }^{\prime \prime }+\left[k^2+a^2{M}^2\pm i(aM{)}^{\prime }\right]u_{\pm }=0.$$

$$i{u}_{-}^{\prime }+rk{u}_{-}-aM{u}_{+}=0,i{u}_{+}^{\prime }-rk{u}_{+}-aM{u}_{-}=0.$$

$$u_{\pm }^{\prime \prime }+(k^2+a^2{M}^2)u_{\pm }-\frac{(aM{)}^{\prime }}{aM}(u_{\pm }^{\prime }\pm irk{u}_{\pm })=0.$$

$$|M(\varphi ){|}^2={\tilde{\mu }}^2+g^2(\varphi -{\varphi }_n{)}^2={\tilde{\mu }}^2+g^2{a}_n^2{\dot{\varphi }}^2(\eta -{\eta }_n{)}^2,$$

$$u_{+}^{\prime \prime }+\left(k^2+a_n^2{\tilde{\mu }}^2+g^2{a}_n^4{\dot{\varphi }}^2(\eta -{\eta }_n{)}^2+i{a}_n\frac{g^2{a}_n^2{\dot{\varphi }}^2(\eta -{\eta }_n)}{\sqrt{{\tilde{\mu }}^2+g^2{a}_n^2{\dot{\varphi }}^2(\eta -{\eta }_n{)}^2}}\right)u_{+}=0$$

$$u_{+}(\eta )$$

$$u_{-}(\eta )$$

$$\Phi (\eta )={\int }_0^{\eta }d{\eta }^{\prime }\sqrt{{\nu }_0^2+{\nu }_1^4{\eta }^{\prime 2}}\simeq \{\begin{array}{ll}\frac{1}{2}{\nu }_1^2{\eta }^2+\frac{{\nu }_0^2}{2{\nu }_1^2}\log \left(\frac{2{\nu }_1^2}{{\nu }_0}\eta \right)+\frac{1}{4}\frac{{\nu }_0^2}{{\nu }_1^2},& \eta \to +\infty \\ -\frac{1}{2}{\nu }_1^2{\eta }^2-\frac{{\nu }_0^2}{2{\nu }_1^2}\log \left(\frac{2{\nu }_1^2}{{\nu }_0}|\eta |\right)-\frac{1}{4}\frac{{\nu }_0^2}{{\nu }_1^2},& \eta \to -\infty \end{array}.$$

$$\omega \simeq \{\begin{array}{ll}-{\nu }_1^2\eta \left(1+\frac{{\nu }_0^2}{2{\nu }_1^4{\eta }^2}\right),& \eta \to -\infty \\ {\nu }_1^2\eta \left(1+\frac{{\nu }_0^2}{2{\nu }_1^4{\eta }^2}\right),& \eta \to +\infty \end{array}.$$

$$a_{n}M\simeq \{\begin{array}{ll}-{\nu }_1^2\eta \left(1+\frac{{\nu }_0^2-k^2}{2{\nu }_1^4{\eta }^2}\right),& \eta \to -\infty \\ {\nu }_1^2\eta \left(1+\frac{{\nu }_0^2-k^2}{2{\nu }_1^4{\eta }^2}\right),& \eta \to +\infty \end{array}$$

$$\eta \to -\infty ,$$

$$\eta \to +\infty ,$$

$$\eta \to -\infty ,$$

$$\eta \to +\infty ,$$

$${\hat{a}}_{\mathbf{k}}(\eta )={\alpha }_{\mathbf{k}}{a}_{\mathbf{k}}+{\beta }_{\mathbf{k}}{b}_{-\mathbf{k}}^{†},{\hat{b}}_{-\mathbf{k}}^{†}(\eta )=-{\beta }_{\mathbf{k}}^{\ast }{a}_{\mathbf{k}}+{\alpha }_{\mathbf{k}}^{\ast }{b}_{-\mathbf{k}}^{†}$$

Full text

SLAC-PUB-16993, SU-ITP-17/07

Heavy Fermion Production and

Primordial $N$-Spectra

Danjie Wenren

Stanford Institute for Theoretical Physics, Stanford University, Stanford, CA 94305, USA and

SLAC National Accelerator Laboratory, 2575 Sand Hill Rd., Menlo Park, CA 94025, USA

Abstract

We compute the non-adiabatic production of heavy fermion during inflation due to its coupling with inflaton. The coupling, partly inspired by axion monodromy, comes from the modulation of the fermion mass by the inflaton field. Even though the fermion mass is always much higher than the Hubble scale and the density of the produced fermions is low, they can still have detectable signatures in the cosmic microwave background. Their contributions to the primordial $N$-spectra are then analyzed in detail and compared with those from the fermion’s bosonic super-partner. At the classical level, where the produced particles are treated as classical sources, the effect on the $N$-spectra is proportional to the density of the produced particles and the fermion and boson cases have the same contribution. Quantum interference, however, leads to distinction between the two cases. Implications of this similarity and distinction are discussed before making general remarks about the limitations of our calculation and possible ways of overcoming them.

1 Introduction

The cosmic microwave background (CMB) observation and potentially the large-scale structure (LSS) observation in the future offer exciting opportunities to test high energy physics whose energy scale is far beyond what other types of experiments can achieve. In particular the oscillatory features of those observations can in principle constrain the possible “ultra-violet complete” theories of gravity. Such theories are needed in order to fully understand inflation, which is currently the most compelling theory framework of explaining the CMB data. With the significant amount of CMB data being collected by Planck [1, 2, 3] it is worth studying from the theory perspective what physics mechanism can generate primordial seeds and providing searching templates for the corresponding signals.

String theory as one of the most prominent candidates of theories of quantum gravity has offered a variety of ideas for inflation, of which axion monodromy [4, 5, 6] is a very natural one. Inspired in part by this mechanism, authors of [7] consider theories where heavy boson fields interact with inflaton through non-derivative couplings. Such interaction comes from the fact that the boson mass is modulated by the inflaton field. As is well known in quantum field theory, such modulation can lead to non-adiabatic production of particles. It is found in [7] that in a well-defined window of parameters as discussed there in, current CMB data is sensitive to such non-adiabatic production, even if the particle mass is always much higher than the Hubble scale $H$. The production’s contribution to the primordial $N$-spectra is studied in detail and it is shown that the mechanism can generate oscillatory non-Gaussianity parametrically larger than that from the previously studied resonant non-Gaussianity.

In this work we consider the effects of non-adiabatic production of heavy fermions which are super-partners of the bosons as in [7]. As shown in the appendix of that paper, the radiative corrections to the inflaton $N$ point functions from these super-partners are suppressed if some degree of supersymmetry is assumed. Therefore in that case one only needs to consider the contributions from the particle production effects. It remains unclear, however, the relative scale between the contributions from boson production and those from fermion production, and if there can be novel feature shapes from fermion production.

In this work we will find that at the classical level, namely treating the produced particles as classical source, the two types of particles have the same contribution. The contribution due to quantum interference, however, is different between the fermion case and the boson case.

The basic setup and particle mass as a function of inflaton field are reviewed in section 2. In order to make a direct comparison with the boson results we work with the same mass function and in the same parameter regime as [7]. The free fermion equation of motion in the inflation background is then solved using WKB approximation in section 3 along with the corresponding Bogoliubov transformations. Due to the relatively more complicated equation of motion and large minimum mass of the fermion field, a different method needs to be used than previous works in literature [8, 9]. We found that in the parameter regime we consider, where the produced particle density is low and gets diluted rapidly by the exponential expansion of the universe, the Bogoliubov transformation is similar as in the boson case. In section 4, we compute two types of contributions to the inflaton $N$ point functions, one from treating the produced particle as classical source and the other from the quantum interference of the heavy particle fields. Comparisons between the fermion case and the boson case are also made in that section. Summary and a brief discussion of future directions are presented in section 5.

2 Setup

It is discussed in [7] that the radiative correction from heavy boson fields can be suppressed by the correction from their fermion partners assuming some degree of supersymmetry. As a result, the only contributions to the $N$-spectra will be those from particle production effects. In order for the radiative corrections to cancel, the boson mass $m_{b}$ and fermion mass $m_{f}$ need to satisfy $|m_{b}|^{2}=|m_{f}|^{2}$ apart from the matching of numerical factors from particle and anti-particle doubling and helicity doubling. This does not, however, immediately determine $m_{f}$ as a function of the inflaton field $\phi$ since in a supersymmetric theory $m_{f}$ can be complex. We will work with an $m_{f}$ in the following that simplifies the analysis.

Denoting the fermion field as $\psi$, its mass $M$ as a function of the inflaton field $\phi$ has two possible forms [7].

1. (a)

This is when the theory has multiple fields that undergo monodromy and each of them reaches its minimum mass when the inflaton traverses an underlying period in the axion field space. Let $n$ label different species of such fields and their masses can be written as

[TABLE]

where we used $a(\phi)=\phi/f$ with $f$ being the decay constant and $g_{a}=\hat{\mu}_{a}/f$. 2. (b)

In this case there is only one field whose mass is sinusoidally modulated by the inflaton field

[TABLE]

The positivity of $|M|^{2}$ requires

[TABLE]

We will use the WKB method to solve wave equations of the $\psi$ field and compute the number density of the produced $\psi$ particles. The production happens near the minimum mass region and for the $n$th such region, the mass can be expanded as

[TABLE]

where $\phi_{n}=(2n+1)\pi f$ and $\phi(t_{n})=\phi_{n}$. In order for the WKB method to be valid, the production region should be within each oscillation period. The production region can be determined by maximizing $\dot{\omega}(p)/\omega^{2}(p)$ where $\omega(p)^{2}=p^{2}+|M|^{2}$ is the physical frequency of $\psi$ modes. A straightforward calculation shows that the maximum is at $|t-t_{n}|=\mu_{b}/(\sqrt{2}g\dot{\phi})$. Therefore the requirement of production happening well inside each oscillation period translates into the following condition

[TABLE]

which can be consistent with condition (2.3).

Notice that in both cases the mass has the same quadratic dependence on $\phi$ in the production region so we will first discuss their single production events and $N$ point correlation functions in a similar way. To explicitly evaluate the contributions of particle production events to the $N$ point functions we focus on case (b) since it provides novel searching templates while the other case overlaps strongly with existing ones.

The action of the fermion field $\psi$ can then be written down as

[TABLE]

where the vierbein $e^{a}_{\mu}$ is the constant $4\times 4$ matrix $\frac{1}{a}I$ with $\nabla_{\mu}$ being the covariant derivative. Substituting their explicit forms in the FRW background into (2.6) we get

[TABLE]

and the gamma matrices satisfy the usual anti-commutation relation $\{\gamma^{\mu},\gamma^{\nu}\}=-2\eta^{\mu\nu}$ where $\eta_{\mu\nu}$ is the metric for flat spacetime $\rm{diag}(-1,+1,+1,+1)$ 111Note that strictly speaking, $\gamma^{\mu}\partial_{\mu}$ should be understood as $\gamma^{a}I_{a}^{\mu}\partial_{\mu}$ and $\gamma^{a}$ are the matrices that satisfy the usual anti-commutation relation in flat spacetime. For notational simplicity we will ignore this detail in this work.. Defining $\psi=a^{-3/2}\chi$ the action simplifies to

[TABLE]

As in [7] we consider the parameter regime where the heavy particles do not back-react on the inflation background $\phi_{0}$ and the $\chi$ equation of motion is solved in this background. More explicitly we will work with the following free and interaction Lagrangians

[TABLE]

3 WKB solution for fermion field and Bogoliubov transformation

The free equation of motion of $\chi$ can be derive from (2.9) as

[TABLE]

It can be solved using Fourier transformation. To begin with, let us write $\chi(\eta,\bf{x})$ as

[TABLE]

where $\int_{\bf{k}}$ is a short-hand notation for $\int\frac{d^{3}\bf{k}}{(2\pi)^{3}}$ and $r=\pm 1$ denotes the helicity of the $\chi$ modes. The mode functions can be further decomposed as

[TABLE]

where the two-component vectors $\psi_{r}$ are helicity eigenstates satisfying ${\bf k}\cdot{\bf\sigma}\psi_{r}=rk\psi_{r}$ with $\sigma^{i}$ being the Pauli matrices.

Before we substitute (3.2) into (3.1) we also need to select a basis for the gamma matrices. While in principle different bases should lead to the same physical result, we will use the Dirac basis where

[TABLE]

It works better than the other commonly used one, namely the Weyl basis where the $\gamma^{i}$ matrices are the same while the other two are

[TABLE]

This is because of the more complicated wave equations the latter leads to. We will first derive the wave equations for $u_{\pm}$ in Dirac basis and then demonstrate why it is not as convenient in Weyl basis. The mode functions $v_{\pm}$ are related to $u_{\pm}$ as

[TABLE]

so we will not explicitly write out the wave equations for $v_{\pm}$ [8].

3.1 Wave equation in Dirac basis

Using (3.1) and other definitions above, one can get the equations of motion for $u_{\pm}$ as

[TABLE]

where prime means taking derivative with respect to $\eta$. These first order equations can then be combined into second order equations similar to the Klein-Gordon equation

[TABLE]

This form is relatively simple and we will solve it using WKB approximation.

One can also easily verify using (3.7) that $(|u_{+}|^{2}+|u_{-}|^{2})^{\prime}=0$ so we will normalize the mode functions as $|u_{+}|^{2}+|u_{-}|^{2}=2$, which can be used as a consistency check for the Bogoliubov coefficients we get. Note that in deriving this normalization condition we used the fact that $M$ is real. In a supersymmetric theory this is not always the case so we will need to be careful about the definition of $M$ in the following. We will prescribe a particular analytic structure for $M(\eta)$ in section 3.3 which ensures that it is real on the real $\eta$ axis.

3.2 Wave equation in Weyl basis

In this basis the first order wave equations become

[TABLE]

Besides the problem of now having the mode functions depending on helicity $r$, we also get more complicated second order equations of motion

[TABLE]

Since this looks much less tractable than the Dirac basis case (3.8), the rest of our calculation will be in Dirac basis.

3.3 Solution in the WKB approximation

In both cases (a) and (b) we can write the mass near the $n$th production region as

[TABLE]

where $\phi=\phi_{n}$ is the point where $\psi$ reaches its minimum mass and dot represents the derivative with respect to physical time $t$. We have also used the relation $d/d\eta=a(d/dt)$. In section 3.1 we have restricted $M$ to be on the real axis but we still need to fix a branch cut when we take the square root of (3.11). This is particularly an issue in the analysis here due to the imaginary term in (3.8). This can be better seen if we plug (3.11) in and expand the derivative term

[TABLE]

Note that since within Hubble the expansion is relatively slower than the changing in mass, we have treated $a$ as a constant $a_{n}$. Also for notational simplicity in this section we will make the substitution $\Delta\eta=\eta-\eta_{n}\to\eta$.

A typical way of using WKB method to solve (3.12) involves getting the exact solution of it and then matching it with the WKB solutions in the $\eta\to-\infty$ limit and the $\eta\to+\infty$ limit222Note that we have made the substitution $\eta-\eta_{n}\to\eta$ for notational simplicity. [8, 9, 10]. However unlike the wave equations in those works, it seems rather intractable to exactly solve the wave equation (3.12) due to the inverse square root term. A simplification one might consider is that since the matching happens in the $|\eta|\to\infty$ region, the wave equation can be simplified in that region. This indeed happens because in that region the imaginary term becomes $iga_{n}^{2}|\dot{\phi}|$ in the $\eta\to+\infty$ limit and $-iga_{n}^{2}|\dot{\phi}|$ in the opposite limit. However this also shows that the solutions to (3.12) cannot be extended to analytic functions on the full $\eta$ plane including $\infty$ since the solution should satisfy different differential equations at $-\infty$ and $+\infty$. We explore this method in the Appendix A and argue that it does not work due to the complicated Stokes phenomenon of the original wave equation (3.12). In this section we instead use a different approach which directly rotates one WKB solution in one limit to the one in the other limit along a semi-circle on the complex $\eta$ plane.

The WKB solution to (3.12) takes the form 333 Note that in our solution the definition for $A$ and $B$ is different from that in [8] because of a sign difference in the free Hamiltonian (3.23).

[TABLE]

with $A$ and $B$ being constants and $\omega^{2}=k^{2}+a_{n}^{2}\tilde{\mu}^{2}+g^{2}a_{n}^{4}\dot{\phi}^{2}\eta^{2}$ the Fourier frequency444Note again that the lower limit [math] in the integral means the production point since we have made the substitution $\eta-\eta_{n}\to\eta$ for the purpose of having less cumbersome equations.. As usual we would have different $A$’s and $B$’s in the $\eta>0$ region and $\eta<0$ region due to the non-adiabaticity in the production region. In general one would need a recurrence relation for $A_{n},B_{n}$ after the $n$th production event in terms of $A_{n-1},B_{n-1}$ before the $n$th production event. However our method will only be able to give result for $B_{n}$ in the region where $|B_{n-1}|\simeq 0$. As argued in [7] this is actually good enough for our purpose since during each production event only a small amount of particles get produced and they are quickly diluted by inflation. We will see later that the density of the produced fermions is proportional to $|B|^{2}$ so we can always assume $|B|$ to be small. Moreover, the normalization condition $|u_{+}|^{2}+|u_{-}|^{2}=2$ enforces $|A|^{2}+|B|^{2}=1$ , therefore for each production event, it would be sufficient to start with only positive frequency modes, that is $A=1$ and $B=0$ in the $\eta\to-\infty$ limit and then solve for $B$ in the opposite limit. Without loss of generality, we will only make the rotation explicitly for $u_{+}$ since $u_{-}$ has the same set of parameters.

Defining $\nu_{0}^{2}=k^{2}+a_{n}^{2}\tilde{\mu}^{2}$ and $\nu_{1}^{4}=g^{2}a_{n}^{4}\dot{\phi}^{2}$ we can immediately write down the following two approximations needed for computing the WKB solutions, namely the phase

[TABLE]

and frequency

[TABLE]

The other one in need is the mass $M$. One needs to be more careful about it because if the usual branch cut is used, the rotation would cross the branch cut at $\arg\eta=\pm\frac{\pi}{2}$ as shown in Fig 1. Because $M$ changes sign at the branch cut, there is no way to keep rotating the same analytic function to the other side of the real axis.

To show this more precisely, we can use the following mass function approximation

[TABLE]

and spell out the WKB solution

[TABLE]

where we have defined $x=\nu_{1}\eta$. It is clear that one cannot rotate the first solution to the second one555Notice that neither of the terms in the $\eta\to+\infty$ limit matches the $x$ dependence in the $\eta\to-\infty$ limit. due to the branch cut along $\arg\eta=-\frac{\pi}{2}$.

To solve this problem a different branch cut needs to be used for the mass function, which is shown in Fig 2. Essentially the branch cut of the square root function is taken from the negative real axis to the positive axis, so that the branch cuts on the $\eta$ plane are taken to the real axis. With this definition of branch cut, one could rotate the WKB solution in the $\eta\to-\infty$ region along a semi-circle to the $\eta\to+\infty$ region, where at both ends of the path the mass is positive. One caveat of this rotation is that in order to keep the mass positive in the end, one would still need to cross the branch from $+\infty(1-i\epsilon)$ to $+\infty(1+i\epsilon)$, where $\epsilon>0$ is an arbitrarily small number. We expect, however, this infinitesimal rotation would not change the Bogoliubov coefficients in a significant way.

With this definition of the mass function, we can write down the WKB solutions at the two ends of the rotation

[TABLE]

where the second limit is taken to be on the lower side of the real axis. It is now clear that (3.20) rotates into the second term of (3.21). To make this rotation, let $x=\rho e^{i\phi}$ and $\phi$ goes from $-\pi+\epsilon$ to $-\epsilon$ where $\epsilon>0$ is an arbitrarily small number. The $(-x)^{\frac{i\nu_{0}^{2}}{2\nu_{1}}}$ factor can be written as $\left(e^{i\pi}e^{-i\pi}\rho\right)^{\frac{i\nu_{0}^{2}}{2\nu_{1}}}$ before rotation and becomes $\left(e^{i\pi}\rho\right)^{\frac{i\nu_{0}^{2}}{2\nu_{1}}}$ afterwards. This should match with the second term in (3.21) and therefore yields $B=e^{-\frac{\pi\nu_{0}^{2}}{2\nu_{1}^{2}}}$.

One might worry about the appearance of the first term in (3.21). To address this we notice that during the rotation, $e^{\pm\frac{i}{2}x^{2}}=e^{\frac{\rho^{2}}{2}(\mp\sin 2\phi+\pm i\cos 2\phi)}$. If one starts with (3.21) and rotates it to (3.20), the first term will be exponentially suppressed compared to the second one and therefore should not be kept in this WKB approximation [11]. In other words, at the level of WKB approximation, (3.20) and (3.21) are the same solution and we cannot fix $A$ using this rotation method. If one is able to solve (3.12) exactly then in principle $A$ can be computed. This possibility is explored in Appendix A.

3.4 Bogoliubov transformation

After each production event the creation and annihilation operators need to be redefined in order to diagonalize the Hamiltonian with a new set of mode functions $u_{+}$ and $u_{-}$. This can fix the Bogoliubov transformation

[TABLE]

in terms of $A$ and $B$. Note that the anti-commutation relation $\{\hat{a}_{\bf k},\hat{a}^{\dagger}_{\bf k^{\prime}}\}$ and $\{\hat{b}_{\bf k},\hat{b}^{\dagger}_{\bf k^{\prime}}\}$ immediately lead to the normalization condition $|\alpha_{\bf k}|^{2}+|\beta_{\bf k}|^{2}=1$. In this subsection we will first solve for $\alpha$ and $\beta$ by explicitly diagonalizing the Hamiltonian and then show that it is equivalent to rewriting the $\chi$ expansion (3.2) in terms of purely positive modes and negative modes.

The free Hamiltonian can be derived from (2.9) as666Despite the fact that it is straightforward to compute the Hamiltonian, we still find a sign difference as compared with [8]. It seems our result agrees with most textbooks. As a result our result (3.4) also has a sign difference with that of [8].

[TABLE]

where we have used the free equation of motion of $\chi$ so that $H_{0}$ has only one time derivative term. Using the fermion field expansion (3.2) it can be expressed in terms of $a$ and $b$ as

[TABLE]

where

[TABLE]

and $E^{2}+|F|^{2}=\omega^{2}$. Using the reverse of (3.22), that is,

[TABLE]

and enforcing the vanishing of off diagonal terms (terms proportional to $\hat{b}\hat{a}$ and $\hat{a}^{\dagger}\hat{b}^{\dagger}$) lead to the following relations

[TABLE]

In order to relate $\alpha$ and $\beta$ to $A$ and $B$ one can plug in the WKB solution (3.13) and (3.14) and get $E=\omega(|A|^{2}-|B|^{2})$ and $F=-2\omega AB$. So (3.27) gives $\alpha/\beta=-A/B^{*}$ and $|\beta|^{2}=|B|^{2}$. This means that the diagonalization of the Hamiltonian does not completely fix $\alpha$ and $\beta$. We pick $\alpha=A$ and $\beta=-B^{*}$ for the requirement that $\chi$ can be written as a sum of purely positive and negative modes in terms of $\hat{a}$ and $\hat{b}$ as shown in the following.

Denote $\hat{u}_{\pm}$ and $\hat{v}_{\pm}$ as purely positive and negative modes, respectively then

[TABLE]

where $\int\omega$ is a short hand for $\int_{0}^{\eta}d\eta^{\prime}\omega(\eta^{\prime})$. Then the WKB solution (3.13), (3.14) can be written as

[TABLE]

With this relation the integrand of the $\chi$ expansion (3.2) becomes

[TABLE]

Comparing (3.30) and (3.22) shows that we need to define $\alpha=A$ and $\beta=-B^{*}$ in order for

[TABLE]

Using this and the result of section 3.3 we have the Bogoliubov coefficients $\beta_{\bf k}=-e^{-\frac{\pi\nu_{0}^{2}}{2\nu_{1}^{2}}}$ and $\alpha_{\bf k}$ a small phase rotation.

4 Power spectrum and non-Gaussianity from fermion production

In this section we compute the contributions to the inflaton correlation function from the fermion production events and compare them with those from the boson production events. As argued earlier, if the bosons and fermions are super partners and satisfy $|M_{b}|^{2}=|M_{f}|^{2}=M^{2}$, their radiative contributions from loop diagrams involving inflaton cancel with each other, therefore we only consider effects from particle production.

The general inflaton $N$-point function can be computed using the $in$-$in$ formalism

[TABLE]

where it is understood that the $\delta\phi$’s on the left hand side are evolved in the full theory and those on the right hand side are the free fields. $H_{I}(\eta)$ is the interaction Hamiltonian, which from (2.10) can be easily derived as

[TABLE]

The $|in\rangle$ state is the initial state which is vacuum in terms of the $a_{r}({\bf k})$ and $b_{r}({\bf k})$ operators. After the $n$th production event, it can be expressed as the squeezed state in the Fock space defined by the $\hat{a}_{r}({\bf k})$ and $\hat{b}_{r}({\bf k})$ operators

[TABLE]

where $\mathcal{N}$ is a normalization constant and $|\hat{0}\rangle$ is the vacuum state in the Fock space defined by $\hat{a}_{r}({\bf k})$ and $\hat{b}_{r}({\bf k})$. The last approximation is due to the fact that $\beta$ is exponentially small and $\alpha$ is a small phase rotation. This can be checked by the following result

[TABLE]

which gives the usual particle density $|\beta_{k}^{2}|$.

To compute the inflaton correlation functions we would also need $\bar{\chi}\chi_{\bf k}(\eta)$. Using the expansion (3.2) and (3.31) we have

[TABLE]

where $\hat{u}_{r}({\bf k})$ and $\hat{v}_{r}({\bf k})$ are defined as

[TABLE]

and $\hat{u}_{\pm},\hat{v}_{\pm}$ are defined in (3.28). As we will show in the following, ${\bf k}^{\prime}$ integrals contain $\mathcal{O}(|\beta|^{2})$ or $\mathcal{O}(|\beta|)$ factors and peak at order $a_{n}\sqrt{g|\dot{\phi}|}$. Since the external momenta ${\bf k}$ are of order $H$, which is parametrically smaller than $a_{n}\sqrt{g|\dot{\phi}|}$, we can approximate ${\bf k}^{\prime}\pm{\bf k}$ by ${\bf k}^{\prime}$. Therefore we have

[TABLE]

where the last approximation is due to the fact that the ${\bf k}^{\prime}$ integral peak location $a_{n}\sqrt{g|\dot{\phi}|}$ is parametrically smaller than $aM$. Similarly

[TABLE]

where again $\int\omega$ is a short hand notation for $\int_{\eta_{n}}^{\eta}d\eta^{\prime}\omega(\eta^{\prime},{\bf k}^{\prime})$. Combining these results we have

[TABLE]

In order to make direct comparisons between the fermion and boson contributions we also compute the relevant quantities for the boson case. The interaction Hamiltonian is

[TABLE]

and the source field is

[TABLE]

where we have dropped the hats on the $a$ operators and from this point on it is understood that they mean the operators defined with respect to purely positive and negative frequency modes.

It is argued in [7] that the oscillating terms in (4.13) only contribute to $N$-point functions of $\delta\phi$ through $(N+2)$-point vertex insertions. We elaborate on this point in section 4.1. Order $|\beta|^{2}$ contributions and order $|\beta|$ contributions are computed in section 4.2 and 4.3, respectively.

4.1 Contributions from various interaction vertex profiles

The oscillating terms in (4) and (4.13) have either two creation or two annihilation operators and can contract with the two annihilation or creation operators from the squeezed state. Therefore it is possible to have only one $\beta$ factor when evaluating the $\delta\phi$ correlators as opposed to the order $|\beta|^{2}$ contributions when only using the $a^{\dagger}a$ terms. It is argued in [7] that such order $\beta$ contribution only shows up when using $(N+2)$ point vertex for computing $N$ point function. The reason is that all the other use of such terms do not resonate and therefore their contribution is negligible. Since such argument applies to both boson and fermion, we will work with boson in this subsection.

A generic term after expanding the $N$ point function (4) is

[TABLE]

where there are $l_{1}$ insertions from the anti-time-ordered operator and $l_{2}$ insertions from the time-ordered operator, with the $i$th vertices being $(o_{i}+2)$ and $(o^{\prime}_{i}+2)$ point vertices, respectively. To contract each of the $\delta\phi_{\bf k}(0)$ operators with a $\delta\phi_{\bf k}(\eta)$ or $\delta\phi_{\bf k}(\zeta)$ operator from $H_{I}$ we also have $\sum_{i=1}^{l_{1}}o_{i}+\sum_{i=1}^{l_{2}}o^{\prime}_{i}=N$. Also notice that since on the right hand side of (4), $\delta\phi_{\bf k}$’s are evolved using the free inflaton Hamiltonian the expectation value factorizes as in (4.1) and the $\delta\phi_{\bf k}$ correlator is computed in the $\delta\phi$ vacuum state while the $\chi^{2}_{\bf k}$ correlator is computed in the $\chi$ squeezed state.

In order to look for the contributions of order $\beta$ or even $\mathcal{O}(1)$, we first consider the case where the first insertion comes from the anti-time-ordered operator, i.e., $l_{1}>0$. Similar to (4.3) we have the squeezed state for the boson case

[TABLE]

If there is no $\frac{\beta^{*}}{\alpha}aa$ on the left of the $\chi^{2}$ operators in (4.1), then the first $\chi^{2}$ factor $\chi^{2}_{-\sum_{i=1}^{o_{1}}{\bf k}^{\prime}_{1,i}}$ must contribute an $aae^{-2i\mu(t_{1}-t_{p})}$ factor777$t_{p}$ denotes the time of minimum $\chi$ mass.. Since $\delta\phi_{{\bf k}^{\prime}_{1,i}}(\eta_{1})$ for $\forall i$ have to contract with some $\delta\phi_{{\bf k}}(0)$ later in the product, they all contribute an oscillating factor $e^{-ik^{\prime}_{1,i}\eta_{1}}$. Using $t_{1}-t_{p}=-\frac{1}{H}\log\frac{\eta_{1}}{\eta_{p}}$, we have the oscillating part of the $\eta_{1}$ integral as $\exp\left(-i\sum_{i}^{o_{1}}k^{\prime}_{1,i}\eta_{1}+i\frac{2\mu\pm\omega}{H}\log\frac{\eta_{1}}{\eta_{p}}\right)$ where $e^{\pm i\frac{\omega}{H}}$ comes from the mass derivatives. It has no resonance in the parameter region we work with where $\mu>\omega$. Therefore there must be a $\frac{\beta^{*}}{\alpha}aa$ factor on the left. In order to get contribution of order $\mathcal{O}(\beta)$ we cannot have a $\frac{\beta}{\alpha^{*}}a^{\dagger}a^{\dagger}$ on the right and the last $\chi^{2}$ insertion must give an $a^{\dagger}a^{\dagger}e^{2i\mu(t-t_{p})}$ factor.

If the last vertex is from the time-ordered operator ($l_{2}>0$), then $\delta\phi_{{\bf p}_{1,i}}(\zeta_{1})$ must contribute an $e^{i\sum_{i=1}^{o^{\prime}_{1}}p_{1,i}\zeta_{1}}$ factor. It does not resonate with the $\exp\left(-i\frac{2\mu}{H}\log\frac{\zeta_{1}}{\eta_{p}}\right)$ factor from $\chi^{2}_{-\sum_{i}{\bf p}_{1,i}}(\zeta_{1})$ for the same reason as in the previous paragraph. If the last vertex is from the anti-time-ordered operator ($l_{2}=0$), then the oscillating part of its time integral is

[TABLE]

which resonates at $-(\sum_{i}k^{\prime}_{l_{1},i})\eta_{l_{1}}=\frac{2\mu\pm\omega}{H}$. All the other time integrals resonate at either $k\eta=-\frac{\omega}{H}>-\frac{2\mu\pm\omega}{H}$ or $-\frac{2\mu\pm\omega}{H}$, which means the resonance point of $\eta_{l_{1}}$ is the earliest. However the anti-time-ordering requires $\eta_{l_{1}}$ to be the latest, so there is still no resonance in the integral region.

Next we consider the case where the first insertion is from the time-ordered operator (i.e., $l_{1}=0$). Similar to the previous case there must be a $\frac{\beta}{\alpha^{*}}a^{\dagger}a^{\dagger}$ factor on the right and no $\frac{\beta^{*}}{\alpha}aa$ on the left. Therefore the first insertion must give $aa\exp\left(i\frac{2\mu}{H}\log\frac{\zeta_{l_{2}}}{\eta_{p}}\right)$ and its time integral resonates at $k\zeta_{l_{2}}=-\frac{2\mu\pm\omega}{H}$. Again due to the time-ordering requirement $\zeta_{l_{2}}>\zeta_{l_{2}-1}>\cdots>\zeta_{1}$ there is no resonance in the integral region.

The analysis above shows that the only way for the $\chi^{2}$’s to contribute order $\beta$ oscillating factor is to have only one $\chi^{2}$, i.e., bring down only one $H_{N+2}$ vertex. In the next subsection we compute the contributions to the $N$ point functions from $N$ three point vertices and show that the fermion case and boson case are the same except for a numerical factor coming from the helicity and anti-particle doubling. Since for all the other ways of using various $H_{I}$ terms (other than using only one $H_{N+2}$), only $a^{\dagger}a$ and $b^{\dagger}b$ contribute, it is clear that contributions to the $N$ point functions should be the same between the fermion case and the boson case in those situations. After that we compute the contributions from only using one $H_{N+2}$ vertex and show that the order $\beta$ contribution in the fermion case is further suppressed than in the boson case.

4.2 Order $|\beta|^{2}$ contributions

The contributions to the inflaton $N$ point function from bringing down $N$ three point vertices can be organized as [12]

[TABLE]

where

[TABLE]

The result in subsection 4.1 shows that only classical source terms $a^{\dagger}a$ and $b^{\dagger}b$ in (4) have contribution in this case and we will see in the following that expectation values of the product of those operators do not rely on their order. Therefore the expectation value of (4.17) factorizes as

[TABLE]

where the subscript “cl” means computing expectation value only using the classical source terms. The expectation value of the $\delta\phi$ commutator can be proven to be

[TABLE]

by mathematical induction. First note that from

[TABLE]

we have

[TABLE]

and it is a $c$-number. Then in the case of $N=2$ it is easy to show directly that

[TABLE]

Next we assume that (4.2) holds for $N-1$. By using the identity

[TABLE]

for any operator $O$ and $P_{i}$ we have

[TABLE]

where in the last equality we used the induction assumption that (4.2) holds up to $N-1$. This completes the proof. With this result (4.2) can be written as

[TABLE]

where the summation $\pi$ is over all permutations of $1,2,\cdots,N$. Each of the permutation can be renamed such that the integrand becomes $G_{{\bf k}_{1}}(0,\eta_{1})\cdots G_{{\bf k}_{N}}(0,\eta_{N})\langle\bar{\chi}\chi_{{\bf k}_{1}}(\eta_{1})\cdots\bar{\chi}\chi_{{\bf k}_{N}}(\eta_{N})\rangle_{\mathrm{cl}}$ and as a result it has a different time ordering from the one in the original integral. Summing over all the permutations covers the entire space spanned by $\eta_{1},\cdots,\eta_{N}$ and therefore the integral simplifies as

[TABLE]

Next we need to compute the $\chi$ correlator

[TABLE]

It is clear that whether each factor contributes an $a^{\dagger}a$ or $b^{\dagger}b$, the expectation value is always

[TABLE]

where $c$ can be either $a$ or $b$ and ${\bf k}_{T}=\sum_{i=1}^{N}{\bf k}_{i}$. Combining the $2^{N}$ of such terms we have

[TABLE]

where we have used

[TABLE]

Also note the extra factor of 2 that comes from the helicity summation, which the boson case does not have. Define

[TABLE]

then summing over different production events gives the contributions from the fermion three point vertices

[TABLE]

To be more concrete, in the rest of this paper we will consider a situation similar to that of [7] with an approximate discrete shift symmetry, with production events evenly spaced in proper time $t$, corresponding to conformal times

[TABLE]

with a constant $\gamma$ and frequency $w=|\dot{\phi}|/f$ derived from (2.2.

Using the boson three point vertex

[TABLE]

and (4.13) we can get the same contributions in the boson case

[TABLE]

except that it is smaller than the fermion case by a factor of $4$, which comes from the helicity doubling and anti-particle doubling.

In the next subsection we will show that

[TABLE]

From this fact and the discussion above one can expect that the contributions from interaction vertex configurations including higher point vertices would also be the same for the fermion case and boson case (except for the case where only one $H_{N+2}$ is brought down for $N$ point function as shown in section 4.3). This is because in both cases the $\delta\phi$ correlators and $\chi^{2}$ (or $\bar{\chi}\chi$) correlators give the same result and they are combined in the same way. While it is worth exploring the possible contributions from various insertion profiles using higher point vertices, we note that if one is working in the same parameter regime as outlined in [7] the contributions from $(N+2)$ three point fermion vertices to $N$ point functions will dominate over other contributions involving higher point vertices.

4.3 Order $\beta$ contributions

In this subsection we consider the contributions from only bringing down one $H_{N+2}$ vertex for $N$ point functions. As discussed in section 4.1 such contributions contain parts that are of order $\beta$. For both the fermion case and the boson case they are

[TABLE]

In the fermion case, the $(N+2)$ point vertices are

[TABLE]

and (4.38) becomes

[TABLE]

To compute the $\delta\phi$ commutator we first notice that in $\langle\delta\phi_{{\bf k}^{\prime}_{1}}(\eta_{1})\cdots\delta\phi_{{\bf k}^{\prime}_{N}}(\eta_{1})\delta\phi_{{\bf k}_{1}}(0)\cdots\delta\phi_{{\bf k}_{N}}(0)\rangle$ each of $\delta\phi_{{\bf k}^{\prime}_{i}}(\eta_{1})$ must contract with one $\delta\phi_{{\bf k}_{j}}(0)$ and therefore the former must be $a_{{\bf k}^{\prime}_{i}}u_{k^{\prime}_{i}}(\eta_{1})$ and the later $a^{\dagger}_{-{\bf k}_{j}}u^{*}_{k_{j}}(0)$. So the correlator is evaluated to be

[TABLE]

and the summation is again over all the possible permutation $\pi$ over $1,2,\cdots,N$. And similarly

[TABLE]

Plugging in

[TABLE]

we have the fermion $(N+2)$ point vertex contribution

[TABLE]

where $k_{T}=\sum_{i=1}^{N}k_{i}$. Similarly the boson $(N+2)$ point vertex

[TABLE]

has the following contribution to the $N$ point correlator

[TABLE]

A comparison of (4.3) and (4.3) shows the difference between the two cases could come from the difference in the two point functions of the produced particle fields. Indeed the $k^{\prime}/aM$ factor in (4) and the fact that ${\bf k}^{\prime}$ integral peaks around $a_{n}\sqrt{g|\dot{\phi}|}\ll aM$ indicate that the order $\beta$ contribution in the fermion case will be further suppressed by a factor of $\sqrt{g|\dot{\phi}|}/\mu$. We confirm this point in the following by evaluating (4.3) and (4.3).

Let us start by computing $\langle\bar{\chi}\chi_{{\bf k}_{T}}(\eta_{1})\rangle$. With $\langle a^{\dagger}_{r}({\bf k}^{\prime})a_{r}({\bf k}_{T}+{\bf k}^{\prime})\rangle=\langle b^{\dagger}_{r}(-{\bf k}_{T}-{\bf k}^{\prime})b_{r}(-{\bf k}^{\prime})\rangle\simeq|\beta_{k^{\prime}}|^{2}\delta({\bf k}_{T})$, $\langle a^{\dagger}_{r}({\bf k}^{\prime})b^{\dagger}_{r}(-{\bf k}_{T}-{\bf k}^{\prime})\rangle\simeq\beta^{*}_{k^{\prime}}\delta({\bf k}_{T})$ and $\langle b_{r}(-{\bf k}^{\prime})a_{r}({\bf k}_{T}+{\bf k}^{\prime})\simeq\beta_{k^{\prime}}\delta({\bf k}_{T})$ the mode expansion (4) becomes

[TABLE]

The integral over $\omega$ can be approximately evaluated as

[TABLE]

where we have defined $A_{n}(\eta)=1-a_{n}^{2}/a(\eta)^{2}$. Using this and the definition (4.31) of particle density $\bar{n}_{\chi}$ the source expectation value is evaluated to be

[TABLE]

Similarly the result for the boson case is

[TABLE]

From (4.49) and (4.50) it is already clear that the fermion order $\beta$ contribution is suppressed by an extra factor of $\frac{a_{n}}{a}\frac{(g|\dot{\phi}|)^{1/2}}{\mu}$.

To proceed we will focus on case (b) of the mass function which produces novel searching templates. For the case of $N=1$, it is clear that the following relation holds

[TABLE]

in the parameter regime outlined in [7], where $\theta_{N}$ is a multiple of $\pi/2$ that gives the right trigonometric function. We can prove this result for general $N$ using mathematical induction. Starting from (4.51), the $(N+1)$th derivative is

[TABLE]

The first term is smaller than the second by a factor of $g^{2}f^{2}/\mu^{2}$ so the left hand side is approximately $\frac{1}{M}g^{2}f^{1-N}\cos\left(\frac{\phi}{f}+\theta_{N+1}\right)$ which is also approximately

[TABLE]

This proves (4.51) holds for general $N$. Combining the result above we can compute the saddle point approximation for the contributions (4.3) and (4.3). The result for the fermion case is

[TABLE]

where we have defined the following quantities

[TABLE]

The summation of $n$ starts from the minimum $n$ (or the latest $\eta_{n}$) that satisfies $\eta_{n_{\mathrm{min}}}<\hat{\eta}_{\pm}$ and ends with the first/earliest production event with $n_{\mathrm{max}}$ 888A different definition of $n$ is used in Section 3 where the first production event is labeled as $n=1$. The definition is slightly changed here for notational simplicity and the meaning of $n$ should be clear within its context.. It is also clear from (4.3) that this result respects the approximate discrete shift symmetry $\log k_{T}\to\log k_{T}+2\pi H/\omega$ given that the shape $(k_{1}/k_{T},\cdots,k_{N}/k_{T})$ does not change, since $k_{T}$ only appears with $\eta_{n}$ as

[TABLE]

Such a shift can be absorbed into a relabeling of $n$. Similarly the result for the boson case is

[TABLE]

A comparison between (4.3) and (4.3) shows that besides the usual suppression comming from $|k_{T}\eta_{n_{\mathrm{min}}}|^{-3}\lesssim\left(\frac{2\mu\pm\omega}{H}\right)^{-3}$, there is an extra suppression factor of $(g|\dot{\phi}|)^{1/2}/\mu$ for the order $\beta$ contributions in the fermion case. The result of the fermion case is also further suppressed by a factor of

[TABLE]

for each production event $n$ in the summation. Since $-\eta_{n}$ gets exponentially larger than $-\eta_{n_{\mathrm{min}}}$ (which is greater than $-\hat{\eta}_{s}$) as $n$ increases, the contributions from the fermion case can be much smaller than those from the boson case, especially those generated by the early production events.

5 Summary and discussion

In this paper we worked out in detail the contributions to the inflaton $N$ point functions from the fermion production events and compared them to those from the boson counterparts. The expectation value of the fermion source operator $\bar{\chi}\chi$ and the boson source operator $\chi^{2}$ have the same classical part, which is of order $|\beta|^{2}$, while for the order $|\beta|$ part that represents quantum interference, the fermion case is relatively suppressed by an extra factor of $k^{\prime}/aM$. It is argued in section 4 that at the level of saddle point approximation the only order $|\beta|$ contributions to inflaton $N$ point functions come from those using one $N+2$ point interaction vertex. At order $|\beta|^{2}$ the contributions of the two types of particle are the same. A detailed evaluation at order $|\beta|$ then shows that besides the usual suppression factor $\left(\frac{2\mu\pm\omega}{H}\right)^{-3}$, the fermion contributions as compared to those of boson are further suppressed by a factor of $(g|\dot{\phi}|)^{1/2}/\mu$ and another coming from the separation between the production time and the resonant saddle point. This distinction between fermion and boson contributions leads to the interesting possibility of getting better understandings of what type of particles and how they interact with inflaton by fitting CMB data with searching templates from these models. Similarly, this could also hint on more difficult versions of this type of models, which consider the non-adiabatic production of strings. We leave these interesting topics to future works.

The successive production of particles is treated approximately as independent from each other in this paper. In regimes where the density of the produced particles is small and quickly gets diluted by inflation this is a valid approximation. In principle, however, there can be other situations where the density is no longer small and the effect of successive productions on each other is not negligible. For the boson case this is not too hard to solve since its wave equation can be solved exactly and the recursive Bogoliubov transformation can be obtained in a similar fashion as in [8, 9]. The fermion case, on the other hand, seems to be less straightforward. First the wave equation (3.12) has an extra imaginary part which makes it much more difficult to solve analytically. Also the fermion production is restricted by the Pauli exclusion principle while in the boson case the production can be enhanced by previously produced particles without limit. In the regime where $|\beta|$ is not parametrically smaller than 1, the Bogoliubov coefficients need to be solved recursively as

[TABLE]

with the initial condition $\alpha_{0}=1,\beta_{0}=0$. $T_{n}$ is the transformation matrix for the $n$th production event that ensures $|\alpha|^{2}+|\beta|^{2}=1$. It would be interesting to explore in this direction in future works, either analytically or numerically.

For general $N$ point functions it seems rather difficult, if possible, to derive generic close form formulae for the contributions from arbitrary interaction vertex profiles. In this work we found simple close form solutions for two extreme cases, where in one case the contributions are from $N$ three-point vertices and in the other from one $(N+2)$-point vertex. For the intermediate cases one needs to consider various possible profiles and potentially combine the methods used for the two extreme cases. Even though there is a well-defined window of parameters where at order $|\beta|^{2}$, contributions from $N$ three-point vertices dominate over those from other profiles, it would still be interesting to study the shapes and the relative scales of those other contributions [13].

acknowledgments

I would like to thank Eva Silverstein for extensive discussion and useful comments on a draft. I am also grateful to Moritz Münchmeyer for initial collaboration. I am supported in part by the National Science Foundation under grant PHY-0756174 and NSF PHY11-25915 and by the Department of Energy under contract DE-AC03-76SF00515.

Appendix A More on Bogoliubov coefficients

In this section we explore the possibility of deriving the exact recurrence relation between the coefficients $A_{n}$, $B_{n}$ after the $n$th production event and those before it. As noted in section 3.3, the wave equation (3.12) can be greatly simplified in the $|\eta|\to\infty$ limit, which is where the matching of WKB solutions happens. One can use the asymptotic behavior of the solution to (3.12) in this region to match with the two WKB solutions. Similar to the rotation method used in the main text, the same problem regarding the branch cut of $M$ arises with this method. In the $|\eta|\to\infty$ limit, the original wave equation (3.12) simplifies to

[TABLE]

where the sign for the imaginary term depends on which branch of $M$ is being considered. In order to make use of the solution to (A.1) for matching, it needs to be on the same branch which covers both of the $\eta\to-\infty$ and $\eta\to+\infty$ limits since otherwise one would have two different approximate differential equations, one with a positive imaginary term and another with a negative one. The solutions then would have to be connected via the small $|\eta|$ region where there is no branch cut. Connecting the two branches in this way is as hard as solving (3.12) exactly. Therefore we pick the same branch cut as shown in Figure 2 so that the imaginary part of (A.1) always has a negative sign.

With this definition, the solution to (A.1) is

[TABLE]

where $C_{1}$, $C_{2}$ are two constants and $D_{p}(z)$ is the parabolic cylinder function. $p_{1}$ and $p_{2}$ are defined as

[TABLE]

while all the other symbols are defined in the main text. In the meantime, the WKB solution is given in (3.13). Note again that there is only one set of $C_{1}$ and $C_{2}$ for both the $\eta\to-\infty$ and $\eta\to+\infty$ limit while $A$, $B$ are really $A_{n-1}$, $B_{n-1}$ before the $n$th production event and $A_{n}$, $B_{n}$ afterwards. The asymptotic behavior of (A.2) in the $\eta\to-\infty$ limit is [14]

[TABLE]

and in the $\eta\to+\infty$ limit

[TABLE]

where “$\sim$” denotes “asymptotic to”. They need to match with those of the WKB solutions, which in the $\eta\to-\infty$ limit is

[TABLE]

and in the $\eta\to+\infty$

[TABLE]

After some algebra, this matching gives the following recurrence relation

[TABLE]

As noted in section 3.3 $A$ and $B$ satisfy the normalization $|A|^{2}+|B|^{2}=1$. This serves as a good consistency check of the above result. A straightforward calculation yields

[TABLE]

This value is not exactly 1 given $|A_{n-1}|^{2}+|B_{n-1}|^{2}=1$. The reason is that the recurrence relation (A.8) and (A.9) are obtained by matching the WKB solutions with (A.2), which is the solution to the approximate differential equation obtained under the assumption $\tilde{\mu}^{2}\ll g^{2}a_{n}^{2}\dot{\phi}^{2}\eta^{2}$ and $g^{2}a_{n}^{4}\dot{\phi}^{2}\eta^{2}\simeq\nu_{0}^{2}=k^{2}+a_{n}^{2}\tilde{\mu}^{2}$. The two approximations imply that $k^{2}\simeq\nu_{0}^{2}$, which along with (A) leads to $|A_{n}|^{2}+|B_{n}|^{2}\simeq 1$. Therefore at this level of approximation, the result (A) is consistent with the normalization condition. We may understand the failure of obtaining the recurrence relation of coefficients $A_{n}$ and $B_{n}$ using this approximation in the following way. The asymptotic behavior of the exact solution to the original differential equation (3.12) should really be

[TABLE]

where now $\tilde{C}_{1}(\eta)$ and $\tilde{C}_{2}(\eta)$ are functions of $\eta$ and they are asymptotically constants. The solution to the approximated differential equation (A.2) assumes that the constants are universal on the lower half plane. However $\tilde{C}_{1}(\eta)$ and $\tilde{C}_{2}(\eta)$ can have more complicated Stokes phenomenon and have different constant values in different sections. We expect that this omission of Stokes phenomenon to be the reason of failing to recover the recurrence relation of $A_{n}$ and $B_{n}$’s.