Contents

Multifield Stochastic Particle Production: Beyond a Maximum Entropy Ansatz

Mustafa A. Amin111mustafa.a.amin@gmail.com, Marcos A. G. Garcia222marcos.garcia@rice.edu, Hong-Yi Xie and Osmond Wen

Department of Physics & Astronomy, Rice University, Houston, Texas 77005-1827, U.S.A.

Abstract
We explore non-adiabatic particle production for 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 and cases. We also present results for the general 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 .

## 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 , when the number of fields and . For , 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 . In general, for simple results, it is essential to have as well as or .

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 () case as a warm-up example. Moving beyond MEA, in section 6, we consider the 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 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 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, . The time variable will be , and overdots will denote derivatives with respect to . We will use a bold sans-serif font for matrices, e.g. . The indices will be used for field indices. We use the metric convention.

## 2 The Model

Consider coupled scalar fields, , with masses and time-dependent stochastic masses and couplings . Here . The quadratic action for these fields is taken to be

 S(2)=∫d4x L=∫d4xNf∑a,b=1(12δab∂μχa∂χb−12Mab(τ)χaχb), (2.1) Mab(τ)=m2aδab+msab(τ).

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 are

 (d2dτ2−∇2+m2a)χa(τ,→x)+Nf∑b=1msab(τ)χb(τ,→x)=0. (2.2)

The equations of motion of the Fourier modes of the field are

 (d2dτ2+ω2a)χa(τ,→k)+Nf∑b=1msab(τ)χb(τ,→k)=0, (2.3) ω2a(k)=k2+m2a,

where we have used the symbol 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 on the wavenumber will be suppressed.

The time dependence of can be quite complex. For example, 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 after the -th non-adiabatic event is

 χaj(τ)≡1√2ωa[βajeiωaτ+αaje−iωaτ], (3.1)

where the overall normalization is chosen for future convenience. The Wronskian of the solutions and is a constant, , which in terms of the coefficient in eq. (3.1) implies .

The coefficients and before and after the -th scattering, are related by the transfer matrix .

 (→βj→αj)=Mj(→βj−1→αj−1). (3.2)

For a single scattering we can always write where contains the necessary information about the scattering (for weak scattering ). We can connect the coefficients before and after scatterings by simply chaining together the transfer matrices

 (→βj→αj)=M(j)(→β0→α0),whereM(j)=MjMj−1…M1. (3.3)

Note that denotes the transfer matrix of scatterings, whereas is the transfer matrix for the -th scattering.

### 3.1 The Occupation Number

We define the occupation number of each field after scatterings by333We note that this definition agrees with the usual one when thinking about 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].

 na(j)=12ωa(|˙χaj|2+ω2a|χaj|2)−12=|βaj|2, (3.4)

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

 n(j)≡Nf∑a=1na(j). (3.5)

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

 n(j)=14Tr[M(j)M†(j)−1]. (3.6)

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

 Rj≡MjM†j,%andcorrespondinglyR(j)≡M(j)M†(j). (3.7)

For the total occupation number, we only care about and not . This simplifies matters considerably. Note that is Hermitian. The total occupation number is related to the sum of the eigenvalues of .

### 3.2 Parametrizing R

How many parameters are needed to describe ? A general parametrization of the matrix is

 R=(u00u∗)U(f~f~ff)F(u†00uT)U†=UFU†, (3.8)

where . In the above expression is a diagonal matrix consisting of eigenvalues and is a unitary matrix parametrized by angular variables. Note that and are matrices with blocks constructed out of and . Thus, is parametrized by total variables. This should be contrasted with , for which a general parametrization has the form [44]

 M=(u00u∗)(√1+n√n√n√1+n)(v00v∗), (3.9)

where , and where is a unitary matrix parametrized by angular variables, for a total of variables characterizing . The parameters of will be denoted by where are the diagonal entries of and are the angular variables for the unitary matrix . In this parametrization, note that

 n(j)=14Tr[R−1]=12Tr[f(j)−1]=12Nf∑a=1(fa−1). (3.10)

### 3.3 The Fokker Planck Equation

Our next goal is to figure out how 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 . This can then be used to get the expectation value of the occupation number (or any other function depending on variables in ).

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

 Pτ+δτ(M)=∫dμ(M2)Pτ(M−12M)Pδτ(M2), (3.11)

where denotes the invariant (Haar) measure with respect to the transfer matrix group, with the property that with 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 -matrix parameters , evolves with time according to the (forward) Kramers-Moyal expansion

 ∂τP =−Nf2+Nf∑b=1∂λb[⟨δλb⟩δτδτP]+12!Nf2+Nf∑b,c=1∂λb∂λc[⟨δλbδλc⟩δτδτP]+… (3.12)

The in the above expression should be thought of as the “small” increment of the parameters (relative to their values before the addition of the interval ) due to a single additional scatterer in the time interval (weak scattering). The expectation value is over the probability distribution describing the properties of the scatterer in the interval . 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 .444In the single field limit it can be explicitly shown that , , and for , where denotes the local mean particle production rate [34].

It is often not necessary to get an explicit solution for . Instead we can integrate eq. (3.12) to obtain equations for the expectation values of quantities of interest. Let us consider a function where is the total occupation number of the fields. Note that only depends on the eigenvalues and not on the angular variables. Multiplying eq. (3.12) with and integrating over all variables we arrive at

 (3.13)

where . In deriving the above expression we have repeatedly used integration by parts, dropped the boundary terms, taken advantage of the fact that the variables of form a compact manifold, and that is a function of the eigenvalues only. We have also used the expression for in terms of the eigenvalues . Note in particular that the sums are over variables only, however the integration is over all variables. Moreover the coefficients as well as might be non-factorizable functions of all variables.

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

 ntyp≡exp[⟨ln(1+n)⟩]−1. (3.14)

The quantity 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 the ratio grows with time (see the discussion in [34]).

With in mind, consider the function . Due to the relation (3.10), eq. (3.13) can be then be written as

 ∂τ⟨ln(1+n)⟩=⟨12(1+n)Nf∑a=1⟨δfa⟩δτδτ−18(1+n)2Nf∑a,b=1⟨δfaδfb⟩δτδτ⟩. (3.15)

To make further progress, let us first turn to more detailed understanding of and other .

### 3.4 Increments in R

Our next task is to find expressions for . To this end, let us consider the change in due to a single additional scattering (see fig. 2)

 R(j+1)=R(j)+δR, (3.16) δR=R(j)m†j+1+h.c.+mj+1R(j)m†j+1.

We will often suppress as an argument of the matrices.

The increment in , , can be expressed in terms of increments in the eigenvalue matrix and the angular matrix (up to second order in the increments) as follows:

 δR =U(δF+[U†δU,F]−[U†δU,F]U†δU+[U†δU,δF])U†. (3.17)

Comparing the expressions for in eqns. (3.16) and (3.17), we get

 δF+[U†δU,F]−[U†δU,F]U†δU+[U†δU,δF]=G, (3.18)

where

 G≡(QF+h.c.)+QFQ†,withQ≡U†mj+1U. (3.19)

Equation (3.18) expresses the increments and in terms of which contains information about the additional scattering via , as well as and . As is evident from the above equation the increments and will each contain - and -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

 δF≡δF(1)+δF(2), δU≡δU(1)+δU(2),G(1)≡(QF+h.c),G(2)≡QFQ†, (3.20)

which then yields

 δF(1)+[U†δU(1),F]=G(1), (3.21) δF(2)+[U†δU(2),F]−[U†δU(1),F]U†δU(1)+[U†δU(1),δF(1)]=G(2).

Recall that and are block matrices, which are constructed out of blocks of functions of and (and ) matrices. The above expressions for and can be unpacked in terms of and directly to yield

 δf(1)+[u†δu(1),f]=g(1), (3.22a) δf(2)+[u†δu(2),f]−[u†δu(1),f]u†δu(1)+[u†δu(1),δf(1)]=g(2), (3.22b) δ~f(1)+u†δu(1)~f−~fuTδu∗(1)=~g(1), (3.22c) δ~f(2)+u†δu(2)~f−~fuTδu∗(2)−(u†δu(1)~f−~fuTδu∗(1))uTδu∗(1)+u†δu(1)δ~f(1)−δ~f(1)uTδu∗(1)=~g(2), (3.22d)

where is the top left block of , and is the top right block of . Although the previous system represents equations for the 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 equations for the components of and , and similarly for the second-order terms.

At the level of the variables that parametrize and , the above equations can be solved for each and in terms of the parameters characterizing the scatterers (elements of ) as well as . In component form the solution can be written explicitly. The first order correction to can be directly obtained from the diagonal elements of eq. (3.22a),

 δf(1)a=g(1)aa. (3.23)

For the second-order change in , eq. (3.22a) yields for the off-diagonal components

 (u†δu(1))ab=g(1)abfb−fa,(a≠b). (3.24)

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

 δf(2)a=g(2)aa+∑c≠ag(1)acg(1)cafa−fc. (3.25)

The first-order corrections for the entries of 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

 (u†δu(1))aa=~g(1)aa2~fa−fa2~f2aδfa. (3.26)

These two equations form a linear system for the perturbations of , which can be solved explicitly as

 δu(1)ab=∑c≠buacg(1)cbfb−fc+uab2~f2b(~fb~g(1)bb−fbg(1)bb). (3.27)

Analogously, one can solve for the second-order corrections for components, which are given by

 δu(2)ab =∑c≠buacfb−fc[g(2)cb+g(1)cd(u†δu(1))db−g(1)bb(u†δu(1))cb] +uab2~fb[~g(2)bb−δ~f(2)b−~g(1)bd(u†δu(1))∗db−δ~f(1)b(u†δu(1))bb−~fb(δu(1)†δu(1))bb]. (3.28)

To make further progress, we need an explicit form of . Below we will choose Dirac-delta scatterers. Once the form of is specified, we can solve, in principle, for the increments in the parameters of for arbitrary . In practice this is a non-trivial calculation made particularly onerous by the large number of parameters () in and () in .

### 3.5 Dirac Delta Scatterers

The general formalism we have presented is independent of the precise form of the scatterers (), 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 , we approximate as follows:

 msab(τ)=2√ωaωbNs∑j=1Λab(τj)δ(τ−τj), (3.29)

where is the number of events, are uniformly distributed, and are Dirac Delta functions. For each , the elements of which characterize the strength of the scatterers are drawn from some distribution. We assume that the distributions are identical and

independent for all . Note that the elements of must be symmetric with respect to and . The scaling with is for future convenience. We will assume that for each

With the parametrization (3.29), the scattering strength variances are functions of the wavenumber and the corresponding scalar field masses (see fig. 3).

By imposing appropriate junction conditions on at , the transfer matrix at for the Dirac delta function scatterers evaluates to

 (3.31)

where

 pj≡⎛⎜ ⎜ ⎜ ⎜ ⎜⎝eiω1τj0…00eiω2τj…0⋮⋮⋱⋮00…eiωNfτj⎞⎟ ⎟ ⎟ ⎟ ⎟⎠,and[Λj]ab=Λab(τj). (3.32)

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

 g(1) ≡[Γ,f]+[Σ,~f]+~f(Σ+Σ∗), (3.33) g(2) ≡−ΓfΓ+Γ~fΣ∗−Σ~fΓ+ΣfΣ∗, ~g(1) ≡Γ~f+Σf+fΣ−~fΓ∗, ~g(2) ≡ΓfΣ−Γ~fΓ∗+Σ~fΣ−ΣfΓ∗,

where we have defined for the -th scatterer

 Σ≡iu†p†j+1Λj+1p∗j+1u∗,andΓ≡iu†p†j+1Λj+1pj+1u. (3.34)

Note that is evaluated at rather than .

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 , we can rewrite eq. (3.23) in terms of the matrix given in eq. (3.33). The following explicit expression for the perturbation is obtained:

 δf(1)a =i~fa∑b,c,d,e(u†abp†bcΛcdp∗deu∗ea−uTabpTbcΛcdpdeuea) =i~fa∑b,c(p∗bp∗cu†abΛbcu†ac−pbpcubaΛbcuca). (3.35)

In the second line we have made use of the fact that is diagonal. Consider now the expectation value . This product contains phase factors due to the presence of the components of , which can be explicitly evaluated, as each is assumed to be uniformly distributed in the interval . For example,

 ⟨papbp∗cp∗d⟩δτδτ=⟨eiΔωτj⟩δτδτ=sin(Δωδτ/2)Δωδτ/2eiΔω(j−1/2)δτ. (3.36)

where . 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, ; this is consistent with the treatment of the dynamics as free in-between scatterings. This assumption implies that, within the product , we can disregard terms that contain unpaired and , as they will vanish after averaging over scatterers; additionally upon averaging. Recalling the statistical properties of in eq. (3.30), we can then write555In eq. (3.37) (and in what follows) we have denoted . Whenever there is no potential for confusion we will also suppress the subindex, to avoid cumbersome expressions.

The second-order change in 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,

 ⟨δf(2)a⟩δτ=⟨g(2)aa+∑b≠ag(1)abg(1)bafa−fb⟩δτ, (3.38)

where

 ⟨g(2)aa⟩δτ=∑b,c,dfdσ2bc[(1+δbc)u†abubdu†dcuca+2u†abubau†dcucd], (3.39)

and

 ⟨∑b≠ag(1)abg(1)bafa−fb⟩δτ =∑b,c,db≠aσ2cd[(fa−fb)(u†acucau†bdudb+δcdu†acudbu†bcuda) (3.40) +~f2b+~f2afa−fb(u†acucau†bdudb+u†acucbu†bduda)].

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

This is an explicit expression for the evolution of the typical occupation number for arbitrary and (assuming Dirac-delta function scatterers). However, note that the angular brackets denote integration with respect to probability distribution . In the absence of an explicit parametrization for 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 -flat probability distribution and for a statistically isotropic (i.e. ). In the subsequent sections, we will evaluate this expression in general for the cases, and for under some restrictive conditions for the scattering strengths .

## 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 with respect to the parameters in . 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 associated with the MEA corresponds to that which maximizes the (Shannon) entropy functional [42],

 S[P]=−⟨lnP⟩δτ, (4.1)

subject to the constraints that the local mean particle production rate is known, and that the evolution of the transfer matrix is under perturbative control, as . With this definition, the MEA recipe is said to provide the least biased estimate of the consistent with the (fixed) local production rate. Under these approximations, it can be shown that the probability density is independent of , . Moreover, for the typical occupation number , the MEA implies the late-time result [34]

 ∂τ⟨ln(1+n)⟩MEA⟶2NfNf+1μj, (4.2)

where denotes the instantaneous mean particle production rate,

 μj≡1Nf⟨nj⟩δτδτ. (4.3)

Here