Multifield Stochastic Particle Production: Beyond a Maximum Entropy Ansatz

TL;DR

This paper develops a framework for calculating particle production rates in multi-field scalar systems with stochastic interactions, extending previous maximum entropy approaches and providing analytic results validated by numerical simulations.

Contribution

It introduces a new analytical approach for multi-field stochastic particle production, connecting it to maximum entropy assumptions and exploring deviations for different parameter regimes.

Findings

Analytic formulas for particle production in 1 and 2 fields

Agreement between analytic results and numerical simulations

Deviations from maximum entropy ansatz when parameters differ significantly

Abstract

We explore non-adiabatic particle production for $N_{\rm f}$ coupled scalar fields in a time-dependent background with stochastically varying effective masses, cross-couplings and intervals between interactions. Under the assumption of weak scattering per interaction, we provide a framework for calculating the typical particle production rates after a large number of interactions. After setting up the framework, for analytic tractability, we consider interactions (effective masses and cross couplings) characterized by series of Dirac-delta functions in time with amplitudes and locations drawn from different distributions. Without assuming that the fields are statistically equivalent, we present closed form results (up to quadratures) for the asymptotic particle production rates for the $N_{\rm f}=1$ and $N_{\rm f}=2$ cases. We also present results for the general $N_{\rm f} >2$ case, but with more restrictive assumptions. We find agreement between our analytic results and direct numerical calculations of the total occupation number of the produced particles, with departures that can be explained in terms of violation of our assumptions. We elucidate the precise connection between the maximum entropy ansatz (MEA) used in Amin & Baumann (2015) and the underlying statistical distribution of the self and cross couplings. We provide and justify a simple to use (MEA-inspired) expression for the particle production rate, which agrees with our more detailed treatment when the parameters characterizing the effective mass and cross-couplings between fields are all comparable to each other. However, deviations are seen when some parameters differ significantly from others. We show that such deviations become negligible for a broad range of parameters when $N_{\rm f}\gg 1$.

Figures

Equations

$$S^{(2)}=\int d^{4}x$$

$${\mathcal{M}}_{ab}(\tau )=m_a^2{\delta }_{ab}+m_{ab}^{\mathrm{s}}(\tau ).$$

$$\left(\frac{d^2}{d{\tau }^2}-{\nabla }^2+m_a^2\right){\chi }^a(\tau ,\vec{x})+\sum _{b=1}^{N_{\mathrm{f}}}{m}_{ab}^{\mathrm{s}}(\tau ){\chi }^b(\tau ,\vec{x})=0.$$

$$\left(\frac{d^2}{d{\tau }^2}+{\omega }_a^2\right){\chi }^a(\tau ,\vec{k})+\sum _{b=1}^{N_{\mathrm{f}}}{m}_{ab}^{\mathrm{s}}(\tau ){\chi }^b(\tau ,\vec{k})=0,$$

$${\omega }_a^2(k)=k^2+m_a^2,$$

$${\chi }_j^a(\tau )\equiv \frac{1}{\sqrt{2{\omega }_a}}\left[{\beta }_j^a{e}^{i{\omega }_a\tau }+{\alpha }_j^a{e}^{-i{\omega }_a\tau }\right],$$

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

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

$$n_a(j)=\frac{1}{2{\omega }_a}\left(|{\dot{\chi }}_j^a{|}^2+{\omega }_a^2|{\chi }_j^a{|}^2\right)-\frac{1}{2}=|{\beta }_j^a{|}^2,$$

$$n(j)\equiv \sum _{a=1}^{N_{\mathrm{f}}}{n}_a(j).$$

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

$${\mathsf{R}}_j\equiv {\mathsf{M}}_j{\mathsf{M}}_j^{†},\text{and correspondingly}\mathsf{R}(j)\equiv \mathsf{M}(j){\mathsf{M}}^{†}(j).$$

$$\mathsf{R}=\underset{\mathsf{U}}{\underbrace{\left(\begin{array}{cc}\mathsf{u}& 0\\ 0& {\mathsf{u}}^{\ast }\end{array}\right)}}\underset{\mathsf{F}}{\underbrace{\left(\begin{array}{cc}\mathsf{f}& \tilde{\mathsf{f}}\\ \tilde{\mathsf{f}}& \mathsf{f}\end{array}\right)}}\underset{{\mathsf{U}}^{†}}{\underbrace{\left(\begin{array}{cc}{\mathsf{u}}^{†}& 0\\ 0& {\mathsf{u}}^{\mathsf{T}}\end{array}\right)}}=\mathsf{U}\mathsf{F}{\mathsf{U}}^{†},$$

$$\mathsf{M}=\left(\begin{array}{cc}\mathsf{u}& 0\\ 0& {\mathsf{u}}^{\ast }\end{array}\right)\left(\begin{array}{cc}\sqrt{1+\mathsf{n}}& \sqrt{\mathsf{n}}\\ \sqrt{\mathsf{n}}& \sqrt{1+\mathsf{n}}\end{array}\right)\left(\begin{array}{cc}\mathsf{v}& 0\\ 0& {\mathsf{v}}^{\ast }\end{array}\right),$$

$$n(j)=\frac{1}{4}\mathrm{Tr}[\mathsf{R}-1]=\frac{1}{2}\mathrm{Tr}[\mathsf{f}(j)-1]=\frac{1}{2}\sum _{a=1}^{N_{\mathrm{f}}}(f_a-1).$$

$$P_{\tau +\delta \tau }(\mathsf{M})=\int d\mu ({\mathsf{M}}_2)P_{\tau }({\mathsf{M}}_2^{-1}\mathsf{M})P_{\delta \tau }({\mathsf{M}}_2),$$

$${\partial }_{\tau }P$$

$${\partial }_{\tau }⟨F(n)⟩=\sum _{a=1}^{N_{\mathrm{f}}}⟨\frac{F^{\prime }(n)}{2}\frac{⟨\delta f_a{⟩}_{\delta \tau }}{\delta \tau }+\sum _{b=1}^{N_{\mathrm{f}}}\frac{F^{\prime \prime }(n)}{8}\frac{⟨\delta f_a\delta f_b{⟩}_{\delta \tau }}{\delta \tau }⟩,$$

$$n_{\mathrm{typ}}\equiv \exp [⟨\ln (1+n)⟩]-1.$$

$$\boxed{{\partial }_{\tau }⟨\ln (1+n)⟩=⟨\frac{1}{2(1+n)}\sum _{a=1}^{N_{\mathrm{f}}}\frac{⟨\delta f_a{⟩}_{\delta \tau }}{\delta \tau }-\frac{1}{8(1+n{)}^2}\sum _{a,b=1}^{N_{\mathrm{f}}}\frac{⟨\delta f_a\delta f_b{⟩}_{\delta \tau }}{\delta \tau }⟩.}$$

$$\mathsf{R}(j+1)=\mathsf{R}(j)+\delta \mathsf{R},$$

$$\delta \mathsf{R}=\mathsf{R}(j){\mathsf{m}}_{j+1}^{†}+h.c.+{\mathsf{m}}_{j+1}\mathsf{R}(j){\mathsf{m}}_{j+1}^{†}.$$

$$\delta \mathsf{R}$$

$$\delta \mathsf{F}+[{\mathsf{U}}^{†}\delta \mathsf{U},\mathsf{F}]-[{\mathsf{U}}^{†}\delta \mathsf{U},\mathsf{F}]{\mathsf{U}}^{†}\delta \mathsf{U}+[{\mathsf{U}}^{†}\delta \mathsf{U},\delta \mathsf{F}]=\mathsf{G},$$

$$\mathsf{G}\equiv \left(\mathsf{Q}\mathsf{F}+h.c.\right)+\mathsf{Q}\mathsf{F}{\mathsf{Q}}^{†},\text{with}\mathsf{Q}\equiv {\mathsf{U}}^{†}{\mathsf{m}}_{j+1}\mathsf{U}.$$

$$\delta \mathsf{F}\equiv \delta {\mathsf{F}}^{(1)}+\delta {\mathsf{F}}^{(2)},$$

$$\delta {\mathsf{F}}^{(1)}+[{\mathsf{U}}^{†}\delta {\mathsf{U}}^{(1)},\mathsf{F}]={\mathsf{G}}^{(1)},$$

$$\delta {\mathsf{F}}^{(2)}+[{\mathsf{U}}^{†}\delta {\mathsf{U}}^{(2)},\mathsf{F}]-[{\mathsf{U}}^{†}\delta {\mathsf{U}}^{(1)},\mathsf{F}]{\mathsf{U}}^{†}\delta {\mathsf{U}}^{(1)}+[{\mathsf{U}}^{†}\delta {\mathsf{U}}^{(1)},\delta {\mathsf{F}}^{(1)}]={\mathsf{G}}^{(2)}.$$

$$\delta {\mathsf{f}}^{(1)}+[{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)},\mathsf{f}]={\mathsf{g}}^{(1)},$$

$$\delta {\mathsf{f}}^{(2)}+[{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(2)},\mathsf{f}]-[{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)},\mathsf{f}]{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)}+[{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)},\delta {\mathsf{f}}^{(1)}]={\mathsf{g}}^{(2)},$$

$$\delta {\tilde{\mathsf{f}}}^{(1)}+{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)}\tilde{\mathsf{f}}-\tilde{\mathsf{f}}{\mathsf{u}}^T\delta {\mathsf{u}}^{\ast (1)}={\tilde{\mathsf{g}}}^{(1)},$$

$$\delta {\tilde{\mathsf{f}}}^{(2)}+{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(2)}\tilde{\mathsf{f}}-\tilde{\mathsf{f}}{\mathsf{u}}^T\delta {\mathsf{u}}^{\ast (2)}-({\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)}\tilde{\mathsf{f}}-\tilde{\mathsf{f}}{\mathsf{u}}^T\delta {\mathsf{u}}^{\ast (1)}){\mathsf{u}}^T\delta {\mathsf{u}}^{\ast (1)}+{\mathsf{u}}^{†}\delta {\mathsf{u}}^{(1)}\delta {\tilde{\mathsf{f}}}^{(1)}-\delta {\tilde{\mathsf{f}}}^{(1)}{\mathsf{u}}^T\delta {\mathsf{u}}^{\ast (1)}={\tilde{\mathsf{g}}}^{(2)},$$

$$\delta f_a^{(1)}=g_{aa}^{(1)}.$$

$$(u^{†}\delta u^{(1)}{)}_{ab}=\frac{g_{ab}^{(1)}}{f_b-f_a},(a\ne b).$$

Full text

1. 1 Introduction

2. 2 The Model

3. 3 The Transfer Matrix Approach

4. 3.1 The Occupation Number

5. 3.2 Parametrizing ${\sf R}$

6. 3.3 The Fokker Planck Equation

7. 3.4 Increments in ${\sf R}$

8. 3.5 Dirac Delta Scatterers

9. 4 The Maximum Entropy Approximation

10. 5 A Warm-up Example: the ${N_{\rm{f}}}=1$ Case

11. 6 Beyond MEA: The ${N_{\rm{f}}}=2$ Case

12. 6.1 Calculation of the Correlators

13. 6.2 The Typical Occupation Number

14. 6.3 Summary

15. 7 Beyond MEA: Generalization to Many Fields

16. 7.1 Calculation of the Correlators

17. 7.2 The Typical Occupation Number

18. 7.3 Summary

19. 8 Summary and Future Outlook

20. A Appendix

21. A.1 The Fokker-Planck Equation

22. A.2 Coefficients of the ${N_{\rm{f}}}=2$ FP Equation

23. A.3 Numerical Method

24. A.4 Haar average of $\mathcal{F}_{\Omega}$

1 Introduction

Repeated bursts of particle production are possible during inflation and reheating after inflation. Such events typically result from a non-adiabatic change in the effective mass or the couplings between fields as the background evolution of the fields passes through special points in field space. For example, sharp turns in field trajectories or certain fields becoming effectively massless can result in significant particle production. The produced particles can impact curvature fluctuations and/or cause a change in the way energy is transferred between the effective inflaton field and the spectator/daughter fields (see for example [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]).

For a fully specified model which includes all the interaction terms between fields, one can always numerically compute the background evolution of the fields. Subsequently, one can calculate particle production in the fields for any finite wavenumber by solving the mode equation in this background. Such fully specified models are rarely available. Even for a fully specified model with many fields there might be many background solutions. Since we might not know a-priori which trajectory is taken, it is not always the most efficient way to go carry our particle production calculations (however, see for example, [32, 33] for some recent numerical approaches in calculating multifield inflationary perturbations).

An alternate approach is to think about particle production in a statistical ensemble of models or trajectories within a model, and focus on calculating the typical particle production rate. This is the approach we take in this paper. The work here is heavy up front, in terms of setting up the framework, but once set up it can be applied to a wide range of models where only statistical information about the parameters is available, and only coarse-grained information about the final observables is needed. Given the complexity of ultraviolet completion of the models of inflation, and the simplicity of the data from the early universe we believe that this statistical approach is a reasonable one to pursue.

An earlier paper [34] by one of us took advantage of a mathematical mapping between current conduction in wires with impurities and particle production in a time dependent background. In that work, a significant simplification in the results is seen due to a maximum entropy ansatz (MEA), which heuristically assumes statistically equivalent fields. In this paper, we take a more detailed approach which allows us to go beyond statistically equivalent fields. However, the cost is an increase in complexity of our derivations. Nevertheless, once obtained, our results are relatively simple to understand and use, and do reduce to the results of [34] in their overlapping domain of validity (which we delineate).

In this paper, we first provide a more explicit connection between a microscopic (but still statistical) description of the scatterers and the MEA. Without relying on this ansatz, we derive a more general Fokker-Planck equation for the time dependent probability distribution of the occupation number of fields. The equation is valid under the assumption that the particle production for each interaction is small. Using this equation, we calculate the typical total occupation number of the fields. For the sake of analytic tractability, each interaction is modeled by a Dirac-delta function in time whose amplitude and location is drawn from different distributions. We are able to derive closed form (up to quadratures) results for the typical total occupation number after a large number of scatterings ${N_{\rm s}}\gg 1$, when the number of fields ${N_{\rm{f}}}=1$ and ${N_{\rm{f}}}=2$. For ${N_{\rm{f}}}>2$, the derivation becomes increasingly complex. By making additional technical assumptions regarding the relative strengths of the effective masses and cross couplings of the fields, we are able to calculate the total occupation number for an arbitrary number of fields. As might be expected, we see a significant simplification in our results for ${N_{\rm{f}}}\gg 1$. In general, for simple results, it is essential to have ${N_{\rm s}}\gg 1$ as well as ${N_{\rm{f}}}=1,2$ or ${N_{\rm{f}}}\gg 1$.

We briefly discuss some important assumptions made in this paper, which we will also reiterate in the conclusions. First, we consider particle production in Minkowski space rather than an expanding universe which would be more appropriate for inflation and reheating applications. Part of the effect of expansion is accounted for by randomizing the location of the events – this mimics the phase scrambling that takes place in an expanding universe [3], but a more careful treatment of the effects due to expansion is left for future work. Second, we treat each Fourier mode independently; it is a linear treatment about a spatially homogeneous background. As significant particle production takes place, this assumption might be broken. Thus, our results are only valid as long as the linearized equations hold. Finally, we assume that the particle production per event is weak, and the number of particle production events is large for the occupation number to build up. More technical assumptions are discussed in the following sections once the appropriate terminology has been introduced.

Without any attempts at being exhaustive, we would like to highlight a few papers that provide some background and context for our work. Non-perturbative particle production in the context of reheating has a long history [1, 2] (see [26] for a review). Similarly, in the context of inflation, see for example [6, 27, 30, 31] (also see the review [25]). In [35, 36, 37], the authors took advantage of the connection between stochastic particle production and random Schrödinger operators in the context of reheating. The implications of weak disorder in the early universe and inflationary perturbations was discussed in [38]. In [34], a statistical framework for calculating non-adiabatic particle production rates in stochastically coupled (but statistically equivalent) multi-field scenario was presented, which in turn took advantage of related existing techniques used to understand current conduction in disordered wires [39, 40, 41, 42].

The rest of the paper is organized as follows. In section 2, we provide a Lagrangian for the perturbations in the fields undergoing non-adiabatic particle production. In section 3, we discuss the Transfer Matrix Approach where we provide the formal expressions for the occupation number of the fields in terms of a product of transfer matrices evaluated at each non-adiabatic event. Under the assumption of weak particle production per event, we derive a general Fokker-Planck equation describing the evolution of the probability density of the parameters characterizing the particle production events. In section 4, we connect the MEA to the statistical description of the underlying model. In section 5, we consider the single field (${N_{\rm{f}}}=1$) case as a warm-up example. Moving beyond MEA, in section 6, we consider the ${N_{\rm{f}}}=2$ case in detail. We first calculate the coefficients of the Fokker-Planck equation by taking averages over the properties of the underlying non-adiabatic events. We then calculate the typical occupation number of the field using the Fokker-Planck equation. In section 7, we generalize to the ${N_{\rm{f}}}>2$ case, and again compare with the corresponding results in [34]. In section 8, we discuss our results, reiterate the assumptions and caveats relevant for the results, and outline our plans for future work. A number of technical details are relegated to the appendices. In Appendix A.1, we give details of the derivation of the multi-field Fokker-Planck equation. In Appendix A.2, we provide some of the explicit calculations for the coefficients in the ${N_{\rm{f}}}=2$ case. In Appendix A.3, we show the explicit numerical algorithm used for calculating the particle production rate. Finally, in Appendix A.4 we discuss technical details of the Haar measure for the transfer matrices.

Notation and Conventions

Throughout, we will use natural units, $c=\hbar\equiv 1$. The time variable will be $\tau$, and overdots will denote derivatives with respect to $\tau$. We will use a bold sans-serif font for matrices, e.g. ${\sf M},{\sf n},{\boldsymbol{\lambda}}$. The indices $a,b,\ldots$ will be used for field indices. We use the $+---$ metric convention.

2 The Model

Consider ${N_{\rm{f}}}$ coupled scalar fields, $\chi^{a}(\tau,{\vec{x}})$, with masses $m_{a}$ and time-dependent stochastic masses and couplings $m^{\rm s}_{ab}(\tau)$. Here $a,b=1,2,\ldots,{N_{\rm{f}}}$. The quadratic action for these fields is taken to be

[TABLE]

The above action is to be interpreted as the action for perturbations around a complicated time dependent background due to the presence of many fields. The equations of motion for the fields $\chi^{a}$ are

[TABLE]

The equations of motion of the Fourier modes of the field $\chi^{a}(\tau,\vec{k})$ are

[TABLE]

where we have used the symbol $\chi^{a}$ for the Fourier transform of the fields as well as the fields themselves. We will be explicitly dealing with Fourier space quantities from now on. To reduce clutter, the dependence of the fields and $\omega_{a}$ on the wavenumber $k$ will be suppressed.

The time dependence of $m^{\rm s}_{ab}(\tau)$ can be quite complex. For example, $m^{\rm s}_{ab}(\tau)$ can consist of a series of well separated “hills” and “valleys”, uniformly distributed in time. Between these hills and valleys, the fields are assumed to be free and uncoupled (see fig. 1).

3 The Transfer Matrix Approach

The Fourier mode of the field $\chi^{a}$ after the $j$-th non-adiabatic event is

[TABLE]

where the overall normalization is chosen for future convenience. The Wronskian of the solutions $\chi^{a}_{j}$ and ${\chi^{a}_{j}}^{*}$ is a constant, $W[\chi^{a}_{j},{\chi^{a}_{j}}^{*}]\equiv 1$, which in terms of the coefficient in eq. (3.1) implies $|\alpha^{a}_{j}|^{2}-|\beta^{a}_{j}|^{2}=1$.

The coefficients $\vec{\beta}_{j}=(\beta^{1}_{j},\beta^{2}_{j}\ldots\beta^{N_{\rm{f}}}_{j})^{T}$ and $\vec{\alpha}_{j}=(\alpha^{1}_{j},\alpha^{2}_{j}\ldots\alpha^{N_{\rm{f}}}_{j})^{T}$ before and after the $j$-th scattering, are related by the transfer matrix ${\sf M}_{j}$.

[TABLE]

For a single scattering we can always write ${\sf M}_{j}=\mathbb{1}+{\sf m}_{j}$ where ${\sf m}_{j}$ contains the necessary information about the scattering (for weak scattering $|{\sf m}_{j}|\ll 1$). We can connect the coefficients before and after $j$ scatterings by simply chaining together the transfer matrices

[TABLE]

Note that ${\sf M}(j)$ denotes the transfer matrix of $j$ scatterings, whereas ${\sf M}_{j}$ is the transfer matrix for the $j$-th scattering.

3.1 The Occupation Number

We define the occupation number of each field after $j$ scatterings by333We note that this definition agrees with the usual one when thinking about $\chi^{a}(\tau,\vec{k})$ as mode functions for a single field undergoing adiabatic evolution. In more general time dependent backgrounds, and with coupled fields, the connection is not always so simple. See for example [43, 26, 24].

[TABLE]

and the sum of the occupation numbers of all the fields is defined as

[TABLE]

For technical reasons which will become obvious towards the end of the calculation, the total occupation number is the quantity we will focus on. A useful identity relating the total occupation number defined above and the transfer matrix is given by

[TABLE]

The above expression for the total occupation number motivates the following definition:

[TABLE]

For the total occupation number, we only care about ${\sf R}$ and not ${\sf M}$. This simplifies matters considerably. Note that ${\sf R}$ is Hermitian. The total occupation number is related to the sum of the eigenvalues of ${\sf R}$.

3.2 Parametrizing ${\sf R}$

How many parameters $\{\lambda_{a}\}$ are needed to describe ${\sf R}$? A general parametrization of the ${\sf R}$ matrix is

[TABLE]

where ${\tilde{\sf f}}\equiv\sqrt{{\sf f}^{2}-1}$. In the above expression ${\sf f}$ is a diagonal matrix consisting of ${N_{\rm{f}}}$ eigenvalues and ${\sf u}$ is a ${N_{\rm{f}}}\times{N_{\rm{f}}}$ unitary matrix parametrized by ${N_{\rm{f}}}^{2}$ angular variables. Note that ${\sf U}$ and ${\sf F}$ are $2{N_{\rm{f}}}\times 2{N_{\rm{f}}}$ matrices with ${N_{\rm{f}}}\times{N_{\rm{f}}}$ blocks constructed out of ${\sf f}$ and ${\sf u}$. Thus, ${\sf R}$ is parametrized by ${N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ total variables. This should be contrasted with ${\sf M}$, for which a general parametrization has the form [44]

[TABLE]

where ${\sf n}\equiv({\sf f}-\mathbb{1})/2$, and where ${\sf v}$ is a ${N_{\rm{f}}}\times{N_{\rm{f}}}$ unitary matrix parametrized by ${N_{\rm{f}}}^{2}$ angular variables, for a total of $2{N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ variables characterizing ${\sf M}$. The parameters of ${\sf R}$ will be denoted by $\{\lambda_{a}\}=\{f_{1},\ldots,f_{{N_{\rm{f}}}},\theta_{1},\ldots,\theta_{{N_{\rm{f}}}^{2}}\}$ where $f_{a}$ are the diagonal entries of ${\sf f}$ and $\theta_{a}$ are the angular variables for the unitary matrix ${\sf u}$. In this parametrization, note that

[TABLE]

3.3 The Fokker Planck Equation

Our next goal is to figure out how ${\sf R}$ evolves as the fields experience scattering events. Assuming that the parameters describing the scatterers are drawn from some distribution, we wish to derive an equation for the probability density $P_{\tau}({\sf R})$. This can then be used to get the expectation value of the occupation number (or any other function depending on variables in ${\sf R}$).

We start by adding a small time interval $\delta\tau$ with a single scatterer to an existing interval $\tau$ with $j$ scattering events. We will consider the case that the scatterer is weak. That is, the change in ${\sf M}$ is small due to the additional interval. The transfer matrix ${\sf M}_{\tau+\delta\tau}$ for the elongated interval can be written in terms of ${\sf M}_{1}\equiv{\sf M}_{\tau}$ and ${\sf M}_{2}\equiv{\sf M}_{\delta\tau}$ as ${\sf M}_{\tau+\delta\tau}={\sf M}_{\delta\tau}{\sf M}_{\tau}$ due to the composition law in eq. (3.3). As the underlying processes occurring in the $\tau$ and $\delta\tau$ strips are assumed to be uncorrelated, the probability densities $P_{\tau}({\sf M}_{1})$ and $P_{\delta\tau}({\sf M}_{2})$ are statistically independent. One can then show that the density $P_{\tau+\delta\tau}({\sf M})$ can be obtained as the convolution [42]

[TABLE]

where $d\mu({\sf M})$ denotes the invariant (Haar) measure with respect to the transfer matrix group, with the property that $d\mu({\sf M})=d\mu({\sf M}_{0}{\sf M})$ with ${\sf M}_{0}$ being any other (fixed) transfer matrix.

Equation (3.11) corresponds to the Smoluchowsky equation for the present (Markovian) process. As it is well known from the theory of Brownian motion, the Smoluchowsky equation implies that the probability distribution of the transfer matrix parameters, and consequently that for the ${\sf R}$-matrix parameters $P_{\tau}(\{\lambda_{a}\})$, evolves with time according to the (forward) Kramers-Moyal expansion

[TABLE]

The $\delta\lambda_{a}$ in the above expression should be thought of as the “small” increment of the parameters $\lambda_{a}$ (relative to their values before the addition of the interval $\delta\tau$) due to a single additional scatterer in the time interval $\delta\tau$ (weak scattering). The expectation value is over the probability distribution describing the properties of the scatterer in the interval $\delta\tau$. For a proof of the expansion in eq. (3.12) we refer the reader to Appendix A.1. In the weak scattering limit, the Kramers-Moyal expansion can be curtailed at second order, and eq. (3.12) may be referred to as the Fokker-Planck (FP) equation for $P_{\tau}(\{\lambda_{a}\})$.444In the single field limit it can be explicitly shown that $\langle\delta n\rangle=(\mu\,\delta\tau)(1+2n)$, $\langle(\delta n)^{2}\rangle=(\mu\,\delta\tau)2n(1+n)+\mathcal{O}[(\mu\,\delta\tau)^{2}]$, $\langle(\delta n)^{2k}\rangle=\mathcal{O}[(\mu\,\delta\tau)^{k}]$ and $\langle(\delta n)^{2k+1}\rangle=\mathcal{O}[(\mu\,\delta\tau)^{k+1}]$ for $k>1$, where $\mu$ denotes the local mean particle production rate [34].

It is often not necessary to get an explicit solution for $P$. Instead we can integrate eq. (3.12) to obtain equations for the expectation values of quantities of interest. Let us consider a function $F(n)$ where $n$ is the total occupation number of the fields. Note that $n=(1/2)\sum_{a=1}^{{N_{\rm{f}}}}(f_{a}-1)$ only depends on the ${N_{\rm{f}}}$ eigenvalues $\{f_{a}\}$ and not on the angular variables. Multiplying eq. (3.12) with $F(n)$ and integrating over all ${N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ variables we arrive at

[TABLE]

where $\langle\ldots\rangle=\int\prod_{a=1}^{{N_{\rm{f}}}^{2}+{N_{\rm{f}}}}d\lambda_{a}P_{\tau}(\{\lambda_{a}\})(\ldots)$. In deriving the above expression we have repeatedly used integration by parts, dropped the boundary terms, taken advantage of the fact that the ${N_{\rm{f}}}^{2}$ variables of ${\sf u}(\{\theta_{\alpha}\})$ form a compact manifold, and that $F(n)$ is a function of the ${N_{\rm{f}}}$ eigenvalues $\{f_{a}\}$ only. We have also used the expression for $n$ in terms of the eigenvalues $\{f_{a}\}$. Note in particular that the sums are over ${N_{\rm{f}}}$ variables only, however the integration is over all ${N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ variables. Moreover the coefficients $\langle\ldots\rangle_{\delta\tau}$ as well as $P$ might be non-factorizable functions of all ${N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ variables.

Our interest will be focused on the so-called typical occupation number, defined as

[TABLE]

The quantity $n_{\rm typ}$ is a good measure of the typical number of particles produced, as its variance grows slower than the square of the mean. For the mean occupation number $\langle n\rangle$ the ratio $(\langle n^{2}\rangle-\langle n\rangle^{2})/\langle n\rangle^{2}$ grows with time (see the discussion in [34]).

With $n_{\rm typ}$ in mind, consider the function $F(n)=\ln(1+n)$. Due to the relation (3.10), eq. (3.13) can be then be written as

[TABLE]

To make further progress, let us first turn to more detailed understanding of $\delta f_{a}$ and other $\delta\lambda_{a}$.

3.4 Increments in ${\sf R}$

Our next task is to find expressions for $\delta\lambda_{a}$. To this end, let us consider the change in ${\sf R}$ due to a single additional scattering (see fig. 2)

[TABLE]

We will often suppress $``(j)"$ as an argument of the matrices.

The increment in ${\sf R}$, $\delta{\sf R}=\delta({\sf U}{\sf F}{\sf U}^{\dagger})$, can be expressed in terms of increments in the eigenvalue matrix ${\sf F}$ and the angular matrix ${\sf U}$ (up to second order in the increments) as follows:

[TABLE]

Comparing the expressions for $\delta{\sf R}$ in eqns. (LABEL:eq:dR1) and (3.17), we get

[TABLE]

where

[TABLE]

Equation (LABEL:eq:dFdU) expresses the increments $\delta{\sf F}$ and $\delta{\sf U}$ in terms of ${\sf G}$ which contains information about the additional scattering via ${\sf m}_{j+1}$, as well as ${\sf F}$ and ${\sf U}$. As is evident from the above equation the increments $\delta{\sf F}$ and $\delta{\sf U}$ will each contain ${\sf F}$- and ${\sf U}$-related variables. As mentioned earlier, each additional scattering is assumed to make a small difference, hence it is convenient to define the following new quantities with the perturbation order in mind

[TABLE]

which then yields

[TABLE]

Recall that ${\sf F}$ and ${\sf U}$ are $2{N_{\rm{f}}}\times 2{N_{\rm{f}}}$ block matrices, which are constructed out of ${N_{\rm{f}}}\times{N_{\rm{f}}}$ blocks of functions of ${\sf u}$ and ${\sf f}$ (and ${\tilde{\sf f}}$) matrices. The above expressions for $\delta{\sf U}^{(i)}$ and $\delta{\sf F}^{(i)}$ can be unpacked in terms of ${\sf f}$ and ${\sf u}$ directly to yield

[TABLE]

where ${\sf g}^{(i)}$ is the top left ${N_{\rm{f}}}\times{N_{\rm{f}}}$ block of ${\sf G}^{(i)}$, and $\tilde{{\sf g}}^{(i)}$ is the top right ${N_{\rm{f}}}\times{N_{\rm{f}}}$ block of ${\sf G}^{(i)}$. Although the previous system represents $4{N_{\rm{f}}}^{2}$ equations for the $2{N_{\rm{f}}}({N_{\rm{f}}}+1)$ unknowns, there is no issue of overdetermination as, for example, the off-diagonal components of the first and second equations can be checked to be degenerate, leading to a system of ${N_{\rm{f}}}({N_{\rm{f}}}+1)$ equations for the components of $\delta{\sf f}^{(1)}$ and $\delta{\sf u}^{(1)}$, and similarly for the second-order terms.

At the level of the variables that parametrize ${\sf f}$ and ${\sf u}$, the above equations can be solved for each $\delta\lambda_{a}^{(1)}$ and $\delta\lambda_{a}^{(2)}$ in terms of the parameters characterizing the scatterers (elements of ${\sf m}_{j}$) as well as $\{\lambda_{b}\}$. In component form the solution can be written explicitly. The first order correction to ${\sf f}$ can be directly obtained from the diagonal elements of eq. (3.22a),

[TABLE]

For the second-order change in ${\bf\sf f}$, eq. (3.22a) yields for the off-diagonal components

[TABLE]

In turn, this expression can be substituted into eq. (3.22b) to obtain the diagonal entries, which correspond to

[TABLE]

The first-order corrections for the entries of ${\sf u}$ can be obtained from the off-diagonal components of eq. (3.22a), and from the diagonal entries of eq. (3.22c). These lead to the relations eq. (3.24) and

[TABLE]

These two equations form a linear system for the perturbations of ${\sf u}$, which can be solved explicitly as

[TABLE]

Analogously, one can solve for the second-order corrections for ${\sf u}$ components, which are given by

[TABLE]

To make further progress, we need an explicit form of ${\sf m}_{j}$. Below we will choose Dirac-delta scatterers. Once the form of ${\sf m}_{j}$ is specified, we can solve, in principle, for the increments in the parameters of $\delta\lambda_{a}^{(i)}$ for arbitrary ${N_{\rm{f}}}$. In practice this is a non-trivial calculation made particularly onerous by the large number of parameters (${N_{\rm{f}}}^{2}$) in ${\sf u}$ and (${N_{\rm{f}}}$) in ${\sf f}$.

3.5 Dirac Delta Scatterers

The general formalism we have presented is independent of the precise form of the scatterers (${\sf m}_{j}$), but having a concrete simple example in mind makes things simpler to present and also allows for explicit calculations. With this in mind, as long as we restrict ourselves to wavenumbers $k\ll w^{-1}$, we approximate $m^{\rm s}_{ab}(\tau)$ as follows:

[TABLE]

where ${N_{\rm s}}$ is the number of events, $\tau_{j}$ are uniformly distributed, and $\delta(\tau-\tau_{j})$ are Dirac Delta functions. For each $\tau_{j}$, the ${N_{\rm{f}}}\times{N_{\rm{f}}}$ elements of $\Lambda_{ab}(\tau_{j})$ which characterize the strength of the scatterers are drawn from some distribution. We assume that the distributions are identical and

independent for all $\tau_{j}$. Note that the elements of $\Lambda_{ab}(\tau_{j})$ must be symmetric with respect to $a$ and $b$. The scaling with $2\sqrt{\omega_{a}\omega_{b}}$ is for future convenience. We will assume that for each $\tau_{j}$

[TABLE]

With the parametrization (3.29), the scattering strength variances $\sigma_{ab}^{2}$ are functions of the wavenumber $k$ and the corresponding scalar field masses (see fig. 3).

By imposing appropriate junction conditions on $\chi^{a}$ at $\tau=\tau_{j}$, the $2{N_{\rm{f}}}\times 2{N_{\rm{f}}}$ transfer matrix at ${\sf M}_{j}$ for the Dirac delta function scatterers evaluates to

[TABLE]

where

[TABLE]

The right hand sides of eqns. (3.22a)-(3.22d) for the delta function scatterers can be written as

[TABLE]

where we have defined for the $j+1$-th scatterer

[TABLE]

Note that ${\sf u}$ is evaluated at $j$ rather than $j+1$.

With the explicit form of the local transfer matrix at hand, we can immediately calculate the disorder-averaged quantities that appear in the expression for the typical occupation number defined in eq. (3.15), for any number of fields. Namely, for the coefficients that depend on the first-order correction to ${\sf f}$, we can rewrite eq. (3.23) in terms of the matrix ${\sf g}^{(1)}$ given in eq. (3.33). The following explicit expression for the perturbation $\delta f_{a}^{(1)}$ is obtained:

[TABLE]

In the second line we have made use of the fact that ${{\sf p}}$ is diagonal. Consider now the expectation value $\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle_{\delta\tau}$. This product contains phase factors due to the presence of the components of ${\sf p}$, which can be explicitly evaluated, as each $\tau_{j}$ is assumed to be uniformly distributed in the interval $[(j-1)\delta\tau,j\delta\tau]$. For example,

[TABLE]

where $\Delta\omega=\omega_{a}+\omega_{b}-\omega_{c}-\omega_{d}$. In what follows, and throughout this paper, we will implicitly assume that the time period of oscillations is much smaller than the mean free path determined by the separation between events, $\omega\delta\tau\gg 1$; this is consistent with the treatment of the dynamics as free in-between scatterings. This assumption implies that, within the product $\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle_{\delta\tau}$, we can disregard terms that contain unpaired ${\bf\sf p}$ and ${\bf\sf p}^{*}$, as they will vanish after averaging over scatterers; additionally $p^{*}_{a}p_{b}\rightarrow\delta_{ab}$ upon averaging. Recalling the statistical properties of $\Lambda_{ab}$ in eq. (3.30), we can then write555In eq. (3.37) (and in what follows) we have denoted $\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle_{\delta\tau}\equiv\lim_{\delta\tau\rightarrow 0}\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle_{\delta\tau}/\delta\tau$. Whenever there is no potential for confusion we will also suppress the $\delta\tau$ subindex, to avoid cumbersome expressions.

[TABLE]

The second-order change in ${\bf\sf f}$ due to a scattering can be computed in an analogous way, directly upon averaging (3.25) over scatterings, a procedure which leads to the following expression for the mean second-order perturbation,

[TABLE]

where

[TABLE]

and

[TABLE]

Substitution of (3.37) and (3.38) into (3.15) then yields

[TABLE]

This is an explicit expression for the evolution of the typical occupation number for arbitrary ${N_{\rm{f}}}$ and $\sigma_{ab}$ (assuming Dirac-delta function scatterers). However, note that the angular brackets denote integration with respect to probability distribution $dP(\{{\sf f},{\sf u}\})$. In the absence of an explicit parametrization for ${\sf u}$ and the solution to the FP equation, equation (3.41) is of limited applicability. In the next section we will evaluate it under the assumption of a ${\sf u}$-flat probability distribution and for a statistically isotropic ${\sf\Lambda}$ (i.e. $\sigma^{2}_{ab}=\sigma^{2}$). In the subsequent sections, we will evaluate this expression in general for the ${N_{\rm{f}}}=1,2$ cases, and for ${N_{\rm{f}}}>2$ under some restrictive conditions for the scattering strengths $\sigma_{ab}^{2}$.

4 The Maximum Entropy Approximation

In the previous section, we have laid the foundations necessary to compute the solution of the FP equation that determines the instantaneous dependence of the probability distribution $P$ with respect to the parameters in ${\sf R}$. However, before proceeding to such analysis, we will review in this section the predictions of the Maximum Entropy approximation, and state our expectations for when we eventually compare them against the analytical result obtained from the integration of the FP equation.

The probability density $P$ associated with the MEA corresponds to that which maximizes the (Shannon) entropy functional [42],

[TABLE]

subject to the constraints that the local mean particle production rate is known, and that the evolution of the transfer matrix ${\sf M}$ is under perturbative control, ${\sf M}_{\tau+\delta\tau}\rightarrow{\sf M}_{\tau}$ as $\delta\tau\rightarrow 0$. With this definition, the MEA recipe is said to provide the least biased estimate of the $P$ consistent with the (fixed) local production rate. Under these approximations, it can be shown that the probability density is independent of ${\sf u}$, $P(\{{\sf f},{\sf u}\})=P(\{{\sf f}\})$. Moreover, for the typical occupation number $n_{\rm typ}$, the MEA implies the late-time result [34]

[TABLE]

where $\mu_{j}$ denotes the instantaneous mean particle production rate,

[TABLE]

Here $\langle n_{j}\rangle_{\delta\tau}$ corresponds to the local change in the occupation number due to a single scattering event. For Dirac-delta scatterers, it can be easily evaluated in terms of the scattering strengths as

[TABLE]

where in the last equality we have used the defining relation (3.30). Thus, the MEA-based result is

[TABLE]

Factorizability, the Haar measure & the MEA: The independence of $P$ on ${\sf u}$ within MEA implies that the differential probability may be written as

[TABLE]

where $d\mu({\sf u})$ denotes the invariant (Haar) measure on the group of unitary matrices $U({N_{\rm{f}}})$. The previous expression is a key prediction of the MEA approximation, and, as we now show, it is the single result responsible for the trace formula in eq. (4.4). For these purposes, let us now assume that the probability density satisfies eq. (4.6), without assuming that it extremizes the Shannon entropy in eq. (4.1). As it was shown in the previous section, in the case of Dirac delta scatterers, the general expression for time-evolution of the function $F(n)=\ln(1+n)$ can be written as eq. (3.41) in terms of an average with respect to the measure $dP(\{{\sf f},{\sf u}\})$. Under the assumption (4.6), the averages for the monomial functions of the components of ${\sf u}$ and ${\sf u}^{\dagger}$ can be immediately computed in terms of the Weingarten relations [42], which for monomials of ranks 2 and 4 correspond to

[TABLE]

and

[TABLE]

By making use of these relations, the general expression (3.41) reduces to

[TABLE]

In the right-hand side we recognize the trace in eq. (4.4). The prefactor is the expectation value of a complicated function of the eigenvalues $f_{a}$. Nevertheless, in the late time limit, it is expected that the occupation number $n\gg 1$. Moreover, in the asymptotic regime, it is expected that the eigenvalues of ${\sf M}$, itself a product of a large number of transfer matrices, separate, and the asymptotic growth is purely exponential, implying that we can assume that a single eigenvalue dominates the evolution at late times [45, 46]. Recalling the relation (3.10) we can therefore write

[TABLE]

which is none other than the MEA result (4.5). This result will be of use in the coming sections.

${\sf\Lambda}$-isotropy and the MEA: The previous calculation demonstrates that factorizability and ${\sf u}$-isotropy (see eq. (4.6)) is sufficient to reproduce the MEA prediction for the evolution of the typical occupation number (a fact already discussed in [34]). However, it tells us nothing about the probability density $P(\{{\sf f}\})$. On one hand, within MEA, its form would be determined by the requirement that it extremizes the entropy functional (4.1) and then solving the resulting FP equation for $P(\{{\sf f}\})$. On the other hand, in terms of the microscopic model built in section 3, $P(\{{\sf f}\})$ would be determined by the particular set of scattering strengths $\sigma_{ab}$ for which the solution of the FP equation (3.12) is ${\sf f}$-${\sf u}$ factorizable and compatible with (4.6). It is a-priori unclear under what conditions on the microscopic model that one recovers the MEA-based $P(\{{\sf f}\})$.

It has been conjectured that the Maximum Entropy approximation is equivalent to the assumption that the interacting fields are statistically equivalent and maximally mixed with a random but uniform distribution of non-adiabatic events [47, 34]. In the context of the present study, it would correspond to $\sigma_{ab}^{2}=\sigma^{2}$ (i.e. a statistically isotropic ${\sf\Lambda}$). Using our general expression in eq. (3.41), we can analytically verify this conjecture in the context of the typical occupation number. Using $\sigma_{ab}^{2}=\sigma^{2}$ in eq. (3.41), we immediately get

[TABLE]

where the second line is valid in the late-time limit under the assumption of a dominant eigenvalue. That is, $n=(1/2)\sum_{a}^{{N_{\rm{f}}}}(f_{a}-1)\rightarrow f_{\rm max}/2$. In the third line we have used a straightforward substitution of $\sigma_{ab}^{2}=\sigma^{2}$ in eqns. (4.4). We have shown that the result in eqns. (4.11) corresponds to the ${\sf\Lambda}$-isotropic limit of the MEA expression (4.5). In the upcoming sections, we will also realize that ${\sf\Lambda}$-isotropy $\rightarrow{\sf u}$-isotropy (in a somewhat restricted sense relevant for the typical occupation number). Fig. (4) is a graphical illustration of which assumption leads to what conclusion.

While we have proven that in the ${\sf\Lambda}$-isotropic case, the MEA result for the occupation number matches the result from explicit calculation, we have yet to see how the exact results deviate (or even if they deviate) from MEA based results when ${\sf\Lambda}$-isotropy is not assumed. In the following sections, we will carry the full analysis assuming neither MEA nor ${\sf u}$ nor ${\sf\Lambda}$-isotropy for the ${N_{\rm{f}}}=1,2$ cases, and in a restricted sense, for the general ${N_{\rm{f}}}>2$ case.

5 A Warm-up Example: the ${N_{\rm{f}}}=1$ Case

Let us consider first the simplest case, that of a single field. The unitary matrix ${\sf u}$ reduces simply to a phase, $u=e^{i\phi}$, meaning that the ${\sf R}$ matrix contains only the two real variables $\{\lambda_{a}\}=\{f,\phi\}$. Straightforward substitution into eqns. (3.22a)-(3.22d) leads to the following expressions for $\{\delta\lambda_{a}\}$,

[TABLE]

where for a single delta function scatterer we have (see eq. (3.33))

[TABLE]

With the results at hand, we can calculate the coefficients of the FP equation by taking averages over the properties of the scatterer as well as the location of the scatterers. For this purpose we will denote the ($\tau_{j}$-independent) variance of $\Lambda_{j}$ by

[TABLE]

As it was discussed around eq. (3.36), under the assumption $\omega\delta\tau\gg 1$, all time-dependent phase factors average to zero. We then find that

[TABLE]

where we have lost the $\delta\tau$ sub-index for notational simplicity. It is worth noting that all expectation values are independent of the angular variable $\phi$. Thus, let us consider the marginal probability distribution $\int d\phi\,P(f,\phi;\tau)$, which somewhat abusing notation we will continue denoting by $P$. The FP equation (3.12) satisfied by $P$ becomes

[TABLE]

Although the exact solution of eq. (5.5) can be found in integral form [48], we are interested in its behavior in the asymptotic regime $\tau\rightarrow\infty$; in this limit we expect that $n\gg 1$ and hence $f\gg 1$. Therefore, in the right-hand side we may replace $(f^{2}-1)\rightarrow f^{2}$. Further switching the independent variable from $f$ to $n$, we find the FP equation

[TABLE]

which for an initial condition $P(n;0)=\delta(n)$ has the log-normal solution,

[TABLE]

Equation (5.7) recovers the single field result derived in [34]; where it was shown that the assumption of maximum entropy implies the independence of the probability density on the angular variable $\phi$, which inevitably leads to the solution in eq. (5.7) after assuming that the local mean particle production rate is known.

Having an explicit expression for the probability density, one can in principle calculate all moments of the distribution. For the typical occupation number in eq. (3.14), we have

[TABLE]

Although the simplicity of the single field case allowed us to calculate a closed form for the probability distribution $P$, the complexity of the analysis increases many-fold even upon the addition of a single additional field, and we will not obtain the full form for $P$ even in the two-field scenario. Nevertheless, as we argued in section 4, the full solution of the FP equation is not needed to obtain an expression for the (total) occupation number. Namely, if we consider the function of the occupation number $F(n)=\ln(1+n)$, we can immediately write the late time behavior of eq. (3.15) as

[TABLE]

which trivially implies (5.8). We will continue to exploit the relation in eq. (3.15) in the two- and multi-field cases discussed below.

To summarize, without considering the MEA, we arrive at a result identical to the one obtained by using the MEA. The equivalence to the MEA arises from the fact the disorder averaged correlators are independent of ${\sf u}$ (parametrized by $\phi$ here); the ${N_{\rm{f}}}=1$ is special is this sense. Let us now consider the more interesting and more involved case with ${N_{\rm{f}}}=2$.

6 Beyond MEA: The ${N_{\rm{f}}}=2$ Case

For the two field case (${N_{\rm{f}}}=2$), the unitary matrix ${\sf u}$ can be parametrized in terms of the $U(2)$ Euler angles as follows,

[TABLE]

where $\phi,\varphi\in[0,\,2\pi)$, $\theta\in[0,\,\pi)$, and $\psi\in[0,\,4\pi)$ [49]. The ${\sf f}$ matrix is diag$\{f_{1},f_{2}\}$. Therefore the ${\sf R}$ matrix contains six real variables,

[TABLE]

6.1 Calculation of the Correlators

In complete analogy with the single field case, and after a significant amount of algebra, we find the following expressions for $\{\delta\lambda_{a}\}$. First, at leading order (see eqns. (3.22a)-(3.22d))

[TABLE]

where $\Delta f=f_{1}-f_{2}$. At next to leading order, we have

[TABLE]

The leading order coefficients should be subsituted in the expressions above to get the second order contributions explicitly. Note that we have only listed a select number of $\delta\lambda_{a}$. The rest are provided in Appendix A.2.

For evaluating the ${\sf g}^{(i)}$ and $\tilde{{\sf g}}^{(i)}$ matrices we need the transfer matrix for a single delta function scatterer with ${N_{\rm{f}}}=2$. This is given by eq. (3.31) with

[TABLE]

Note that $\Lambda_{12}(\tau_{j})=\Lambda_{21}(\tau_{j})$. In order to calculate the coefficients of the FP equation we will need the variances of the entries of ${\sf\Lambda}_{j}$. These $j$-independent variances will be denoted as follows,

[TABLE]

Care should be taken regarding the notation $\sigma_{i}^{2}$ and $\sigma_{ab}^{2}$. In particular

[TABLE]

We apologize for what might seem like confusing notation here; it is a direct consequence of the way we defined $\sigma^{2}_{ab}$ in eq. (3.30).666This notation will also be useful for ${N_{\rm{f}}}>2$ case. Moreover, we assume that $\langle e^{i\omega_{a}\tau_{j}}\rangle_{\delta\tau}=0$ for $a=1,2$ as well as $\langle e^{i(n\omega_{1}-m\omega_{2})\tau_{j}}\rangle_{\delta\tau}=0$ as long as $n\omega_{1}\neq m\omega_{2}$ with $n,m\in$ Integers.777We have found that our numerical results match the analytical predictions even when $\omega_{1}=\omega_{2}$. In addition, numerical evidence in the multifield case also indicates that the equal $\omega_{a}$ case is not special.

With these assumptions at hand, the computation of the coefficients of the FP equation proceeds in a straightforward way, albeit requiring a significant amount of algebra. The full list of correlators is provided in Appendix A.2. Here, we simply list the ones that are directly needed to complete the calculation:

[TABLE]

where the $\gamma$- and $l$-functions, which depend only on the polar angle $\theta$ and the scattering strengths $\sigma_{1,2,\perp}^{2}$, are defined in the Appendix A.2 in eqns. (A.9a)-(A.10d). We list $\gamma_{1}(\theta)$ and $l_{1}(\theta)$ here to provide an idea of their structure.

[TABLE]

6.2 The Typical Occupation Number

Although we have now at hand all coefficients of the FP equation (3.12), we will not attempt its solution. Instead, we will focus on the calculation of the typical particle production rate at late times (see eq. (3.15)), when it is assumed to have saturated to a constant value. In the limit of $\tau\rightarrow\infty$ (or equivalently $N_{s}\rightarrow\infty$) we expect that $n\gg 1$ and that a single eigenvalue dominates the evolution at late times, $f_{1}\gg f_{2}$ or $f_{2}\gg f_{1}$. In the former case, using the coefficients we found in the previous section, we obtain

[TABLE]

(the case $f_{2}\gg f_{1}$ is recovered by substituting $1\rightarrow 2$). The previous expression implies that, remarkably, the late-time particle production rate depends only on the probability distribution of the angular variable $\theta$. Therefore, to evaluate the right hand side, we need only the marginal distribution

[TABLE]

To calculate $w(\theta)$ we need to solve the Fokker-Planck equation (3.12), which in the $\tau\rightarrow\infty$ limit can be approximated as stationary. Integrating both sides of eq. (3.12) over $\{f_{1},f_{2},\phi,\varphi,\psi\}$ one finds the differential equation satisfied by the marginal probability density $w(\theta)$,

[TABLE]

where the coefficient functions are given by

[TABLE]

and

[TABLE]

In the step going from eq. (6.19) to eq. (6.20) we have assumed level repulsion ${\rm sgn}(\Delta f)={\rm sgn}(\Delta\sigma^{2})$. This assumption is justified at the perturbative level as long as $\sigma_{\perp}^{2}\leq\sigma^{2}_{\rm max}$ where $\sigma_{\rm max}^{2}=\max\{\sigma_{1}^{2},\sigma_{2}^{2}\}$. With the coefficient functions at hand, one can solve the differential equation (6.18). Integration yields the equation

[TABLE]

where $j$ denotes the current in the $\theta$-direction. Expansion around $\theta=0$ reveals that the marginal distribution would behave as $w(\theta)\sim\theta-(j/2\sigma_{\perp}^{2})\,\theta\ln\theta$, meaning that for $j\neq 0$ the probability distribution becomes negative in $(0,\pi)$. Thus, the normalizable, positive definite marginal probability $w(\theta)$ is found by straightforward integration of eq. (6.22) with $j=0$:

[TABLE]

where for general values of the exponents $\nu$ the normalization constant $\mathcal{N}$ must be determined numerically. In the above expression

[TABLE]

The $w(\theta)$ calculated above can now be used to evaluate the right hand side of eq. (6.16). In the right panel of figure 5, we show the two-field theoretical prediction for $\langle\ln(1+n)\rangle$ for different values of the ratios $\sigma_{a}/\sigma_{1}$ (solid lines), and compare them to the numerical results (points). For a description of the numerical method see Appendix A.3. Since $w(\theta)$ only depends on the ratio of $\sigma_{a}/\sigma_{1}$, in the vertical axis the particle production rate has been normalized with respect to $\sigma_{1}^{2}$.

The left panel of fig. 5 shows a numerical solution for the evolution of the occupation number $n$ for an ensemble of fifty different realizations of randomly generated scatterers, verifying that $n$ executes a drifting random walk. The analytic result with $\sigma_{2}^{2}=\sigma_{1}^{2}/10$ and $\sigma_{\perp}^{2}=\sigma_{1}^{2}/4$ is shown as a orange line, and is clearly a good fit for the instantaneous average occupation number. Even more remarkable is the agreement between the numerical and analytic results in the right panel of fig. 5. For comparison, we have also provided the simple MEA inspired predictions, shown as the diagonal dashed gray curves. From eq. (4.5) with ${N_{\rm{f}}}=2$, the gray curves are given by

[TABLE]

It is clear from fig. 5 that, for $\sigma_{1}^{2}=\sigma_{2}^{2}=2\sigma_{\perp}^{2}$ (equivalently all $\sigma^{2}_{ab}$ being the same), our calculation and the MEA based agree exactly. This is to be expected based on our proof in section 4. Indeed, in the ${\sf\Lambda}$-isotropic limit $\sigma_{1}^{2}=\sigma_{2}^{2}=2\sigma_{\perp}^{2}$ , the typical occupation number (6.16) trivially reduces to the MEA result,

[TABLE]

Moreover, in this all $\sigma^{2}_{ab}$ equal limit, the angular dependence of the late-time probability density is given simply by

[TABLE]

where $d\mu(U(2))$ denotes the invariant Haar measure for $U(2)$. This $U(2)$-‘flat’ probability density is also consistent with the MEA assumption of isotropy in the fundamental strip. Therefore, the MEA result is fully recovered assuming only ${\sf\Lambda}$-isotropy.

It is worth discussing the agreement/difference between the MEA based results and the current calculation further. We note that the MEA based (top) dashed gray line is ‘tangential’ to the actual analytical result, implying that its range of validity extends significantly beyond the values $2\sigma_{\perp}^{2}=\sigma_{1}^{2}=\sigma_{2}^{2}$. Also note that looking at the second from the top gray dashed line, the deviation of the actual result from the MEA based is small unless $\sigma_{\perp}$ is significantly different from $\sigma_{1},\sigma_{2}$. These observations show that MEA-inspired ansatz (4.5) is a good approximation to the exact result in the limit when all scattering strengths are comparable, but when they differ significantly care must be taken in using the approximation (at least in the two field case).

We also observe that, although there exist other special values of $\sigma_{a}^{2}$ (for example, large $\sigma^{2}_{\perp}$) for which the MEA based results intersect with the results of this section, they do not correspond to special points (i.e. points of enhanced symmetry). We also note that the asymptotic behavior of the analytical curves in fig. 5, which for $\sigma_{\perp}^{2}\gg\sigma_{1}^{2},\sigma_{2}^{2}$ results in $\partial_{\tau}\langle\ln(1+n)\rangle\sim 4\sigma_{\perp}^{2}/\pi$, is different from the MEA prediction $\partial_{\tau}\langle\ln(1+n)\rangle\sim 4\sigma_{\perp}^{2}/3$. So the MEA lines should not necessarily be interpreted as asymptotic curves in the large $\sigma_{\perp}$ limit.888It is unclear to us whether this regime is allowed by our approximations. Numerically we found that when $\sigma_{\perp}^{2}\gtrsim 10\sigma_{1}^{2}$, we start seeing deviations from numerical simulations. This large $\sigma_{\perp}$ regime might violate our assumption that ${\rm sgn}(\Delta f)={\rm sgn}(\Delta\sigma^{2})$.

6.3 Summary

Under the weak scattering approximation, we have solved the two field case in general without using the MEA, or ${\sf\Lambda}$-isotropy. As was to be expected from our proof in section 4, we found precise agreement with MEA results when we have ${\sf\Lambda}$-isotropy and deviations from it otherwise. We verified that the MEA based trace formula for the typical occupation number (eq. (4.5)) provides a good approximation to the more accurate calculation when the strengths of the self and cross couplings are all comparable to each other. Our analytic expressions (up to quadratures) here match the results from numerical simulations, and correctly capture the deviations from the MEA result.

7 Beyond MEA: Generalization to Many Fields

The analysis introduced in section 6 for two fields can be extended to an arbitrary number of fields in a straightforward way, if the parametrization of the ${N_{\rm{f}}}\times{N_{\rm{f}}}$ unitary matrix ${\sf u}$ is known explicitly. However, for large ${N_{\rm{f}}}$ the computational expense increases dramatically, as now $({N_{\rm{f}}}^{2}+{N_{\rm{f}}})({N_{\rm{f}}}^{2}+{N_{\rm{f}}}+3)/2$ correlators are necessary to solve the Fokker-Planck equation (3.12) for the probability density $P$. Therefore, in order to obtain a result that will allow us to draw generic conclusions for any number of fields, we will make the following simplifying assumptions: all off-diagonal couplings are equally distributed, and only one of the ${N_{\rm{f}}}$ fields has a statistically different self-coupling; that is

[TABLE]

Here we have associated $\sigma_{1}$ with the ${N_{\rm{f}}}$-th field for convenience; recall that the extra factors of $1/2$ for the diagonal entries stem from the definition in eq. (3.30). In what follows we will show that under the previous assumptions, and some seemingly reasonable conditions, only four correlators are necessary to calculate the typical occupation number in the general case.

7.1 Calculation of the Correlators 999In this section we denote the total number of fields by $N$ instead of ${N_{\rm{f}}}$ for notational simplicity.

Let us explicitly calculate first the coefficients of the FP equation that depend on the first-order correction to the eigenvalues $f_{i}$. Luckily, the hardest part of the work has already been carried out in section 4, which led to the expression (3.37) for the correlator $\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle$. Let us then consider the late time limit $N_{s}\rightarrow\infty$, in which, as in the two-field scenario, we expect one eigenvalue $f_{a}$ to be overwhelmingly dominant. We will identify this dominant eigenvalue with the $N$-th one, $f_{N}\gg\{f_{a},1\}$. In this limit the sum $\sum_{a,b=1}^{N}\langle\delta f_{a}^{(1)}\delta f_{b}^{(1)}\rangle_{\delta\tau}/\delta\tau$ will be dominated by its last term, which in turn, under the assumption (7.1), can be approximated as

[TABLE]

In the previous expression $\stackrel{{\scriptstyle\infty}}{{=}}$ denotes an equality valid in the $N_{s}\rightarrow\infty$ regime.

To go further, an explicit parametrization of a generic $N\times N$ unitary matrix is needed. In analogy with the two field case, it is convenient to express the unitary matrix ${\sf u}$ in terms of the (generalized) Euler angular variables. As outlined in [50], an Euler parametrization of an element of $SU(N)$ can be iteratively constructed as follows:101010It is sufficient to restrict ourselves to the $SU(N)$ subgroup of $U(N)$ as the extra phase is canceled in the products appearing on (7.2).

[TABLE]

where the relevant generators of the $su(N)$ Lie algebra are given by

[TABLE]

This parametrization implies that the magnitudes of the elements of the last column of ${\sf u}$ can be written as

[TABLE]

where, in order to facilitate the comparison with the two-field result, we have denoted $\alpha_{2(N-1)}\equiv\theta$. Therefore, the dominant correlator in eq. (7.2) evaluates to

[TABLE]

Here $\mathcal{F}(\Omega_{N})$ is a function of the Euler angles other than $\theta$ appearing in eq. (7.6). It can be calculated iteratively as

[TABLE]

with $\mathcal{F}(\Omega_{2})=1$. In what follows we will denote $\mathcal{F}(\Omega_{N})$ by $\mathcal{F}_{\Omega}$.

The general expression for the second-order change $\delta f_{a}^{(2)}$ is given by eq. (3.38). Specializing to the particular case with eq. (7.1), in the late time limit it is straightforward to show that the perturbations to the non-dominant eigenvalues are negligible, $\langle\delta f_{a}^{(2)}\rangle\stackrel{{\scriptstyle\infty}}{{=}}0$ for $a\neq N$, while for the dominant eigenvalue we have

[TABLE]

The first- and second-order corrections for the entries of ${\sf u}$ have been written in eqns. (3.27) and (3.28). In principle then, the angular correlators can be evaluated by solving these two equations in terms of the angular perturbations; namely

[TABLE]

In practice, the solution for this system of equations and the subsequent evaluation of the first and second order correlators cannot be done for a generic $N$ and must be performed on a case by case basis (c.f. the two-field calculation). Nevertheless, for the $\theta$ correlators, eqns. (3.27), (3.28) and (7.6) can be used to relate $\delta u_{NN}$ and $\delta\theta$. In complete analogy with the previous calculation for the $\delta f$ expectation values, in the $f_{N}\gg f_{a}$ limit one obtains

[TABLE]

and

[TABLE]

7.2 The Typical Occupation Number

In the previous section we have shown that under the assumptions stated in eqns. (7.1) we can explicitly calculate the eigenvalue correlators in the $\tau\rightarrow\infty$ limit, as well as the set of correlators for one angular variable. Following the same argument for the two-field case discussed in section (6.2), we can calculate the typical occupation number in the ${N_{\rm{f}}}$-field case by considering the late-time limit of (3.15). Substitution in eq. (3.15) of the expressions (7.7) and (7.9) we obtain

[TABLE]

where the functions $l(\theta)$ and $\gamma(\theta,\Omega_{N})$ are the generalization of their two-field counterparts,

[TABLE]

Note first that, as expected, in the isotropic-${\sf\Lambda}$ limit, the left-hand side of eq. (7.13) is angle-independent, and the typical occupation number takes the MEA value consistent with eq. (4.5):

[TABLE]

As a second particular example, in the large ${N_{\rm{f}}}$-limit, the angular dependence of the expression inside the angular brackets in eq. (7.13) becomes subdominant. Therefore, one obtains the universal limit

[TABLE]

In the general case, evaluation of eq. (7.13) requires the knowledge of the late-time marginal probability density which is a function of all the angles that parametrize (7.6), due to the explicit angular dependence of $\mathcal{F}_{\Omega}$. Analogously to the two-field case, such function could in principle be found for a given ${N_{\rm{f}}}$ after the inversion of eq. (7.10) and the computation of all angular correlators. We have not attempted to take on such a task, and instead we will go beyond MEA by calculating the right-hand side of eq. (7.13) in complete analogy with the two-field case. Namely, we will assume that in the stationary $N_{s}\rightarrow\infty$ limit, the dependence on angles other than $\theta$ of the probability density is factorizable. The compact nature of the angular variables and the previous assumption imply that, similarly to the two-field case, the required marginal distribution $w(\theta)$ may be obtained from the solution of the single-variable Fokker-Planck equation (6.18), where the coefficient functions $c_{1}$ and $c_{11}$ are now determined from the expectation values noted in eqns. (7.11) and (7.12), and correspond to

[TABLE]

In the previous expressions, we have introduced the notation $\langle\mathcal{F}_{\Omega}\rangle=\mathcal{F}$. To be able to go further, we need the explicit form of this expectation value. Although in general this can only be achieved by solving the full FP equation and averaging, remarkably, in the ${\sf u}$-isotropic limit, a closed-form for the expectation value of the angular function $\mathcal{F}_{\Omega}$ can be calculated; for the interested reader the details are shown in appendix A.4. It results in the following,

[TABLE]

We will then consider the additional simplifying approximation $\mathcal{F}=\langle\mathcal{F}_{\Omega}\rangle_{\rm Haar}$. Despite the fact that such approximation is in general unjustified, we expect it to lead to reasonable results in the regime $\sigma_{1}^{2}\sim\sigma_{2}^{2}\sim\sigma_{\perp}^{2}$. With these assumptions, the vanishing-current solution of eq. (6.18) is given by

[TABLE]

where now

[TABLE]

Here $\mathcal{N}$ denotes the normalization constant, which must be calculated numerically in general. For negative values of $Q$ the expression (7.22) is to be understood in the analytically continued sense. Note that in the ${N_{\rm{f}}}\rightarrow 2$ limit, the solution reduces to the one in eq. (6.24), as expected. The $w(\theta)$ calculated above can now be used in the right hand side of eq. (7.13) to calculate the typical occupation number.

The numerical solution for the evolution of the occupation number $n$ for fifty different realizations of randomly generated scatterers and ${N_{\rm{f}}}=3,10$ is shown in the left panels of figure 6. Similarly to the two-field case, the analytic result (shown as the orange line) accurately approximates the typical occupation number.

In the right panels of figure 6 we show the three- and ten-field theoretical predictions for different values of the ratios $\sigma_{a}/\sigma_{1}$, and compare them to the numerical results. It must be noted that the vertical axis has been rescaled by ${N_{\rm{f}}}$ compared to the two-field case for a better reading of the results. In both cases, but more evidently in the ten-field case, the deviation between the theoretical prediction and the numerical results is significant in the regime $\sigma_{\perp}^{2}\ll\sigma_{1}^{2}$, $\sigma_{2}^{2}\sim\sigma_{1}^{2}$. We believe that this deviation is related to the breakdown of our assumptions regarding $\mathcal{F}$ and $w(\theta)$; namely, the factorizability of the angular dependence of the probability density and/or the naive assumption of a scattering strength-independence of $\mathcal{F}$.

The MEA ansatz by itself provides a good fit to the data in the $\sigma_{1}\sim\sigma_{2}\sim\sigma_{\perp}$ regime. Here the MEA solution can be written as

[TABLE]

The corresponding curves are shown in fig. 6 (right) as the dashed gray lines.

Unlike the two-field case, for ${N_{\rm{f}}}>2$ fields we lack the knowledge of the full set of coefficients of the FP equation, and thus we cannot claim that, in general, ${\sf\Lambda}$-isotropy implies ${\sf u}$-isotropy, although we can prove that the distribution of the angle relevant for the typical occupation number takes the ${\sf u}$-isotropic form. In the ${\sf\Lambda}$-isotropic limit, the analytical solution in eq. (7.21) evaluates to

[TABLE]

which is nothing but the $U({N_{\rm{f}}})$-flat result w(\theta)\big{|}_{\rm Haar}, see Appendix A.4.

At any rate, for ${N_{\rm{f}}}\gg 1$, based on eqns. (7.17) and (7.25) we expect the full analytical result and the result based on the MEA to be coincident in most of the parameter space for a large number of fields, save for $\sigma_{\perp}^{2}\ll\sigma_{1}^{2}/{N_{\rm{f}}}$ with $\sim{N_{\rm{f}}}^{2}$ contributions from $\sigma_{\perp}^{2}$ dominating the answer. In terms of fig. 6, the separation between the curves, and the magnified region are observed to shrink as $\sim 1/{N_{\rm{f}}}$.

For generic scattering strengths, not necessarily of the form in eq. (7.1), it has been argued that, due to the concentration of measure property of the unitary group, the parameter space region for which the dependence of $\delta f$ deviates significantly from that predicted by MEA, becomes negligibly small as ${N_{\rm{f}}}\rightarrow\infty$ [51]. We therefore expect that the validity of the trace formula (4.10) carries over beyond the assumption (7.1), as long as the disorder strengths do not greatly deviate from ${\sf\Lambda}$-isotropy.

7.3 Summary

We have derived the typical occupation number in the ${N_{\rm{f}}}>2$ case without appealing to the maximum entropy ansatz, but with some restrictions (see eq. (7.1)) on the form of the $\sigma_{ab}^{2}$. We had to make a couple of ad-hoc assumptions along the way: (1) $P(\theta,\alpha_{2}\ldots)=w(\theta)p(\alpha_{2},\ldots)$. (2) $\mathcal{F}=\langle\mathcal{F}_{\Omega}\rangle_{\rm Haar}$. These were motivated by their validity within the MEA based approach. We found that the results then match the numerical simulations exceptionally well, apart from the region $\sigma_{\perp}^{2}\ll\sigma_{1}^{2}/{N_{\rm{f}}}$. The MEA based result (see eq. (7.25)) and the result from this section match quite well apart from a small region $\sigma_{\perp}^{2}\ll\sigma_{1}^{2}/{N_{\rm{f}}}$ which scales as $1/{N_{\rm{f}}}$. Thus in the large ${N_{\rm{f}}}$ limit, it is safe to use the simple MEA based formula (see eq. (4.5)) for the typical occupation number.

8 Summary and Future Outlook

Non-adiabatic particle production during and after inflation can impact the generation of curvature perturbations as well the efficiency of energy transfer during reheating. In this work, we derived a general Fokker-Planck equation for the probability distribution of the total occupation number of ${N_{\rm{f}}}$ coupled fields under the assumption that the typical particle production rate per interaction is small. Our analytic calculation of the typical total occupation number based on the FP equation agreed exceptionally well with direct numerical computations in the limit of large number of interactions.

For analytic tractability, we restricted our attention to cases where the effective stochastic masses and couplings of the fields can be described by Dirac-delta function with amplitudes and locations drawn from different distributions. Furthermore, we ignored effects of expansion and backreaction on the perturbations. Physically this corresponds to non-adiabatic interactions whose temporal width is small compared to the wavelength of the mode-functions of interest, and that we are restricting our attention to sub-horizon, linearized perturbations.

In the present manuscript, we have extended the framework of [34] for calculating the typical occupation number of ${N_{\rm{f}}}$ coupled fields. The key difference from [34] was that we did not rely on the maximum entropy ansatz (MEA). We elucidated the domain of validity of the MEA in terms of the relative strengths of the effective masses and cross couplings of the fields ($\sigma_{ab}^{2}$) and the number of fields (${N_{\rm{f}}}$). We proved that MEA based results coincide with our present results when all $\sigma^{2}_{ab}$ are equal. We verify that a simple and useful trace formula based on MEA agrees well with detailed numerical results when the $\sigma^{2}_{ab}$ are comparable (not necessarily equal), or more generally if ${N_{\rm{f}}}\gg 1$ even when some of the $\sigma^{2}_{ab}$ differ significantly. The deviations from this trace formula are captured by our more detailed framework here, albeit at an increased level of complexity in the derivation.

In detail, we note that the expressions for the typical occupation number in the ${N_{\rm{f}}}=1$ and ${N_{\rm{f}}}=2$ are given in closed form (up to quadratures). To do the same for ${N_{\rm{f}}}>2$ fields, additional simplifying assumptions are needed. These assumptions, however, mattered less when the number of fields is large (${N_{\rm{f}}}\gg 1$) where we find the MEA based trace formula to be a good approximation for a very wide range of effective masses and couplings.

The overall approach advocated here is to avoid relying on detailed model building, and rather take a coarse-grained approach to the particle production in the early universe. We believe that this is reasonable way forward given (i) the complexity of fundamental physics models of the early universe, and (ii) the simplicity of the existing data from the early universe. While the technical approach is different here, the EFT of inflation and reheating [52, 53, 54, 55], as well as disorder during inflation [38], is based on a similar underlying philosophy. There is also similarity in the underlying view taken by the authors of [14, 13, 56], where simple power spectra of the curvature perturbations are shown to exist in spite of the complexity of the underlying models.

A possible means of moving away from the weak scattering limit is to use Random Matrix Theory (RMT), taking advantage of the large number of fields ${N_{\rm{f}}}$ and the large number of scatterings ${N_{\rm s}}$ [45, 46, 57]. In the case of rapid particle production, a detailed investigation of back-reaction and nonlinear mode effects might be needed after a few interactions. Such effects, however, are beyond the domain of applicability of the present framework.

For applications of this framework to inflation, the impact on the curvature fluctuations as well as backreaction of curvature fluctuations on particle production is needed. The results from [58, 38, 30, 56] should provide useful guidance for such calculations. For reheating applications, the impact of different expansion histories on particle production and vice versa is needed. We hope that the technical results presented here and in [34] will be useful as starting points for calculations related to observables.111111In the present work we have focused on determining the total typical occupation number which is related to the sum of the eigenvalues of the transfer matrices. This is a “basis independent” quantity, and was chosen for reasons of calculational convenience as well as usefulness in potential applications (for example, reheating studies on sub-horizon scales). In some applications, especially those related to superhorizon scales, such an occupation number need not have a useful physical meaning. Moreover, even on subhorizon scales, the quantities of interest may correspond to the most probable branching ratio $(n_{a}/n_{b})_{\rm typ}$ of the individual occupation numbers $n_{a}$ for different fields or to other combinations of terms in the transfer matrices. Such quantities can in principle be evaluated, using this present work as a guide. We leave these calculations for future work.

Finally, we note that while we were motivated by the desire to tackle particle production in the early universe, the new results presented here can be re-used for the problem of current conduction in disordered wires (where the original motivation for the statistical framework came from). More generally, within the limits of the assumptions stated, our framework and results might be useful for analyzing systems described by a set of ${N_{\rm{f}}}$ coupled oscillators with ${N_{\rm s}}\gg 1$ localized and stochastic variations in their frequencies and couplings.

Acknowledgements

We would like to thank Anthony Chan, Matthew Foster, Daniel Baumann, Daniel Green, Mafalda Dias and Jonathan Frazer for helpful conversations. M.A.A. and M.A.G.G. acknowledge support from the US Department of Energy Grant DE-SC0018216. This work was completed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293.

Appendix A Appendix

A.1 The Fokker-Planck Equation

In this section we present a derivation of the Fokker-Planck equation (3.12), where the coefficient correlators are evaluated as functions of the matrices ${\sf R}(j)$ and ${\sf R}_{j+1}$. For this purpose we rewrite the Smoluchowski equation in the form

[TABLE]

or, equivalently, in terms of the $2{N_{\rm{f}}}^{2}+{N_{\rm{f}}}$ parameters ${\boldsymbol{\lambda}}$ that characterize each transfer matrix,

[TABLE]

where the measure $d{\boldsymbol{\lambda}}$ is suitably defined, and where $\tilde{{\boldsymbol{\lambda}}}$ denotes the functional dependence of the parameters of a product of transfer matrix, i.e. ${\sf M}(\tilde{{\boldsymbol{\lambda}}})={\sf M}_{1}({\boldsymbol{\lambda}}_{1}){\sf M}_{2}({\boldsymbol{\lambda}}_{2})$. Upon Fourier-transforming the previous equation, and introducing the change of variables $\tilde{{\boldsymbol{\lambda}}}={\boldsymbol{\lambda}}_{1}+\delta{\boldsymbol{\lambda}}$, we can write

[TABLE]

Using the inverse Fourier transform,

[TABLE]

where

[TABLE]

Upon division by $\delta\tau$, and integration with respect to the set of parameters that characterize ${\sf v}$ (which form a compact manifold, c.f. (3.9)), the resulting equation corresponds to (3.12), where abusing notation we have also denoted by $P$ the probability density marginalized with respect to ${\sf v}$.

A.2 Coefficients of the ${N_{\rm{f}}}=2$ FP Equation

The list of increments $\delta\lambda_{a}$, needed to calculate the coefficients appearing in the ${N_{\rm{f}}}=2$ FP equation is provided in this appendix. Using eqns. (3.22a)-(3.22d)), we have

[TABLE]

where $\Delta f=f_{1}-f_{2}$. At next to leading order, we have

[TABLE]

The correlators $\langle\delta\lambda_{a}\rangle$ and $\langle\delta\lambda_{a}\delta\lambda_{b}\rangle$ that appear in the FP equation are listed next. The dependence of coefficients on the polar angle $\theta$ and the scatterings strengths can (almost) be encoded via the $\gamma,l$-functions defined below,

[TABLE]

and

[TABLE]

The correlations $\langle\delta\lambda_{a}\delta\lambda_{b}\rangle$ are given by

[TABLE]

while the expectation values $\langle\delta\lambda_{a}\rangle$ correspond to

[TABLE]

A.3 Numerical Method

The transfer matrix approach described in section 3 can be translated in a straightforward way to a simple numerical method to compute the occupation number after a certain number of scatterings; the occupation numbers $n_{a}(j)$ are obtained by keeping track of the iterated multiplication of matrices and evaluating its trace at each step. However, as it is shown in figure 7, in order to reach the stationary late-time limit for the particle production rate, for which our analytical results are valid, an enormous amount of calculations need to be considered as ${N_{\rm{f}}}$ increases. Since such a computation using the naive approach would almost inevitably lead to numerical overflow, a better algorithm is needed.

Recall that the solution for the Fourier mode of the field $\chi^{a}$ after $j$ scattering events is given by (3.1), where the amplitudes of the exponential terms for all fields can be arranged in the coefficient vector

[TABLE]

Note first that

[TABLE]

and $|\boldsymbol{\Psi}(0)|^{2}={N_{\rm{f}}}$. Therefore, if we split the chain of ${N_{\rm s}}$ scatterings into $m$ blocks containing each $N_{m}$ interactions, ${N_{\rm s}}=m\times N_{m}$, the typical particle number rate in the ${N_{\rm s}}\gg 1$ limit can be approximated as

[TABLE]

i.e. we compartmentalize the matrix multiplication; the argument of the logarithm is identified with the amount that the arbitrary unit vector is scaled during each block. At the start of each $j^{\rm th}$ block we re-normalize the outcome of the previous block, $|\boldsymbol{\Psi}((j-1)N_{m})|^{2}=1$, and apply the transfer matrix $N_{m}$ times to obtain $|\boldsymbol{\Psi}(jN_{m})|^{2}$; the procedure is then repeated for the $(j+1)^{\rm th}$ block and the results are added. Upon finding the particle production rate for a particular realization of the scattering locations and strengths, the procedure is repeated for several realizations, and the expectation value returns the desired mean particle production rate.

A.4 Haar average of $\mathcal{F}_{\Omega}$

In the notation of section 7.1, the Haar measure of the group $SU(N)$ can be explicitly written as [50]

[TABLE]

where the kernel of the measure is given by

[TABLE]

with

[TABLE]

and

[TABLE]

For the subset of angles that parametrize $\mathcal{F}_{\Omega}$, one can immediately check that

[TABLE]

Therefore, the expectation value $\langle\mathcal{F}_{\Omega}\rangle_{\rm Haar}$ can be calculated inductively using the defining expression (7.8), as follows:

[TABLE]

As for the last step we have $\alpha_{2(N-2)}$, we obtain $\langle\mathcal{F}_{\Omega}\rangle_{\rm Haar}=2/N$.