Aspects of wave turbulence in preheating II: Rebirth of the nonminimal coupled models

TL;DR

This paper investigates wave turbulence during preheating in a nonminimally coupled inflaton model with a quartic potential, highlighting energy transfer mechanisms and thermalization effects through numerical simulations.

Contribution

It demonstrates the role of wave turbulence in energy transfer during preheating in nonminimally coupled models, including the impact of backreaction on the universe's dynamics.

Findings

Turbulence influences the power spectra of inflaton fluctuations.

Energy transfer from inflaton to inhomogeneous modes is characterized.

Effective equation of state is determined considering backreaction.

Abstract

We study the nonlinear stage of preheating in a model consisting of a single inflaton field $\phi$ nonminimally coupled to the spacetime curvature and considering a self-coupling quartic potential V ($\phi$) = $\lambda \phi^4/4$. As the first motivation, we mention that this nonminimally coupled model agrees with the observational data. The second and the central issue of the present work is to exhibit aspects of wave turbulence associated with the mechanism of energy transfer from the inflaton field to the inhomogeneous fluctuations towards a state of thermalization. The imprints of the turbulence phase are mainly shown in the power spectra in time and wavenumber of relevant quantities such as the variance and the energy density of the inflaton field. We performed the simulations for several values of the nonminimally coupled models and taking into account the backreaction of the inflaton fluctuations in the dynamics of the Universe, and that allowed to determine the effective equation of state.

Figures

Equations

$$\mathcal{L}=\sqrt{-g}\left[\frac{R}{2m_{pl}^2}-\frac{1}{2}{\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi -\frac{1}{2}\xi R{\varphi }^2-V\left(\varphi \right)\right],$$

$$(1-\xi m_{pl}^{-2}{\varphi }^2)G_{\mu \nu }=\frac{1}{m_{pl}^2}\left\{(1-2\xi ){\partial }_{\mu }\varphi {\partial }_{\nu }\varphi -g_{\mu \nu }\left[\left(\frac{1}{2}-2\xi \right)(\partial \varphi {)}^2-V(\varphi )+2\xi \varphi \square \varphi \right]\right\}$$

$$\square \varphi -\xi R\varphi -\frac{\partial V}{\partial \varphi }=0,$$

$$V(\varphi )=\frac{1}{4}\lambda {\varphi }^4.$$

$$\mathrm{d}{s}^2=-\mathrm{d}{t}^2+a^2\left(t\right)\mathrm{d}\mathbf{x}.\mathrm{d}\mathbf{x},$$

$$G_{\mu \nu }=\frac{1}{m_{pl}^2}⟨T_{\mu \nu }^{\mathrm{eff}}⟩,$$

$${\mathbf{x}}_p=\sqrt{\lambda }{\varphi }_e\mathbf{x},\mathrm{d}{t}_p=\sqrt{\lambda }{\varphi }_e{a}^{-1}\mathrm{d}t,{\varphi }_p=a{\varphi }_e^{-1}\varphi .$$

$${\varphi }_p^{\prime \prime }-\frac{a^{\prime \prime }}{a}{\varphi }_p-{\nabla }_p^2{\varphi }_p+{\varphi }_p^3+\frac{\xi a^2}{\lambda {\varphi }_e^2}⟨R⟩=0,$$

$$⟨R⟩=6\left(\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}\right)=6\lambda {\varphi }_e^2\frac{a^{\prime \prime }}{a^3},$$

$${\varphi }_p^{\prime \prime }-\left(1-6\xi \right)\frac{a^{\prime \prime }}{a}{\varphi }_p-{\nabla }_p^2{\varphi }_p+{\varphi }_p^3=0.$$

$$H_p^2$$

$$\frac{a^{\prime \prime }}{a}$$

$$H_p^2=\frac{{\varphi }_e^2}{3m_{pl}^2}⟨\frac{...}{1-\xi m_{pl}^{-2}{a}^{-2}{\varphi }_e^2{\varphi }_p^2}⟩\approx \frac{{\varphi }_e^2}{3m_{pl}^2}\frac{⟨...⟩}{\left(1-\xi m_{pl}^{-2}{a}^{-2}{\varphi }_e^2⟨{\varphi }_p^2⟩\right)}.$$

$${\varphi }_N({\mathbf{x}}_p,\tau )=\sum _{k_x,k_y,k_z=-N/2}^{N/2-1}{\hat{\varphi }}_{\mathbf{k}}(\tau ){\psi }_{\mathbf{k}}({\mathbf{x}}_p),$$

$${\psi }_{\mathbf{k}}({\mathbf{x}}_p)={\mathrm{e}}^{\frac{2\pi i}{L_p}\mathbf{k}\cdot {\mathbf{x}}_p},$$

$${\varphi }_N({\mathbf{x}}_p,\tau )={\alpha }_0(\tau )+\sum _{\mathbf{k}\ne 0}{\hat{\varphi }}_{\mathbf{k}}(\tau ){\psi }_{\mathbf{k}}({\mathbf{x}}_p),$$

$$\mathrm{Res}({\mathbf{x}}_{p\mathbf{n}},\tau )={\varphi }_{p\mathbf{n}}^{\prime \prime }-(1-6\xi )\frac{a^{\prime \prime }}{a}{\varphi }_{p\mathbf{n}}+\frac{4\pi }{L_p^2}\sum |\mathbf{k}{|}^2{\hat{\varphi }}_{\mathbf{k}}(\tau ){\psi }_{\mathbf{k}}({\mathbf{x}}_{p\mathbf{n}})+{\varphi }_{p\mathbf{n}}^3=0.$$

$${\varphi }_{p\mathbf{n}}={\varphi }_N({\mathbf{x}}_{p\mathbf{n}},\tau )=\sum {\hat{\varphi }}_{\mathbf{k}}(\tau ){\psi }_{\mathbf{k}}({\mathbf{x}}_{p\mathbf{n}}).$$

$$⟨(...)⟩=\frac{1}{L_p^3}{\int }_{\mathcal{D}}(...){\mathrm{d}}^3{\mathbf{x}}_p\approx \frac{1}{L_p^3}\sum _{\mathbf{n}}(...{)}_{\mathbf{n}}{w}_{\mathbf{n}},$$

$${\sigma }_{\varphi }^2=⟨{\left({\varphi }_p-⟨{\varphi }_p⟩\right)}^2⟩=\sum _{\mathbf{k}}{\left|{\hat{\varphi }}_{\mathbf{k}}\right|}^2-{\alpha }_0^2.$$

$$w\left(\tau \right)=\frac{⟨p⟩}{⟨\rho ⟩},$$

$$\rho$$

$$p$$

$$P\left(k\right)\propto k^{-\gamma }.$$

$$\text{var}\left({\varphi }_p\right)={\left({\varphi }_p-⟨{\varphi }_p⟩\right)}^2.$$

$$\bar{\rho }=\frac{a^4\rho }{\lambda {\varphi }_0^4}.$$

$$\text{var}\left({\varphi }_p\right)$$

$$\bar{\rho }$$

$$P_{\sigma }\left(k\right)\propto k^{-5/3}{\exp }^{-7.0\times 10^{-4}{k}^{2.15}},$$

$$P_{\rho }\left(k\right)\propto \frac{k^3}{{\mathrm{e}}^{bk}-1},$$

$${\sigma }_{\varphi }^2\left(\tau \right)=\sum _{j=-N/2+1}^{N/2}{\hat{s}}_j{\mathrm{e}}^{{\omega }_j\tau },$$

$$P\left(\omega \right)\propto {\omega }^{-\gamma },$$

Full text

††: \JCAP

1 Introduction

The inflationary phase is characterized by accelerated expansion of the universe in its initial stages. It is remarkable that although cosmic inflation was formulated almost 40 years ago, it has resisted to all observational tests [1, 2]. We highlight that the transition from an empty and cold universe to radiation dominated and hot universe constitutes a challenging enterprise to connect the early stages to the description provided by the standard cosmological model.

Preheating is an initial phase of the process of reheating after inflation characterized by a rapid transfer of energy from the inflaton field into the other matter fields [3]. The parametric resonance governs the initial stage of this enegery transfer from the inflaton field coupled to other field fluctuations and can be studied using a linear approximation. Since we can only use the linear approximation for a limited interval of time once certain field fluctuations grow exponentially, a fully nonlinear evolution must be considered. Several works in this direction indicated the role played by the turbulence in the nonlinear stages of preheating towards the thermalization [4, 5, 6, 7, 8, 9]. We point out that the turbulence here must be understood as the wave turbulence since we have no usual fluids but only fields.

The first simulations of the nonlinear phase of preheating considered the simplest models as the single field inflation with quartic self-interaction potential $V(\phi)=\lambda\phi^{4}/4$ [10], or the inflaton with potential $V(\phi)=m^{2}\phi^{2}/2$ but coupled with another field [11]. However, according to the last survey of observational data, the oldest single field chaotic inflationary models are ruled out in particular due to the high level of the tensor-to-scalar ratio these models predict [1]. Remarkably, the model $V(\phi)=\lambda\phi^{4}/4$ including the nonminimal coupling with gravity described by $\xi\phi^{2}R/2$, where $R$ is the scalar of curvature, is in perfect agreement with the observational just assuming a small absolute value of the coupling constant $\xi$ [1, 12, 13].

In this paper, we studied the dynamics of a scalar field nonminimally coupled with gravity mentioned above taking into consideration distinct values of the coupling parameter within a range dictated by the observational data. In this case, after a successful inflationary phase, we start the evolution of the inflaton field together with its fluctuations, or inhomogeneous modes. We mention that Tsujikawa et al. [14] first studied the nonlinear preheating stage in this model taking into account the Hartree approximation in the field equation. Here, we have followed a distinct numerical approach including the backreaction inflaton fluctuations into the metric that allows us to determine the evolution of the effective equation of state. Another important aspect is that we focused on exhibiting the features of the wave turbulence through the power spectra of relevant quantities whose power-law behavior in specific regions signalizes the energy cascades from the inflaton to the inhomogeneous fluctuations. By assuming the form of the Planck distribution for the power spectrum of the inflaton energy density, we estimate the temperature at the thermalization phase.

The present work is organized as follows. In Sections I and II we present the field equations of the model and the numerical scheme, respectively. In Section IV, we explore the relevant aspects of the dynamics and the development of turbulence by presenting the power spectra of the variance and the energy density of the scalar field. We show in Section V the reheating temperatures for several values of the coupling parameter, and in the last Section, we conclude.

2 The model

Let us consider a scalar field $\phi(\mathbf{x},t)$ nonminimally coupled with the space-time geometry and represented by the Lagrangian [14]

[TABLE]

where $R$ is the scalar curvature, $\sqrt{-g}$ is the determinant of the space-time metric, $m_{pl}$ is the reduced Planck mass and $\xi$ is the coupling constant. The field equations reads

[TABLE]

[TABLE]

where $\square\phi=g^{\mu\nu}\nabla_{\mu}(\partial_{\nu}\phi)$. We consider the following self-coupling potential

[TABLE]

This model can be reconciled with the observational data for $\xi\leqslant-10^{-3}$ [2, 12, 13].

The spacetime is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) flat spacetime with the line element

[TABLE]

where $a(t)$ is the scale factor. The scalar field $\phi(\mathbf{x},t)$ has a homogeneous component, $\phi_{0}(t)$ that plays a relevant role during the inflationary phase, and an inhomogeneous component, $\delta\phi(\mathbf{x},t)$ due to the quantum fluctuations of the scalar field. The dymanmics of the spacetime has the contribution of the small inhomogeneites according with the following form of the field equations

[TABLE]

where $\left\langle...\right\rangle$ denotes the spatial average in the physical domain and $T_{\mu\nu}^{\mathrm{eff}}$ is the effective energy-momentum tensor obtained after writting the Eq. (3) in the form shown above. We impose the condition $\left\langle\delta\phi(\mathbf{x},t)\right\rangle=0$ for the fluctuations.

We introduce dimensionless variables to write the field equations for the numerical integration. For the model under consideration, we have

[TABLE]

The subscript $p$ indicates the dimensionless version of the variable in the computational domain, $\phi_{e}=\phi_{0}(t_{\mathrm{end}})$ is the amplitude of the homogeneous inflaton field at the end of inflation, and $t_{p}\equiv\tau$ is the conformal dimensionless time. With these new definitions, the Klein-Gordon equation (3) becomes

[TABLE]

where $\phi^{\prime}\equiv\partial\phi/\partial\tau$ and $\nabla^{2}_{p}$ is the flat three-dimensional Laplacian operator. Since

[TABLE]

it follows that

[TABLE]

The remaining relevant field equations are

[TABLE]

where $H_{p}=a^{\prime}/a$.

We integrate numerically the system of field equations (10) and (12) starting at the end of inflation and using the numerical scheme described in the next Section. We point out that Tsujikawa, Maeda and Torii [14] studied the preheating in this model in the context of the Hartree approximation and assuming a further approximation in the Hamiltonian constraint (11). Using the present notation, they have taken

[TABLE]

The authors of Ref. [14] pointed out the expectation that the above approximation did not affect significantly the dynamics of preheating, but it was not clear under which conditions it is valid. Then, we have opted to include in the numerical code the calculation of the spatial average present in Eqs. (11) and (12) without any approximation other than the provided by the numerical procedure.

3 Numerical approach

We integrate the field equations (10) and (12) numerically in a three-dimensional square box of comoving size $L$ with periodic boundary conditions as a suitable model of the Universe. The algorithm based on the collocation or pseudospectral method generalizes the previous two-dimensional version [9] devised to study the dynamics of several minimally coupled scalar field models in the nonlinear phase of preheating. For the sake of completeness, we present briefly the main aspects of the numerical procedure.

The scalar field is approximated by a series expansion with respect to a set of basis functions identified as the Fourier functions due to the periodic boundary conditions [15]. Then, we have

[TABLE]

where $\mathbf{k}=(k_{x},k_{y},k_{z})$ is the comoving momentum, $N$ is the truncation order, $\widehat{\phi}_{\mathbf{k}}(\tau)$ represents the unknown modes. The Fourier functions $\psi_{\mathbf{k}}(\mathbf{x}_{p})$ are

[TABLE]

with $L_{p}$ being the dimensionless size. The modes are complex functions, i. e., $\widehat{\phi}_{\mathbf{k}}(\tau)=\alpha_{\mathbf{k}}(t_{p})+i\beta_{\mathbf{k}}(\tau)$, but satisfy the relations $\widehat{\phi}_{\mathbf{k}}(\tau)=\widehat{\phi}^{*}_{\mathbf{-k}}(\tau)$ to produce a real scalar field in Eq. (14). In addition, the spectral approximation (14) comprises naturally the homogeneous and the inhomogeneous pieces of the scalar field:

[TABLE]

where $\alpha_{\mathbf{0}}(\tau)$ is the homogeneous component $\phi_{0}$ of the inflaton field.

The next step is to substitute the approximation (14) into the Klein-Gordon Eq. (10) to form the corresponding residual equation. According to the pseudospectral method, the unknown modes are such that the residual vanishes at the collocation, or grid, points $\mathbf{x}_{p\mathbf{n}}=L_{p}\mathbf{n}/N$. Therefore, we arrive at

[TABLE]

In the above expression, we express the Laplacian term in function of the modes; $\phi_{p\mathbf{n}}$ are the values of the scalar field at the collocation points connected with the modes through

[TABLE]

We use the fast Fourier transform (FFT) to invert this relation whenever necessary.

Concerning Eq. (12) we calculate at each time step the spatial average using quadrature formulae schematically indicated by

[TABLE]

where $(...)_{\mathbf{n}}$ indicates that the corresponding term is calculated at the quadrature points and $w_{\mathbf{n}}$ are the weights. For the sake of computational economy, we made the quadrature points to coincide with the collocation points. At this point, it is important to mention other codes for exploring the nonlinear phase of preheating like the LATTICEEASY [17], the DEFROST [18] and the PSpectre [19]. Although they are very efficient, we have implemented an alternative and valid pseudospectral algorithm to taken into account the backreaction of the field fluctuations besides the power spectra in the time and space domains for the selected quantities.

The integration starts at the end of inflation with $\alpha_{\mathbf{0}}(0)=1$ and $\alpha^{\prime}_{\mathbf{0}}(0)=0$ for the homogeneous inflaton field, $a(0)=1$ and the initial modes and their first derivatives, $\widehat{\phi}_{\mathbf{k}}(0),\widehat{\phi}^{\prime}_{\mathbf{k}}(0)$, have typical amplitudes of order of $10^{-4}$ [10]. We employ the FFT to obtain the values $\phi_{p\mathbf{n}}(0)$ and $\phi^{\prime}_{p\mathbf{n}}(0)$ associated with the inhomogeneous fluctuations and their first time derivatives necessary for the equations of motion (17). Taking these initial values we are able to calculate $a^{\prime}(0)$ from the Friedman equation. Therefore, with these intial values we can determine $a^{\prime},\phi_{p\mathbf{n}},\phi^{\prime}_{p\mathbf{n}}$ at the next time level (here the FFT is used to update the modes $\widehat{\phi}_{\mathbf{k}},\widehat{\phi}^{\prime}_{\mathbf{k}}$) and the whole process repeats resulting in the time march. We adopted the fourth order Runge-Kutta integrator using a stepsize $h=0.001$.

4 Wave turbulence towards the thermalization

We investigate the influence of the parameter $\xi$ on the aspects of wave turbulence developed in the nonlinear stage of preheating. We considered distinct values of $\xi$: $\xi=0$ representing the minimal coupling whose dynamics is already known; $\xi=-0.1$ that is in good agreement with the observational data [2, 13]; and $\xi=-1$ and $-10$ as examples of strong nonminimal couplings.

In all numerical experiments, we set $N=32$ that corresponds to a three-dimensional grid with $32^{3}$ points. For the sake of numerical resolution, we evolved the system in a total of $2^{24}$ time steps, or a time interval $\Delta\tau=16777.216$. Also, we choose the self-coupling constant $\lambda=10^{-8}$ that is consistent with the observational bounds $10^{-12}\lesssim\lambda\leqslant 10^{-4}$ imposed by the WMAP measurements to the tensor-to-scalar ratio and the spectral index [2, 12]. It represents an alleviation of the quite small value $\approx 10^{-14}$ in the case of the minimally coupled inflaton [16].

In the single field models with the quartic potential (4), the linear phase of preheating is governed by a Lamé-like equation with a narrow parametric resonant window. The number of resonant modes depends on the size of the box $L_{p}$ and the wavenumber magnitude $|\mathbf{k}|$. We choose $L_{p}\approx 42.7$ such that the resonant modes in the linear stage have $|\mathbf{k}|^{2}=74$.

4.1 Interaction between the modes

To study the energy transfer through different scales and the establishment of turbulent stages of the system, we analyze the decay of the homogeneous mode $\alpha_{\mathbf{0}}(\tau)$ and the variance of the scalar field $\sigma^{2}_{\phi}(\tau)$ given by

[TABLE]

This quantity gives us a relevant information about the contribution of the inhomogeneous modes of the inflaton field throughout the time evolution.

We notice that for the selected values of $\xi$, the time evolution of $\alpha_{\mathbf{0}}(\tau)$ and $\sigma^{2}_{\phi}(\tau)$ have a similar structure: the homogeneous mode oscillates with a constant amplitude until $\tau\approx 100$, and the variance (20) experiences exponential growth with a resonant peak approximately at the same instant. These aspects characterize the linear stage of the system when some of the inhomogeneous modes growth by the parametric resonance mechanism. From this instant, the nonlinearities of the system play a nonnegligible role leading the system to a stage of energy transfer through different scales called “driven turbulence”, where $\alpha_{\mathbf{0}}$ plays the “energy source” role. Besides, the amplitude of $\alpha_{\mathbf{0}}^{2}$ decays approximately with a power law $\tau^{-2/3}$ and the variance grows with power-law $\tau^{1/5}$, as shown in [4, 5]. However, we did not observe the decay of the variance with power-law $\tau^{-2/5}$ as found in [5], but instead, it seems that the variance reaches an approximate plateau.

According to Figs. 1 - 6, the dynamical behavior of these models is qualitatively similar. We noticed that the influence of the nonminimal coupling in the linear stage is to decrease the natural frequency of the homogeneous model slightly as illustrated in Fig. 7. The periods for $\xi=0,-0.1,-1.0$ and $-10.0$ are approximately $\Pi_{\phi}\approx$ 7.19, 7.20, 7.36 and 7.597, respectively. Therefore, for a small value of $|\xi|$ the behavior of the variance is not altered substantially, contrary to the results presented by Tsujikawa et al. [14].

4.2 Back reaction effects

We proceed with the influence of the nonminimal coupling parameter $\xi$ on the backreaction of the created particles in the effective equation of state and the expansion of the Universe described by the scale factor $a(\tau)$. The effective equation of state is denoted by the parameter $w$ given by

[TABLE]

where $\left<p\right>$ and $\left<\rho\right>$ are the effective pressure and energy density, respectively with the corresponding expressions derived from the field equations (6):

[TABLE]

We show in Fig. 8 the time evolutions of the effective equation of state $w$ for $\xi=0$ (left) and $\xi=-1$ (right). The mean behavior of $w$ tends to $w\approx 1/3$ (black line in the right side of Fig. 9) for all values $\xi$ that we have considered. This value would be expected for the radiation era. The noticeable effect of the parameter $\xi$ is to change the speed of growth of the scale factor $a(\tau)$ (left side of Fig. 9) but since the oscillations of the homogeneous component is not drastically affected, the energy transfer mechanism seems is almost the same for small values of $\xi$. The same occurs with the effective equation of state whose describing a radiation fluid in average.

4.3 Power spectra in spatial domain

Power spectra in the spatial domain (wave vectors) provide essential aspects about the energy transfer through different scales of the system in the turbulent phase with the homogeneous mode $(k\equiv\left|\mathbf{k}\right|=0$) acting as the energy source. The energy flows from larger scales, represented by the modes with small $k$, to smaller scales represented by modes with greater $k$, through constant fluxes named energy cascades. The presence of these cascades is signalized by scale laws on the power spectra in the spatial domain. In particular, for the hydrodynamic turbulence, some intervals of these power spectra should be described by power laws like

[TABLE]

To analyze the features of the wave turbulence, we have constructed the power spectra of two quantities related to the system. The first one is the function $\textrm{var}\left(\phi_{p}\right)$ given by

[TABLE]

And the second function is the rescaled energy density $\bar{\rho}$

[TABLE]

In both cases we evaluated the expansion of these quantities in terms of Fourier basis:

[TABLE]

where we set $N_{e}=8N$ and evaluated these expansions with the same FFT algorithm employed for the spectral expansion of $\phi_{p}$. The power spectra $P_{\sigma}\left(k\right)$ and $P_{\rho}\left(k\right)$ were obtained calculating the the root mean squares of the modes $\left|\widehat{v}_{\mathbf{k}}\right|$ and $\left|\widehat{e}_{\mathbf{k}}\right|$, respectively, having the same $k\equiv\left|\mathbf{k}\right|$. In particular, $\widehat{v}_{\mathbf{0}}=\sigma^{2}_{\phi}$ (cf. Eq. (20)).

We show in Figs. 10 - 12 the power spectra of $\textrm{var}\left(\phi_{p}\right)$ and $\bar{\rho}$ evaluated at $\tau\approx 8380$ (panels on the left) and at $\tau\approx 16770$ (panels on the right) for $\xi=0,\xi=-0.1$ and $-1$, respectively. We noticed that depending on the range of $k$ they share the same scale law. In the wavenumber interval $0\leqslant k\apprle 29$ we found $P_{\sigma}\simeq k^{-\gamma}$ with $\gamma\approx 1/3$ at $\tau=8380$ for all values of the parameter $\xi$; at $\tau=16770$, the exponent $\gamma$ changes according to $\xi$ as: $\gamma\approx 1$ for $\xi=0$ and $-0.1$, $\gamma\approx 1/2$ for $\xi=-1$ and $\gamma\approx 3/4$ for $\xi=-10$. More interestingly, in the interval $29\lesssim k\leqslant k_{\mathrm{max}}=80$, the scale law

[TABLE]

is valid for all values of the parameter $\xi$ and in both instants under consideration. In particular, this type of spectrum decay in wave numbers occurs in magnetohydrodynamic turbulence [20] and reflects a Kolmogorov-like law due to the turbulent cascade.

The power spectrum of the energy exhibits a more complex structure than the corresponding to the variance. For $k\lesssim 10$, all the spectra are described by the power-law (24) with $\gamma\simeq 1/7$. For $k>10$ we were able to identify two distinct components described by the following one-parameter function

[TABLE]

where $b\approx 0.18$ is independent of $\xi$ as well the instants of time under consideration. We identify this scaling law as the Planck’s law that describes the radiation spectrum of a black body. It might be remarked that the value of $b$ found in the present three-dimensional model is the same we have obtained in a two-dimensional model [9] of preheating after inflation with a minimally coupled inflaton field.

The scaling law (30) indicates the presence of the energy cascade that drives the system to the typical radiation energy distribution. As a consequence, the wave turbulence is the primary process driving the Universe to the radiation era. This feature is not altered by the values of the parameter $\xi$ despite the details of the power spectrum are slightly altered as it can be verified in Figs. 13 - 15. Another aspect denoting the cascade of energy the presence of an intermediate range that can be described approximately by the power law (cf. eq. 24) with $\gamma\approx 5/3$ (blue lines) that reveals the so-called inertial range. In the present simulations, we found the inertial range in the around the interval $10\lesssim k\lesssim 50$.

According to Nazarenko [21], 3D numerical simulations of Bose-Einstein condensation (BEC) phenomenon within the nonlinear Schrödinger (NLS) equation presents a gradual transition between a KZ power-law energy cascade to a Bose-Einstein distribution with chemical potential $\mu=0$, that is a Planck’s distribution. This author cites [22, 23, 24, 25, 26, 27], where was found evidences of KZ power law spectra with $\gamma\approx 5/3$. These authors argue that the presence of KZ power laws in 3D NLS models is natural, because of the similarity between these models and Euler’s hydrodynamics systems observed by Madelung transformations [21].

This kind of structure is very similar that we have found in our simulations. More specifically, in the top component of $P_{\rho}\left(k\right)$. This fact seems to be very reasonable from a physical perspective. Micha and Tkachev [5] identify the homogeneous component of the inflaton scalar field, that is the energy source in the driven turbulent stage, as a Bose-Einstein condensed. Besides representing waves, the Fourier modes of the inflaton spectral expansion are related to creation and annihilation operators. Naturally, the major part of the created particles should be bosons. Thus, the turbulent process in the preheating of the Universe should have a direct connection with the wave turbulence in the BEC phenomenon. In the next subsection, we show other aspects related to this fact and observed in the power spectra in the time domain.

We noticed that the power spectrum of the energy density, $P_{\rho}\left(k\right)$, exhibits two components and the division between these components occurs in $k\approx 30$. In the present simulations, we used the grid resolution $N=32$ and the maximum value for $k$ is $k_{max}=16\sqrt{3}\approx 27.7$. Nazarenko [21] comments about the consequences of a finite Fourier expansion on the numerical simulations of the type we performed here, and it might be an explanation for the split of the energy power spectrum near $k=k_{max}$. We interpret the second component lying at the bottom of the spectrum as describing a distribution of created particles that tend to thermalization, where the energy distribution is associated to the Planck’s law with the same temperature of the reheated Universe, $T_{reh}$.

4.4 Power spectra in time domain

In addition to the analysis of the power spectra in the spatial domain, we can also obtain insights about the energy flux from the inflaton to the particle production after analyzing the power spectra in time or the frequency domain. We have evaluated tha Fourier expansion of the variance $\sigma_{\phi}^{2}$, given by

[TABLE]

where $\omega_{j}=2\pi j/\left[\left(N-1\right)h_{e}\right]$ and $h_{e}$ is related to the discretization in the time domain. To obtain the best possible resolution in these spectra, we considered $h_{e}=h$. The graphics of the power spectra $P_{\sigma}\left(\omega\right)$ were built from the independent modes of the Fourier expansion of the variance expansions that are in the range $0\leqslant j\leqslant N/2$ due to their symmetries.

The Figs. 16, 17 e 18 show the power spectra corresponding to the intervals $\Delta\tau=2^{21}h=2097.152$ (panels on the left) and $\Delta\tau=2^{24}h$ (panels on the right), for $\xi=0$, $\xi=-0.1$ e $\xi=-1$, respectively. The feature shared by all power spectra is the presence of the Kolmogorov-Zakharov power law described by

[TABLE]

with $\gamma$ assuming distinct values such as $\gamma\approx 1,5/3,5/2$ depending on the ranges of $\omega$. The presence of the above power-law is a signature of any chaotic system; however by trying to understand the origin of these exponents can reveal valuable physical information.

In general, in all physical wave systems characterized by a unidimensional parameter (for example, the wave number $k=\left|\mathbf{k}\right|$) there exists a dispersion relation given by

[TABLE]

where the physical aspects of the system determine the constants $\lambda_{0}$ and $\alpha$. This feature allows the determination of the scaling law exponents for the energy, $\gamma_{E}$, and particles, $\gamma_{N}$, fluxes by the expressions (cf. [21])

[TABLE]

where $d$ is the number of the spatial domain dimensions in the models and $N_{w}$ is the number of the nonlinear interacting waves. In the preheating case, we have $d=3$ and $\alpha=1$ [4, 5] and the Eqs. 34 and 35 give us

[TABLE]

If $\gamma=\gamma_{E}$, the power spectra represents an inverse KZ energy cascade; if $\gamma=\gamma_{N}$, we have a double cascade, where there are a direct energy cascade and an inverse wave cascade. It means there is an energy flux, but the particle number is conserved, and this kind of behavior can drive the system to a BEC phenomenon [21].

Nazarenko [21] mentions the double cascade has physical meaning that are occurring energy transfer in direction of the greater energy levels and a far from equilibrium particle transfer to the lower energy levels. We noticed that the exponent $\gamma=1$ is independent of the kind of nonlinear interaction for 3D spatial models. In fact, this exponent is present in all power spectra, mainly in the ranges $\omega<\omega_{min}$ and $\omega>\omega_{max}$. This behavior takes place in generalized quantum models of Bose-Einstein gas, where it represents a direct energy cascade that originated in the condensates with lower $\omega$ values in direction to a range described by a thermal cloud for higher $\omega$ values. The description for lower $\omega$ values is given by a classical wave model and for higher $\omega$ values is given by classical particle models. This effect is called “wave-particle crossover” [21]. In particular for these models, the homogeneous mode frequency is $\omega_{0}=2\pi/\Pi_{\phi}$, the minimum and maximum frequencies, respectively $\omega_{min}$ and $\omega_{max}$, of the spectral expasion are $\omega_{min.}=2\pi/L_{p}$, $\omega_{max.}=32\pi\sqrt{3}/L_{p}$ and the resonant frequencies in the linear stage are $\omega_{res.}=c_{p}2\pi\sqrt{74}/L_{p}$, with $c_{p}\in\mathbb{N}$.

Another interesting point mentioned by Nazarenko [21] is the case of a system without dissipation on the lowest energy levels in which the inverse cascade represents a transient state and is characterized by an auto-similar solution. This solution is represented by a power law with $\gamma\approx 2.48\approx 5/2$, for $d=3$, that is different of the expected value for a KZ inverse cascade ($\gamma=7/6$ for $d=3$). It is a typical situation for power spectra with finite capacity in magneto-hydrodynamics models. Based on these facts and knowing that $\gamma\approx 5/3$ is a universal exponent in the Kolmogorov’s theory, we can conclude that the KZ energy cascades dominate in the interval $\omega_{min}<\omega<\omega_{max}$. The ranges in frequency with $\gamma\approx 5/3$ are related to direct cascades, and those with $\gamma\approx 5/2$ are related to inverse cascades that emerge from the auto-similar evolution cited by Micha and Tkachev in [4]. These cascades are concentrated around the $\omega_{min}$ and $\omega_{max}$, $\omega_{0}$ and $\omega_{res}$ .

When we compare the power spectra for different $\Delta\tau$, we can see that these cascade patterns are changing and the laws with $\gamma\approx 5/2$ and $\gamma\approx 1$ become more relevant that $\gamma\approx 5/3$. We believe this fact has a direct relation with the stage of turbulent transition as described in [5]. For the lowest $\Delta\tau$ we noted that stronger energy fluxes between the characteristics frequencies of the spectral expansion and major power laws with $\gamma\approx 1$ are out of this range. This aspect should be related to typical energy cascades of the turbulent process in the BEC phenomenon described by Nazarenko in [21]. For greater $\Delta\tau$, we have observed a tendency to power laws with $\gamma\approx 1$ and $\gamma\approx 5/2$. This fact can indicate that the driven turbulence stage is finishing, since the number of ranges with stationary fluxes is decreasing, and the free turbulence stage is establishing when the BEC effects are dominants, and the Universe will be driven to a thermalized state.

Again, we shall mention the possible influence of the nonminimal coupling on this process. We noted that all power spectra in the time domain, with the same $\Delta\tau$, present a very similar structure for $\xi=0$ and $\xi=-0.1$. For $\xi=-1$ there is subtle change of the spectrum structure, with the more predominance of power laws with $\gamma\approx 5/3$ in $\omega_{min.}<\omega<\omega_{max.}$. For greater $\Delta\tau$, all power spectra tend to have a similar power law structure, independently of $\xi$ value.

5 Energy distribution temperature and Universe thermalization

Following the same approach that we have done for two spatial dimensions [9], we can obtain the temperature of the energy distribution for differents values of $\xi$. In order to estimate this value, we shall restore the physical variables in the exponential argument, where we have $bk=\hbar ck_{phys}/k_{B}T$, $k_{phys}$ is the physical momentum, $k_{B}$ é is the Boltzmann’s constant and $T$ is the temperature of the distribution. For nonminimal coupling models, we have $k=L_{p}k_{phys}/\sqrt{\lambda}\phi_{0}$ and

[TABLE]

We note a decrease of the temperature as the parameter $\xi$ increases in absolute value. For $\xi=0$ the temperature is about $T\sim 10^{12}\,\mathrm{GeV}$ that is smaller than $T\sim 10^{14}\,\mathrm{GeV}$ we obtained for the model in two-dimensions [9]. Two factors explain this discrepancy: the difference of $L_{p}$ from $11.12$ (2d model) to $42.70$ in the present 3d model (the parameters $\phi_{0}$ and $b$ are unchanged) and a forgotten factor $8\pi$ due to the use of reduced Planck mass. Then, the present value of $T(\xi=0)\sim 10^{12}\,\mathrm{GeV}$ is consistent.

The obtained temperatures of the energy distribution might agree with the reheating temperature estimates, in short, $T_{\mathrm{reh}}$. The precise value of $T_{\mathrm{reh}}$ depends on the underlying microphysical processes; nevertheless, an upper bound of $T_{\mathrm{reh}}<10^{15}\,\mathrm{GeV}$ is estimated in Ref. [28], but distinct upper and lower bounds are possible such as $10^{8}\,\mathrm{GeV}\lesssim T_{\mathrm{reh}}\lesssim 10^{10}\,\mathrm{GeV}$ [29] or $10^{8}\,\mathrm{GeV}\lesssim T_{\mathrm{reh}}\lesssim 10^{11}\,\mathrm{GeV}$ [30]. As we have mentioned, the observation rule out the minimally coupled one field models, but left enough room to validate the nonminimally coupled fields with sufficient small $\xi$. In this case, the obtained temperatures $T_{\mathrm{reh}}$ are in agreement with the upper bounds of the reheating temperature.

6 Discussion

We presented the consequences of the nonlinear stages of preheating after inflation considering a single nonminimally coupled scalar field with quartic potential. The recent study of inflationary models [2] showed that even with a small coupling parameter $|\xi|>10^{-3}$ this single field model is not ruled out by the observations. We integrated the field equations starting from the end of inflation and taken into account the backreaction of the inhomogeneous fluctuations of the scalar field. For this aim, we have extended the code of Ref. [9] a three-dimensional spatial domain modeled by a box of periodic boundary conditions.

We explored the wave turbulence features through the power spectra of the variance and the energy density of the scalar field either in time and in function of the wave number. As a consequence of the energy transfer from the homogeneous inflaton field to the inhomogeneous fluctuations in the turbulent phase, some intervals of the power spectra are described by power-laws $k^{-\gamma}$ and $\omega^{-\gamma}$, with $k$ and $\omega$ representing the wavenumber and the frequency, respectively.

There are two other results worth of mentioning. The first is the equation of state arising from the average spatial procedure of the energy-momentum tensor that received the contribution of the inhomogeneous fluctuations of the inflaton field. We found that after the turbulent regime has been established the matter content behaves as a radiation fluid with $w(\tau)\approx 1/3$ for several values of the parameter $\xi$. The second result refers to the estimate of the reheating temperature from the power spectrum of the energy density in the wavenumber identified as the Planck law describing the radiation spectrum of a black body. We estimated the temperature associated with this spectrum the value $T\approx 10^{11}\,\mathrm{GeV}$ for $\xi$ inside the acceptable observational constraints, and is acceptable as a reheating temperature.

We continue the present investigation by studying the aspects of wave turbulence in the nonlinear stages of preheating in a distinct class of two-fields models in the next paper.

Acknowledgments

J. A. Crespo acknowledges the financial support of the Brazilian agency Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ). H. P. de Oliveira thanks to Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) (Grant No. E-26/202.998/518 2016 Bolsas de Bancada de Projetos (BBP)).