Stochastic Particle Production in a de Sitter Background

TL;DR

This paper investigates stochastic particle production in a de Sitter universe, modeling the effects of random mass variations and non-adiabatic events on scalar fields, revealing log-normal distributions and a Wiener process behavior.

Contribution

It introduces a transfer matrix approach to analyze stochastic particle production, deriving a Fokker-Planck equation for the evolution of system parameters in a de Sitter background.

Findings

Field amplitudes follow a log-normal distribution.

Logarithm of the field amplitude behaves like a Wiener process.

Analytical solutions for distributions in the weak-scattering limit.

Abstract

We explore non-adiabatic particle production in a de Sitter universe for a scalar spectator field, by allowing the effective mass $m^2(t)$ of this field and the cosmic time interval between non-adiabatic events to vary stochastically. Two main scenarios are considered depending on the (non-stochastic) mass $M$ of the spectator field: the conformal case with $M^2=2H^2$, and the case of a massless field. We make use of the transfer matrix formalism to parametrize the evolution of the system in terms of the "occupation number", and two phases associated with the transfer matrix; these are used to construct the evolution of the spectator field. Assuming short-time interactions approximated by Dirac-delta functions, we numerically track the change of these parameters and the field in all regimes: sub- and super-horizon with weak and strong scattering. In all cases a log-normally distributed field amplitude is observed, and the logarithm of the field amplitude approximately satisfies the properties of a Wiener process outside the horizon. We derive a Fokker-Planck equation for the evolution of the transfer matrix parameters, which allows us to calculate analytically non-trivial distributions and moments in the weak-scattering limit.

Equations

$${\ddot{\chi }}_k+3H{\dot{\chi }}_k+\left[\frac{k^2}{a^2}+M^2+m^2(t)\right]{\chi }_k=0,$$

$${\partial }_{Ht}⟨\ln |{\chi }_k{|}^2⟩={\mu }_1-2\mathrm{and}{\partial }_{Ht}\mathrm{Var}[\ln |{\chi }_k{|}^2]={\mu }_2.$$

$$|{\chi }_k(t){|}_{\mathrm{typ}}^2\equiv e^{⟨\ln |{\chi }_k(t){|}^2⟩},$$

$$⟨Z_k(t)Z_{k^{\prime }}(t^{\prime })⟩\simeq {\mu }_{2}H\text{min}[t-t_k,t-t_{k^{\prime }},t^{\prime }-t_k,t^{\prime }-t_{k^{\prime }}],$$

$$⟨|{\chi }_{k_1}(t_1){|}^2\cdots |{\chi }_{k_n}(t_n){|}^2⟩=\exp \left[\sum _{i=1}^n⟨\ln |{\chi }_{k_i}(t_i){|}^2⟩+\frac{1}{2}\sum _{i,j=1}^n⟨Z_{k_i}(t_i)Z_{k_j}(t_j)⟩\right].$$

$$\mathcal{S}$$

$$\displaystyle=\;\frac{1}{2}\int d^{3}{\bf x}\,d\tau\,a^{2}\Big{[}(\partial_{\tau}\hat{\chi})^{2}-(\nabla\hat{\chi})^{2}-\left(M^{2}+m^{2}(t)\right)a^{2}\hat{\chi}^{2}\Big{]}\,,$$

$$\displaystyle=\;\frac{1}{2}\int d^{3}{\bf x}\,d\tau\,\Big{[}(\partial_{\tau}\hat{X})^{2}-(\nabla\hat{X})^{2}-\left(a^{2}M^{2}+a^{2}m^{2}(\tau)-\frac{\partial^{2}_{\tau}a}{a}\right)\hat{X}^{2}\Big{]}\,,$$

$$\hat{X}\equiv a\hat{\chi },$$

$$\left[{\partial }_{\tau }^2-{\nabla }^2-\frac{a^{\prime \prime }}{a}+a^2(M^2+m^2(\tau ))\right]\hat{X}(\tau ,\mathbf{x})=0.$$

$$\hat{X}({\bf x},\tau)\;=\;\int\frac{d^{3}{\bf k}}{(2\pi)^{3/2}}\,e^{-i{\bf k}\cdot{\bf x}}\Big{[}X_{k}(\tau)\hat{a}_{{\bf k}}+X^{*}_{k}(\tau)\hat{a}^{\dagger}_{-{\bf k}}\Big{]}\,,$$

$$X_k^{\prime \prime }(\tau )+\left[k^2-\frac{a^{\prime \prime }}{a}+a^2(M^2+m^2(\tau ))\right]X_k(\tau )=0,$$

$$X_k(\tau )={\alpha }_{k,j}{f}_k(\tau )+{\beta }_{k,j}{f}_k^{\ast }(\tau ),$$

$$f_k^{\prime \prime }(\tau )+\left[k^2-\frac{a^{\prime \prime }}{a}+a^2{M}^2\right]f_k(\tau )=0.$$

$$|{\alpha }_{k,j}{|}^2-|{\beta }_{k,j}{|}^2=1.$$

$$H(t)=\dot{a}(t)/a(t)=const.\text{and}\tau =-1/aH<0.$$

$$f_k^{\prime \prime }(\tau )+k^2{f}_k(\tau )=0,$$

$$f_k(\tau )=\frac{e^{-ik\tau }}{\sqrt{2k}}.$$

$$f_k(\tau )=\frac{e^{-ik\tau }}{\sqrt{2k}}\left(1-\frac{i}{k\tau }\right).$$

$$\left(\begin{array}{c}{\beta }_j\\ {\alpha }_j\end{array}\right)={\mathsf{M}}_j\left(\begin{array}{c}{\beta }_{j-1}\\ {\alpha }_{j-1}\end{array}\right).$$

$$\left(\begin{array}{c}{\beta }_j\\ {\alpha }_j\end{array}\right)=\mathsf{M}(j)\left(\begin{array}{c}{\beta }_0\\ {\alpha }_0\end{array}\right),\text{where}\mathsf{M}(j)={\mathsf{M}}_j{\mathsf{M}}_{j-1}\cdots {\mathsf{M}}_1.$$

$$\mathsf{M}(j)=\left(\begin{array}{cc}{e}^{i\varphi }& 0\\ 0& e^{-i\varphi }\end{array}\right)\left(\begin{array}{cc}\sqrt{1+n}& \sqrt{n}\\ \sqrt{n}& \sqrt{1+n}\end{array}\right)\left(\begin{array}{cc}{e}^{i\psi }& 0\\ 0& e^{-i\psi }\end{array}\right),$$

$$\varphi (j)=\frac{1}{2}\arg ({\beta }_j{\alpha }_j^{\ast }),\psi (j)=-\frac{1}{2}\arg ({\beta }_j{\alpha }_j),$$

$$n(j)$$

$$=\frac{1}{4}\mathrm{Tr}\left[\mathsf{M}(j){\mathsf{M}}^{†}(j)-1\right].$$

$$m^2(t)$$

$$=\sum _j\frac{m_j}{a({\tau }_j)}\delta (\tau -{\tau }_j)=m^2(\tau ).$$

$$X_k^{\prime \prime }(\tau )+\left[k^2-\frac{a^{\prime \prime }}{a}+a^2{M}^2+\sum _j{m}_{j}a({\tau }_j)\delta (\tau -{\tau }_j)\right]X_k(\tau )=0.$$

$$X_{k,j+1}({\tau }_j)$$

$$X_{k,j+1}^{\prime }({\tau }_j)$$

$${\mathsf{M}}_j$$

$$m_{j}a({\tau }_j)|f({\tau }_j){|}^2\ll 1.$$

$$P_{t+\delta t}(\mathsf{M})=\int d{\mathsf{M}}_2{P}_t({\mathsf{M}}_2^{-1}\mathsf{M})P_{\delta t}({\mathsf{M}}_2).$$

Full text

1. 1 Introduction

2. 2 Summary of the main results

3. 3 Mathematical preliminaries

4. 3.1 Spectator field in an expanding Universe

5. 3.1.1 Mode functions in de Sitter spacetime

6. 3.2 The Transfer Matrix formalism

7. 3.3 The Fokker-Planck equation

8. 4 Conformally massive field in de Sitter background

9. 4.1 Numerical results

10. 4.1.1 Individual realizations

11. 4.1.2 Means and variances

12. 4.1.3 Probability densities

13. 4.2 The field two-point function

14. 4.3 Analytical results

15. 4.3.1 Sufficiently sub-horizon

16. 4.3.2 Outside the horizon

17. 5 Massless field in de Sitter background

18. 5.1 Numerical results

19. 5.1.1 Individual realizations

20. 5.1.2 Means and variances

21. 5.1.3 Probability densities

22. 5.2 The field two-point function

23. 5.3 Analytical results

24. 5.3.1 Sufficiently sub-horizon

25. 5.3.2 Outside the horizon

26. 6 Conclusions

27. A Appendix

28. A.1 Typical vs. average

29. A.2 Convergence tests

30. A.3 Backreaction

31. A.3.1 Conformally massive field

32. A.3.2 Massless field

1 Introduction

The embedding of the inflationary paradigm within ultraviolet completions of particle theories often involves many fields with potentially complicated interactions that may lead to a chaotic evolution as a function of the initial conditions and values of the model parameters. The presence of such a large number of degrees of freedom can also dramatically complicate the dynamics of post-inflationary reheating.

Although the full deterministic description of such models can be highly model-dependent, one might expect that in the limit of many fields/interactions, emergent universal properties may arise. Moreover, the coarse-grained nature of the available cosmological observations is unlikely to shed light on all the microscopic details of the fundamental theory. A theoretical framework that advocates a coarse-grained and approximately model-independent approach, with a focus on universal features, is at the heart of our efforts here and in [1, 2, 3] (also see [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]).

It is plausible that the complex dynamics of fields can lead to repeated non-adiabatic particle production in the inflationary [19, 20, 21, 22, 23] and post-inflationary [24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38] universe. With sufficient complexity, the strength of the interactions and intervals between them can be treated stochastically; statistical tools can then be invoked without having to rely on detailed model building. In earlier work [1, 2], we were motivated by the connection between particle production in cosmology and current conduction in disordered wires  (also see [39, 40]). Apart from the elegant mathematical correspondence, the primary drive there was that certain universal features, such as Anderson Localization [41] in one dimension, arise independent of the details of the systems – motivating a search for similar universality in particle production.

In previous works [1, 2], the problem of stochastic, non-adiabatic particle production has been formulated exclusively in a non-expanding Minkowski background, for simplicity. In [1] the evolution of the occupation number for a single scalar degree of freedom was studied in detail, in the limit of narrowly localized interactions in time; this allows for a quasi-discrete description of the dynamics by means of the Transfer Matrix formalism. Under the assumption that each scattering can be treated as a perturbation of the transfer matrix, the authors derived a Fokker-Planck equation describing the dynamical evolution of the probability distribution for the occupation number of the scalar field. The results were then generalized for multiple fields by imposing a maximality constraint on the Shannon entropy of the probability distribution; this constraint is known as the Maximum Entropy Ansatz [42] (MEA), and results in a dramatic reduction of the effective degrees of freedom that describe the average behavior of the system.

In [2] the Fokker-Planck formalism was extended to the case of multiple statistically inequivalent fields, with stochastically varying effective masses, cross couplings and intervals between interactions, which, as we will demonstrate later, somewhat mimics the phase scrambling that takes place in an expanding universe. The main results therein were (1) a practical demonstration of the equivalence between the MEA and statistically equivalent interacting fields, and (2) the convergence to the MEA in the limit of large number of (possibly statistically inequivalent) fields.

Particle production in an expanding universe is distinct from the flat space case. An expanding universe introduces a competition between particle production from the interactions and dilution. More importantly, the existence of the Hubble horizon introduces an additional scale into the problem, with qualitatively different behavior of particle production expected in the spectator fields on super-horizon scales compared to the sub-horizon case. In spite of these complications, we find a surprisingly simple and universal behavior of the non-adiabatically excited fields on sub-horizon (which is expected from earlier work) and on super-horizon scales (which is new to this work). In upcoming work, the results from this manuscript will be used to calculate the curvature fluctuations resulting from the particle production during inflation, and to estimate the efficiency of reheating after inflation.

For the sake of simplicity we will mostly restrict ourselves to the single spectator field case in de Sitter space. Most of our mathematical framework is valid for a general expansion history and a general mass of the spectator field. However, to contain this already long paper to a manageable size we have limited the detailed discussion to (1) conformal mass ($M^{2}=2H^{2}$) and (2) massless spectator fields ($M^{2}=0$) in de Sitter space ($H={\rm const.}$). These choices are phenomenologically interesting. As an example, in supergravity models with minimal kinetic terms, the large vacuum density during inflation $V\sim H^{2}M_{P}^{2}$, where $M_{P}$ denotes the Planck mass, typically leads to an induced mass for all scalar fields of order $H$ [43, 44]. As we shall see, setting the constant of proportionality to $\sqrt{2}$ in the conformally massive case greatly simplifies our calculations. Similarly, the massless case can approximate light fields ($M\ll H$), which may be easily perturbed during inflation and source curvature fluctuations [45]. While $M=0$ and $M=\sqrt{2}H$ are special in terms of their calculational convenience, they are not special in terms of the physical implications of our results.

Our analysis will be restricted to the linear regime of the spectator field, for which each Fourier mode can be treated independently and the backreaction on the homogeneous expanding background can be ignored; these assumptions can break down when the energy density of the spectator fields becomes sufficiently large. We will show that significant amount of scattering is allowed for a sufficient number of e-folds to make this analysis worthwhile, and relevant for calculating observables. In addition to the above simplifying assumptions, we will also consider for the sake of analytic and numerical tractability that each interaction can be modeled as a Dirac-delta function333ie. we assume that the physical wavelengths are large compared to the duration of the non-adiabatic interactions. in time whose amplitude and location are drawn from different distributions.

The narrow-width interactions allow us to use the Transfer Matrix Formalism quite efficiently since the evolution between scatterings is that of free fields. We will go beyond the assumption of small changes per scattering in our numerical explorations, although our analytical understanding based on a Fokker-Planck equation will be robust in the weak scattering case only. In distinction with earlier papers, we prefer to follow Fourier modes of the spectator field rather than the occupation number density. This is natural since the occupation number density is ill-defined on super-horizon scales. For an application to inflation which we will pursue in an upcoming paper, appropriate combinations of these Fourier modes of the spectator field will serve as a source for curvature perturbations. We can also use these to calculate gravitational wave production from this period.

The rest of the paper is organized as follows:

Section 2 provides a bird’s eye view of the most important, and simplest to state, results of our analysis. We caution that a lot is left out here; we intend this section to be more of an invitation to explore the analysis in the rest of the text.

Section 3 contains the formalism necessary to study the dynamics of a spectator field excited by a non-adiabatic, stochastic mass term in an expanding background. In Section 3.1 we introduce the effective single-field model in an expanding Universe that we will study. In particular we discuss a conformally massive and a massless field in a de Sitter background. In Section 3.2 we describe the transfer matrix formalism that will allow us to track the evolution of the scalar field and its number density after each consecutive scattering. Section 3.3 contains a brief summary of the Fokker-Planck formalism, with emphasis on single-field models. The general results of this section are applied to the specific cases of the conformally massive (and massless) spectator field in the subsequent sections.

Section 4 contains the results for the evolution of a conformally massive scalar field, its occupation number and other transfer matrix parameters in a de Sitter background. Section 4.1 shows numerical results in the weak- and strong-scattering limits for the field amplitude and the transfer matrix parameters, including their values given individual realizations of an ensemble of location and scattering amplitudes, as well as their probability densities and their lowest moments. In Section 4.3 we describe the analytical results obtained from the application of the Fokker-Planck formalism, which are valid in the weak-scattering regime, for which the instantaneous change in the transfer matrix can be treated perturbatively. We also demonstrate how these results (for sub-horizon modes) correspond to a natural generalization of the Minkowski result discussed in earlier papers.

Section 5 provides a discussion of the corresponding numerical (Section 5.1) and analytical (Section 5.3) results for a massless scalar field in a de Sitter background.

Section 6 contains a summary of our results and our conclusions.

In Appendix A we provide some essential checks for our results and some justification for our focus on certain variables in the main text. In the Appendix A.1 we discuss the difference between typical and average quantities when the distributions are not normal. In Appendix A.2, we verify the approximate independence of our results from the details of the distribution from which the effective mass and time of scattering are drawn. Finally, in Appendix A.3 we determine the regime in which the excitation of the field is strong enough to backreact on the expanding background. We show the domain of validity of our results in terms of the strength of the non-adiabatic events and the number of $e$-folds of inflation during which the scalar field is excited.

2 Summary of the main results

Consider a Fourier mode $\chi_{k}(t)$ of the spectator field $\chi(t,\bf x)$ in de Sitter space satisfying the equation of motion444While we discuss $\chi_{k}$ here, we find it more convenient to use the scaled field $X_{k}=a\chi_{k}$ in the main text. We also drop the subscript $k$ (ie. the momentum dependence) in denoting $X$ and other related quantities in much of the main text (though we do analyze the behavior with $k$). (see Section 3 for details):

[TABLE]

where $H$ is the expansion rate, $M$ is the mass of the field and we include a stochastic mass term $m^{2}(t)=\sum_{j=1}^{N_{s}}m_{j}\delta(t-t_{j})$ to capture the complicated interaction that this field is undergoing with other fields/background. The masses $m_{j}$ and locations $t_{j}$ are drawn from independent distributions. We assume that each $m_{j}$ is independent and identically distributed (not necessarily Gaussian), with $\langle m_{i}m_{j}\rangle=\sigma^{2}\delta_{ij}$ where $\langle\ldots\rangle$ represents an average over the ensemble. For convenience, we define $\mathcal{N}_{s}$ as the number of non-adiabatic events per Hubble time and assume that $\mathcal{N}_{s}\gg 1$.555In this limit the specific form of the distribution of $m_{i}$ is irrelevant, see Appendix A.2. The single dimensionless scattering parameter $\mathcal{N}_{s}(\sigma/H)^{2}$ is sufficient to determine the statistical behavior of the fields outside the horizon.

The behavior of $\ln|\chi_{k}(t)|^{2}$ for different realizations of $m^{2}(t)$ are shown in Fig. 1. We provide the most important and simple to understand takeaways from our analysis below:

The $\chi$ field is approximately in its vacuum state sufficiently inside the horizon (i.e. $|\chi_{k}|^{2}\simeq 1/(2ka^{2})$).666This statement has caveats, in terms of the magnitude of $\mathcal{N}_{s}$ and $\mathcal{N}_{s}(\sigma/H)^{2}$ as well as initial conditions. A large magnitude of these parameters can lead to deviations from the vacuum behavior as one would expect. We explore these caveats and details further in the main text. Outside the horizon, $\ln|\chi_{k}|^{2}$ evolves linearly with cosmic time (in an ensemble averaged sense), with

[TABLE]

where the variance and mean are over different realizations of the effective mass $m^{2}(t)$. The rates $(\mu_{1},\mu_{2})$ are functions of $\mathcal{N}_{s}(\sigma/H)^{2}$. The values of $\mu_{1}$ and $\mu_{2}$ as a function of $\mathcal{N}_{s}(\sigma/H)^{2}$ are shown in Fig. 8 (conformal mass) and Fig. 26 (massless). 2. 2.

Importantly, $\ln|\chi_{k}(t)|^{2}$ is normally distributed on super and sub-horizon scales at all times (as an ensemble over realizations of $m^{2}(t)$). Equivalently, $|\chi_{k}(t)|^{2}$ is log-normally distributed. This means that

[TABLE]

is a better representative of the ensemble rather than $\langle|\chi_{k}(t)|^{2}\rangle$, which will be dominated by the largest values of $\chi_{k}(t)$ in the ensemble. 3. 3.

On super-horizon scales, $\ln|\chi_{k}|^{2}$ satisfies the properties of a drifted random walk. In particular, as mentioned above, the mean and variance of $\ln|\chi_{k}|^{2}$ grow linearly with time and for the drift-less variable $Z_{k}(t)\equiv\ln|\chi_{k}|^{2}-\langle\ln|\chi_{k}|^{2}\rangle$, we find

[TABLE]

where $t_{k}$ is the time when the $k$-mode exits the horizon (see Fig. 2). Note that the above condition contains within it the statement that the increments $Z_{k}(t_{4})-Z_{k}(t_{3})$ and $Z_{k}(t_{2})-Z_{k}(t_{1})$ are uncorrelated when the interval $(t_{3},t_{4})$ does not overlap with $(t_{1},t_{2})$ on super-horizon scales.

The behavior of $\ln|\chi_{k}|^{2}$ implies that $|\chi_{k}(t)|$ performs a geometric random walk. With this understanding, all $n$-point correlation functions for the field magnitude can be computed in terms of the $Z_{k}$ two-point functions,

[TABLE] 4. 4.

The phase $\arg\chi_{k}$ is randomly distributed inside the horizon and converges to an $\mathcal{N}_{s}(\sigma/H)^{2}$-dependent value outside the horizon. 5. 5.

In the weak scattering limit, using a Fokker-Planck equation, we analytically derive the log-normal probability distribution for $|\chi_{k}(t)|$ on sub-horizon scales (though we cannot do so yet on super-horizon scales). On super-horizon scales, we can derive the time rate of evolution for $\langle\ln|\chi_{k}|^{2}\rangle$ in the weak scattering limit. These results are consistent with our numerical simulations.

While the behavior of the field is easiest to discuss, and perhaps the most useful for future calculations, we found it useful and at times necessary to understand the behavior of the transfer matrix parameters individually (a combination of which yields the field amplitude and phase). These parameters include an “occupation number density”($n_{k}$) and two phases $\phi_{k}$ and $\psi_{k}$ (see Section 3 for definitions):

The occupation number density grows exponentially inside and outside the horizon. The growth rates of the ensemble mean and variance of $\ln(1+n_{k})$ are linear in cosmic time outside the horizon, and are determined by $\mathcal{N}_{s}(\sigma/H)^{2}$. 2. 2.

The phases are uniformly distributed on subhorizon scales, but their distribution will in general depend on $\mathcal{N}_{s}(\sigma/H)^{2}$ outside the horizon. 3. 3.

The means and variances of $\phi_{k}$ and $\psi_{k}$ undergo non-trivial evolution at horizon crossing, but converge to constants (which can be $\mathcal{N}_{s}(\sigma/H)^{2}$ dependent) outside the horizon. 4. 4.

In the weak scattering limit, we derive the behavior of the occupation number and phases, including their distribution and evolution rates of the lowest moments, inside the horizon using a Fokker-Planck equation. Outside the horizon, we derive a highly non-trivial distribution of the phase $\phi_{k}$ as well the evolution rates of the mean of $\ln(1+n_{k})$.

We re-iterate that this is a rather cursory summary. Details, caveats and many other relevant results, which are not included here, are discussed for the main text.

3 Mathematical preliminaries

3.1 Spectator field in an expanding Universe

Consider a spectator field777$\chi$ is assumed to contribute a negligible amount to the total energy density of the universe. $\chi(t,{\bf x})$ of mass $M$ in a homogeneous and isotropic expanding universe. For this field, we use an additional effective mass $m(t)$ to parametrize the coupling of this field to a time-dependent background, especially including random non-adiabatic events arising from complicated interactions with other fields. The quadratic action for such a (quantum) field is taken to be

[TABLE]

where $a$ is the scale factor, the conformal time $\tau$ is related to cosmic time $t$ via $dt/d\tau=a$, and where we have defined

[TABLE]

in going from the second to the third line. In slight abuse of notation, $m^{2}(\tau)=m^{2}(t(\tau))$ and $a(\tau)=a(t(\tau))$ above. While it is most convenient to write the necessary equations and formalism in terms of $X$ and $\tau$, some of our results are most naturally written in terms of $\chi$ and $t$.

The equation of motion for the $\hat{X}$ field is then given by

[TABLE]

The mode expansion for this field can be written as

[TABLE]

where $[\hat{a}_{{\bf k}},\hat{a}^{\dagger}_{{\bf k}^{\prime}}]=\delta({\bf k}-{\bf k}^{\prime})$, $[\hat{a}_{{\bf k}},\hat{a}_{{\bf k}^{\prime}}]=[\hat{a}^{\dagger}_{{\bf k}},\hat{a}^{\dagger}_{{\bf k}^{\prime}}]=0$. The mode functions $X_{k}(\tau)$ satisfy

[TABLE]

and are normalized888For consistency with the canonical commutation relations between $\hat{X}(t,{\bf x})$ and its conjugate momentum. by means of the Wronskian condition $X_{k}(\tau)X_{k}^{*^{\prime}}(\tau)-X_{k}^{\prime}(\tau)X_{k}^{*}(\tau)\;=\;i\,$. In addition, the physical mode functions are chosen so that in the infinite past the vacuum is of Bunch-Davies type.

In general, the time dependence of the stochastic mass $m(\tau)$ can be complicated. To track the evolution of the field, we will assume that this effective mass consists of localized, non-adiabatic events. In between these non-adiabatic events, the free field solutions to (3.7) will have the form

[TABLE]

where $f_{k}(\tau)$ satisfies (3.7) in absence of $m^{2}(\tau)$:

[TABLE]

The coefficients $\alpha_{k,j}$ and $\beta_{k,j}$ are the Bogoliubov coefficients after the $j$-th non-adiabatic event, with

[TABLE]

This constraint forces $f_{k}(\tau)$ to satisfy the same normalization condition as $X_{k}(\tau)$ above. $f_{k}(\tau)$ is then completely specified (up to an irrelevant phase) provided it also satisfies the Bunch-Davies vacuum initial condition in the infinite past. Note that before any non-adiabatic interactions, $\alpha_{k,0}=1$ and $\beta_{k,0}=0$. The quantity $|\beta_{k,j}|^{2}$ can be interpreted as the occupation number density of particles of the field $\hat{X}$ with momentum $k$ after the $j$th non-adiabatic event.999The total number density would be $\int d^{3}k\,|\beta_{j,k}|^{2}$. We caution that this interpretation of $|\beta_{k,j}|^{2}$ as an occupation number density, however, does not carry over easily on super-horizon scales [46].

3.1.1 Mode functions in de Sitter spacetime

The focus of our discussion from Section 4 onwards will be the study of spectator fields in a de Sitter background where

[TABLE]

General expressions for $f_{k}(\tau)$ are available in this case in terms of Hankel functions (see e.g. [46]). There are two cases where the form of $f_{k}(\tau)$ is even simpler.

Conformally Massive Fields : When the mass of the field $M^{2}=2H^{2}$ in a de Sitter background, Eq. (3.9) becomes the equation of motion for a free field in a non-expanding universe:

[TABLE]

with a solution

[TABLE]

Note that while $f_{k}(\tau)$ does not see the effects of expansion, the full solution $X_{k}(\tau)$ can depend on expansion through the non-adiabatic term $a^{2}m^{2}(\tau)$ in Eq. (3.7).

Massless Fields : In this case $M^{2}=0$, and the mode functions have the form

[TABLE]

3.2 The Transfer Matrix formalism

The assumption of localized interactions in the effective mass $m^{2}(\tau)$ allows for a transfer matrix approach for the determination of the coefficients $\alpha_{k}$ and $\beta_{k}$. We define the transfer matrix at the location of the $j$-th scattering, ${\sf M}_{j}$, to be such that

[TABLE]

Here we have ignored the momentum dependence (ie. subscript $k$) for notational simplicity. By chaining together ${\sf M}_{j}$ with all transfer matrices ${\sf M}_{i<j}$ we can construct $\alpha_{j},\beta_{j}$ from their initial values,

[TABLE]

Note that starting with fields in the Bunch-Davies vacuum is consistent with choosing $\{\beta_{0},\alpha_{0}\}=\{0,1\}$. A general parametrization of the matrix ${\sf M}(j)$ can be written as [47]

[TABLE]

where the angular parameters can be identified as follows,

[TABLE]

and $n(j)$ is the occupation number density:

[TABLE]

Note that two complex numbers $\alpha$ and $\beta$ keep track of the field evolution which we are interested in (see Eq. (3.8)). These numbers have one constraint (3.10) which leaves 3 independent real numbers. These three numbers are conveniently parametrized by $n$, $\phi$ and $\psi$ in the transfer matrix ${\sf M}$.

Clearly, the functional dependence of ${\sf M}_{j}$ will be determined by the explicit form for $m^{2}(t)$. Following the analysis of [1, 2], in order to obtain analytically tractable expressions we will restrict ourselves to the assumption of Dirac-delta scatterers in cosmic time,

[TABLE]

Physically, we are assuming that the temporal width of the scatterers is much smaller than the characteristic period of $f_{k}(\tau)$. With such Dirac-Delta scatterers, (3.7) takes the form

[TABLE]

The junction conditions at each $\tau_{j}$ correspond to

[TABLE]

Translated to (3.15), it implies the following general form for the transfer matrix:

[TABLE]

By multiplying transfer matrices with the form derived above, we can numerically compute the evolution of the occupation number density, phases and the field amplitude.

Note that the addition of an extra scattering event to ${\sf M}(j)$ may be treated as a perturbation if the instantaneous scattering amplitude is such that

[TABLE]

For our general investigation (including our numerical simulations), we do not assume a small perturbation. However, the Fokker-Planck formalism described immediately below relies on this small perturbation assumption.

3.3 The Fokker-Planck equation

The stochastic nature of the effective mass $m^{2}(t)$ implies that the non-adiabatic event amplitudes $m_{j}$ and locations $t_{j}$ are assumed to be drawn from some distribution. This in turn implies that the transfer matrices ${\sf M}_{j}$ and ${\sf M}(j)$ will also be stochastic in nature, and will take different values for different realizations of the $m_{j},t_{j}$. We will therefore consider an ensemble of realizations for the amplitudes and locations of the scatterings, over which we can define a probability density $P_{t}({\sf M})$. Any physically meaningful quantity can then be obtained from expectation values with respect to this density.

The evolution equation for $P_{t}({\sf M})$ can be constructed by considering the addition of a small time interval $\delta t$ with a single weak scatterer (c.f. 3.26) to an existing interval $t$ with $j$ scattering events. It can then be shown that the probability density of the enlarged time interval $P_{t+\delta t}({\sf M}={\sf M}_{2}{\sf M}_{1})$, where ${\sf M}_{1}\equiv{\sf M}(j)$ and ${\sf M}_{2}\equiv{\sf M}_{j+1}$, corresponds to the convolution of the density for the transfer matrix of the extra strip of width $\delta t$: $P_{\delta t}({\sf M}_{2})$, with the density related to $j$ prior scatterings scatterings $P_{t}({\sf M}_{1}={\sf M}_{2}^{-1}{\sf M})$101010For a detailed derivation of the Smoluchowski and Fokker-Planck equations, see [42, 2, 1]

[TABLE]

This integral equation (known as the Smoluchowski or Chapman-Kolmogorov equation) can be equivalently written as the Fokker-Planck equation

[TABLE]

where $\{\lambda_{a}\}$ denote the set of $2{N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ parameters that characterize a general ${N_{\rm{f}}}$-field transfer matrix ${\sf M}$; in the single-field case $\{\lambda_{a}\}=\{n,\phi,\psi\}$. The $\delta\lambda_{a}$ denote the small increment in the parameters $\lambda_{a}$ due to the addition of an extra scattering. The expectation value is over the probability distribution describing the properties of the scatterer in the interval $\delta t$ (which includes location within this interval as well as strength/shape of the scatterers).

In the present work we will focus on the evolution of the occupation number density and the magnitude of the scalar field mode functions. When this is the case, it is convenient to not track the full transfer matrix ${\sf M}$ but its square ${\sf R}(j)\equiv{\sf M}(j){\sf M}^{\dagger}(j)$. This matrix is Hermitian and depends only on ${N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ variables in the general case. For a single field, we can write

[TABLE]

where, somewhat abusing notation, we have defined $\lambda\equiv 2n+1$ and $\tilde{\lambda}\equiv\sqrt{\lambda^{2}-1}$. The Fokker-Planck equation (3.28) would then be written as a three-variable PDE for the probability density $P(\lambda,\phi;t)$.

The re-parametrization in terms of ${\sf R}$ is particularly convenient for our purposes since the perturbations in the parameters $\lambda,\phi$ upon the addition of an extra scattering are known. They were derived in [2] under the assumption of Dirac-delta scatterers with zero-mean uncorrelated amplitudes,

[TABLE]

Since ${\sf R}$ is quadratic in the ${\sf M}_{j}$ matrices, its instantaneous change contains first- and second-order corrections in the order parameter ${m}_{j}a(\tau_{j})|f(\tau_{j})|^{2}$. In terms of $\lambda$ and $\phi$, these are given by

[TABLE]

where the $g,\tilde{g}$ are themselves functions of the parameters $\lambda,\phi$,

[TABLE]

To go further, the explicit expression for the mode-functions are necessary. We will show that, in spite of the apparent complexity of the equations, the equations predict certain universal results, which we turn to after we have discussed our numerical results first. For the moment, note that simplified, general expressions for expectation values of functions of $\lambda,\phi$ can be obtained by integration of the Fokker-Planck equation (as discussed in [2, 1]).

We will return back to this Fokker-Planck equation in Section 4.3 and Section 5.3 and use it to explain aspects of our numerical results for the conformal mass and massless fields.

4 Conformally massive field in de Sitter background

4.1 Numerical results

In this section we focus on numerical results for the evolution of the occupation number, the scalar field amplitude and its phase in the conformal mass case using the transfer matrix approach discussed in Section 3.2. We separate the discussion of our numerical results into four regimes, namely the sub- and super-horizon regimes with weak and strong scattering.

The sub- and super-horizon regimes correspond to physical wavelengths smaller and larger than the horizon scale, $k/aH\gg 1$ and $k/aH\ll 1$, respectively. In an expanding de Sitter background $a(t)=e^{H(t-t_{i})}$, any given comoving wavelength that starts inside the horizon will eventually cross outside the horizon at a time $t_{k}$ satisfying:

[TABLE]

Numerically, this allows us to explore both the sub- and super-horizon regimes for a given Fourier mode by starting the computation for some time $t_{i}<t_{k}$ and finishing it at $t_{f}>t_{k}$; for definiteness, we have considered a total range of 40 Hubble times between the initial and final times, centered at horizon crossing. We will also set $a(t_{i})=1$.

We will talk of weak or strong scattering depending on whether the parameter

[TABLE]

is much smaller or larger than unity. A more careful delineation will be provided later. Here, $N_{s}$ is the number of scatterers in the interval $t_{f}-t_{i}$, $\mathcal{N}_{s}$ denotes the number of scatterers per Hubble time, and $\sigma^{2}=\langle m_{i}^{2}\rangle$ characterizes the strength of the scatterers. As we will demonstrate numerically and analytically below, it is the combination $\mathcal{N}_{s}(\sigma^{2}/H^{2})$ that really determines the growth rate for the occupation number and the field amplitude.

Note that our delineation of strong and weak scattering is different from our perturbativity condition (3.26). For the conformal de Sitter scenario ($M^{2}=2H^{2}$), the perturbativity condition can be rewritten as $(\sigma/k_{\rm phys})^{2}\ll 1$ where $k_{\rm phys}=k/a$. In the subhorizon regime, this condition is always satisfied if $(\sigma/H)^{2}\leq 1$, while outside the horizon the scattering amplitudes must satisfy the much more restrictive constraint $(\sigma/H)^{2}\ll(k/aH)^{2}\ll 1$. It is clear then that, for any given $(\sigma/H)^{2}<1$, perturbativity will be eventually lost outside the horizon. Because of this restriction, we naively expect the Fokker-Planck approach will properly account for the evolution inside the horizon for scattering amplitudes not greater than the Hubble scale, while outside the horizon it will fail unless $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$. Our numerical approach does not require any such restriction.

As we will discuss below, the universality of our results relies on the assumption that the number of scatterers per Hubble time $\mathcal{N}_{s}$ is large (see Fig. 41). This implies that many thousands of operations over hundreds or even thousands of realizations are necessary to reach a stationary regime. This complexity, coupled with exponentially increasing or decreasing quantities, cries out for a numerical code capable of handling the extremely high precision required. To achieve this, we have built our (Fortran) code making extensive use of the thread-safe arbitrary precision package MPFUN-For written by David H. Bailey [48]. We have confirmed that the precision used in our numerical simulations (500 digits) is adequate by ensuring that the constraint on the Bogoliubov coefficients, $|\alpha_{j}|^{2}-|\beta_{j}|^{2}=1$, holds up to the chosen precision for all realizations.

4.1.1 Individual realizations

Fig. 3 shows the evolution of the the field amplitude and its phase, as well as transfer matrix parameters $\{n,\phi,\psi\}$ as functions of time. Here we have assumed that both the amplitudes and the locations of the non-adiabatic events are uniformly distributed in the intervals $m_{j}\in(-\sqrt{3}\sigma,\sqrt{3}\sigma)$ and $\delta t_{j}\in(0,1/H\mathcal{N}_{s})$, where $\delta t_{j}$ denotes the interval between event locations, $\delta t_{j}\equiv t_{j}-t_{j-1}$. For definiteness we have taken $N_{s}=300$ and $k=e^{20}H$. The factor of $e^{20}$ allows for $20$ e-folds before horizon crossing.111111In order to allow a simpler reading of our numerical results, in this figure and all other figures that follow, we present the canonically normalized field re-scaled by its magnitude in the Bunch-Davies vacuum. That is, in all figures the physical value of $X$ can be recovered by taking

$X\;\longrightarrow\;\sqrt{2k}\,X\,.$

(4.3) ,121212Our choice of 20 e-folds after horizon crossing runs afoul of backreaction constraints for strong scattering (see Appendix A.3). Nevertheless, we display our results for ease of comparison with weak scattering.

Each column in Fig. 3 corresponds to a single realization of the disorder $\{(m_{1},t_{1}),(m_{2},t_{2}),\ldots\}$. The different columns correspond to disorder realizations drawn from distributions which correspond to different values of the parameter $\mathcal{N}_{s}(\sigma/H)^{2}$, one in the weak regime (left), one for a “moderate” value (center), and one in the strong scattering regime (right). Despite the fact that the results in Fig. 3 correspond to a single realization of the amplitudes and locations of the non-adiabatic events, we can still read off the main features of the evolution of the parameters of interest, namely

The magnitude of the (re-scaled) canonically normalized field, $|X|$, remains very close to one, with virtually no influence from scattering on subhorizon scales. In other words, $|\chi|=|X|/a$ decreases exponentially with its decay rate determined by the scale factor. Outside the horizon, the magnitude of $X$ remains $\mathcal{O}(1)$ for weak scattering, which implies that the decaying trend for $|\chi|$ is continued after horizon crossing. For moderate scattering, the effect of the stochastic non-adiabaticity is capable of exciting $|X|$ by a couple of orders of magnitude away from its vacuum value, although no clear increasing or decreasing trend is noticeable. Finally, in the case of strong scattering, $|X|$ grows exponentially outside the horizon, with a rate dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$. In the case shown in Fig. 3, this growth is sufficiently large to overcome the decay of $|\chi|$ and to make it grow for $|k\tau|\ll 1$. 2. 2.

The scalar field phase, $\arg X$, is uniformly distributed in the interval $(0,2\pi)$ before horizon crossing. After horizon crossing, if scattering is weak, this phase freezes asymptotically to a small value, $\arg X\ll 1$. If scattering is moderate or strong, the phase becomes almost frozen along a random direction, with $X$ evolving along a ray in the complex plane. Also see Fig. 4 and the discussion after this list of observations. 3. 3.

The occupation number density, $n=|\beta|^{2}$, grows exponentially (for the given Fourier mode). For weak scattering the exponential growth rate is constant throughout the evolution (with $n\propto a^{2}$), while for strong scattering, the rate increases shortly before horizon crossing. The numerical value of these rates depends on $\mathcal{N}_{s}(\sigma/H)^{2}$. Note that for very weak scattering (or when $n\ll 1$ more generally) $n\propto a^{2}$ seems counter-intuitive at first glance. However, this result follows from the observation that each scatterer (as seen in Eq. (3.22)) comes with an increasing strength $a(\tau_{j})m_{j}$ as a function of time. 4. 4.

The transfer matrix phase $\phi$ (c.f. (3.17)) is naturally defined on the domain $(0,\pi)$. The phase varies randomly over this domain inside the horizon, and freezes asymptotically to $\phi\sim\pi/2$ far outside the horizon. 5. 5.

From the statistical point of view, as we will discuss below, the natural range for the second transfer matrix phase $\psi$ (c.f. (3.17)) corresponds to $(-\pi/2,\pi/2)$, where it varies randomly for $|k\tau|\gg 1$. For weak scattering, it freezes to $|\psi|\ll 1$ in super-horizon scales, while for strong scattering it freezes to a seemingly random value.

The curious behavior of $\arg X$ for moderate and strong scattering is displayed in a clearer fashion in Fig. 4. There, the evolution of $X$ is shown in the complex plane for $\mathcal{N}_{s}(\sigma/H)^{2}=1$ and $10^{2}$. It is immediately clear that, as discussed above, the field amplitude is constant in time in subhorizon scales, and the random, but uniformly distributed phase, results in a random walk of the (re-scaled) field $X$ on the unit circle. As $k\sim aH$, the phase of $X$ locks along a random line in the complex plane, and $X$ evolves along this ray; for strong scattering it grows exponentially, jumping between diametrically opposite directions. These diametrical jumps can be understood as follows. In terms of the transfer matrix ${\sf M}$-parameters $\{n,\phi,\psi\}$ defined in (3.17), the spectator field can be written in general as

[TABLE]

Inside the horizon, $n\ll 1$, and $\arg X\simeq-(\phi+\psi+k\tau)\pmod{2\pi}$, which is randomly distributed, not only because $\phi$ and $\psi$ themselves are, but most importantly because $|k\tau|\gg 1$, completely scrambling the phase. Outside the horizon, however, $n\gg 1$, and

[TABLE]

which implies that $\arg\chi$ is mostly determined by $\psi$. However, a curious behavior regarding the sign of $X$ arises due to the asymptotic behavior of $\phi$. As we discussed above, Fig. 3 shows that $\phi\rightarrow\pi/2$ as $|k\tau|\rightarrow 0$. Moreover, Fig. 5 demonstrates that the argument of the cosine in (4.5) is driven exponentially fast in cosmic time towards $\pi/2$, alternating signs randomly. Straightforward expansion implies then that

[TABLE]

with $\zeta=0\,\,{\rm or}\,\,1$ after each scattering.

4.1.2 Means and variances

In the previous subsection we discussed the evolution of the transfer matrix parameters $\{n,\phi,\psi\}$ and the scalar field amplitude $|X|^{2}$ and its argument for particular realizations of the locations and amplitudes of the scattering events. We now turn to the description of the dynamics of the system given an ensemble of realizations of the scatterers. In this section we will discuss the evolution of the lowest moments of the angles $\phi$ and $\psi$, as well as those for $\ln(1+n)$ and $\ln|X|^{2}$; in Section 4.1.3 we study the form of their probability distributions. Note the focus on logarithms of $n$ and $|X|$ is related to the observation that both $n$ and $|X|$ show an exponential behavior with cosmic time. We will also find that $\ln|X|^{2}$ is normally distributed both inside and outside the horizon (for any strength of scattering), making it a simpler variable to work with.

Figs. 6 and 7 show the dependence of the mean and variance of $\ln(1+n)$, $\ln|X|^{2}$, $\phi$ and $\psi$ on the scattering-strength parameter $\mathcal{N}_{s}(\sigma/H)^{2}$ and the wavenumber $k$, respectively. For Fig. 6, we have set $k/H=e^{20}$ for definiteness, while for Fig. 7 we have fixed $\mathcal{N}_{s}(\sigma/H)^{2}=60$. In both cases we consider $N_{s}=2000$ non-adiabatic events, drawn from an ensemble of 2000 members for the scattering locations and strengths $m_{j}$ and $\delta t_{j}$. For simplicity we have assumed that the amplitudes and locations of the scatterings are drawn from uniform distributions, as in Fig. 3. Nevertheless, we will show in Appendix A.2 that the results discussed here are not sensitive to the ensemble distributions, provided that both the $m_{j}$ and $\delta t_{j}$ are random.

Moments of $\ln|X|^{2}$ on sub-horizon scales: In the first row of Figs. 6 and 7, the panels show the evolution of the moments for the (logarithmic) scalar field amplitude. The mean $\langle\ln|X|^{2}\rangle$ is approximately constant on sub-horizon scales, and independent of the value of $\mathcal{N}_{s}(\sigma/H)^{2}$. In light of (4.3), we can then simply write

[TABLE]

Although it is not obvious from the figure, unlike the mean, the variance of $\ln|X|^{2}$ increases exponentially with cosmic time, with a rate that is independent of $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$. Its functional dependence can be approximated as

[TABLE]

where $k_{\rm phys}=k/a$. Note that while the mean of $\ln|X|^{2}$ remains similar to its value in the vacuum, the variance is growing $\propto a^{2}$.

Moments of $\ln|X|^{2}$ on super-horizon scales: In the case of very weak scattering, $\ln|X|^{2}$ continues to be approximately constant as it is inside the horizon; it is only for strong scattering that the field can overcome the expansion and either decay or grow. Moreover, this rate is also independent of $k$. The variance of $\ln|X|^{2}$ is also a linear function of cosmic time for $|k\tau|\ll 1$, with a rate that is independent of the wavenumber $k$. We can then write

[TABLE]

where $\mu_{1}$ and $\mu_{2}$ are functions of $\mathcal{N}_{s}(\sigma/H)^{2}$ and are shown in Fig. 8. The curves displayed in the figure are the result of a linear fit to the averaged moments over 400 realizations in the super-horizon regime.

For $\mathcal{N}_{s}(\sigma/H)^{2}\ll 10^{-2}$, $\partial_{Ht}\langle\ln|X|^{2}\rangle\sim 0$. Along with Fig. 8, this feature is also evident in the top left panel of Fig. 6. For $10^{-2}\lesssim\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 2$, the rate of change of $\langle\ln|X|^{2}\rangle$ is negative (notice the dip in the top panel of Fig. 8). This peculiar behavior is also demonstrated by the green curve in the third panel on the right of Fig. 6. In other words, $\langle\ln|\chi|^{2}\rangle$ decays faster than in the vacuum (recall that $X=a\chi$) in this regime of $\mathcal{N}_{s}(\sigma/H)^{2}$. If scattering is stronger, this rate of decay is smaller, until it vanishes for $\mathcal{N}_{s}(\sigma/H)^{2}\simeq 60$. When this is the case, $\langle\ln|\chi|^{2}\rangle\sim{\rm const.}$, or equivalently $\langle\ln|X|^{2}\rangle\sim a^{2}$ while scatterings continue taking place. For $\mathcal{N}_{s}(\sigma/H)^{2}\gg 1$, the field grows exponentially, with rate $\mu_{1}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.3}-2$ (top panel of see Fig. 8).

The dependence of ${\mu}_{2}=\partial_{Ht}{\rm Var}[\ln|X|^{2}]$ on $\mathcal{N}_{s}(\sigma^{2}/H^{2})$ is shown in the lower panel of Fig. 8. A remarkable feature of this variation is the sharp distinction in the evolution for weak and strong scattering. In the former case, the variance grows with $\mu_{2}\sim\mathcal{N}_{s}(\sigma/H)^{2}$, while in the later it grows as $\mu_{2}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.33}$.

Moments of $\ln(1+n)$ on sub-horizon scales: The second row of Fig. 6 shows the evolution of $\langle\ln(1+n)\rangle$ and Var$[\ln(1+n)]$ in a log-scale, both showing an exponential growth ($\propto a^{2}$) in the sub-horizon regime. The growth rates are independent of $\mathcal{N}_{s}(\sigma/H)^{2}$, but the absolute magnitude of the moments depends on it. Notice that this trend is maintained until $\langle\ln(1+n)\rangle\sim\mathcal{O}(1)$, which occurs shortly after horizon crossing for weak scattering and before for strong scattering. The second row of Fig. 7 further shows that while the mean and the variance on sub-horizon scales depend on the value of $k$, their growth rates are independent of $k$. These observations lead us to write:

[TABLE]

which we derive analytically in Section 4.3.1. Note that the exponential growth of $\langle\ln(1+n)\rangle$ and Var$[\ln(1+n)]$ comes from $k_{\rm phys}=k/a$ where $a=e^{H(t-t_{i})}$. At first sight this seems inconsistent with the results from earlier papers by some of us [1, 2] which calculated particle production in a non-expanding universe and found that $\ln(1+n)$ grows linearly with time. On sufficiently sub-horizon scales, one could expect the above result and the non-expanding case to agree. A closer look reveals that there is no inconsistency. The linear growth of $\ln(1+n)$ was really true for $n\gg 1$, whereas above we have $n\ll 1$. In detail, for $n\ll 1$, since $\ln(1+n)\sim n$, both $\ln(1+n)$ and $n$ grow exponentially with cosmic time. When $n\gg 1$, find that $\ln(1+n)\sim\ln(n)$ will grow linearly with cosmic time, whereas $n$ will grow exponentially.

Moments of $\ln(1+n)$ on super-horizon scales: Third row from the top in Figs. 6 and 7, the left and right panels show the evolution of the moments of the occupation number, but now in a linear scale to demonstrate the linear increase of $\langle\ln(1+n)\rangle\simeq\langle\ln(n)\rangle$ for $|k\tau|\ll 1$ (ie. on super-horizon scales). In this regime the growth rate is clearly dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$, being steeper for strong scattering, and it is seemingly independent of the wavenumber $k$. In analogy with $\ln|X|^{2}$, we define these super-horizon growth rates of $\langle\ln(1+n)\rangle$ and Var$[\ln(1+n)]$ as follows:

[TABLE]

The functions $\tilde{\mu}_{1}$ and $\tilde{\mu}_{2}$ are also shown in Fig. 8. As seen in this figure, these rates are closely connected with the corresponding rates for $\ln|X|^{2}$. Explicitly, $\tilde{\mu}_{1}\simeq\mu_{1}+2$ everywhere, whereas $\tilde{\mu}_{2}\simeq\mu_{2}$ for $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 10^{-1}$.

For $\mathcal{N}_{s}(\sigma/H)^{2}\leq 10^{-2}$, the growth rate of the mean is constant, $\partial_{Ht}\langle\ln(1+n)\rangle\simeq 2$. For the occupation number, this implies that $n\sim a^{2}$ at all times when scattering is weak. We will revisit this result with the Fokker-Planck formalism in Section 4.3. For $10^{-2}\lesssim\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 2$, the typical occupation number grows at a slower rate. In the strong scattering regime, the mean grows in a power-like fashion with the scattering strength parameter, $\tilde{\mu}_{1}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.3}$.

The rate of growth with time for the variance of $\ln(1+n)$ is shown in the lower panel of Fig. 8. For $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$, it is approximately constant, $\tilde{\mu}_{2}\sim 0.025$, but it rises sharply as the scattering parameter increases. For $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 10^{-1}$ the rate follows a power-law dependence, $\tilde{\mu}_{2}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.33}$, slightly steeper than that of the mean.

It is important to note that the rapid growth of the variance in sub- and super-horizon scales (for both $\ln(1+n)$ and $\ln|X|^{2}$) sheds some doubt on our characterization of their means corresponding to the “most probable” member of the ensemble of realizations. We dismiss these concerns in detail in Appendix A.1, by constructing ratios of means and standard deviations of these quantities and showing that while the standard deviations grow, the means grow even faster.

Moments of $\phi$ on sub-horizon scales: The evolution of the angular parameter $\phi$ is shown in the fourth row in Figs. 6 and 7. Deep inside the horizon we find

[TABLE]

values which are consistent with a uniformly distributed random variable in $(0,\pi)$.

Moments of $\phi$ on super-horizon scales: As the mode leaves the horizon, both the mean and variance oscillate about these values, with the amplitude of these oscillations being dependent on the scattering strength parameter. Once the mode is far outside the horizon, the oscillations stop, and the moments settle down to

[TABLE]

with an exponentially decreasing variance. Numerically we find that the final value (at $H(t-t_{k})=20$) of the mean of $\phi$ is equal to $\pi/2$ for all $\mathcal{N}_{s}(\sigma/H)^{2}$ up to a numerical error smaller than one part in $10^{-6}$. We also find for the time rate of the log of the variance that

[TABLE]

for any scattering strength, up to a $\lesssim 8\%$ deviation that lacks a simple dependence on $\mathcal{N}_{s}(\sigma/H)^{2}$.

Moments of $\psi$ on sub-horizon scales: Finally, the time-dependence of the moments of $\psi$ are shown in the bottom left and right panels of Figs. 6 and 7. Similarly to $\phi$, these results are consistent with a uniformly distributed random variable inside the horizon (on the interval $(-\pi/2,\pi/2)$):

[TABLE]

Also similar to the $\phi$ case, the mean and variance of $\psi$ appears to be perturbed away from these values during horizon crossing.

Moments of $\psi$ on super-horizon scales: In the super-horizon regime, both moments of $\psi$ appear to asymptote to values dependent on the scattering strength parameter. From Fig. 6 it is apparent that $\psi$ retains its uniform distribution for strong scattering, while for weak scattering the moments are consistent with a narrow probability density centered at $\psi\simeq 0$. Numerically we find that a good approximation to the lowest moments of $\psi$ is given by

[TABLE]

To arrive to the previous expressions we have ignored a mild but complicated dependence on the scattering parameter for $\langle\psi\rangle$. The maximum deviation is found at $\mathcal{N}_{s}(\sigma/H)^{2}\sim 1$, for which $|\langle\psi\rangle|\lesssim 0.1$. Note that the variance freezes at small values in the case of weak scattering, while for strong scattering it freezes with the same value as in (4.15), consistent with a non-evolving probability distribution.

In our previous discussion we have mostly focused on the super-horizon behavior of the transfer matrix parameters under the assumption that this late-time evolution is controlled only by the scattering strength parameter $\mathcal{N}_{s}(\sigma/H)^{2}$. Furthermore, all the results previously presented assume an underlying uniform distribution for the strength and location of the non-adiabatic events that drive particle creation. Nevertheless, it can be shown that these results are unchanged if the previously mentioned assumptions are broken, as long as the density of scatterers $\mathcal{N}_{s}$ is sufficiently high. See Appendix A.2 for details.

4.1.3 Probability densities

In the previous section we have described the sub- and super-horizon evolution of the lowest moments of $\ln|X|^{2}$ and the transfer matrix parameters. In this section we now study in a mostly qualitative fashion the form and dynamics of the full probability density functions (pdf) for these random variables.

Figs. 10-12 display snapshots of the time evolution of the instantaneous normalized pdfs for the field and ${\sf M}$-parameters for selected values of $\mathcal{N}_{s}(\sigma/H)^{2}$. In all cases we have considered $k/H=e^{20}$ and $N_{s}=2000$ with $N_{r}=5000$ realizations. As before, the total evaluation time interval corresponds to $H(t_{i}-t_{f})=40$. The pdfs are built using a Gaussian kernel density estimator of variable bin size. For the angular variables $\phi$ and $\psi$, the data has been extended periodically in the cases where the pdf support is of size $\pi$, to minimize edge effects.

Pdf for $\ln|X|^{2}$: Fig. 9 shows the pdfs for the logarithmic field amplitude $\ln|X|^{2}$. The description of the distribution is exceptionally simple: for all times and values of $\mathcal{N}_{s}(\sigma/H)^{2}$ a normal pdf is a good fit for the data. We prove this fact analytically for sub-horizon modes in Section 4.3.1 (see Fig. 18). Far outside the horizon, at $H(t-t_{k})=20$, the normal form of the pdf is demonstrated in Fig. 13, where a Gaussian fit is superimposed for small and large values of the scattering strength. We can therefore conclude that:

The squared-field amplitude, $|X|^{2}$, follows a log-normal distribution both inside and outside the horizon, for weak and strong scattering.

Pdf for $\ln(1+n)$: The instantaneous pdfs for $\ln(1+n)$ are shown in Fig. 10. In the first two columns from the left, turning a blind eye to the axis tick values, it is clear that this pdf “flows” in a way almost independent of the strength of the scattering. At very early times it starts with a highly right-skewed, almost exponential shape.131313This is consistent with a coefficient of variation $\tau_{\ln(1+n)}\simeq 1$ (see Appendix A.1). In Section 4.3.1 we will confirm this fact analytically (see Fig. 17). As time increases, the position of the maximum increases, together with the width of the distribution, maintaining its shape in all cases save for the strongest scattering case, where the shape is now slightly distorted. It is not until the mode is stretched to super-horizon scales that the difference between weak and strong scattering is evident.

The middle column of Fig. 10 shows the transition regime, and demonstrates the delay in evolution of the weak scattering case compared to the strong scattering ones. As it is clearly exhibited by the $\mathcal{N}_{s}(\sigma/H)^{2}=10^{-2}$ case, the pdf shifts from the left-lobed exponential-like distribution to a right- or center-lobed normal-like distribution. The last two columns of the figure show finally the pdf outside the horizon. In all cases the distribution has a marked peak, with the weak cases retaining a significant tail of realizations with low occupation numbers, while the strong cases show symmetry with respect to the peak of the distribution.

In Fig. 14 we show four of the five panels of the last column of Fig. 10, with a skew-normal fit to the data shown as a black dashed curve. As a reminder, a random variable $x$ is skew-normal distributed if its pdf is given by

[TABLE]

where $x_{0}$, $\omega$ and $\alpha$ denote the location, scale and shape parameters, respectively [49, 50]. As it is clear in the top panels of the figure, the distributions for low values of $\mathcal{N}_{s}(\sigma/H)^{2}$ show significant skewness, deviating noticeably from normality. In contrast, large values of the scattering strength parameter exhibit a Gaussian shape, with $\alpha=0$ and $x_{0}$ and $\omega$ given by the mean and variance of the discrete data for the fit. Hence, the occupation number $n$ will be log-skew-normally distributed in general, with a decreasing skewness for increasing $\mathcal{N}_{s}(\sigma/H)^{2}$.

Pdf for $\phi$: The $\phi$-distribution is shown in Fig. 11. This pdf presents a non-trivial evolution, strongly dependent on time and $\mathcal{N}_{s}(\sigma/H)^{2}$. As anticipated in Sections 4.1.1 and 4.1.2, the distribution of $\phi$ is approximately uniform in sub-horizon scales. Although some structure is visible for some pdfs (at $H(t-t_{k})=-12$ with $\mathcal{N}_{s}(\sigma/H)^{2}=10^{-2}$, for example), we believe that upon increasing the number of realizations any features will be mostly smoothed out. Note that a uniform distribution is obtained independent of the strength of scattering.

The central column of Fig. 11 shows the transition forms of $\mathcal{P}[\phi]$. As the mode leaves the horizon the uniformity of the pdf is lost, and a two-lobed distribution arises, with the lobes being of approximately the same size and located symmetrically with respect to $\phi=\pi/2$. As time increases, these lobes approach, but never fully merge. This is expected from the “jumping” behavior of $\phi$ with respect to $\pi/2$ shown earlier in Fig. 5.

The super-horizon form of the $\phi$ pdf is shown in the two rightmost columns of Fig. 11. In these cases, the maxima of the two lobes approach $\pi/2$ exponentially fast, while their widths also decrease exponentially. Clearly the rates of approach and narrowing are dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$; for strong scattering the rates are so high that our pdf estimator is not capable of showing clearly the structure of the distribution. In Section 4.3.2 we discuss an analytical approximation to the super-horizon evolution of $\mathcal{P}[\phi]$ in terms of a Fokker-Planck equation (see Fig. 20) which captures this behaviour.

Pdf for $\psi$: Fig. 12 shows the evolution of the probability distribution for the transfer matrix phase $\psi$. In analogy with the pdf for $\phi$, for $|k\tau|\gg 1$ the distribution of $\psi$ is uniform for any scattering strength, albeit over the interval $(-\pi/2,\pi/2)$. We also observe some features on the pdfs, but we believe that they are mostly an artifact of our finite ensemble of realizations. For strong scattering, the uniformity of the distribution is preserved into super-horizon scales, where a frozen pdf is evident in the last two rows of the figure in question. This is consistent with (4.16). For weak scattering, the distribution becomes two-lobed as the mode leaves the horizon. However, unlike the $\phi$ case, these two lobes merge in a finite time around $\psi\simeq 0$, and lead to a peaked, frozen distribution for $|k\tau|\ll 1$. Fig. 15 shows the four upper right panels of Fig. 12 compared to a normal distribution of zero mean and variance $(\pi/2)\mathcal{N}_{s}(\sigma/H)^{2}$, as per (4.16). This clearly shows that for weak scattering, $\psi$ is normally distributed outside the horizon. Finally, the central row of Fig. 12 shows a super-horizon pdf intermediate between a uniform and a normal distribution, clearly dependent on the value $\mathcal{N}_{s}(\sigma/H)^{2}$.

4.2 The field two-point function

Arguably, the most remarkable result from our previous numerical explorations consists in the fact that the spectator field amplitude $|X|^{2}$ is lognormally distributed at all times for any scattering strength. Moreover, outside the horizon, the one-point pdf of $\ln|X|^{2}$ possesses a mean and a variance that increases linearly with cosmic time; eqs. (4.9) may be rewritten as

[TABLE]

where $t_{k}$ is the time of horizon-crossing for the given mode, and we have restored for convenience the momentum-dependence of the mode function. A normal one-point pdf with mean and variance that linearly grow with time are characteristic features of Brownian motion (random walk, Wiener processes) with drift [51]. Also characteristic of Wiener processes is the property that the unequal time two point correlation function is linearly proportional to the smaller of the two times. We check this property below.

In order to compute the two-point function of the logarithm of the field amplitude, let us define the driftless (zero mean) variable

[TABLE]

We then define the expectation value

[TABLE]

Fig. 16 shows the time-dependence ($t$) of the two-point function in two cases.

Unequal time: The left column corresponds to considering the equal momenta $(k^{\prime}=k)$, unequal time ($t\neq t^{\prime}$) scenario for discrete values of $t^{\prime}$ and $\mathcal{N}_{s}(\sigma/H)^{2}$. For all three values of $\mathcal{N}_{s}(\sigma/H)^{2}$, the qualitative behavior of the curves shown in similar. For $t<t_{k}$, the magnitude of the two-point function is negligible for any $t^{\prime}$. If $t^{\prime}\lesssim t_{k}$, this non-growing trend is preserved after the mode in question crosses outside the horizon, as demonstrated by the purple and blue curves. Assuming now that $t^{\prime}>t_{k}$ (green, orange and red curves), we observe that the two point function grows at the same rate as the variance of $\ln|X_{k}|^{2}$ does (c.f. Section 4.1.2) for $t_{k}<t<t^{\prime}$, and is frozen at its value at $t=t^{\prime}$ for $t>t^{\prime}$. In summary,

[TABLE]

indicating that we have an approximately Wiener process on super-horizon scales.

It is worth noting that for fixed $k$, $Z_{k}(t)$ describes a Gaussian process and therefore all its higher point correlation functions may be computed in terms of its two-point function (4.21),

[TABLE]

Unequal Momenta: The right column of Fig. 16 corresponds to the equal time but unequal momenta case. Note here that for the purple and blue curves, for which $k^{\prime}\leq k$, the two-point function grows linearly with time for $t>t_{k}$. For the green, orange and red curves, which correspond to $k^{\prime}>k$, the two-point function grows only for $t>t_{k^{\prime}}$. We therefore conclude that

[TABLE]

which again confirms our expectation for a Wiener process when both modes are super-horizon.

Unequal Momenta and Time: For the general case of unequal time and unequal momenta correlators, we have thus found that the equation

[TABLE]

is a good approximation.

The above results show that $\ln|X_{k}|^{2}$ satisfies the auto-correlation properties of a Brownian motion with drift for super-horizon $k$. Hence, the field amplitude $|X_{k}|^{2}$ (and consequently $|\chi_{k}|^{2}$) describes a geometric Brownian motion with drift in cosmic time as long as modes are outside the horizon. For unequal momenta, when either of the modes is inside the horizon, the correlation is vanishing for unequal times.

It is worth noting that, in terms of the two-point function (4.24), the $n$-point correlation function for the squared field magnitude can be written in general as follows,

[TABLE]

This result follows trivially from the lognormality of $|X_{k}|^{2}$.

4.3 Analytical results

Up to this point we have only discussed the numerically-obtained trends and values for the transfer matrix parameters and the scalar field magnitude, without referring to analytical expectations. We have decided to follow this “inverted” program because the numerical results compose an almost complete picture that will not be attainable with the analytical tools at our disposal. In particular, our discussion following Eq. (4.2) suggests that only the very weak scattering regime $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$ can be reliably probed analytically using the formalism laid out in Section 3.3. As an example, obtaining the precise functional form for the functions $\mu_{1,2}$ and $\tilde{\mu}_{1,2}$ is beyond our present study. Nevertheless, as we will show, it is possible to confirm the functional dependence of some moments on $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$ for weak scattering, as well as that of the probability densities. In the next two sections we derive the form of the Fokker-Planck equation (3.28) which corresponds to a conformally-massive scalar field in a de Sitter expanding background in the sub- and super-horizon regimes and we use it to derive analytical and semi-analytical expressions for moments and pdfs. The main results in this section include: (i) the full pdf of $\phi$, $\ln(1+n)$ and $\ln|X|^{2}$ on sub-horizon scales, and (ii) the rather non-trivial pdf of $\phi$ and the time-evolution rates for $\langle\ln(1+n)\rangle$ and $\langle\ln|X|^{2}\rangle$ on super-horizon scales. These results are consistent with our numerical investigations.

4.3.1 Sufficiently sub-horizon

The general form for the coefficients for the two-parameter FP equation was derived in Section 3.3, under the assumption of Dirac-delta scatterers with uncorrelated amplitudes of vanishing mean. For a conformally massive scalar field in an expanding de Sitter background, the coefficient functions (3.32a)-(3.32d) can be written as follows,

[TABLE]

where $\lambda=2n+1$, $\tilde{\lambda}=\sqrt{\lambda^{2}-1}$. With these expressions at hand we can compute the coefficients of the FP equation (3.28). By means of an example we will be able to find a pattern that will allow us to bypass the need to compute the disorder averages for all coefficients. From (3.31a) it follows that

[TABLE]

Let us evaluate the first term of the previous expression in full detail,

[TABLE]

Here in the second line we have made the variable change $t_{j}\rightarrow y=k\tau_{j}$, we have recalled the definition for $\sigma^{2}=\langle m_{j}^{2}\rangle$ in (3.30), and we used $\delta t=a(\tau)\delta\tau$. Note that the time interval over which the disorder average is taken should be at most of the order of the separation between the non-adiabatic events, $\delta t\lesssim\delta t_{i}$. Assuming for simplicity that the scattering locations are uniformly distributed in cosmic time, we can identify

[TABLE]

where we again recall that $\mathcal{N}_{s}$ is the number of scatterers per Hubble time. This implies that the parameter delineating different regimes for the coefficients of the FP equation in the previous calculation corresponds roughly to the ratio of the physical wavenumber to the Hubble scale, weighed by the density of scatterers per Hubble time,

[TABLE]

In the super-horizon regime in (4.28), $|k\delta\tau|\ll 1$, the disorder average is equivalent to simply evaluating the coefficient function at $\tau=\tau_{j}$; this is to be expected as the time period of oscillations of the mode function is larger than the mean free path determined by the separation between scattering events. On the other side, deep inside the horizon, the mode function oscillates a large number of times in between events, resulting in a vanishing expectation value. We can generalize this result for any non-oscillatory (e.g. polynomial) function $g(\tau)$ as follows:

[TABLE]

With the previous result at hand, we can immediately write the full set of correlators for the FP equation in the deep sub-horizon regime,

[TABLE]

In a completely analogous manner to the non-expanding scenario, all the expectation values are independent of the angular variable [1, 2]. Therefore, we can immediately conclude that the probability density $P$ which is a solution to the FP equation (3.28) is independent of $\phi$, or equivalently,

$\phi$ is uniformly distributed deep inside the horizon.

Moments of $\ln(1+n)$: Before attempting to solve the FP equation, let us consider the expectation value $\langle\ln(1+n)\rangle$. The multiplication of (3.28) by $\ln(1+n)$ and integration with respect to both $\lambda$ and $\phi$ leads to the expression

[TABLE]

Using (4.2), (4.29) and (4.33a)-(4.33e) we can immediately rewrite the above equation as follows:141414Note that this expression is consistent with the expectation from a non-expanding universe [1, 2] where the right hand side of the above equation was $\sigma^{2}/k^{2}$). The extra factor of $4$ is explained by a slight change in the definition of $\sigma^{2}$, whereas the appearance of $\mathcal{N}_{s}$ is related to the choice of time variable.

[TABLE]

Integration with respect to time leads to $\langle\ln(1+n)\rangle\simeq\mathcal{N}_{s}(\sigma/k_{\rm phys})^{2}/8$ (c.f. 4.10), as verified via numerical simulations. Note that no assumptions regarding the magnitude of $n$ have been made to derive this result.

We can repeat this exercise for the variance if we multiply the general FP Eq. (3.28) with $[\ln(1+n)]^{2}$ and integrate over $\phi$, and $\lambda$, obtaining

[TABLE]

where in the third line we have approximated $\langle\ln(1+n)\rangle\simeq\langle n\rangle$ deep inside the horizon, and we have used (4.35) for its value. The previous expression can be integrated to give

[TABLE]

which reproduces the numerically obtained result $\langle\ln(1+n)\rangle\simeq{\rm Var}[\ln(1+n)]^{1/2}$ (c.f. 4.10).

Pdf for $\ln(1+n)$: Let us now consider the $\phi$-independent FP equation. Upon substitution of (4.2), (4.29) and (4.33a)-(4.33e) in the general FP Eq. (3.28), we have

[TABLE]

In terms of the new variables

[TABLE]

it can be rewritten as

[TABLE]

which has the integral-form solution [52]

[TABLE]

We can find an approximate expression for this probability density in the deep sub-horizon regime ($\xi\ll 1)$, if we re-write it in terms of $\ln\rho=\ln(1+n)\ll 1$,

[TABLE]

where $\xi=\mathcal{N}_{s}(\sigma^{2}/k_{\rm phys})^{2}/8\propto a^{2}$. From this distribution it is straightforward to verify that $\langle\ln(1+n)\rangle=\left({\rm Var}[\ln(1+n)]\right)^{1/2}=\xi$, consistent with the previously derived and numerically verified results. Moreover, the exponential form of the pdf is compatible with the correlation coefficient $\tau_{\ln(1+n)}\simeq 1$, computed numerically in Appendix A.1. Fig. 17 shows the agreement between the expression (4.42) and four selected panels from Fig. 10.

Pdf for $\ln|X|^{2}$: With the pdf for $\phi$ and $n$ at hand, we can now compute the corresponding pdf for the scalar field amplitude inside the horizon. Starting from (4.4) we can write

[TABLE]

As $n\ll 1$, we can approximate the logarithm of the amplitude as

[TABLE]

and therefore,

[TABLE]

where recall that $\xi=\mathcal{N}_{s}(\sigma^{2}/k_{\rm phys})^{2}/8\propto a^{2}$. Consistent with our numerical exploration of section 4.1.3, we have found that $\ln|X|^{2}$ is normally distributed, or equivalently, $|X|^{2}$ is log-normally distributed inside the horizon. The pdf (4.45) immediately implies that $\langle\ln|X|^{2}\rangle\simeq-\ln(2k)$ and ${\rm Var}\,[\ln|\chi|^{2}]\simeq 2\xi$, in agreement with the numerical fits (4.7) and (4.8). Fig. 18 further shows the agreement between (4.45) and four selected weak scattering panels from Fig. 9.

4.3.2 Outside the horizon

Let us now find the form of the FP equation far outside the horizon. In this case, the disorder average for the FP coefficients must be taken as in (4.32), where averages are replaced essentially by their instantaneous values. In this late-time limit, we will assume for simplicity that the occupation number has grown sufficiently so that the approximation $\tilde{\lambda}\approx\lambda\approx 2n$ is valid (recall that $\lambda\equiv 2n+1$ and $\tilde{\lambda}\equiv\sqrt{\lambda^{2}-1}$). After some algebra we obtain the following set of correlators,

[TABLE]

Pdf for $\phi$: In terms of the shifted variable

[TABLE]

the FP equation (3.28) takes then the form

[TABLE]

We do not attempt to find a closed-form solution to this equation. Nevertheless, we can find an approximate expression for the time-dependent marginal probability distribution

[TABLE]

which in turn will allow us to calculate the mean particle production rate. Integrating both sides of (4.48) with respect to $\lambda$, and re-parametrizing the time-dependence in terms of $\xi$ defined in (4.39), we obtain the following expression for the equation of motion of $w$,

[TABLE]

In the very-late time limit one could naively expect that the temporal dependence is negligible, and the marginal distribution tends to a limiting pdf. Were this the case, the FP equation for this limit distribution would have the form

[TABLE]

which has the solution $w(\phi)\propto\sec^{2}\phi$. However, this function is divergent at $\phi=\pi/2$ and it is not normalizable, implying that the time dependence in (4.50) cannot be outright disregarded. Nevertheless, as it turns out, this solution correctly describes the qualitative behavior of the marginal distribution at late times, save for a time-dependent cutoff of the divergence at $\pi/2$. Fig. 19 shows the numerical solution of the FP equation (4.50), in solid curves, for different values of the temporal parameter $\xi$, assuming the initial condition $w(\phi,0)=1/\pi$, i.e. the sub-horizon uniform distribution. The dashed black curve is given by the approximation

[TABLE]

where we find the cutoff $\delta$ to be approximately given by

[TABLE]

Fig. 20 shows a comparison between the numerical solution of the marginal FP equation (4.50) and the fully numerically calculated pdf for $\phi$ for selected time slices, with weak scattering $\mathcal{N}_{s}(\sigma/H)^{2}=10^{-4}$ (c.f. Section 4.1.3). The agreement between both results is clear, and it improves as we move farther outside the horizon. From the approximation (4.52) we also obtain

[TABLE]

which agree with the numerical results (4.13), (4.14).

Rates $\langle\ln(1+n)\rangle$ and $\langle\ln|X|^{2}\rangle$: Given the list of correlators (4.46a)-(4.46e), we can calculate the expectation value for the particle production rate for $|k\tau|\ll 1$. Substitution into (4.34) gives

[TABLE]

At very late times $\delta\ll 1$, hence we recover the result $\partial_{Ht}\langle\ln(1+n)\rangle\simeq 2$ outside the horizon for weak scattering (c.f. Fig. 8). Moreover, from $|X|=\sqrt{2n/k}|\varphi-\pi/2|$ (see Eq. (4.6)), we can write

[TABLE]

Integration using the approximation (4.52) yields $\langle\ln|\varphi-\pi/2|^{2}\rangle\simeq 2\ln\delta+2$. Therefore, in the super-horizon regime (where $\delta\ll 1$) from Eq. (4.55) and (4.56), we have

[TABLE]

in agreement with the numerical result shown in Fig. 8 in the weak scattering limit.

5 Massless field in de Sitter background

5.1 Numerical results

We now turn to the discussion of the numerical results for the massless case. This analysis will mirror our previous study for the conformal case: we consider the sub- and super-horizon regimes with weak and strong scattering, where the former are explored by computing the evolution of a Fourier mode over a range of 40 Hubble times centered at horizon crossing, while the later are defined depending on the magnitude of the scattering strength parameter $\mathcal{N}_{s}(\sigma/H)^{2}$ defined in (4.2). A straightforward substitution of the mode functions for the free massless field (3.14) into (3.25) provides the instantaneous transfer matrix used to the derive the results discussed below.

5.1.1 Individual realizations

Fig. 21 shows the evolution of the field amplitude and its phase, as well as transfer matrix parameters $\{n,\phi,\psi\}$, as functions of time. All assumptions on the scattering parameters coincide with those for the conformal mass case discussed in Section 4.1.1: the amplitudes and the locations of the non-adiabatic events are uniformly distributed in the intervals $m_{j}\in(-\sqrt{3}\sigma,\sqrt{3}\sigma)$ and $\delta t_{j}\in(0,1/H\mathcal{N}_{s})$, with $N_{s}=300$ and $k=e^{20}H$, and to avoid cumbersome notation we have also adopted the re-scaling convention (4.3). Each plot corresponds to a single realization of the disorder for the same three scattering strength parameters as in Fig. 3. The features that one can read from these results are similar to those for the conformal case, except for a few key differences,

The magnitude of the canonically normalized field $X$ is constant with virtually no spread in subhorizon scales. For $|\chi|$, this implies an exponential decrease with the rate determined by the inverse of the scale factor. Outside the horizon, $|X|$ grows exponentially. In the case of weak scattering, the growth rate is to a good approximation exactly that given by the scale factor; equivalently, $|\chi|$ is frozen to a constant value after horizon crossing, as expected from the mode function (3.14). For moderate scattering, $|X|$ clearly grows at a slightly slower rate than $a$, signifying an exponential decrease for $|\chi|$. In turn, for strong scattering, the rate of growth of $|X|$ is significantly larger than that for the scale factor, leading to the exponential increase of $|\chi|$. 2. 2.

The behavior of the field phase is analogous to that of a conformally massive field, it fluctuates uniformly in $(0,2\pi)$ for $|k\tau|>1$, and it is frozen at $\pi/2$ with $|k\tau|<1$ for weak scattering, or along a random direction for strong scattering. 3. 3.

The occupation number grows exponentially. In all cases the growth rate increases approximately at horizon crossing, with a value dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$. 4. 4.

Unlike the conformal case, the angular parameter $\phi$ is defined here on the domain $(-\pi/2,\pi/2)$, over which it fluctuates randomly inside the horizon. It freezes asymptotically to $\phi\simeq 0$ far outside the horizon. 5. 5.

In the massless scenario, the natural domain for the phase $\psi$ is $(0,\pi)$. This angular parameter varies randomly over its whole range for $|k\tau|\gg 1$. When scattering is weak, $\psi\rightarrow\pi/2$ outside the horizon, while for strong scattering $\psi$ freezes to a random value.

Note that, qualitatively, in the strong scattering limit $\arg X$ behaves in the same manner as in the conformal case, being locked along a ray in the complex plane and jumping between diametrically opposite directions after the mode in question has left the horizon. This result is more clearly seen in Fig. 22, which shows the evolution of $X$ in the complex plane with a color coded time dependence. Analytically, in terms of the transfer matrix parameters, the massless scalar can be written as

[TABLE]

Deep inside the horizon, with $|k\tau|\gg 1$ and $n\ll 1$, $\arg X\simeq-(\phi+\psi+k\tau)\pmod{2\pi}$, implying a random uniform distribution for the field phase. Conversely, outside the horizon,

[TABLE]

showing that $\arg X$ is determined by $\psi$ up to a sign, determined by the asymptotic value of $\phi$. Fig. 23 shows an enhancement around $\phi=0$ of the $\phi$-panels in Fig. 21. Is is clear that, with $|k\tau|\rightarrow 0$, the argument of the sine in (5.2) will change signs as it is driven to zero. Expanding around this value we can then write

[TABLE]

with $\zeta=\{0,1\}$ randomly, leading to the diametrical flip of the field phase.

5.1.2 Means and variances

Let us now discuss the dynamics of the moments of the field and transfer matrix parameters given an ensemble of realizations. As we did in Section 4.1.2 in the conformal case, we will focus on the lowest moments for $\ln|X|^{2}$, $\ln(1+n)$, $\phi$ and $\psi$, leaving the discussion of their probability distributions to the next section.

The dependence of the means and variances of $\ln|X|^{2}$, $\ln(1+n)$, $\phi$ and $\psi$ on $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$ are shown in Figs. 24 and 25, respectively. To allow a straightforward comparison, all parameters have been chosen as their conformal counterparts shown in Figs. 6 and 7. The same is true for the distribution of the local disorder parameters, namely a uniform distribution for amplitudes and locations of scatterings. This assumption may be broken to allow for different disorder distributions, but it always leads to the same results provided that both the $m_{j}$ and $\delta t_{j}$ are random. We state this fact without explicitly showing our checks, as they are qualitatively indistinguishable from those discussed in Appendix A.2. A similar argument follows for the convergence test for the particle production rate at large $N_{s}$, also addressed in Appendix A.2 for the conformal case.

Moments of $\ln|X|^{2}$ on sub-horizon scales: The time-evolution of the massless scalar field amplitude is shown in the top of Figs. 24 and 25. Inside the horizon, $\langle\ln|X|^{2}\rangle$ is constant, independently of both the scattering strength and the wavenumber of the corresponding mode. From (5.1), for $n\ll 1$, we can parametrize this behavior using the conformal expression (4.7). Conversely, ${\rm Var}\left[\ln|X|^{2}\right]$ grows exponentially, with a rate that is independent of $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$, and approximately equal to (4.8). Note that this implies that the variance grows $\propto a^{2}$.

Moments of $\ln|X|^{2}$ on super-horizon scales: Outside the horizon, the difference between the massless and the conformal scenarios becomes evident. Unlike the conformal case, in which the magnitude of $X$ remains approximately constant for weak scattering, for a massless spectator field the magnitude grows with the scale factor in super-horizon scales if $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$; this is of course the expected behavior for a massless adiabatic mode, for which $|\chi|\sim{\rm const}$. In order to break from the adiabatic limit, the strong scattering regime must be considered. Fig. 25 suggests that the rate of this growth is independent of the mode wavenumber. In analogy with the conformal case we therefore write,

[TABLE]

where the functions $\mu_{1,2}$ are shown in Fig. 26; therein the number of realizations as well as the scattering parameters are chosen as in its conformal counterpart, Fig. 8. As it is clear, $\langle\ln|X|^{2}\rangle$ grows as $a^{2}$ outside the horizon for $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 0.3$, it increases at a reduced rate for $0.3\lesssim\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 50$, and grows exponentially with rate $\mu_{1}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.25}-4>2$ for $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 50$.

The time-rate dependence on the scattering parameter for the variance of $\ln|X|^{2}$ is shown in the lower panel of Fig. 26. This dependence is very similar to that for the conformal mass case, with the variance growing as $\mu_{2}\sim\mathcal{N}_{s}(\sigma/H)^{2}$ in the weak scattering limit, while for strong scattering we have $\mu_{2}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.32}$. In analogy with the conformal case, one can verify that $\langle\ln|X|^{2}\rangle$ characterizes the typical member of the ensemble of realization, despite the rapidly growing variance (see Appendix A.1 for details).

Moments of $\ln(1+n)$ on sub-horizon scales: The evolution of the mean and variance of $\ln(1+n)$ is shown in a log-scale in the second row of Fig. 24. Both moments grow exponentially ($\propto a^{2}$) inside the horizon, with values dependent on the scattering strength parameter $\mathcal{N}_{s}(\sigma/H)^{2}$, but with rates independent of it. The exponential growth continues until the mode leaves the horizon, for weak scattering, or shortly before horizon crossing, for strong scattering. Fig. 25 further demonstrates that the occupation number is dependent on the wavenumber, but the growth rate is independent of it. It is straightforward to check that the expressions (4.10), valid in the conformal case, also correctly describe the sub-horizon evolution of the occupation number in the massless case. This is consistent with the fact that, in the $|k\tau|\gg 1$ limit, the mode functions in both cases have the same plane-wave form (3.13).

Moments of $\ln(1+n)$ on super-horizon scales: The third row of Figs. 24 and 25 displays the time evolution of the occupation number moments, but in this case in a linear scale. The linear growth of $\langle\ln(1+n)\rangle$ is evident, with a rate that is seemingly independent of $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$ for weak scattering, and independent only of $k$ for strong scattering,

[TABLE]

The functional dependence parametrized by the $\tilde{\mu}_{1,2}$ functions in the massless case is shown in Fig. 26. Note that for the massless case, $\tilde{\mu}_{1}\simeq\mu_{1}+4$ for any scattering strength, while $\tilde{\mu}_{2}\simeq\mu_{2}$ for $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 10^{-1}$. The growth rate for the mean is approximately constant, $\partial_{Ht}\langle\ln(1+n)\rangle\simeq 6$, in the weak scattering regime, with $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 0.4$. When $0.4\lesssim\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 40$, the typical occupation number grows at a slower rate. Finally, for strong scattering, the mean grows with the scattering strength parameter, $\tilde{\mu}_{1}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.25}$.

The time rate for ${\rm Var}\left[\ln(1+n)\right]$ is shown in the lower panel of Fig. 26. For very low values of the scattering strength parameter, $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 10^{-2}$, the rate appears approximately constant, $\tilde{\mu}_{2}\sim 4\times 10^{-3}$. For larger values, $10^{-2}\lesssim\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 1$, the growth rate becomes scattering strength-dependent and relatively steep, $\tilde{\mu}_{2}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{1.4}$. In the strong scattering regime the power-law dependence is gentler with $\tilde{\mu}_{2}\sim\left[\mathcal{N}_{s}(\sigma/H)^{2}\right]^{0.32}$. Although this rate is steeper than that of the mean, we expect $\exp[\langle\ln(1+n)\rangle]$ to be a good measure of the typical number of particles produced, in analogy with the conformal scenario. We confirm this fact in detail in Appendix A.1.

Moments of $\phi$ on sub-horizon scales: The fourth row of Figs. 24 and 25 show the time evolution of the angular parameter $\phi$. For all values of $\mathcal{N}_{s}(\sigma/H)^{2}$ and $k$ we find that, inside the horizon,

[TABLE]

the expected results for a uniformly distributed random variable in $(-\pi/2,\pi/2)$. Similarly to the conformal case, both moments oscillate about these values as the mode leaves the horizon, with amplitudes dependent of the scattering strength parameter and independent of the wavenumber of the mode.

Moments of $\phi$ on super-horizon scales: Far outside the horizon, the mean and the variance of $\phi$ are driven to

[TABLE]

Numerically we find that the final value of $\langle\phi\rangle$ deviates from zero for any $\mathcal{N}_{s}(\sigma/H)^{2}$ up to a numerical error $\lesssim\mathcal{O}(10^{-20})$. We also find that the decay of the variance may be parametrized as

[TABLE]

for any scattering strength.

Moments of $\psi$ on sub-horizon scales: The time-dependence of the moments of $\psi$ is shown in the bottom panels of Figs. 24 and 25. In this case, the mean and the variance correspond to those of a uniformly distributed random variable in $(0,\pi)$ inside the horizon,

[TABLE]

Moments of $\psi$ on super-horizon scales: Outside the horizon, both the mean and the variance of $\psi$ are dependent on the scattering strength, but independent of $k$. We find for their asymptotic values

[TABLE]

Here we have ignored a mild dependence on the scattering parameter for $\langle\psi\rangle$, which corresponds to a $\lesssim 6\%$ variation. Note that the variance of $\psi$ depends on the scattering parameter only for weak scattering, and it is frozen at the same constant value as in (5.9) for strong scattering, consistent with a non-evolving probability distribution.

5.1.3 Probability densities

We now turn to the description of the evolution of the probability density functions for the random variables $\ln|X|^{2}$, $\ln(1+n)$, $\phi$ and $\psi$. We base our analysis on Figs. 28-30, which show snapshots of the instantaneous normalized pdfs for the field and transfer matrix parameters for selected values of $\mathcal{N}_{s}(\sigma/H)^{2}$. The scattering parameters are taken as those for the conformal case discussed in Section 4.1.3, to allow for a simple comparison. As in that case, the pdfs are built using a Gaussian kernel density estimator of variable bin size, with a periodic extension of the data for the angular variables in the case where the pdf support is of width $\pi$.

Pdf for $\ln|X|^{2}$: The evolution of the numerical probability distribution for $\ln|X|^{2}$ is shown in Fig. 27. Similar to the conformal case, the description is straightforward, as the pdf exhibits the bell-shape for all times and values of scattering strength. For super-horizon modes, confirmation is provided in Fig. 31, where Gaussian fits are superimposed over the numerical pdf for four different values of the scattering parameter. In the case of sub-horizon modes, we prove the normality of $\ln|X|^{2}$ in Section 5.3.1 (see Fig. 36). We therefore conclude in this case as well that

 

The squared-field amplitude, $|X|^{2}$, follows a log-normal distribution both inside and outside the horizon, for weak and strong scattering.

Pdf for $\ln(1+n)$: Fig. 28 shows the evolving probability distribution for $\ln(1+n)$. The qualitative resemblance with the conformal case of Fig. 10 is evident. The leftmost two columns show pdfs that flow in a manner almost independent of the value of $\mathcal{N}_{s}(\sigma/H)^{2}$, albeit with different numerical values. As anticipated in Section 5.1.2, the early time pdf is of a highly skewed, almost exponential form, with a coefficient of variation $\tau_{\ln(1+n)}\approx 1$ (see Appendix A.1). In the next section we will confirm this result analytically (see Fig. 35). The distributions grow in width and mean, maintaining their quasi-exponential shape until horizon crossing, albeit some deformation is clear for larger values of the scattering strength parameter.

In the middle column of Fig. 28 we can see the form of the pdfs at the horizon exit transition. Unlike the conformal case, the transitional form for the pdfs is not evident here, as they have all clearly evolved toward their right- or center-lobed forms, depending on the value of $\mathcal{N}_{s}(\sigma/H)^{2}$. The last two columns of the figure in question show the deep super-horizon form of the distributions. Similarly to the conformally massive field, all pdfs are notably peaked about their means, with the weak scattering distributions retaining a tail of low $n$ realizations. Fig. 32 shows four of the five panels at $H(t-t_{k})=20$ together with skew-normal fits to the data (c.f. Eq. (4.17)). For $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$ the distributions present significant skewness. In turn, $\mathcal{N}_{s}(\sigma/H)^{2}\gg 1$ leads to a normally distributed $\ln(1+n)$, implying a log-normally distributed occupation number $n$.

Pdf for $\phi$: Fig. 29 displays the distribution of the transfer matrix phase $\phi$ at several times and scattering strengths. Comparison with Fig. 11 reveals that the massless and conformal forms for the $\phi$ pdf have similar structures and evolution, although the central values and the narrowing rates differ. The two leftmost columns show approximately uniform distributions over $(-\pi/2,\pi/2)$ inside the horizon for all $\mathcal{N}_{s}(\sigma/H)^{2}$. Some structure is visible, although it is likely due to the finite number of realizations considered for the calculation.

The central column of Fig. 29 shows the form of $\mathcal{P}[\phi]$ at the horizon exit transition. It is clear here that the density evolves to a two-lobed distribution with a rate that increases with $\mathcal{N}_{s}(\sigma/H)^{2}$. In the present case the two lobes are approximately symmetrically located with respect to $\phi=0$. This structure is clearly visible for weak scattering, but it is lost in the strong scattering case due to the rapid sharpening of the distribution. Unlike the conformal scenario, it is not clear that the height of both lobes is equal for $\mathcal{N}_{s}(\sigma/H)^{2}\gg 1$, as some asymmetry is clearly visible, consistent with the distribution of points in Fig. 23.

The two rightmost columns of Fig. 29 attempt to show the form for the $\phi$ pdf with $|k\tau|\ll 1$. In the top panels we can see the two lobes of the distribution approaching each other exponentially fast, with decreasing widths. As the rates for these processes are dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$, this evolution cannot clearly be seen in the panels for $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 1$; similarly to the conformal case our pdf estimator is not capable of showing clearly the structure of the distribution. We further discuss the super-horizon dynamics of the $\phi$ distribution in Section 5.3.2, where we compute an analytical approximation by means of a reduced Fokker-Planck equation (see Fig. 37).

Pdf for $\psi$: Fig. 30 shows the instantaneous form of the pdf for $\psi$. Deep inside the horizon the distribution is approximately uniform at any time and for any scattering strength. Some fluctuating features are visible, likely an artifact of our finite sized ensemble. As in the conformal case, the uniformity of the distribution is preserved after horizon exit for $\mathcal{N}_{s}(\sigma/H)^{2}\gg 1$, and its freeze-out is clear from the rightmost three columns. For weak scattering the distribution evolves toward a bell shape centered at $\psi\simeq\pi/2$ outside the horizon. Fig. 33 shows the four upper right panels of Fig. 30 compared to a normal distribution with mean and variance (5.10); $\psi$ is therefore normally distributed in this limit. Finally, for $\mathcal{N}_{s}(\sigma/H)^{2}=1$ the limiting distribution appears intermediate between a uniform and a normal distribution.

5.2 The field two-point function

We now proceed to analyze numerically the two-point function of the logarithm of the field amplitude defined in (4.20), in terms of the driftless variable $Z_{k}(t)$, defined in (4.19). Fig. 34 shows the time-dependence of the two point function in two cases.

Unequal time: The left column of Fig. 34 corresponds to the unequal time ($t^{\prime}\neq t$) but equal momenta ($k^{\prime}=k$) case, for selected values of $t^{\prime}$ and the scattering parameter. For all three values of $\mathcal{N}_{s}(\sigma/H)^{2}$ the behavior of the curves is very similar to that in the conformal case. Namely, for $t<t_{k}$ the two point function has negligible value for any $t^{\prime}$. If $t^{\prime}\lesssim t_{k}$, the two point function does not grow as the mode crosses outside the horizon (purple and blue curves). However, if $t^{\prime}>t_{k}$, $\langle Z_{k}(t)Z_{k}(t^{\prime})\rangle$ increases linearly with time after horizon crossing, with rate $\mu_{2}$, to posteriorly freeze at its value at $t=t^{\prime}$ for $t>t^{\prime}$ (green, orange and red curves). We can therefore summarize our results by Eq. (4.21), where $\mu_{2}$ denotes in this case the rate of growth of the variance of the logarithm of the massless field magnitude.

Unequal momenta: The right column of Fig. 34 shows the evolution of the field two-point function for equal time but unequal momenta. Here, if $k^{\prime}\leq k$, the two-point function grows linearly with time after the mode leaves the horizon (purple, blue). If instead $k^{\prime}>k$, the two-point function grows only after the mode with momentum $k^{\prime}$ crosses the horizon (green, orange, red). This behavior is analogous to that in the conformal case, and therefore it is described by Eq. (4.23).

Unequal time and momenta: For the general case of unequal time and unequal momenta correlators, we have found that Eq. (4.24) correctly approximates the super-horizon behavior of the field two-point function. This result suggests that $|X_{k}|^{2}$ describes a geometric Brownian motion with drift in cosmic time, in full analogy with the conformal case. Also analogous is the observation that a non vanishing correlation between unequal modes exists only if both modes are outside the horizon. Finally, note that Eq. (4.25) for the $n$-point correlation function for the squared field magnitude is equally valid for the massless case.

5.3 Analytical results

In the previous sections we have discussed purely numerical results concerning the expectation values and distributions for the field amplitude and the transfer matrix parameters. Following the program established by our study of the conformally massive field, we now turn to the analytical results that are derivable for a massless spectator field using the Fokker-Planck formalism. As we previously discussed in Sections 4.1 and 4.3, we will only be able to probe analytically the very weak scattering regime $\mathcal{N}_{s}(\sigma/H)^{2}\ll 1$. However, as we have shown numerically, these results are universal in the sub-horizon limit, while in the super-horizon case we will be able to derive the form of the lowest moments of the occupation number, the squared field magnitude and the ${\sf M}$-parameter $\phi$, as well as the probability density of the later. We therefore split our discussion depending on the magnitude of $\mathcal{N}_{s}^{-1}|k\tau|$ (c.f. Eq. (4.30)).

5.3.1 Sufficiently sub-horizon

Section 3.3 contains the general form for the coefficients of the FP equation given Dirac-delta scatterers with uncorrelated amplitudes of vanishing mean. In the present scenario with a massless spectator field, the coefficient functions (3.32a)-(3.32d) take the form

[TABLE]

It is clear that, in the deep sub-horizon regime $|k\tau|\ll 1$, the previous expressions reduce simply to their conformal counterparts (4.26a)-(4.26d). Therefore, all the analysis carried out in Section 4.3.1 is also valid in the massless case. This is a consequence of the Bunch-Davies (BD) initial condition imposed on any free-field mode function, which looks precisely like the free conformal mode function (3.13) [53]. We can then immediately conclude that, in the massless scenario, $\phi$ is uniformly distributed inside the horizon, and

[TABLE]

[TABLE]

with $\xi=\frac{\mathcal{N}_{s}}{2}\left(\frac{\sigma}{2k_{\rm phys}}\right)^{2}$, as defined in (4.39). These expressions are in agreement with the numerical results for the moments of $\ln(1+n)$ discussed in Section 5.1.2, and, as Fig. 35 shows, they also agree with the numerically obtained pdfs at early times in the weak scattering regime. Similarly, for the scalar field amplitude we obtain a log-normal distribution inside the horizon,

[TABLE]

This distribution agrees with the numerical results shown in Fig. 36, which correspond to four selected weak scattering panels from Fig. 27. Note that the deviation between the numerical and the analytical results is maximal for the largest values of $\mathcal{N}_{s}(\sigma/H)^{2}$ and time, indicating that $X$ experiences some sub-horizon drift as the strength of the scattering is increased.

5.3.2 Outside the horizon

We now turn to the study of the FP equation outside the horizon. Following the recipe (4.32), and assuming that $n\gg 1$, we obtain the following ${\sf R}$-matrix correlators,

[TABLE]

Pdf for $\phi$: Substitution of the correlators (5.15a)-(5.15e) into (3.28) leads to the FP equation

[TABLE]

where the shifted angular variable $\varphi$ was defined in (4.47). Integrating both sides of (5.16) with respect to $\lambda$, we obtain the following equation for the marginal probability distribution $w$ defined in (4.49),

[TABLE]

where

[TABLE]

In analogy with the conformal mass case, one would expect the temporal dependence to be negligible in the very-late time limit, with a definite limiting pdf. If this was the case, this distribution would correspond to the solution of the time-independent equation

[TABLE]

which has the solution $w(\phi)\propto\csc^{2}\phi$, that is divergent at $\phi=0$ and non-normalizable. Nevertheless, as the conformal case showed us, we can use this expression as an approximation to the time-dependent solution, provided that we choose a suitable time-dependent cutoff around $\phi=0$. The numerical solution of equation of (5.17) is shown in Fig. 37, for different values of $\tilde{\xi}$, under the assumption of an initially uniform distribution. Also in this figure the approximation

[TABLE]

is shown as the dashed black curve. By comparison with the numerical result, we find the cutoff to be approximately equal to

[TABLE]

This approximation implies that

[TABLE]

which coincide with the numerical results (5.7), (5.8). Fig. 38 shows the agreement between the numerical solution of (5.17) and the fully numerical pdf computed at selected time slices, in the weak scattering case $\mathcal{N}_{s}(\sigma/H)^{2}=10^{-4}$.

Rates for $\langle\ln(1+n)\rangle$ and $\langle\ln|X|^{2}\rangle$: We now calculate the super-horizon particle production rate in the weak scattering regime. Substituting the correlators (5.15a)-(5.15e) into the general expression (4.34) gives

[TABLE]

Far outside the horizon $\delta\ll 1$. Additionally, Eq. (5.3) implies that

[TABLE]

The approximation (5.20) yields $\langle\ln|\varphi|^{2}\rangle\simeq 2\ln\delta+2$, which in turn implies that $\partial_{Ht}\langle\ln|\varphi|^{2}\rangle\simeq-6$. Together with (5.23), this signifies that for a massless spectator field, in the super-horizon and weak-scattering limits,

[TABLE]

in agreement with the numerical result shown in Fig. 26.

6 Conclusions

In this paper we have studied the non-adiabatic, stochastic production of particles due to the excitation of a spectator field in a de Sitter (inflating) background. To accomplish this, we have extended the framework of [1, 2], which was valid in Minkowski space, to an expanding universe.

To simplify the analysis and focus on the impact of expansion, we considered the case with spectator fields being non-adiabatically excited by of a series of delta-function scatterers with random strengths.151515We note that the formalism is general enough to accommodate any expansion history and different types of (localized in time) non-adiabatic interactions. We carried out detailed numerical calculations using the transfer matrix approach for a wide range of parameters (including strong and weak scattering), and we explored both the individual and the statistical properties of quantities related to particle production. We developed a Fokker-Planck equation to analytically understand the results of our numerical simulations in the weak scattering limit. In the limits tested we found excellent agreement between the analytical and numerical results, both in the sub- and super-horizon limits. While the sub-horizon behavior is consistent with previous work, the results for super-horizon behavior are new.

While we have already summarized our main results in Section 2, we re-emphasize three general results again: (1) In the limit of large number of scatterings per Hubble time $\mathcal{N}_{s}\gg 1$, a single parameter $\mathcal{N}_{s}(\sigma/H)^{2}$ determines most of the statistical properties of the spectator field outside the horizon. Here $\sigma$ characterizes the strength of the individual scatterers. (2) The field amplitude is log-normally distributed independent of the strength of scattering, the size of the physical momentum relative to the Hubble parameter, or the bare mass of the field. (3) The logarithm of the field amplitude satisfies the necessary properties for an approximately Wiener process outside the horizon.

We note that in Minkowski space, we had found a universal result that the occupation number is log-normally distributed (in the late time limit). Occupation numbers, however, lose their physical meaning outside the horizon. It is the field amplitude that connects more directly to observables, especially on super-horizon scales.

In more detail, in the sufficiently sub-horizon regime the evolution of the occupation number mimics the growth that would be observed in a Minkowski background, apart from the time-dependent rescaling $k/H\rightarrow k/aH$ (as expected). It is worth noting that this result is independent of the bare mass of the field when physical wavelength is sufficiently small inside the horizon. Since the growth rate for $\ln(1+n)$ is now inversely proportional to the square of the physical momentum, we observe the expected exponential growth in the occupation number, albeit with $n\ll 1$ when $|k\tau|^{2}\gg\mathcal{N}_{s}(\sigma/H)^{2}$. Moreover, as $n$ is so small, the field decays on average as it would do in the absence of non-adiabatic events.

Outside the horizon, the mass of the spectator scalar determines the form of the mode function, and therefore the growth rate for the occupation number and the field magnitude. Not surprisingly, when scattering is very weak, the field amplitude $|\chi|^{2}$ is impervious to the non-adiabatic excitation, decaying with the inverse of the scale factor in the conformal case or remaining frozen in the massless scenario. Nevertheless, it is not immediately obvious why this immunity to growth extends all the way to scattering strength parameters that are $\mathcal{O}(1)$, a domain far outside the reach of our Fokker-Planck formalism.

For strong scattering, field and occupation numbers grow exponentially with a rate dependent on $\mathcal{N}_{s}(\sigma/H)^{2}$. The dependence on $\mathcal{N}_{s}(\sigma/H)^{2}$ extends to ${\sf M}$-parameters $\phi$ and $\psi$. The phase $\phi$ approaches a double-delta function distribution with exponential speed, while the form of the distribution of $\psi$ depends strongly on $\mathcal{N}_{s}(\sigma/H)^{2}$, interpolating between a normal and a uniform distribution depending on the magnitude of this parameter.

Our present work has been focused on describing stochastic particle production in an expanding universe; we have not addressed or calculated the observational implications here. We plan to make use of the results presented here to study the effects of a stochastically excited field acting as a source for curvature perturbations during inflation. Such sourcing (which can be taken as an indicator of a complex inflationary sector) could lead to features in the power spectrum [54], as well as potentially universal scaling relations between higher point correlations functions [55] of the curvature perturbations. The sourcing can also generate gravitational waves [56, 57]. The stochastic framework can also be applied to the early stages of non-perturbative reheating before non-linearity and thermalization take over. Finally, extending this framework to particle production in higher spin fields would be interesting to pursue. We postpone these studies for future publications.

Acknowledgements

We would like to thank Daniel Baumann, Horng-Sheng Chia (Amsterdam U.) for many stimulating and helpful discussions, Eva Silverstein (Stanford) for insightful discussions, and Jia-Liang Shen (Rice U.) for his involvement in calculating particle production in an expanding universe during the early stages of this project. Numerical results were obtained from a custom Fortran code utilizing the thread-safe arbitrary precision package MPFUN-For written by David H. Bailey. MA and MG are supported by the US Dept. of Energy grant DE-SC0018216. DG is supported by the US Dept. of Energy grant DE-SC0019035.

Appendix A Appendix

A.1 Typical vs. average

As it is discussed in detail in Sections 4.1.2 and 5.1.2, the logarithms of the occupation number $n$ and the scalar field magnitude $|X|$ are characterized by rapidly growing variances both inside and outside the horizon. Such a large amount of dispersion could signify that the mean does not provide a good approximation to the behavior of a typical member of the ensemble. We address below these concerns.

In order to address the suitability of $\ln(1+n)$ as our variable of choice to describe the evolution of the occupation number over the ensemble of realizations, we consider the so-called coefficient of variation (or noise-to-signal ratio)

[TABLE]

which determines the extent of variability in relation to the mean, assuming that the random variable in question is measured on a ratio scale (in this case $\ln(1+n)\geq 0$). The time dependence of $\tau_{\ln(1+n)}$ is shown in Fig. 39 for selected values of $\mathcal{N}_{s}(\sigma/H)^{2}$ for both the conformal (top) and massless (bottom) cases. In both scenarios, it is clear that, for any scattering strength, the coefficient of variation is approximately equal to one inside the horizon. This suggests an exponential distribution for $\ln(1+n)$, as confirmed in Sections 4.1.3 and 5.1.3. It is also worth noting a transient growth of $\tau_{\ln(1+n)}$ during horizon crossing for weak scattering, which is eventually overcome by the growth of the mean. Most importantly, in all cases we find that $\tau_{\ln(1+n)}$ decreases monotonically as a function of time far outside the horizon. The mean of $\ln(1+n)$ is therefore a good measure of the number of particles produced in most regimes.

In the case of the scalar field magnitude, the coefficient of variation is not a suitable measure of the quality of $\langle\ln|X|^{2}\rangle$ as a descriptor of the behavior of a typical member of the ensemble, given that the value of this mean depends on the wavenumber $k$ and can be positive or negative. Nevertheless, a qualitative argument can be built by comparing the evolution in time of $\langle\ln|X|^{2}\rangle$ and several individual trajectories, for a given ensemble of amplitudes and locations of the non-adiabatic events. This comparison is made in Fig. 40, where 50 individual trajectories, shown in gray, are displayed together with $\langle\ln|X|^{2}\rangle$, shown in black, for weak, moderate and strong scattering, and for both the conformal and massless cases. Also shown therein in blue is the logarithm of the expectation value of the squared field magnitude, $\ln\,\langle|X|^{2}\rangle$. As it is clear, the mean-of-the-log is in all cases a better descriptor of the behavior of a typical member of the ensemble.

The previous result is consistent with the fact that $\ln|X|^{2}$ is normally distributed, as we have determined in sections 4.1.3 and 5.1.3. The normal distribution has zero skewness, and for it all average quantifiers (mean, median and mode) coincide. In contrast, $|X|^{2}$ has a lognormal distribution, which becomes heavily skewed as ${\rm Var}\,[\ln|X|^{2}]$ grows with time (it grows a “fat tail” of improbable trajectories). Its median coincides with the typical value, $|X|^{2}_{\rm typ}=e^{\langle\ln|X|^{2}\rangle}$, while its mean is dominated by unlikely events for which $|X|^{2}$ is large, as it is clear in Fig. 40. In the limit of an infinite number of realizations, this mean can be related to the moments of $\ln|X|^{2}$ as follows,

[TABLE]

The right hand side of the previous equation is shown as the red, dashed curves in Fig. 40. Clearly, as the variance grows, $\langle|X|^{2}\rangle$ becomes less reliable as an estimate for the most probable ensemble member. Nevertheless, for a finite-sized ensemble, it is a reliable bound for the value of $|X|^{2}$, up to very improbable outliers.

A.2 Convergence tests

We demonstrate here that, given a sufficiently high density of scatterers $\mathcal{N}_{s}$, the super-horizon evolution of the transfer matrix parameters and the scalar field magnitude is controlled solely by the scattering strength parameter $\mathcal{N}_{s}(\sigma/H)^{2}$. We also show that this evolution is independent of the underlying distribution of scatterer locations and amplitudes, assuming that they are randomly drawn from their corresponding ensembles. Our discussion here is restricted to the case of a conformally massive spectator field, but the conclusions are equally applicable to the massless field case.

Fig. 41 shows the dependence of the particle production rate $\tilde{\mu}_{1}=\partial_{Ht}\langle\ln(1+n)\rangle$ on the density of scatterers $\mathcal{N}_{s}$ for fixed values of $\mathcal{N}_{s}(\sigma/H)^{2}$. To construct this plot we have effectively varied $N_{s}$ from $10^{2}$ to $10^{4}$, which corresponds to varying $\mathcal{N}_{s}$ from 2.5 to 250, and we have plotted the ratio of the production rate to its asymptotic value corresponding to that with the largest $N_{s}$, $\tilde{\mu}_{1}^{\infty}\equiv\tilde{\mu}_{1}^{N_{s}=10^{4}}$. It is clear then that, for weak scattering, the rate is independent of both $\mathcal{N}_{s}$ and $\sigma^{2}$, since we found that in this regime $\tilde{\mu}_{1}\simeq 2$. In the case of strong scattering, we observe a mild dependence on $\mathcal{N}_{s}$, which is enhanced as scattering becomes stronger. Nevertheless, for the values considered here, this deviation is at most of $\sim 33\%$, and it clearly decreases for $\mathcal{N}_{s}\gg 1$, implying a universal value of the rate for fixed $\mathcal{N}_{s}(\sigma/H)^{2}$. This justifies our parametrization.

The dependence of the particle production rate on the distribution of the scatterer locations and amplitudes is shown in Fig. 42. To construct these results we have considered $N_{r}=200$ realizations per value of $\mathcal{N}_{s}(\sigma/H)^{2}$ for each of the following scenarios:

1. (a)

Neither $m_{i}$ or $\delta t_{i}$ are random (purple curve); their values are fixed to be equal to $+\sigma$ and $\langle\delta t_{i}\rangle/2$, respectively, of the corresponding random cases. 2. (b)

The amplitudes $m_{i}$ are fixed, while the locations $\delta t_{i}$ are random and uniformly distributed (blue curve). The set $\{\delta t_{i}\}$ is different between realizations for a given value of $\mathcal{N}_{s}(\sigma/H)^{2}$. 3. (c)

The amplitudes are random and uniformly distributed, while the locations are fixed and equispaced (teal curve). In this and all following cases, the ensemble of realizations is built out of $N_{r}$ distinct sets of $m_{i}$. 4. (d)

Both $m_{i}$ and $\delta t_{i}$ are uniformly distributed, but the ensemble of realizations is built out of a single set $\{\delta t_{i}\}$, denoted as $\delta t_{i}^{(1)}$, for any given $\mathcal{N}_{s}(\sigma/H)^{2}$ (light green curve). This scenario is advantageous, as it guarantees that all realizations can be averaged at any given time step. 5. (e)

Both $m_{i}$ and $\delta t_{i}$ are uniformly distributed, and for each realization a different, random set $\{\delta t_{i}\}$ is considered; we denote this ensemble of $N_{r}$ different sets by $\delta t_{i}^{(N_{r})}$ (orange curve). This scenario represents more faithfully the random ensemble of realizations, at the cost of an underlying coarse graining in time to allow for the calculation of expectation values. 6. (f)

The amplitudes $m_{i}$ follow a normal distribution centered at zero with standard deviation $\sigma$. The $\delta t_{i}$ are also normally distributed, with means and variances that coincide with those for case (e). Here the set of locations $\{\delta t_{i}\}$ is different between realizations (red curve).

From Fig. 42, it is clear that the stochasticity of the scattering amplitudes $m_{i}$ is what determines the time evolution of the transfer matrix parameters, here illustrated for the production rate $\tilde{\mu}_{1}$. In the absence of stochasticity in the amplitudes (cases (a) and (b)), a structure reminiscent of a band is visible. This band structure disappears if the $m_{i}$ are drawn from probability distributions, and the functional dependence of $\tilde{\mu}_{1}$ on the scattering strength parameter is independent of the details of these underlying distributions (cases (c), (d), (e) and (f)). We have verified that these results are independent of the value of the wavenumber $k$.

A.3 Backreaction

In Sections 4.1.2 and 5.1.2 we have found that the rate of growth of the mean and the variance of the scalar field magnitude can be $\mathcal{O}(1)$ or larger for $\mathcal{N}_{s}(\sigma/H)^{2}\gtrsim 1$. This rapid growth will inevitably lead to a $\chi$-dominated universe at sufficiently late times. That is, the mean energy density of the field would be comparable or larger than that of the background. In this Appendix we determine when the backreaction from the excitation of $\chi$ cannot be ignored.

In order to avoid the discontinuous jumps in the time derivative of the scalar field due to the junction condition (3.24), we consider the energy density $\rho_{\chi}$ in-between scatterings. We can then write in general the mean energy density between scatterings as follows,

[TABLE]

Note that this quantity depends on the mean squared-magnitude of the field, $\langle|X_{k}|^{2}\rangle$, and not on its “typical” (median) value, $\exp\langle\ln|X_{k}|^{2}\rangle$. In the present context, this is a sensible result: for a given ensemble of realizations, the average energy density will be dominated by those for which $|X_{k}|^{2}$ is large. Although the median would still provide a valuable measure of the energy density of the most probable member of the ensemble, for half of the members of the ensemble it will underestimate the onset of the backreaction regime.

Note that sufficiently deep inside the horizon, (A.3) will approximately lead to the (UV-divergent) energy density of the free scalar field. We can therefore restrict our calculation to modes that have left the horizon at the given (conformal) time $\tau$. In order to enforce this distinction, we will consider the “adiabatically regulated” energy density $\Delta\langle\rho_{\chi}\rangle$, which we define as

[TABLE]

We now proceed to evaluate this regulated energy density for the conformally massive and massless cases. The backreaction regime will then be defined as that for which this mean energy density becomes comparable to the background, $\Delta\langle\rho_{\chi}\rangle\simeq 3M_{P}^{2}H^{2}$.

A.3.1 Conformally massive field

In light of our definition (A.4), we can immediately disregard the contribution to the mean energy density from sufficiently deep sub-horizon modes. More specifically, Eq. (4.7) shows that the mean of $\ln|X_{k}|^{2}$ is always well approximated by its free, vacuum value, while (4.8) implies that the variance is $\ll 1$ provided that

[TABLE]

In order to simplify our following estimate of $\langle\rho_{\chi}\rangle$, we will assume that the deviation of the mean from its free value and the magnitude of the variance can be disregarded until horizon crossing, $|k\tau|\simeq 1$. Note that this in general will somewhat overestimate the computed energy density in the case of weak scattering, while it may underestimate it for strong scattering.

Let us now note that, in between scatterings, Eq. (4.4) implies that for super-horizon modes with $n\gg 1$,

[TABLE]

up to a $\mathcal{O}(|k\tau|^{2})\ll 1$ correction. This allows us to write the regulated energy density as

[TABLE]

In the strong scattering regime, the scalar field grows exponentially fast outside the horizon with a rate that can even exceed that determined by the scale factor. It is therefore natural to expect that the backreaction regime will be reached within a few $e$-folds of inflation if $\mathcal{N}_{s}(\sigma/H)^{2}$ is sufficiently large. Hence, we now assume that scatterings have not taken place for an arbitrarily long period of time but instead started when $\tau=\tau_{0}$. When this is the case, the form of $\langle|X_{k}|^{2}\rangle$ will depend on the moment in which the mode $k$ leaves the horizon. If the mode leaves the horizon during the non-adiabatic epoch, while scatterings are active ($|k\tau_{0}|>1$), we can write (see Eq. (A.2))

[TABLE]

where we have used the lognormality of $|X_{k}|^{2}$, $\langle|X_{k}|^{2}\rangle=e^{\langle\ln|X_{k}|^{2}\rangle+\frac{1}{2}\langle Z_{k}^{2}\rangle}$. In passing, note that this implies that the “typical” energy density can be recovered from our expressions by formally taking the limit $\mu_{2}\rightarrow 0$. If instead the mode with wavenumber $k$ crosses outside the horizon before scatterings are active, we have

[TABLE]

Combining these results we can rewrite (A.7) as

[TABLE]

Therefore, the mean energy density in $\chi$ will become comparable to the background if, during scatterings, $\Delta\langle\rho_{\chi}\rangle\simeq 3M_{P}^{2}H^{2}$, or equivalently, if

[TABLE]

where $N_{e}(\tau)$ denotes the number of e-folds between $\tau_{0}$ and $\tau$.

Fig. 43 shows the region of parameter space excluded by the backreaction constraint $\langle\rho_{\chi}\rangle\lesssim 3H^{2}M_{P}^{2}$, for $H=10^{13}\,{\rm GeV}$ (light gray). The boundary curve has been constructed from the numerical results shown in Fig. 8. Note that for $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 6.8$ the duration of the non-adiabatic epoch is not bounded by backreaction. For strong scattering, the bound depends on the scattering strength parameter, decreasing from $N_{e}\lesssim 10^{2}$ for $\mathcal{N}_{s}(\sigma/H)^{2}\sim 10$ to $N_{e}\lesssim 1$ for $\mathcal{N}_{s}(\sigma/H)^{2}\sim 10^{4}$.

Also shown shaded in the figure in dark gray is the backreaction constraint on the typical member of the ensemble. As we discussed above, it may be obtained as the $\mu_{2}\rightarrow 0$ limit of the mean energy density, and we have denoted it by $\rho_{\chi}^{\rm typ}$. Note that for it, as expected, the constraint is milder. It is inexistent for $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 67$ and becomes dependent on the scattering strength parameter in an almost parallel way to the constraint on $\langle\rho_{\chi}\rangle$, for strong scattering. Also note in this case that $N_{e}\lesssim 2$ for $\mathcal{N}_{s}(\sigma/H)^{2}\sim 10^{4}$.

A.3.2 Massless field

Let us now consider the massless spectator field. Similarly to the conformal case, we will assume that the free field solution is a valid approximation to the dynamics of the excited field up to horizon crossing. From Eq. (5.1) we obtain that outside the horizon, with $n\gg 1$,161616Note that (A.6) and (A.13) imply that

$\frac{d\chi_{k}}{d\tau}\;\simeq\;\frac{\chi_{k}}{\tau}\times\begin{cases}1\,,&M^{2}=2H^{2}\,,\\ |k\tau|^{2}\,,\quad&M^{2}=0\,.\end{cases}$

(A.12)

[TABLE]

The mean energy density takes then the approximate form

[TABLE]

where we have approximated the free mode function outside the horizon.

Following the same arguments as for the conformal case, for a mode $k$ that leaves the horizon when scatterings are active, we can immediately write

[TABLE]

where $\mu_{1,2}$ are now given by their massless values, shown in Fig. 26. If instead the mode crosses the horizon before scatterings are active, we now have

[TABLE]

Substitution of (A.15) and (A.16) into (A.14) leads to

[TABLE]

Note the similarity between this result and (A.10), with the identification $\mu_{1}\rightarrow\mu_{1}-2$. We can therefore immediately write the number of $e$-folds that the excitation of the massless field $\chi$ can last before its mean energy density becomes comparable to that of the background:

[TABLE]

Fig. 44 shows the exclusion regions determined by $\langle\rho_{\chi}\rangle$ (light gray) and $\rho_{\chi}^{\rm typ}$ (dark gray) as functions of $N_{e}$ and the scattering strength parameter. Notice that, unlike the conformal case, the number of $e$-folds is unconstrained by $\langle\rho_{\chi}\rangle$ for $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 16$, and by $\rho_{\chi}^{\rm typ}$ for $\mathcal{N}_{s}(\sigma/H)^{2}\lesssim 287$. Nevertheless, the boundary contours are steep functions of the scattering strength parameter and, similarly to the conformal case, only $N_{e}\sim\mathcal{O}(1)$ is allowed for the largest value of the scattering strength considered, $\mathcal{N}_{s}(\sigma/H)^{2}\sim 10^{4}$.