Numerical study of inflationary preheating with arbitrary power-law potential and a realization of curvaton mechanism

TL;DR

This paper develops a numerical simulation to study inflationary preheating with arbitrary power-law potentials, demonstrating how entropy field fluctuations can dominate curvature perturbations, thus realizing a curvaton-like mechanism.

Contribution

It introduces a numerical method for analyzing preheating in models with power-law inflaton potentials and demonstrates a curvaton mechanism via entropy fluctuations during preheating.

Findings

Numerical analysis of preheating with arbitrary power-law potentials.

Entropy field fluctuations can dominate curvature perturbations.

Validates the approach with axion-monodromy inflation model.

Abstract

During inflationary preheating, the energy stored in the inflaton field is rapidly converted into excitations of other entropy fields. This stage is characterized by exponential particle production due to parametric resonance and is notoriously difficult to analyze using analytic methods. We develop a detailed numerical simulation to investigate inflationary preheating when the potential of the inflaton is a power-law function with arbitrary power index. To achieve a successful graceful exit from a primordial inflationary phase to a smooth, oscillatory phase during preheating, we assume the inflaton potential reduces to a quadratic function in the infrared regime, which may be regarded as a mass term at low energy scales. With this simplification, our numerical method may be applied to unconventional chaotic inflation models. To demonstrate its validity, we numerically analyze the preheating stage of axion-monodromy inflation which is inspired by string theory. By performing perturbation analyses, our result shows that the power spectrum of primordial curvature perturbation can be dominated by fluctuations of entropy field rather than those of inflaton, which can be regarded as a particular realization of the curvaton mechanism through a preheating process.

Figures

Tables

Table 1. Table 1: Model parameters and initial conditions during inflation.

\toprule$p$	${\varphi }_{\text{I}}\left(M_p\right)$	$\lambda$	${\mu }_{\text{max}}$	$g$	$q$
$1/3$	-7.9	$3.7\times {10}^{-10}$	0.043	0.0868	1.48968
$1/2$	-9.2	$3.7\times {10}^{-10}$	0.05	0.075864	2.17388
$2/3$	-10.2	$3.2\times {10}^{-10}$	0.058	0.057	2.8855
$1.0$	-12.1	$2.0\times {10}^{-10}$	0.078	0.0288	4.92016
\toprule$2.0$	-17.4	$1.7\times {10}^{-11}$	0	$0.003291$	$\setminus$

Equations

$$d{s}^2=-d{t}^2+a^{2}d{\mathrm{x}}^2,$$

$$H^2=\frac{1}{3M_p^2}\rho ,\frac{\ddot{a}}{a}=-\frac{1}{6M_p^2}(\rho +3p),$$

$$\mathcal{L}=-\frac{1}{2}{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi -\frac{1}{2}{\partial }_{\mu }\chi {\partial }^{\mu }\chi -V(\varphi ,\chi ),$$

$$\displaystyle V(\phi,\chi)=\lambda M_{p}^{4-p}\Big{[}\big{(}\phi^{2}+\mu^{2}\big{)}^{p/2}-\mu^{p}\Big{]}+\frac{1}{2}g^{2}\phi^{2}\chi^{2}\,.$$

$$\rho =\frac{1}{2}\sum _I{\dot{\phi }}_I^2+V,P=\frac{1}{2}\sum _I{\dot{\phi }}_I^2-V,$$

$$\ddot{\varphi }+3H\dot{\varphi }-\frac{{\nabla }^2}{a^2}\varphi +V_{,\varphi }=0,$$

$$\ddot{\chi }+3H\dot{\chi }-\frac{{\nabla }^2}{a^2}\chi +V_{,\chi }=0.$$

$$Q_I\equiv \delta {\phi }_I+\frac{{\bar{\phi }}_I^{\prime }}{\mathcal{H}}\psi ,$$

$${\ddot{Q}}_I$$

$$\displaystyle\sum_{J}\bigg{[}V_{,IJ}-\frac{1}{M_{p}^{2}a^{3}}\frac{d}{dt}\Big{(}\frac{a^{3}}{H}\dot{\varphi}_{I}\dot{\varphi}_{J}\Big{)}\bigg{]}Q_{J}=0~{},$$

$$\epsilon=\frac{M_{p}^{2}}{2}\Big{[}\frac{V(\phi,0)^{\prime}}{V(\phi,0)}\Big{]}^{2}~{},$$

$$V\approx \lambda M_p^{4-p}[{\varphi }^p-{\mu }^p],$$

$${\varphi }_{ϵ=1}\ge 10\mu .$$

$$m_{{\varphi }_{\text{eff}}}^2=2\lambda M_p^{4-p}{\varphi }^{p-2}.$$

$$V(\varphi ,\chi )=\lambda \frac{p}{2}\frac{{\mu }^{p-2}{\varphi }^2}{M_p^{p-4}}+\frac{1}{2}{g}^2{\chi }^2{\varphi }^2+O({\varphi }^4).$$

$$m_{{\varphi }_{\text{eff}}}^2=p\lambda M_p^{4-p}{\mu }^{p-2},$$

$$\ddot{\varphi }+\frac{2}{t}\dot{\varphi }+m_{{\varphi }_{\text{eff}}}^2\varphi =0.$$

$$\varphi =\frac{{\Phi }_0}{m_{{\varphi }_{\text{eff}}}t}\sin (m_{{\varphi }_{\text{eff}}}t),$$

$$m_{{\varphi }_{\text{eff}}}^2=p\lambda M_p^{4-p}{\mu }^{p-2}>m_{{\chi }_{\text{eff}}}^2=g^2{\Phi }^2,$$

$$g<\sqrt{100p\lambda M_p^{4-p}{\mu }^{p-4}}.$$

$${\ddot{X}}_k+{\omega }_k^2{X}_k=0,$$

$${\omega }^2=\frac{k^2}{a^2}+\frac{g^2{\Phi }_0^2}{m_{{\varphi }_{\text{eff}}}^2{t}^2}-\frac{g^2{\Phi }_0^2}{m_{{\varphi }_{\text{eff}}}^2{t}^2}\cos (2m_{{\varphi }_{\text{eff}}}t).$$

$$X_k^{\prime \prime }+[A_k-2q\cos (2z)]X_k=0,$$

$$n_k=\frac{{\omega }_k}{2}\left(\frac{|\dot{X_k}{|}^2}{{\omega }_k^2}+|X_k{|}^2\right)-\frac{1}{2},$$

$$\mathcal{R}=H\frac{Q_{\varphi }\dot{\varphi }+Q_{\chi }\dot{\chi }}{{\dot{\varphi }}^2+{\dot{\chi }}^2},$$

$$y_{n+1}=y_n+\frac{h}{6}[K_1+2K_2+2K_3+K_4]$$

$$K_1=f(x_n,y_n)$$

$$K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_1)$$

$$K_3=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_2)$$

$$K_4=f(x_n+h,y_n+h{K}_3)$$

$$N_{\text{inf}}$$

$$P_R$$

$$(\varphi ,\dot{\varphi },Q_{\varphi },{\dot{Q}}_{\varphi },Q_{\chi },{\dot{Q}}_{\chi }{)}_{|\text{end}}$$

$$\ddot{X_k}$$

Full text

cosmological parameters, inflation, theory

1 Introduction

Reheating is a hypothetical mechanism used to generate a thermal phase of hot big bang expansion from an initially cold and empty primordial Universe. During the reheating phase, energy is quickly transferred from a primordial matter field, namely, the inflaton scalar field, to other matter fields through interaction terms in the Lagrangian. For inflationary cosmology, this process occurs after the stage of slow-roll inflation, but earlier than the time of thermal equilibrium. The energy transfer from the inflaton to matter field was initially analyzed using first-order perturbation theory and discussed in terms of decaying inflaton quanta into Standard Model particles (Abbott, 1982; Dolgov, 1982; Albrecht, 1982).

It was realized that the coherence of the inflaton condensate at the beginning of the reheating phase can play an important role in the analyses with nonlinear effects (Traschen, 1990). Matter generation in the inflaton condensate was studied in a semi-classical analysis involving non-perturbative effects in (Traschen, 1990). It was found that a parametric resonance instability plays a crucial role. Parametric resonance after inflation was further studied in (Kofman, 1994; Shtanov, 1994) and the particle production process was analyzed in (Kofman, 1997). In recent years, this topic has been extensively studied in the literature in the framework of inflationary cosmology (Boyanovsky, 1994; Baacke, 1996; Greene, 1997; Felder, 2000; Cormier, 2001; Shuhmaher, 2005; Dufaux, 2006; Abolhasani, 2009) as well as bounce cosmology (Cai, 2011a, 2013; Quintin, 2014; de Haro, 2015). We refer to (Bassett, 2005; Allahverdi, 2010; Amin, 2014) for comprehensive reviews.

While the reheating process can provide a dynamical mechanism for the origin of entropy and particles observed in our universe, details remain somewhat heuristic due to the existence of various uncertainties. For instance, the evolution of the equation-of-state (EoS) parameter strongly depends on the specific construction (Podolsky, 2005). Consequently, one way of probing information about reheating is to track the expansion of the Universe between the moment CMB scales cross the Hubble radius during inflation and the time they re-enter, which connects the energy scale of inflation with the reheating temperature directly. This method was extensively studied in the literature, see for instance, (Liddle, 2003; Martin, 2006, 2010; Adshead, 2010; Mielczarek, 2010; Easther, 2011; Dai, 2014; Martin, 2014; Cai, 2015; Cook, 2015; Domcke, 2015; Drewes, 2015; Lozanov, 2016).

Nevertheless, such constraints on the reheating temperature strongly rely on measurements of primordial curvature perturbations and relic gravitational waves, from various high precision experiments of cosmic microwave background (CMB) photons and their polarizations (Aghanim, 2018; Ade, 2015a). So far these studies were based on the theoretical predictions of particularly chosen inflation models (Akrami, 2018; Ade, 2015b). Moreover, it was observed that some nonlinear processes, namely parametric resonance effects, which could yield a potential amplification on the amplitudes of primordial density perturbations (Taruya, 1997; Bassett, 1998, 1999a, 1999b; Finelli, 2000; Tsujikawa, 2002; Chambers, 2007; Brandenberger, 2008; Bond, 2009; Bethke, 2013; Moghaddam, 2014; McDonough, 2016; Svendsen, 2016; Moghaddam, 2016; Graef, 2017; Gu, 2018), have not yet been taken into account in a self consistent way. In particular, if these amplifications have an impact on curvature perturbations due to a conversion from iso-curvature modes, the theoretical predictions for inflation models may be strongly altered or even completely destroyed (Finelli, 2000; Brandenberger, 2008; Moghaddam, 2014; McDonough, 2016). Therefore, it is necessary to numerically examine whether these amplifications occur in specific inflationary models and to further revisit the corresponding predictions for CMB observations.

For those lattice codes of preheating such as LATTICEEASY (Felder, 2000), DEFROST (Frolov, 2008), PSpectRe (Easther, 2010) and HLattice (Huang, 2011), the cosmological system has been treated as a fixed comoving box. However, in the present study we would like to trace the evolutions of cosmological perturbations at super-Hubble scales so that we could examine whether these modes would continue the evolutions after the Hubble radius crossing. Accordingly, in this article we will use the C++ computer language to develop a new program code to perform the numerical simulation of the inflationary preheating phase, in which the potential of the inflaton is taken to be a power-law function with an arbitrary power index.

To ensure that the oscillatory behavior of the inflaton field around the bottom is continuous and smooth, we introduce a cutoff for the potential in the infrared regime, which behaves as a quadratic potential at low energy scales.The entire background evolution then becomes numerically trackable, even if the underlying inflation model has some unconventional potential forms with singular field derivatives at high energy scales. This treatment is particularly appropriate for a class of the so-called axion-monodromy inflation models inspired by string theory. Here we consider the ordinary chaotic inflation models with quadratic and quartic potentials, and then study particle production in models of power-law potentials with power index being $1$, $2/3$ and $1/2$, respectively.

The present article is organized as follows. In Section 2, we introduce the construction of the inflaton potential and introduce the background equations applied in the numerical simulations. In Section 3, we give a detailed description of the numerical method we developed to calculate the preheating process. In Section 4, we apply the code to several inflationary models to demonstrate the efficiency of the numerical method. In the same Section we firstly consider well-known cases with power index of the potential $p=2$, to demonstrate the validity of the method; and then, investigate some nonconventional models with $p=1$, $2/3$, $1/2$ and $1/3$, which are inspired by models of axion monodromy inflation. Finally, we summarize our results with a discussion in Section 6.

Throughout the article we adopt metric signature $(-,+,+,+)$, use natural units with $c=\hbar=1$ and use the reduced Planck mass $M_{p}=1/{\sqrt{8\pi G}}$.

2 Model construction

We begin with a homogeneous and isotropic universe, described by the (observationally favored, spatially flat) Friedmann-Lemaître-Robertson-Walker (FLRW) metric:

[TABLE]

where $a$ is the scale factor. Accordingly, the field equations of General Relativity yield the Friedmann equations for the background

[TABLE]

where $H\equiv\dot{a}/a$ is the Hubble parameter characterizing the expansion rate of the Universe, and, the dot represents differentiation with respect to cosmic time $t$. Here $\rho$ and $p$ denote energy density and pressure, respectively.

As matter sources, we take two canonical scalar fields, one responsible for driving a period of slow roll inflation at high energy scale, the other as an entropy field. The Lagrangian is

[TABLE]

with potential including interaction term,

[TABLE]

In the above, $\phi$ denotes the inflaton field and $\chi$ is an entropy field which is assumed to be much lighter than the inflaton, and only excited during the preheating phase. The coefficients $\lambda$ and $g$ are the coupling constants. Here, $p$ is the model-dependent power-law index we discuss in later sections (see also Figure 1 for a sketch). In addition, we have introduced a very small mass term $\mu$ in order to ensure the potential behaves smoothly in the infrared regime (in particular, for circumstances when the values of $p$ are smaller than unity).

The background energy density and pressure are,

[TABLE]

where we have ignored spatial gradient terms, $\varphi_{I}$ represent the inflaton field $\phi$ and the entropy field $\chi$, respectively. The equations of motion for two fields can be derived by varying the Lagrangian with respect to the field variables, respectively:

[TABLE]

The field perturbation in spatially flat gauge is the Mukhanov-Sasaki variable

[TABLE]

which is used here for analyzing the perturbation, including the metric perturbation $\psi$, where $\mathcal{H}$ is the comoving Hubble parameter. The equations of motion of the perturbation fields are

[TABLE]

where the comma denotes the derivative with respect to the field $\phi$ or $\chi$. After operating the Fourier transform to Equation (2), the partial derivative with respect to space coordinate, $-\nabla^{2}$, changes to the comoving wavenumber $k^{2}$ in Fourier space. The Partial Differential Equations (PDE) then become Ordinary Differential Equations (ODE), simplifying computations.

During inflation, the contribution of the $\chi$ field is negligible since the $\chi$ energy density is completely diluted by the exponential expansion of the universe without fine-tuning parameters. Consequently, it is natural to take $\chi_{I}=\dot{\chi}_{I}=0$ as the initial conditions at the beginning of inflation. Following the traditional slow-roll paradigm, we introduce the slow-roll parameter

[TABLE]

and inflation ends when $\epsilon=1$. We set the initial conditions of the perturbations to be the Bunch-Davies vacuum. As the power index $p$ of the inflaton potential changes, the initial conditions of $\phi$ and $\dot{\phi}$ must change in order to achieve a sufficient number of e-foldings, $N\approx 50\sim 60$. Furthermore, in order to ensure the potential behaves like $\phi^{p}$

[TABLE]

during the inflation, we require that $|\phi|\gg 10\mu$, restricting

[TABLE]

Thus, for different $p$, we can yield different upper limits for $\mu$. Note that, during inflation, the effective mass of the inflaton is given by

[TABLE]

After inflation, the universe is dynamically driven by the oscillations of $\phi$ about the vacuum, which is called preheating. The energy transfers from the inflaton to the matter field, exciting relativistic particles. In order to transfer the energy, the effective mass of the inflaton must on average continually decrease. However, for those string theory inspired axion-monodromy models where $p<2$, the inflaton becomes heavier and heavier, yielding an inefficient preheating process. In our model, when $|\phi|$ is smaller than $0.1\mu$, the potential $V$ becomes

[TABLE]

The effective mass of inflaton is

[TABLE]

is constant instead of becoming heavier. The graceful exit from $\phi^{p}$ to $\phi^{2}$ results in the constant effective mass of the $\phi$ field, and efficient particle production during the preheating phase. During this time, the EoS parameter is as during matter-domination, so that $a\propto t^{2/3}$, $H=\frac{\dot{a}}{a}=\frac{2}{3t}$, and the equation of motion of the inflaton becomes

[TABLE]

This equation of motion is solved by

[TABLE]

where $\Phi_{0}$ is the coefficient when the inflaton enters the $V\approx\phi^{2}$ regime. Since the coupling is $1/2g^{2}\phi^{2}\chi^{2}$ in the model, we have a kind of scattering process, so the matter field effective mass $m_{\chi_{\text{eff}}}^{2}$ should be lighter than the effective mass of the inflaton $m_{\phi_{\text{eff}}}^{2}$. Since during the $\phi^{2}$ preheating period, $\phi<0.1\mu$, so $\Phi<0.1\mu$, and

[TABLE]

we choose the upper limit of $\Phi$, so that

[TABLE]

The coupling parameter $g$ should be carefully chosen for parametrically resonance. Since the potential is proportional to $\phi^{2}$ during preheating, the condition for parametric resonance is quite similar to that of pioneer work on $\phi^{2}$ inflation (Kofman, 1997). We rescale the matter field $\chi_{k}$ with $X_{k}=a^{3/2}\chi_{k}$, and rewrite the equation of motion:

[TABLE]

where $\omega^{2}_{k}=\frac{k^{2}}{a^{2}}+g^{2}\phi^{2}-\frac{9}{4}H^{2}-\frac{3}{2}\dot{H}$, the last two terms vanish as the EoS parameter approaches zero, and

[TABLE]

In the standard Mathieu function form (McLachlan, 1947) we have,

[TABLE]

where $z=m_{\phi_{\text{eff}}}t$, a prime denotes the derivate with respect to $z$, and $A_{k}=\frac{k^{2}}{m_{\phi_{\text{eff}}}^{2}a^{2}}+\frac{g^{2}\Phi_{0}^{2}}{m_{\phi_{\text{eff}}}^{4}t^{2}}$, $q=\frac{g^{2}\Phi_{0}^{2}}{4m_{\phi_{\text{eff}}}^{4}t^{2}}$. We use $q=\frac{g^{2}\Phi^{2}}{4m_{\phi_{\text{eff}}}^{2}}>1$ to determine the broad resonance, where $\Phi$ is the initial amplitude. We see that when $\mu$ decreases, the effective resonance area decreases, thus the lower limit for $g$, $g>2\Phi m_{\phi_{\text{eff}}}$, increases, indicating stronger coupling.

In order to numerically investigate the evolution of the matter field $\chi$ in this phase, we introduce a small, homogenous and isotropic perturbation $\delta\chi$ as the initial value for $\chi$ and $d\chi$ at the beginning of preheating. The background homogeneous part of the $\chi$ field, though very small, can be interpreted as the condensate matter field. The inhomogeneous part, $\chi_{k}$, can, in turn, be interpreted as the matter particles with different comoving wavenumbers excited by the scattering with the inflaton field. It is important to distinguish $\chi_{k}$, and the spatially flat field perturbation, the Mukhanov-Sasaki variable $Q_{\chi}$, which includes the metric perturbation in our code. The occupying number $n_{k}$ of the $\chi$ field is defined by $X_{k}=a^{3/2}\chi_{k}$ instead of $Q_{\chi}$,

[TABLE]

and the comoving curvature perturbation

[TABLE]

where the curvature perturbation mainly depends on $Q_{\phi}$ and $Q_{\chi}$. Previously, people didn’t include the metric perturbation when they try to consider the change of curvature perturbation caused by matter field. Instead, they introduce the back reaction afterward. However, due to the improving ability to calculate, we are able to consider the metric perturbation at the beginning, which highly improves the precision.

3 Description of the numerical method

In this section, we present a detailed description of the numerical method we apply to the preheating analysis. Previous codes only considered the field fluctuations instead of including metric perturbations, which leads to an ambiguous definition between the curvature perturbation and the occupation number. In our code, we separate the field perturbation and the Mukhanov-Sasaki variable, making the definition clear. In addition, the complicated non-linear process in lattice simulations is unnecessary, since the linear perturbation theory of two fields would, in general, modify the power spectrum of the curvature perturbation. We write an MPI-version to run the evolution for different wave vector in parallelization and get the power-spectrum all at once. We write the code for inflationary preheating in the C++ computer language and hope to revisit the inflationary predictions of the $(n_{s}-r)$ contour plot. The numerical algorithm adopted in our code is based on the 4th order Runge-Kutta approach, which is comparably precise and fast.

We introduce the application of 4-th order Runge-Kutta method in Section 3.1 and discuss the structure of the code together with the initial values in Section 3.2.

3.1 4-th Order Runge-Kutta Method

Here we give a brief introduction to 4-th order Runge-Kutta method for our ODE set, (see Equations (2), (2) and (2)). For a differential equation $\frac{dy}{dx}=f(x,y)$, the next step of $y_{n+1}$ is

[TABLE]

where $h$ is the time step. Using this method the numerical computation can converge very fast and remain high precision.

3.2 The structure of the code

We first calculate the inflationary part. The matter field $\chi$ has almost diluted away so it is appropriate to set $\chi=\dot{\chi}\equiv 0$ during inflation. $\phi_{I}$ and $\dot{\phi}_{I}$ is determined by $N_{inf}$, the e-folding number, to ensure a theoretically viable phase of inflation. The Mukhanov-Sasaki variable $Q_{\phi}$ and $Q_{\chi}$ are introduced at the beginning of the inflation, where the Bunch-Davies vacuum gives the initial values for these two variables. We refer to Appendix B for the detailed description.

With the strong assumption that cosmological perturbations are not altered throughout the whole preheating phase, one can directly restrict the value of the model parameter $\lambda$ by applying observational data from CMB experiments. Namely, one can (roughly) set the e-folding number and the amplitude of power spectrum of primordial curvature perturbations to be,

[TABLE]

where $\epsilon$ is the slow-roll parameter. When $\epsilon\sim\mathcal{O}(1)$, inflation ends, and $\phi_{\text{end}}$ is used for determining $\mu$ by setting the upper limit $\mu_{\text{max}}=0.1|\phi_{\text{end}}|$.

After inflation, we take the ending values of

[TABLE]

as the initial conditions of the preheating phase $(\phi~{},\dot{\phi}~{},Q_{\phi}~{},\dot{Q}_{\phi}~{},Q_{\chi}~{},\dot{Q}_{\chi})_{|\text{p}}$. For the $\chi$ field, we introduce the Bunch-Davies Vacuum as the initial value $\chi_{\text{p}}(k)$ during the preheating period (See Appendix A). The reason we do not let $\chi$ evolve through the entire process is due to its complete dilution due to the hierarchy between $\chi$ and $\phi$ during inflation. We determine the coupling constant $g$ by the semi-analytical method.

During preheating, the inflaton $\phi$ oscillates near $\phi=0$ with the potential $V\sim\phi^{2}$, thus we can approximate $\phi\sim\Phi\sin{(m_{\phi_{\text{eff}}}t)}$, where $\Phi$ is the amplitude and $m_{\phi_{\text{eff}}}=\sqrt{p\lambda M_{p}^{4-p}\mu^{p-2}}$ is the effective mass during preheating. We take $\Phi\approx 0.1\mu$ for simplicity, where $0.1\mu$ represents the resonance area. The equation of motion of $X_{k}(t)=a^{3/2}(t)\chi_{k}(t)$ is,

[TABLE]

where we have neglected the $(-\frac{9}{4}H^{2}-\frac{3}{2}\dot{H})$ term since the background is matter-dominated. One can straightforwardly rewrite Equation (25) in the form of the Mathieu equation as follows,

[TABLE]

where the prime denotes the derivative with respect to $z=m_{\phi_{\text{eff}}}t$. The parameter $q$ is often used to distinguish the narrow resonance ($q\ll 1$) and the broad one ($q\gg 1$). The narrow resonance is inefficient and only suitable for narrowband near $k\approx m_{\phi_{\text{eff}}}$. In the present article, we analyze $q\sim 1.5$ which induces broad resonance for all long wavelength modes. These parameters are capable of efficiently producing matter particles, changing $P_{\mathcal{R}}$, the power spectrum of curvature perturbation, by adding entropy perturbation and also restricting $P_{\mathcal{R}}$ within the observational limit.

The moment when the energy densities of the two fields are equivalent, i.e., $\rho_{\phi}\simeq\rho_{\chi}$ denotes the end of preheating. By comparing $\frac{Q_{\chi|\text{end}}}{Q_{\chi|\text{p}}}$, one can determine whether efficient preheating occurs. Practically speaking, we run the code until the order of $Q_{\chi}$ stops rising and calculate the curvature perturbation there.

Since we focus on the linear perturbation, we use MPI to run parallel computation for different comoving wavenumbers. Choosing $10^{-4}~{}\text{Mpc}^{-1}<k<10^{-1}~{}\text{Mpc}^{-1}$ and converting them into Planck units, we get the range for comoving wavenumber as follows,

[TABLE]

However, the observational range of $k$ is based on the fact that our current scale factor $a_{0}=1$. In our code, we set $a_{\text{I}}=1$ at the beginning of inflation for calculation convenience. To make $k/a$ in the same order with the observation, we need to adjust the range of $k$ in our code. Supposing the preheating temperature $T_{\text{re}}\simeq 10^{13}~{}\text{GeV}$ and from the relation \frac{T_{\text{re}}}{T_{0}}=\big{(}\frac{43}{11g_{s,\text{re}}}\big{)}^{1/3}\frac{a_{0}}{a_{\text{re}}} (Dai, 2014), we can estimate the ratio of $a_{0}$ to $a_{\text{re}}$ by $10^{13}~{}\text{GeV}/10^{-4}~{}\text{eV}\sim 10^{26}$. The initial condition should be chosen to ensure $k/aH\gg 1$ at the beginning of inflation so that every wave vector is inside the Hubble radius. That is,

[TABLE]

where $N_{k}$ is the e-folding number during inflation, $N_{\text{re}}$ is the e-folding number from the start of preheating to today and $H_{\text{ini}}$ is the initial Hubble parameter. Since we have already have $N_{\text{re}}=\ln(10^{26})$, we can use the above conditions to get the right range of the wavenumbers for different case of $p$.

Generally speaking, the scale factor will bump $e^{70}\sim 10^{30}$ during inflation no matter which $p$ we choose. Together with the additional $26$ orders after the inflation, the total e-folding number is about $55$. It is proper to choose the wavenumber range for our code to be $(10^{-5}~{}l_{p}^{-1},10^{-2}~{}l_{p}^{-1})$. In the following sections, we show the pivot scale, $k=10^{-3}~{}l_{p}^{-1}$, for illustration.

To test precision, we use two different scale factor evolution methods. One is $\dot{a}/a=H=\sqrt{\rho/(3M_{p}^{2})}$, and the other $\ddot{a}/a=-(\rho+3p)/(6M_{p}^{2})$. They give equivalent results for the following sections.

4 Applications to specific examples of the background field evolution

We now apply the numerical method introduced above to our potentials with different polynomial power laws. We begin with the classic preheating model (eg, $p=2$), which has been analyzed thoroughly (Kofman, 1997). From the plot, we will show that our method agrees with previous results. We then consider other cases (e.g, $p=1,~{}2/3,~{}1/2,~{}1/3$ ) to analyze the process of preheating for small field inflation.

4.1 Regular chaotic inflation with $p=2$

Following (Bassett, 2005), we plot the production of particles in for $p=2$. When choosing the proper initial conditions, the preheating period is efficient, producing enough matter particle. Our result matches the work of (Kofman, 1997) while they chose the mode $k=0$, which can be understood as producing the condensate $\chi$ particles.

From Figure 2, we see that for the case of $p=2$, the $\chi$ field can experience the resonance excitation in a chosen parameter band, and the physical particle number increases when $\chi$ is resonating. However, potentials with $\phi^{2}$ or $\phi^{4}$ are now significantly disfavored by the CMB data measured by the Planck satellite (Akrami, 2018), motivating consideration of potentials with $p<1$. The key point is to see whether there will be enough particle production and whether the curvature perturbation predictions are significantly changed.

4.2 Nonconventional models with $p=1,2/3,1/2,1/3$

Now we apply the code to analyze the preheating process when $p<1$. As has been discussed previously, we ensure a smooth transition from $\phi^{p}$ to $\phi^{2}$ by fixing the upper limit of $\mu$. As a result, we can make the approximation $\phi\sim\sin{m_{\phi_{\text{eff}}}t}$, which simplifies the analytical calculation of $\chi$ fields. It can be seen in the evolutions of EoS, which are shown in Figure 3. Different potentials ultimately reduce to the case of ($p=2$, $w=0$) at the end.

It is well known that for $\phi^{2}$, the particle production is sufficient. As our potentials always come back to $\phi^{2}$, the particle production should also be sufficient. Indeed, our numerical results indicate the $p<1$ cases efficiently produce particles.

Figure 4 shows the $\chi$ field and occupying number in the case of $p=1/2$. We can see that at the early time, the comoving particle number density increases while $\chi$ is resonating, and then remains unchanged in the comoving universe. The analysis recovers the results found in the $\phi^{2}$ case, as expected.

5 Applications to specific examples of the perturbation and power spectrum

One of the most powerful observational results of the CMB data is the $(n_{s}-r)$ plot which almost rules out $\phi^{2}$ and $\phi^{4}$ models. However, it is unfair to avoid considering the significant preheating process when approximating the spectral index $n_{s}$. The effect that the $\chi$ field acting as an entropy field to unfreeze the curvature perturbation on super-Hubble scales cannot be neglected. In this section, we analyze how the second field influences the curvature perturbation and how it may correct the predictions on the power spectrum of single field inflation.

5.1 Regular chaotic inflation with $p=2$

We begin with an analysis of the standard quadratic model, $V\sim\phi^{2}$.

Figure 5 shows that the power spectrum of the double field curvature perturbation is also nearly scale-invariant, like that of the single field curvature perturbation; however, due to particle production, the magnitude of the curvature perturbation rises by nearly 3 orders of magnitude. In other parameter bands, this amplification will be even more significant and uncontrolled. This provides a new way to constrain the parameter space. Moreover, the analysis demonstrates that the curvature perturbation is not frozen after exiting the Hubble horizon and will keep evolving driven by entropy perturbations. It is therefore inappropriate to use the single field curvature perturbation when extra matter fields are introduced to make efficient particle production during preheating. One ought to consider a combination of both fields when analyzing the CMB information from inflationary predictions, which will be addressed in the follow-up project.

It is interesting to note that, for a sizeable parameter space as illustrated in the above numerical estimate, the power spectrum of primordial curvature perturbation remains nearly scale-invariant, but the amplitude can be altered significantly due to the contribution of the fluctuations of the entropy field. This scenario coincides with the so-called curvaton mechanism (Mollerach, 1989; Linde, 1996; Lyth, 2001; Moroi, 2001; Enqvist, 2001), which requires at least two primordial fields with one being the inflaton and the other the curvaton field. When the curvaton is subdominant it only provides entropy perturbations during inflation (Lyth, 2002), and afterward, these entropy perturbations can be converted to curvature perturbations through various dynamical processes (Sasaki, 2006; Easson, 2007; Li, 2008; Huang, 2008a, b; Cai, 2009b, a; Kobayashi, 2009; Zhang, 2009; Gong, 2009; Cai, 2010, 2011b; Alexander, 2014; Addazi, 2016) so that the curvaton field can start to dominate the universe. For traditional curvaton models, the universe has to enter the standard thermal history after the curvaton decays. However, in our scenario, the curvaton mechanism itself is achieved by the preheating process.

5.2 Nonconventional models with $p=1,2/3,1/2,1/3$

Figure 6 shows that the double fields curvature perturbation is not frozen after exiting the Hubble horizon, unlike the single field curvature perturbation. The peaks correspond to the particle production. The curvature perturbation has an overall increase on average.

Figure 7 is the case of $p=2/3$ showing that the second field makes such a significant contribution to lift the power spectrum that it cannot simply be neglected. At the late time of preheating, $\phi^{2}$ dominates so the power spectrum has similar results as the $p=2$ model. The matter field perturbation contributes significantly to the overall curvature perturbation. In general, the power spectrum of the double field curvature perturbation for $p=1$, $2/3$, $1/2$ and $1/3$ all have an amplification of several magnitudes compared to that of single field inflationary models. In this regard, our present study provides a concrete realization of the curvaton mechanism via the preheating process for a generic inflation model with an arbitrary power-law potential.

Our results demonstrate the significance of possible matter field contributions to the curvature perturbation. The traditional $(n_{s}-r)$ plot used to constrain inflationary models do not take into account the preheating process (largely due to uncertainties in the details of preheating mechanics). However, our analysis demonstrates that preheating may strongly affect the curvature perturbation, radically altering the predictions of single field inflation models.

6 Conclusion and Discussions

In the present article, we developed a code based on the C++ computer language and MPI parallelization, aiming at the numerical simulation of the preheating process after inflation. Users can easily change the potential for the inflaton $\phi$ and matter $\chi$ fields. Here we have mainly focused on the non-perturbative particle producing process and its influences on the curvature perturbation.

We showed that the case of $p=2$ corresponds to the $\phi^{2}$ model studied in the literature, as a viability test for the code. We then extended the analysis to the small $p$ case where the extra term $\mu$ in our potential smooths the singularity. We find these $\phi^{p}$ models ultimately transition to the familiar $\phi^{2}$ model so that after choosing the proper parameters, particle production remains efficient. Because of this, our method can be used to study preheating in axion-monodromy inflationary models, by eliminating the singularity at the minimum of the potential, providing a smooth transition from the axion-monodromy like potential inflation to an ordinary $\phi^{2}$ potential, making the preheating process as efficient as the quadratic case.

Taking into account the effect of the matter field fluctuations, the curvature perturbation can keep evolving at super-Hubble scales, which conflicts with the single field curvature prediction. The extra field can introduce an entropy perturbation that acts as a source of the curvature perturbation, amplifying the power spectrum. The amplification can be large, leading to additional constraints on the parameter space. This can be regarded as a particular realization of the curvaton mechanism, in which the primordial power spectrum is dominated by the fluctuations of the curvaton field instead of those of the inflaton field. However, unlike the traditional curvaton paradigm which requires a process of curvaton decay, in our model, the preheating phase itself serves this purpose. Our study demonstrates that the ordinary predictions for the $(n_{s}-r)$ plot which are based on the single field inflationary phase may fail to capture crucial information due to the unavoidable preheating phase. We conclude that the preheating phase needs to ultimately be taken into account when we confront various inflationary models with cosmological observations in high accuracy. This fact calls for significant improvement in our fundamental understanding of the production of particles in the early universe. To gather more information connecting inflationary cosmology with cosmological observations, the detailed investigation on primordial gravitational waves during preheating process must be included, which has been extensively studied in the literature (Khlebnikov, 1997; Easther, 2006a, b, 2008; Garcia-Bellido, 2007, 2008; Dufaux, 2007; Price, 2008; Martin, 2014). We would like to leave the extension of this analysis based on our numerical code for a follow-up study.

Acknowledgments

We are grateful to J. Chen, J. Delabrouille, Z. Huang, J. Liu, X. Tong, L. Xue, D.-G. Wang, Z. Wang and S.-Y. Zhou for valuable discussions. YFC and JJ are supported in part by the NSFC (Nos. 11653002, 11722327, 11421303, J1310021), by National Youth Thousand Talents Program of China, by the CAST Young Elite Scientists Sponsorship Program (2016QNRC001), and by the Fundamental Research Funds for the Central Universities. QYL is supported in part by the Department of Physics & Astronomy at the University of Pennsylvania. YZ is supported in part by the NSFC (Nos. 11675165, 11633001), the CAS Strategic Priority Research Program “The Emergence of Cosmological Structures” (No. XDB09000000) and SRFDP. All numerical computations are operated on the computer clusters LINDA & JUDY in the particle cosmology group at USTC.

Appendix A Initial condition for the matter field $\chi$

We introduce the initial condition for the matter field $\chi$ at the end of the inflationary period, when $\epsilon=1$. We choose the Bunch-Davies vacuum to be the initial condition for cosmological perturbations.

We rewrite the equation of motion of the matter field by eliminating the first order derivative,

[TABLE]

where $X_{k}(t)=a^{\frac{3}{2}}(t)\chi_{k}(t)$ the $\omega_{k}$ is

[TABLE]

so the Bunch-Davies vacuum in the rescaled coordinate is

[TABLE]

the derivative is

[TABLE]

so in the physical coordinate (real part) the corresponding vacuum is

[TABLE]

the derivative is

[TABLE]

Appendix B Initial conditions for cosmological perturbations $Q_{\phi}$ and $Q_{\chi}$

Similar to the former subsection, we rescale the equation of cosmological perturbations in comoving time

[TABLE]

where the prime denotes the derivative with respect to the comoving time, $\mathcal{H}$ is the Hubble parameter in comoving time, $M_{\phi\phi}$ is given by V_{,\phi\phi}-\frac{8\pi G}{a^{3}}\frac{d}{dt}\bigg{(}\frac{a^{3}}{H}\dot{\phi}\dot{\phi}\bigg{)}, $M_{\phi\chi}$ is V_{,\phi\chi}-\frac{8\pi G}{a^{3}}\frac{d}{dt}\bigg{(}\frac{a^{3}}{H}\dot{\phi}\dot{\chi}\bigg{)}. Initially the matter field $\chi$ is absent, and hence $M_{\phi\chi}=0$.

We introduce the canonical variable $\nu_{\phi}=aQ_{\phi}$, and then derive the equation without first order derivative

[TABLE]

and thus, the Bunch-Davies vacuum for the perturbation of the Mukhannov-Sasaki variable takes $\nu_{\phi}=aQ_{\phi}\rightarrow\frac{e^{-i\omega_{k}\tau}}{\sqrt{2\omega_{k}}}$, where $\omega_{k}=k^{2}-\mathcal{H}^{\prime}-\mathcal{H}^{2}+M_{\phi\phi}$, and the derivative is

[TABLE]

Accordingly, in physical scale, the initial condition for the perturbation takes

[TABLE]

and the associated derivative term is expressed as

[TABLE]

Both real and imaginary part should be taken into consideration.