Quantum Out-of-Equilibrium Cosmology

In this work, our prime focus is to study the one to one correspondence between the conduction phenomena in electrical wires with impurity and the scattering events responsible for particle production during stochastic inflation and reheating implemented under a closed quantum mechanical system in early universe cosmology. In this connection, we also present a derivation of fourth order corrected version of the Fokker Planck equation and its analytical solution for studying the dynamical features of the particle creation events in the stochastic inflation and reheating stage of the universe. It is explicitly shown from our computation that quantum corrected Fokker Planck equation describe the particle creation phenomena better for Dirac delta type of scatterer. In this connection, we additionally discuss It$\hat{o}$, Stratonovich prescription and the explicit role of finite temperature effective potential for solving the probability distribution profile. Furthermore, we extend our discussion to describe the quantum description of randomness involved in the dynamics. We also present a computation to derive the expression for the measure of the stochastic non-linearity arising in the stochastic inflation and reheating epoch of the universe, often described by Lyapunov Exponent. Apart from that, we quantify the quantum chaos arising in a closed system by a more strong measure, commonly known as Spectral Form Factor using the principles of Random Matrix Theory (RMT). Additionally, we discuss the role of out of time order correlation (OTOC) function to describe quantum chaos in the present non-equilibrium field theoretic setup. Finally, for completeness, we also provide a bound on the measure of quantum chaos arising due to the presence of stochastic non-linear dynamical interactions into the closed quantum system of the early universe in a completely model-independent way.


Introduction
Quantum fields in an inflationary background  or during reheating [26][27][28][29][30][31][32] gives rise to the burst of particle production, which has been extensively studied in ref. [33][34][35]. This has been studied to a great extent in the background of the inflationary scenario of the universe in ref. [36][37][38]. Such phenomena has been compared to that of the scattering problem in quantum mechanics with a specific effective potential arising due to the impurity in the conduction wire, which can approximately be solved using the well known WKB technique [34,36] 1 . It is important to note that such particle production events are completely random (or chaotic) when the evolution is non-adiabatic or tachyonic in nature.

Slow-rolling inflaton
Inflation Figure 1. This schematic diagram shows the correspondence between the conduction phenomena in electrical wires with impurity to that of the cosmological random particle-creation events during the non-adiabatic stage of the early universe.
A non-adiabatic change in the time dependent effective mass profiles of the fields (which is actually coming from integrating out the heavy degrees of freedom from the UV complete theory and after path integration finally one gets the time dependent effective coupling parameters between fields) as the background evolution of the fields passes through special points in field space produces these burst of particle creation in (quasi) de Sitter space time. 1 In the context of cosmology conformal time dependent effective mass profile exactly mimics the role of impurity potential in electrical conduction wire. Due to such one to one correspondence the time evolution equation (i.e. Klien Gordon equation) of the Fourier modes corresponding to the quantum fluctuation in the context of primordial cosmology can be described in terms of the Schrodinger equation in electrical conduction wire with specific impurity potential. We have investigated this possibility in detail in this paper. Additionally, it is important to note that such time dependent effective mass profiles are also important to study the role of quantum critical quench and eigen state thermalization [] during the reheating epoch of universe.
There lies a physical and mathematical equivalence between such cosmological events to that of the stochastic random phenomena occurring in mesoscopic systems where fluctuations in physical quantities play a significant role of producing stochastic randomness in the system under consideration. We also discuss the cosmological systems which have been considered to be rather non-linear and dissipative due to the significant amounts of quantum fluctuations in the effective coupling terms (or in the time dependent effective mass profile) of the interactions between the fields. Important reviews on the non-linear and dissipative effects arising in the context of cosmology were put forward in the refs. [9,[39][40][41]. In this paper we explicitly discuss bout the various non-linear and dissipative effects in cosmological set up that arises in (quasi) de-Sitter space with m 2 > 0, where the term m 2 represents the effective mass squared of the created particle in (quasi) de-Sitter background. In this connection it is important to note that, the massless scalar field gets "thermalize" due to the effective time dependent interaction in the (quasi) de-Sitter background. The cosmological events that we talk about in this paper are identified with those of the particle production stochastic random events. In this paper, we present the dynamical features of inherent chaos (stochastic randomness) in the physical system and its connection with the quantum mechanics in detail. The model is exactly similar to that of "massless scalar field" interacting with a scatterer in the background which are treated to be the heavy fields and are mainly responsible for cosmological particle production in (qusi) de-Sitter space (see [36]). In this context, when the free massless scalar field interacts with the the heavy field in the background space time, it mimics the role of thermalization phenomena of the field which occurs during the epoch of reheating of the universe. The specific problem we will discuss here is similar to one presented in ref. [42]. This problem is similar to that of a scattering problem in presence of impurity in quantum me-chanics where the Schrödinger equation yields approximate solutions to the wave-function of the particle which encounters a effective impurity potential barrier V (x) of a given strength. The similarity in the following model is drawn between the current carrying electrons responsible for conduction in electrical wires to that of the particle creation in cosmology as a result of the non-adiabatic random events occurring in the early (inflation and reheating) stage of the universe. In this present problem for the sake of simplicity we consider an one dimensional conducting electrical wire, which implies that the current carrying electrons in the electrical wire has only a single propagating degree of freedom. As mentioned earlier, this has been considered to reduce clutter in our computation. But the similar problem can be generalized to more complicated situation 2 . Since, a current carrying wire consists of a large number of impurities, these act like the potential barriers V (x), which are randomly distributed across the wire. Therefore, the motion of the electrons while confronting these scatterers gets hindered due to the presence of these randomly placed scatterers. One of the most important outcome of such an event is known as Anderson Localization as appearing in the context of condensed matter systems. Usually this is characterized by probability density of the localized wave-function: with ξ being the localization length of the quantum mechanical wave-function ψ(x). This phenomena of Anderson Localization usually occurs due to the interference of the waves scattered from the impurities present in the conduction wire. By formulating cosmological particle production as a random scattering problem, it has been shown in [42] that Anderson localization maps to a problem of estimating exponential particle production, as given by: where µ k is the mean particle production rate which is characterized by the conformal time dependent scalar field φ k (τ ). A striking similarity has been observed between such scattering problems in conducting wires to that of the burst of particle production in cosmological random events shown in ref. [42]. In such cases, it has been observed that the solving a scattering problem in quantum mechanics using Schrödinger equation is similar to solving a Klein-Gordon equation for a massless scalar field in presence of a conformal time-dependent effective mass squared coupling parameter m 2 (τ ). In this context the scalar field with timedependent mass m 2 (τ ) mimics the role of coupling strength parameter which characterizes the scattering to the massless scalar field in (quasi) de Sitter background. For more details see refs. [43][44][45]. Moreover, such stochasticity in a cosmological set up arises due to the stochastic time evolution of Hubble parameter H(t), so that the inflaton (or the field participating in reheating) evolves with time stochastically due to the quantum fluctuations in the FLRW background. In the similar context the role of interacting scalar field has been studied to a great deal in ref. [43].
In this context,we have presented the the amount by which the quantum mechanical system deviates with respect to the initial conditions. This means that more the value of this exponent, more is the chaos or stochastic randomness in the system under consideration in this paper. This exponent plays a significant role in our scenario as the number of particles produces in a given scattering event per unit time is random in nature. In a system of randomly spaced scatterers chaos emerges out of the random scattering events that an electron encounters while drifting across the wire with some drift velocity v within the conducting wire. The number that quantifies this increase in stochastic randomness or chaos in the system is the Lyapunov Exponent. In refs. [42,43], particle production phenomena in cosmological non-adiabatic events has been exclusively studied which yields the fact that the particle occupation number depends on Floquet indices µ k , which finally control the number of produced particles with the following number density: as well as the variances in the field fluctuation. The quantum fluctuations in the inflationary state of the universe results in the randomization of these bursts of particle production. The number density has been a random variable which is rendered stochastic due to the scattering events in the context of early universe cosmology. Our main objective in this paper to quantify this characteristic number for the massless scalar field having a conformal time-dependent mass coupling with it. One of the prime reasons for finding a signature of chaos in such a system is the well known thermalization phenomena, which means that the FLRW background which embeds the massless scalar field into it is being thermalized by the massive field in interaction with the FLRW set up, which constantly being giving rise to a burst of particle production in the context of early universe cosmology. The scalar fields that we considering in our paper are said to be massive or heavy fields (m ≥ H) which mimics the role of the scatterers in the Schrödinger problem in quantum mechanics where the strength of the effective potential or the scatterer is given by the probability distribution function of the effective potential function. We draw a picturesque landscape by considering three distinct mass profiles: which exactly mimics the role of cosmological scatterers in early universe. We thereby investigate the momentum scale dependent behaviour of the Lyapunov exponent. In this context, the incoming momenta of the mode functions of the quantized massless scalar field having random interactions with the scatterer. In the following class of model, the Bogoliubov coefficients arise due to the interaction between massless scalar field with the heavy field. These Bogoliubov coefficients gives the information about the transmission coefficient viz.a.viz in similar problem to that of a scattering problem in quantum mechanics, that we solve using the well known WKB approximation technique 3 . These WKB solutions are ex-tremely useful as it tell us the dynamical feature of the particle production in cosmological scattering events.
In continuation with this, we discuss about the epoch of reheating which occurs after the end of inflationary stage of the universe which finally results in the stochastic random burst of particle production. The dynamics of these stochastic random bursts of particle production can be well understood by using a Fokker-Planck equation, which gives us a statistical interpretation of the number density of particles created per scattering event.
Since, the number of particles created in a given non-adiabatic event is not discrete in nature but rather its random, which means that there must be a probability distribution function associated with the particle number. The various dynamical features of this type of probability distribution and its physical consequences has been studied in ref. [42]. It has been phenomenologically proposed in ref. [42] that such probability density function would necessarily is Gaussian one. The occupation number of the produced particles, n k , executes a drifting Brownian motion and a Fokker-Planck (FP) equation that evolves the probability distribution, P (n k ; τ ), emerging out of this Brownian motion has been studied in ref. [42]. We further compute the analytical expressions for the mean, variance and other higher order moments which are commonly known as, skewness and kurtosis and such additional statistical higher order moments are very useful to study the exact mathematical form and asymptotic limits of the probability distribution function. The evolution of mean, variance, skewness and kurtosis finally gives a coarse-grained analysis of the Fokker-Planck dynamics to more corrected orders of magnitude in quantum regime. We show in this paper explicitly that though Gaussianity is an inherent part of the probability density function, but the consideration of the higher order moments in the Fokker-Planck equation tells us that the density function may not be a Gaussian one but with some higher-order corrections entailed into it due to the quantum mechanical origin. Therefore, to a greater extent we extend the more corrected quantum version of the Fokker-Planck equation used to describe the dynamics of the probability distribution function used in ref. [42] that tells us the dynamics of the bursts of particle production in these random scattering events. The more quantum corrected version tells us that the probability amplitude of the particle production in the scattering events is more than a Log normal distribution. The distribution profile of the probability distribution function depends largely on the profile of the scatterer, i.e., the effective potential V (x) in the Schrödinger-like equation. While calculating the Fokker-Planck dynamics we observe that the skewness gives us a clue about the rate at which the particle production occurs meaning that longer the trailing part of the profile more is the number density of particles in the scattering event for a given time in the frame of the observer, whereas, kurtosis tells us the width of the probability distribution function which is essentially the amplitude with which the particle production phenomena occurs, which more suggestively tells us about the standard deviation of the density function from Gaussianity. This may be a signature of non-Gaussianity that arises in various models in early universe cosmology.
In this connection it is important to note that, such stochastic approaches to the early universe scenario have been studied in details in [48,49], where the authors give an account of how chaos arises in the context of eternal inflation. As any rapidly oscillating classical field looses its energy by creating pairs of elementary particles, these particles interact with each other and comes to a state of showing thermal behaviour at some temperature T . This implies that we must eliminate the necessary assumption of the universe being in thermal equilibrium. This means that the inflating universe is rather thermal in the sense that the particle creation events that occurs during the quantum fluctuation in the randomly distributed scalar fields φ which results in a chaotic model of the inflationary scenario of the universe thereby leading to a generation of stochastic idea of the particle creation events during the thermalization of the quantum states of the field randomly distributed over the space-time. These particle creation events are more phenomenologically associated with one of the fundamental ideas in out-of-equilibrium statistical mechanics known as Fokker-Planck equation which gives the rate of the particle production during theses random events in stochastically emerging space-time along with the distribution function that this rate charts out. In ref. [42], such a phenomenon of particle creation events by the randomly spaced scatterers in due context of cosmology has been shown where the the statistics of the produced particles as a function of time which is the probability distribution function P (n k , τ ) has been predicted to be following a Log-Normal distribution. The entire process have been carried out with the delta-scatterers which are localized in space-time.
Following ref. [42], in this paper we give a more improved quantum corrected version of the same approach to the probability density function of the particle production events and our prediction from the results show that the higher order quantum correction terms being included into the Fokker-planck equation introduces an approximation to the theory. This tells us that the number of particles produced in a given non-adiabatic event during the reheating stage of the universe is quantized, which would mean that the rate of particle production in a given event gives rise to a discrete set of occupation number n k . Furthermore, the quantum corrected terms obtained by deriving the Fokker-Planck equation takes the general form, which is linear in n k being the first order in τ . Using this information we calculate further the leading order, second and third order terms in the Fokker-Planck equation. Hence, we derive the analytic expression of the quantum corrected version of Fokker-Planck equation. We also calculate the various higher moments in order to get an overview of the nature of the solution of the quantum corrected Fokker-Planck equation which are -standard deviation, skewness and kurtosis which gives the hint of how the probability density function deviates from its Gaussian nature when the higher order quantum corrections are taken into account in the computation. This in turn may will be another indirect signature of the primordial non-Gaussianity in cosmology other than obtaining the signatures provided by the 3-point functions from scalar fluctuations.
Apart from that, we discuss about spectral form factor (SFF), which measures the random distribution of eigen values of the energy hamiltonian of a chaotic system. For this computation of SFF we use the principles of random matrix theory (RMT) in this paper. In the present context an upper bound on SFF denotes the saturation of eigen value distribution hence supports the ref. [50] for quantum chaotic system. Within the framework of quantum physics, chaotic systems can be characterized using only some additional constraints. This theoretical approach is discussed in refs. [51][52][53][54][55] and the authors use the theory of random matrices to characterize quantum mechanical system. In this method, any arbitrarily complicated many-body Hamiltonian can be replaced by matrix of random numbers drawn from a Gaussian statistical ensemble. This random matrix approach towards quantum mechanics help to characterize and understand the underlying features of the chaotic random system. After studying the behaviour of SFF with time one can further comment that whether it is valid for a cosmological particle production event (semi-classical) or not. For our purpose we discuss generalized version of SFF for different even order polynomial structure of random potential and then extend that result to describe the cosmological particle production events [56,57]. For any random potential we can use this method of SFF and we can deal with scatterer of any arbitrary type. For any such scatterer we can get a bound on randomness in the chaotic system characterised by SFF. Also using specific transfer matrix for different conformal time dependent effective mass profiles which are precisely known in this paper, we can finally compute Lyapunov exponent which also measure stochastic randomness.
Also it is important to note that in ref. [42], the scatterers were considered to be some localized potential functions in space-time. On the contrary the choice of our specific time dependent mass profiles mimics the role of thermalized fields or effective potential functions, which are playing the role of scatterers in this context. We see that the choice of these time dependent mass profiles leads to particle production which is chaotic in nature and therefore, to determine the rise of chaos in such a system we quantify as well as analyse chaos by a well known quantities known as the , Lyapunov exponent [58] and Spectral Form Factor (SFF) [59]. Here fusing the principles of random matrix theory (RMT) we provide a generalized bound on randomness (or stochasticity) for any general random scatterer whose potential can be expressed in terms of an even polynomial. More precisely, we provide a possible method to compute the degree of randomness in a chaotic system and from that one can check the bound on chaos.
The plan of the paper is as follows -In section 2 we discuss about the model which is responsible for the quantum description of chaos during the cosmological particle production and have similarities with the quantum mechanical problem of electrical conducting wire with impurities. In section 3, we have presented the analytical expressions for the Bogoliubov coefficients, transmisson and reflection coefficients, Lyapunov exponent, conductance, and resistance for different time dependent mass profile. We have discussed the correspondence between In section 5 the specific role of Spectral Form Factor (SFF) to quantify chaos in the context of particle production rate is discussed. In section 6 the particle production event with quantum corrected Fokker-Planck equation is discussed by taking contribution upto fourth order and also different higher order moments from the quantum corrected probability density function are explicitly computed. Finally, in section 7 we conclude with the future prospect and physical impacts of our work.
Additionally it is important to note that, throughout this paper, we use natural system of units, = c = 1.

Modelling randomness in cosmology
The background model which we consider in this section to quantify quantum chaos in cosmology consists of a massless scalar field interacting with coupled with a background scalar field with conformal time dependent mass profile which in principle have heavier or comparable to the Hubble scale (m ≥ H) [46,60,61]. It is important to note that such heavy mass profiles play significant role in finding various cosmological correlation functions and also can be treated as an additional probe to break the degeneracy between various models of inflation from the perspective of implementing cosmological perturbation theory in (quasi) de Sitter background. We know that in usual set up of primordial cosmological perturbation such heavy fields are not appearing in the low energy effective field theory action. For that case in the simplest situation we actually start with an one field set up where the kinetic term is canonical in nature and the field is minimally coupled with the background gravity which is treated to be classical usually. Also such field has an effective structure of the interaction potential which play crucial role to study the time dynamics in FLRW cosmological background. Here specifically the field is treated to be massless compared to the Hubble scale (m << H). However, this is not the complete story yet. To explain this let us start with a Ultra Violet (UV) complete set up of quantum field theory (QFT) such as string theory in higher dimensions. There are various examples of string theory from which one can start the computation, which are -Type II A, Type II B, Heterotic, Mtheory etc. Also the low energy extension of such theories (supergravity) are also useful for the computation in the context of cosmology. Here it is important to note such all such theories contain massive (m >> H), intermediate mass (m ≈ H) and massless (m << H) fields in the matter multiplet. To write down an effective field theory (EFT) one need to integrate out all such heavy degrees of freedom from the UV complete version of the action 4 .
1. In this context, one can construct an EFT by utilizing the underlying symmetries appearing in the field theoretic set up. In such a generalized description where EFT is constructed by following the top down approach, we really don't care about the exact UV completion of the parent theory i.e. detailed quantum field theory origin at high energy scale of such effective constructions are not important in this case. See ref. [4,24] for more technical details.
2. In a more generalized prescription of EFT one can construct the set up which requires to correctly account for all relevant self interactions of adiabatic modes around and after the cosmological horizon crossing. Specifically the adiabatic mode contains all types of EFT relevant operators, including transient reductions in the effective sound speed c S each time the background field undertakes nongeodesic motion in background target space. In an EFT framework with single field setting, where heavy directions are such that the mass of the field under consideration is heavy compared to the Hubble scale i.e. m >> H, one gets transient drops in the effective sound speed c S during slow roll if the potential is such that the field traverses a bend even if the parent theory consists of canonically normalized scalar fields. So for general consideration one can allow many more possibilities without following any restriction to time dependent mass profile. However, these three specific types of time dependent mass profiles are very popular in the context of the study of quantum critical quench in a analytical fashion. We have considered these profiles particularly as our future objective is to apply the idea of quantum quench in the context of De Sitter space to quantify randomness during reheating phase. It is important to note that, using a simple field redefinition at the level of quantum fluctuation to the Mukhanov-Sasaki variable which results in a time dependent mass for the rescaled variable appearing with additional contributions of the mathematical form,ċ S /c S ∼ SH. Here S =ċ S /Hc S is the associated slow roll parameter with the effective time dependent sound speed c S and H is the Hubble parameter and it is associated with the changes in the radius of curvature of the inflaton trajectory. In this case the effective sound speed is given by, c −2 S ≡ 1 + 4φ 2 /κ 2 M 2 , where κ is the radius of curvature of the background inflaton (φ) trajectory and M is the effective cut-off scale of the EFT at high energy (UV) regime. Equivalently, it refers to the degree to which effective sound speed c S is reduced, which actually quantify the distance from the adiabatic minimum of the potential in the background inflaton trajectory is forced by its evolution. Each of these possibilities has different applications in the low energy limiting region of EFT. When the effective sound speed c S << 1 andċ S ∼ 0 is fixed over few e-folds of expansion then it is extremely difficult to maintain a meaningful derivative expansion without considering other types of special symmetries appearing in After doing dimensional reduction along with applying various compactification techniques one can derive various types of UV complete effective field theories at cosmological scale where the effective couplings of various relevant and irrelevant Wilsonian operators have time dependent profile in FLRW background and in such a case from the relevant quadratic operator one can also get the time dependent effective mass which is in general heavy (m ≥ H). It is further important to mention here that, such heavy fields can give rise to non vanishing one point function for scalar (curvature) perturbation in cosmology, which carries the signature of Bell's inequality violation in primordial universe [46,[61][62][63][64][65][66]. Also it is important to note that such Bell violating set up can be explained using the theory of quantum entanglement in (quasi) de Sitter background and can give rise to non-vanishing quantum information theoretic measure i.e. Von Neumann entropy, Rényi entropy, quantum discord, logarithmic entangled negativity [67][68][69][70] etc. Additionally, one can get correct expression for two point function and also the three point function from scalar (curvature) perturbation, which will show significant effect in estimating primordial non-Gaussianity from single field models of inflation. Apart from this one can consider a simplest situation in four space-time dimensions where the cosmological dynamics is explained in terms of two interacting scalar fields. The light field (m << H) is participating in inflation and the other heavy field (m >> H) is participating to explain the dynamics of reheating. If we path integrate out the reheating degrees of freedom then we get an effective field theory of inflation which is exactly same as we have explained earlier. But here one can consider the other possibility as well in which one can path integrate out the light inflaton degrees of freedom and write down an effective field theory to describe reheating in terms of the heavy fields (m ≥ H). In such a description this reheating field have mass and in the effective field theory description one can write down some time dependent coupling in terms of the integrated inflaton degrees of freedom and the mass of the reheating field appearing in the coefficient of the relevant quadratic operator. In this description such time dependent coupling is treated to be the time dependent effective mass parameter profile which is considered in the present discussion. So it is evident from this discussion that using both the effective field theory of inflation and reheating one can actually explain the the set up. However, as mentioned earlier, within certain limits one can consider an adiabatic region where the effective sound speed c S << 1 andċ S ∼ c S H and S =ċ S /Hc S ∼ 1 is fixed over a very small e-fold of expansion and this in turn generate all possible consistent transient strong coupling parameters without violating perturbative uniterity and these terms are explicitly appearing in the derivative expansion in the EFT. Consequently, the nature of these two types of features in the effective sound speed c S give rise to distinctive contributions to the physical observables studied in the EFT set up. The positive detection of these physical observables in different experiments allow to extract the underlying non-trivial physics from the EFT set up. In the technical ground the adiabatic mode is identified with the Goldstone boson, which is appearing due to spontaneously broken time translational symmetry prior to the path integration of the heavy fields. In this context, the invariance of the parent theory completely fixes the entire non-perturbative structure of all possible Wilsonian EFT operators and the associated coupling parameters can be expressed entirely by the effective sound speed c S of adiabatic perturbations, where the adiabaticity conditions c S << 1 andċ S ∼ c S H are respected. In principle, c S can be computed in terms of the parameters of the parent theory. Thus the additional contributions appearing in the adiabatic limit c S << 1 andċ S ∼ c S H directly justifies the validity of our treatment in this paper. For further technical details of this EFT set up see ref. [71][72][73]. origin of such time dependent effective mass profiles in four dimensions. However, in this paper since our objective is to study the cosmological particle production phenomena, we will mostly focus on the reheating epoch of the universe.
The dynamics of this fluctuating scalar field 5 in FLRW cosmological background with a time-dependent coupling obeys the following Klein-Gordon equation 6 : where m 2 (τ ) is the time dependent mass of the scalar field with which is originating from the effective field theory (EFT) of massless scalar field coupled with other heavy degrees of freedom by following the two possibilities: 1. In EFT time dependent couplings are appearing after path integrating out the massive degrees of freedom. This prescription is usually used to construct a most generic EFT of inflation.
2. In EFT time dependent couplings are appearing after path integrating out the massless degrees of freedom. This prescription is usually used to construct a most generic EFT of reheating.
Here φ k (τ ) is the associated Fourier mode of the fluctuating scalar field with momentum k, where it plays the role of wave number in the present context. In this paper, our prime objective is to find a precise equivalence between the dynamics of this scalar field resulting in stochastic particle production in cosmological events during reheating and the similarity with the dynamics of the electron transport in conduction wires. To establish this equivalence we start with the fact that the above mentioned Klein-Gordon equation for the fluctuating scalar field in (quasi) de Sitter background shows a striking similarity with the time-independent one dimensional Schrödinger equation appearing in the context of quantum mechanical system which describes the space evolution of electron inside a wire in presence of impurity as given by: where, V (x) corresponds to the time-dependent potential which is appearing as an outcome of impurity in the electrical wire and plays the similar role of negative of the square of timedependent mass profile as appearing in the context of cosmology i.e. −m 2 (τ ). Also E 5 Here it is important to note that, for inflation this scalar field is actually massless and in the effective field theory description one can construct the time dependent effective mass profile. On the other hand, in the context of reheating the scalar field is massive and in the effective field theory description one can construct time dependent effective mass in terms of the original mass of the reheating field and other degrees of freedom which are integrated out from the original theory. 6 Here we have assumed that the effective sound speed parameter, c S = 1, which indirectly implies the fact that for background time evolution we are considering a single scalar field with canonical kinetic term minimally coupled to the gravity. Effective mass of the scalar field is m(τ ), which has time dependent profile. However, one can generalize this prescription for any general non-canonical single field (i.e.P (X, φ) theory) theoretic framework where the effective sound speed parameter c S = 1.  Table 1. A brief overview of the connection between the scattering problem in quantum mechanics to that of cosmological particle creation events.

Scattering in conduction wire
represents the energy eigen value which mimics the role of the wave number squared i.e. k 2 . Finally, ψ(x) represents the wave function of the quantum mechanical system under consideration which is similar to the Fourier modes of the time dependent fluctuating scalar field in the context of cosmology i.e. φ k (τ ). The above set up can be re-expressed in terms of solving a transfer matrix problem since the scatterers can be thought as potential profiles in Schrödinger problem in quantum mechanics with the incoming and outgoing modes of the scalar field related to each other with the Bogoliubov coefficients. A complete overview of the connection between the variables that quantify the scattering problem in the context of quantum mechanics to the one in the cosmological particle production problem has been shown in table (1).
It is very well known fact that the conductance of the electrical wire is related to the transmission probability of electrons across the wire and this can be obtained by explicitly solving the time-independent Schrödinger equation (see Eq (2.3)) in the presence of the impurities. Before going to the further details of the computation here we begin by reviewing the scattering problem by a single impurity localized at the position x = x j . To the left (L) and the right (R) of the impurity potential, the wave-function can be written as a linear combination of right-propagating waves (exp (ikx)) and the left-propagating waves (exp (−ikx)) as: This is essentially a scattering problem in the context of quantum mechanics in which the impurities act as interaction potentials or scatterers across which the electrons get transmitted within the conduction wire. The map between the Bogoliubov coefficients (β R , α R ) from the right (R) side and the Bogoliubov coefficients (β L , α L ) from the left (L) side can be expressed in terms of the following Bogoliubov transformation equation as: where we define: 5) and in this context the transfer matrix for the j-th scatterer M j is given by the following expression: which is essentially an unitary matrix related the incoming and the outgoing wave functions and their normalization coefficients. Ultimately, using this methodology our objective is to connect several impurities together. This is particularly very easy to describe in terms of the transfer matrix approach, since the total transfer matrix across N s number of scatterers is simply given by the simple matrix multiplication of the individual transfer matrices as given by the following expression: For our choice of convenience of symbols we will drop the term N s for the N s number of scatterers and hence we will be considering this to be equal to M. In Fig:3 we show the electron(wave) encounter a potential(impurity or scatterer).It transmit and reflect through it.From simalirty of Klein-Gordon equation and the time-independent one dimensional Schrödinger equation we calculate R and T for particle production event. Further, let us consider the simplest possibility of having two (N s = 2) scatterers across which the transmission probability can be written as:  where the transmission and reflection coefficients for the j − th scatterer can be expressed as: and additionally e iθ is the overall phase factor which describes the shift in phase between the reflecting waves across the scatterers due to the presence of impurities. If the distance between the two impurities is random in nature and uniformly distributed over a region with the assumption, k∆x >> 1 (where ∆x = x 2 − x 1 is the distance between the scatterers), then the phase θ is also uniformly distributed over the interval 0 < θ < 2π. Using this fact explicitly we take logarithm on both sides of the above equation and further doing average over the phase within the interval 0 < θ < 2π we finally get 7 : The phase-averaged logarithm of the total transmission probability across N s number of 7 Following this discussion, one can generalize this statement for N s number of scatterers as: scatterers then further simply can be written as: where γ is known as the Lyapunov exponent, which is defined as: This actually determines the rise of chaos in the system. Using this information the typical transmission probability is defined as: which corresponds to the most probable transmission probability in the ensemble of random potentials. Also it is important to note that, represents the total length of the conduction wire. Here the localization length is defined as: In one spacial dimension, the localization length is of the same order as the transport mean free path as pointed in ref. [74,75]. If the mean distance between scatterers, ∆x, and the average logarithm of the transmission probability per scattering, γ, are fixed, then the total transmission probability decays exponentially with the length L of the conduction wire 8 . This is commonly known as Anderson localization [76]. Naturally, it is well known that the resistance of the conduction wire scales inversely with the total transmission probability. At zero temperature, all one-dimensional conduction wire are therefore can be treated as an insulator, which is independent of the strength of the impurities appearing in the wire. However, the mathematical structure of the total transmission probability T is preserved for N s number of such scatterers and this can be shown as: For a one-dimensional non relativistic electron in conduction wire under the influence of certain potential V (x) evolution of the wave function ψ(x, t) is given by: with the Hamiltonian for this particle is given by: If we consider the particle is initially prepared in presence of potential V (x) wave-packet take the specific form of ψ(x, t). The final stationary density distribution |ψ(x, t)| 2 at long time carries important information both in their average and fluctuations. The quantum mechanical wave can tunnel through potential hills and reflect for by small fluctuations. So the initial wave packet split on each potential fluctuations into a transmitted and a reflected part. After huge number of scattering instances this reduce to a random walk problem and on average the motion at long times will have the diffusion constant in it. This is exactly the case of electron is propagating in a conduction wire. At long times average dynamics [77] of the wave packet freeze and it takes the shape as given by the equation Here, ξ is the localization length as discussed in Eq:-2.16. An electron in random potential is normally studied using statistical ensemble of random one-electron matrix Hamiltonians. Using Tight Binding approximation in orthogonalized lattice-site basis representation. The diagonal matrix elements are chosen from a flat probability distribution of width W,the strength of disorder. The off-diagonal hopping matrix elements for every pair of nearestneighbour sites and represented by 2x2 matrix where potential take the form, where, I, σ x , σ y , σ z are identity and Pauli spin matrices forming the complete basis set.
Here µ is the random spin-orbit coupling strength, t x , t y , t z are independent random variables taken from uniform distribution on interval [−1/2, 1/2]. The metal-insulator transformation occur at specific values. Below that mobility edges appear in band separating localised states near edges from extended states near band center. The tight-binding random matrix ensembles (TBME) classified scheme is possible on symmetry. Orthogonal Ensemble in random potential..Localisation and mobility occur in all 3-D tight-binding ensembles and in 2-D for symplectic and unitary classes. From this distinction one found striking similarities with symmetry classification of Gaussian random matrix ensemble. The Gaussian ensemble belongs to high dimensionality limit of TBME and always metallic. So the metallic phase is well approximated by Gaussian random matrix theory. From our discussion on RMT we use the Nearest Neighbour Spacing Distribution function [P(ω)see-Eq5.1] and measure it in units of mean level spacing ∆. Around the mobility edge and intermediate law from P (ω) can be obtained. Anderson localization in this context mimics quantum chaotic transition. The fluctuations for the density of states are partially responsible for the conductance fluctuations. Although average density of states is insensitive to Anderson transition its higher order moments are sensitive to it. On an other approach we can relate Anderson localization to RMT using Lyapunov Exponents. Equation 2.13 and 2.16 relates Anderson localization to Lyapunov exponents. Now statistical property of the Lyapunov spectrum with large number of degrees of freedom can be described universally by RMT. [78] As described in [78], the spectrum of Lyapunov exponents is well approximated by the following expression: Here λ max is the time independent parameter which approximately equals to bound of Lyapunov exponent. This equation shows striking similarity with Wigner law [Eq:-5.146]. In this approach we can also show the connection between Anderson localization and RMT. But there is a striking difference too. Random matrix theory takes all its entries from Gaussian random variables but for electronic models [Scattering matrix theory] matrix ensemble have short-ranged and sparse random matrix with most of the matrix elements having main diagonal non-zero.
3 Randomness from conduction wire to cosmology: Dynamical study with time dependent protocols In this section, our objective is to explain the various features from the time dependent effective mass profiles which are related to the quantum mechanical scattering problem in conduction wire as mentioned earlier. These features are appended bellow: 1. Lyapunov exponent: It actually quantify the amount of chaos appearing in the quantum mechanical systems that we are studying in the context of early universe cosmology. In our discussion it tells us the degree of randomness in the stochastic particle production. In our case, the chaos emerges due to the random scattering events which are non adiabatic and we call these as cosmological scattering events leading to particle production. In this section, we discuss about Lyapunov exponent and try to discuss their behaviour for the different time dependent mass profiles. In thus context, Lyapunov exponent is defined as [42,79]: where, T is the transmission coefficient given by the following expression: with t and t * being the transmission amplitude of the incoming and the outgoing wave.
In the present discussion, the transmission coefficient can be expressed as: where β and α are the Bogoliubov coefficients. Also, it is important to note that, in the present context one can define the reflection coefficient as: wherer andr * being the reflection amplitude of the incoming and the outgoing wave. Finally from Eq (6.176) and Eq (3.4), we get the following conservation equation: where we have used the following normalization condition for the Bogoliubov coefficients, as given by: 2. Conductance: It quantify the degree of support of the flow of electron inside an electrical conduction wire. this is exactly reciprocal of resistance. In the present context, conductance refers to the ability of the massless scalar fields to transmit through the massive fields which are the specific heavy mass profiles that we have discussed above. This may be more suggestive in telling us about the interaction of the massless scalar field with the massive fields. More value of conductance refers to the larger transmitivity of the background fields through the scatterers and vice-versa. Thus, conductance also carries a valuable information about the transmission coefficient of the scalar field interacting with the scatterer. In this context, the conductance can be expressed as: where λ is the Lyapunov exponent, T is the transmission coefficient, |t| is the transmission amplitude of the incoming/outgoing wave and β, α are the Bogoliubov coefficients as mentioned above.
3. Resistance: It quantify the degree of oppose of the flow of electrons inside an electrical conduction wire. It is the property by the virtue of which the scatterers (which are the time dependent mass profiles in our case) resist the massless scalar field to tunnel through them. In other words, it is the same Schrödinger formulation in quantum mechanics where the incoming wave interacts with a potential barrier and the strength of the barrier is the measurement of resistance to the tunneling of the incoming particle through it. This means that more the resistance to the incoming wave, more is the lower is the transmission probability across the barrier. Resistance is defined as the reciprocal of conductance G(k), which gives: (3.8) We will discuss details of these features for three different mass profiles as mentioned in Eq (1.4). All of these mass profiles that we choose here mimics the role of scatterers inside the conduction wire. Such scatterers provide the way for scattering events to occur resulting in random particle production in cosmological space-time. To study the cosmological particle creation problem during early epoch of universe (specifically during reheating) we use the analogy with the quantum mechanical scattering problem inside an electrical conduction wire in presence of time dependent effective mass profile we will perform the computation in (quasi) de Sitter space using FLRW spatially flat metric.
Here we consider a massive free scalar field with time-dependent mass 9 : where the scalar field satisfies the following constraint: (3.10) and the Fourier transform of the field is defined as: Also in the (quasi) de Sitter background the scale factor a(τ ) can be expressed in terms of conformal time as 10 : where is the slow-roll parameter in quasi de Sitter space, which is defined as: (3.14) 9 Here it is important to note that our approach is similar to that of used in refs. [58,80] to explain the time dynamics of quantum quench. 10 In de Sitter and quasi de Sitter space one can compute the relation between the conformal time (τ ) and the physical time (t) as given by the following expressions: Here, we define a new slow-roll parameter with respect to the conformal time:

De Sitter
Now we use the following field redefinition in Fourier space: Consequently, the scalar field action as stated in Eq (3.17) can be recast in terms of the newly defined field φ k (τ ) as: Further, varying the above action with respect to the redefined field φ * k (τ ) we get the following equation of motion: Further, Eq (3.18) can be simplified for de Sitter and quasi de Sitter space as: Quasi De Sitter : It is important to note that, the main contribution to particle production is originating from the excitations of the field with k/a >> m >> H, at the stage of oscillations. Therefore, in the first approximation we can neglect the expansion of the Universe, taking the scale factor a(τ ) as a constant during reheating. We call it reheating approximation 11 Consequently, 11 Important note: In the present context, the analysis is perfectly valid for the highly localized particle production events after neglecting the cosmological expansion during reheating approximation. But this approximation fail for the events that are sufficiently spaced out. If we don't neglect the cosmological expansion in this computation then the conformal time dependent mass term of the form 2τ 2 is restored from the background cosmological background. This actually implies that the scattering problem is being performed on a conformal time dependent potential of the form 1/r 2 (inverse square), which makes the analytic computations of the Bogoliubov coefficients and all the other derived physical quantities to quantify quantum randomness from the present set up extremely difficult. Here it is important to note that, for long wavelength cosmological observables particle production appears more than an e-fold apart and consequently the corrections appearing due to the cosmological expansion seem certainly relevant in the computation as the incoming and outgoing wave functions depart from plane waves. Although for localized particle production events, the reheating approximation considered in this paper perfectly holds good. In the present context this approximation breaks when we consider the particle production events for a sustained period of time or may be separated by times approaching an e-fold expansion.
one can approximately write Eq (3.18) in the following simplified form 12 : The Fourier modes of the scalar field follow the equation of motion in as stated in Eq 3.22, with every Fourier mode satisfying the Schrödinger equation where −m 2 (τ ) playing the role of a potential. In Fig:4 the particle produced show fluctuation from ground state and from Figure 4. This diagram shows that ground state fluctuations from the past can in future be amplified which can be measured by the coefficient α whereas particle excitation from ground state can be measured by β.
calculating the Bogoliubov coefficients we predicted all its properties. For the solution we refer to ref. [36], for the field φ k (τ ) can be expressed in two distinctive ways, as given by: where u in,in (k, τ ) and u in,out (k, τ ) are the 'ingoing' and 'outgoing' wave-functions. Also, the in-and out-oscillators are related to each other through the Bogoliubov coefficients α(k) and β(k) Here it is important to note that, since the scale factor a(τ ) is approximately a constant during reheating (reheating approximation), then conformal time (τ ) and the physical time (t) is related through the following coordinate rescaling transformation: (3.21) Now, we calculate the various electrical properties and also the expression for the Lyapunov exponent to quantify quantum chaos for the various time dependent effective mass profiles which are equivalent to the impurity potential term in the time Independent Schrödinger Equation describing a scattering problem inside a conduction wire.  Here we start with the following mass profile: The corresponding Schrödinger problem for this potential function can be solved by using the potential function as given bellow: We can find the following explicit solutions for u in (k, t) and u out (k, t), as given by: where we define ω ± , ω in and ω out in the following: (3.29)

Bogoliubov coefficients
For this specific mass profile the Bogoliubov coefficients can be expressed as:

Optical properties: Reflection and transmission coefficients
For this specific mass profile the transmission and reflection coefficients can be expressed as: (3.31)   coefficients with wave number k.

Chaotic property: Lyapunov exponent
For this specific mass profile the Lyapunov exponent can be expressed as: (3.33) In fig. 8, we observe that with increase in wave number k the Lyapunov exponent decreases. This shows that the Lyapunov exponent is dependent on the momenta values of the fields interacting with the massive field acting as a scatterer. Furthermore, we discover that for the mass profile I, the chaos in the event reduces with increase in the wave number. This suggests that lesser the number of fields interacting with the massive field more is the chaos in the quantum system considered in this paper. Since, a negative value of Lyapunov Exponent pulls a system out of chaos, this further tells us that the Lyapunov exponent is inversely related to the number of background fields interacting with the scatterer or the massive field. This may be interpreted in the following way in the context of Schrödinger problem in quantum mechanics that a higher value of wave number k of the incoming wave would be able to cross a potential barrier of a given strength and would be able to get transmitted through the barrier and the pulse won't damp easily than that of a wave with lower k value. This means that the scatterer acts as a definitive medium which allows only certain wave numbers to pass through thus reducing the chaos in the system.

Conduction properties: Conductance and Resistance
For the given mass profile the expression for conductance and resistance can be expressed as: In fig. 9(a) we have shown the variation of conductance with wave number k. This figure shows that with increase in the momenta value of the massless scalar field, the conductance also increases. Now, accounting for m 0 values, we see that for m 0 = 1 the conductance shoots up at a much lower k value than that of m 0 = 2 and m 0 = 3. This suggests that for m 0 = 1 the field has a much higher transmission probability than that of m 0 = 2 and m 0 = 3. An increase in transmission probability gives a direct evidence of the conductance value. Therefore, we conclude that for m 0 = 1 the field has more conductance value in comparison to m 0 = 2 and m 0 = 3. We also conclude that larger the momenta value, more is the transmission coefficient and thereby shoots up the conductance of the system. This means that an incoming wave with large momenta value would eventually cross a barrier potential field thereby increasing the conductance of the system as the transmission probability would be much higher than an incoming wave with lower momenta value.
In fig. 9(b) we have shown the variation of resistance with wave number. We observe that with an increase in the value of k the resistance starts decreasing which suggests that   with an increase in momenta value the transmission probability across the scatterer. This may be viewed in accordance with the potential barrier in the Schrödinger equation in quantum mechanics also starts increasing thereby allowing the incoming wave to tunnel through the barrier thereby increasing the transmission probability and hence,reducing the resistance. We also observe that with an increase in k value the resistance reduces less rapidly for m 0 = 1 than that of m 0 = 3 and m 0 = 2. Whereas, it reduces more rapidly for m 0 = 3 suggesting that higher the value of the constant m 0 lower is the value of resistance offered.
3.2 Protocol II: m 2 (τ ) = m 2 0 sech 2 (ρτ ) Here we consider the following mass profile: The corresponding Schrödinger problem for this potential function can be solved by using the potential function as given bellow:   Now using the coordinate transformation y = e 2ρτ [58,80] one can recast the equation of motion, analogous to the time independent Schrödinger equation takes the following form: The solution of this equation is given by: where we define a parameter α as:

Bogoliubov coefficients
Now we fix C 1 = 1 and C 2 = 0, which gives the incoming solution u in (k). Further taking the t → +∞ limit and using Bogoliubov transformation we can express incoming solution in terms of the outgoing solution as given by: where α(k) and β(k) are the Bogoliubov coefficients, which are defined as: , β(k) = i sin(πα)cosech πk ρ .

Optical properties: Reflection and transmission coefficients
For this specific mass profile the transmission and the reflection coefficients can be expressed as: In fig. (12(a)) and fig. (12(b)), we have shown the variation of the transmission and reflection coefficients with wave number k.

Chaotic property: Lyapunov exponent
The Lyapunov exponent for this case may be given as:   In fig. 13, we have shown the wave number dependence of Lyapunov exponent. Here we observe that with increase in k value the Lyapunov exponent decreases. This implies that the Lyapunov exponent is dependent on the momenta values of the fields interacting with the massive field acting as a scatterer. Furthermore, we also discover that for this mass profile II, the chaos in the event reduces with increase in the k value.

Conduction properties: Conductance and Resistance
For this specific mass profile the expression for the conductance and resistance can be computed as:  In fig. 14(a) we have shown the wave number dependence of conductance. We observe that for m 0 = 1 conductance starts increasing at a larger value of k than that of m 0 = 2 and m 0 = 3. But, in contrast to the variation of conductance with momenta k in the above figure, here the conductance starts increasing rapidly for m 0 = 3 than that for m 0 = 1 which suggests that the transmission probability for m 0 = 3 is much higher than m 0 = 1 and m 0 = 2, thereby making it more conductive than the other two.
In fig. 14(b), we have shown the wave number dependence of resistance. Here like the first mass profile the resistance for m 0 = 3 falls more rapidly than that of m 0 = 2 and m 0 = 1. This suggests that for the given mass profile II, as the value of m 0 increases, the value of resistance also decreases. But unlike the first mass profile, the resistance for m 0 = 3 falls more rapidly suggesting that for m 0 = 3 this specific mass profile offers more resistance than the first one. Therefore, we conclude that for the same values of m 0 this mass profile offers less resistance in comparison to the first mass profile.   Here we consider the following time dependent mass profile: This Θ function in τ makes the mass profile a quenched one. The corresponding Schrödinger problem for this potential function can be solved by using the potential function as given bellow:

Bogoliubov coefficients
For this specific mass profile the Bogoliubov coefficients can be expressed as: with the solution of the incoming and the outgoing waves are given by the following expressions:   coefficients with wave number k.

Optical properties: Refeclection and transmission coefficients
For this specific mass profile the transmission and the reflection coefficients can be computed as: (3.51) In fig. (17(a)) and fig. (17(b)), we have shown the variation of the transmission and reflection coefficients with wave number k.

Chaotic property: Lyapunov exponent
The Lyapunov in this case is written as:    From fig. 18 we observe that with increase in wave number the Lyapunov exponent decreases more like a rectangular hyperbolic fashion. In comparison to the other two mass profiles where the reduction in the value of the Lyapunov exponent is much less rapid in comparison to this mass profile discussed here. This suggests that since, the mass profile is a heavy side theta function, which is a quenched mass protocol, the Lyapunov exponent also gives a similar like profile. This shows that the Lyapunov exponent is dependent on the wave number of the fields interacting with the massive field acting as a scatterer. Furthermore, we discover that for this given mass profile, the chaos in the event reduces with increase in the k value. So, in this case the Lyapunov exponent decays much rapidly than the first two mass profiles. Next, we will try to find an upper bound of Lyapunov exponent using the definition of [81]

Conduction properties: Conductance and Resistance
For this specific mass profile the expression for the conductance and resistance can be written as:   In fig. 19(a), we have shown the wave number dependence of conductance. This figure shows that with increase in the wave number of the massless scalar field, the conductance also increases. Now, accounting for m 0 values, we see that for m 0 = 1 the conductance shoots up at a much lower k value than that of m 0 = 2 and m 0 = 3. This suggests that for m 0 = 1 the field has a much higher transmission probability than that of m 0 = 2 and m 0 = 3. An increase in transmission probability gives a direct evidence of the conductance value. Therefore, we conclude that for m 0 = 1 the field has more conductance value in comparison to m 0 = 2 and m 0 = 3.
In fig. 19(b), unlike the mass profile I the resistance for m 0 = 1 falls more rapidly than that of m 0 = 2 and m 0 = 3. This suggest that for the given mass profile, as the value of m 0 increases, the value of resistance also increases suggesting that heavier the field gets lesser is the transmission probability of the incoming wave to tunnel through it thereby reducing the value of conductance for this specific mass profile.
4 Quantum chaos from out of time ordered correlators (OTOC) 4

.1 Chaos bound in out-of-equilibrium quantum field theory (OEQFT) and its application to cosmology
We know that in the context of quantum field theory it is possible to achieve the following universal bound on the Lyapunov exponent [50] 13 14 : Universal chaos bound in OEQFT : where k B is the Boltzmann constant and T is the temperature associated with the dynamical system. This upper bound of the Lyapunov exponent is treated as the saturation bound of chaos 15 . Our aim is to establish this bound in the context of cosmology and study its further consequences. This bound was first proposed in the context of quantum information theory of black hole [82][83][84][85]. Additionally it is important that, the bound on the Lyapunov exponent saturates in the context of Sachdev-Ye-Kitaev (SYK) model [56,57,[86][87][88][89][90][91][92], which describes the quantum features of Majorana fermions in presence of infinitely long range disorder. Saturation of the Lyapunov exponent implies that SYK model mimics a quantum 13 For this specific discussion only we keep the Planck's constant and the Boltzmann constant k B in our computation. But for the rest of the paper we fix = 1 and k B = 1 for which the parameter β can be written as, β = 1/T . In such a situation chaos bound is given by, λ < 2π/β. 14 In the context of weakly coupled gauge theory one can introduce 't Hooft coupling λ T which is independent of N and in such a theory the Lyapunov exponent is given by the following expression: Lyapunov exponent in gauge theory : 15 Considering the bulk contribution weakly coupled with string theory with large radius of curvature one can show that the perturbative stringy correction to the Einstein gravity computation of the scrambling can give rise to the following first order corrected expression for the Lyapunov exponent []: where L s is the stringy length scale and µ 2 is a specific constant which is appearing in the shock wave equation propagating along the horizon.
description of black hole via AdS/CFT correspondence. In the strict classical limit → 0 the Lyapunov exponent take any values, which is consistent with the requirement.
To give an explicit derivation of the chaos bound on Lyapunov exponent in the context of cosmology let us follow the steps appended below: 1. Let us start with a completely mathematical problem described by a time dependent function g(τ ), which satisfy the following set of properties: (a) In the complex plane g(τ + iT ) is analytic in the half strip described within τ > 0 and − β 4 ≤ T ≤ β 4 . In this context, τ and T represent the real and imaginary part of the complex number τ + iT after analytical continuation in complex plane. (b) The function g(τ ) is completely real at T = 0. (c) After analytical continuation the function in the complex plane satisfy the following constraint: which is perfectly valid in the complete half strip.
2. Next, we actually conformally map the entire half strip to a unit thermal circle in the complex plane, which can be done using the following Möbius transformation: where ∆ β (τ + iT ) is the temperature dependent function in the complex plane, described by the following expression: 3. Further using Eq (4.4) one can further say that the complex function g(z) is an analytic function from one to one conformal map from unit disk to unit disk.

A variant of the Schwarz lemma can be represented as a invariant contribution under
analytic automorphisms on the unit disk, which implies the bijective holomorphic mappings of the unit disc to itself. This specific variant is known as the SchwarzPick theorem.
5. Now the hyperbolic metric in complex plane is defined as: (4.7) Further using this metric and applying SchwarzPick theorem one can write: .  6. Further applying the fact that the function g(τ ) is real at T = 0 and using Eq (4.8) we get the following simplified result: 7. Further, rearranging Eq (4.9) we get the following final result: which is the outcome of Schwarz-Pick inequality and very very useful to prove the universal chaos bound in OEQFT.
Now it is important to note that in this context, This further implies that: Now at very large time scale (τ → ∞) or at very high temperature (β = 1/T → 0) one can neglect the contribution from the second sub-leading term. As a result we get the following inequality: 8. Further, we take the following phenomenological function: where k is constant and λ is the Lyapunov exponent. This function satisfy all the requirements that we have mentioned earlier explicitly. Further substituting this function in the result obtained in Eq (4.10) we get the following simplified result 16 : which proves the Universal chaos bound in OEQFT.
9. This bound on the Lyapunov exponent is an unique feature of all classes of OEQFT set up. It has a very strong impact in the context of early universe cosmology, specifically during reheating epoch. By knowing specific time dependent couplings in the context of effective field theory (EFT) it is possible to give an estimate of Lyapunov exponent in such OEQFT set up. We will show this feature in the next section for three known model of interactions appearing in EFT. In such a situation one can give an estimate of the upper bound on reheating temperature using this bound, which is again obviously an universal bound itself. The earlier study in the context of reheating actually predicts a very crude estimate of reheating temperature which is based on the assumption that reheating is extremely model dependent. It actually means that to write an EFT of reheating we need to know the all interacting relativistic degrees of freedom in a specific model. In this framework the total energy density during reheating can be expressed in terms of total number of relativistic degrees of freedom by the following expression: Using this expression of energy density during reheating epoch one can able to express the reheating temperature as: where g * (T reh ) is the effective number of total relativistic degrees of freedom present in the thermal bath at temperature T = T reh and V reh is the scale of reheating which can be obtained by fixing the field value at φ = φ reh for a specific model. Counting all the degrees of freedom in the particle physics model one can fix g * (T reh ) in the present context. To find the reheating constraint from the prescribed set up let us further introduce the number of e-foldings at the epoch of reheating, which is defined as: 16 Henceforth we set = 1 for which the bound is translated to λ ≤ 2π β , which we will use for the further application purposes.
where N total is the total number of e-foldings which is defined as: Here t e , t i and t reh are the representative time to specify end of inflation, starting of inflation and time scale at the end of reheating respectively. Similarly φ e and φ reh are the field values at the end of inflation and reheating respectively, which can be computed for a given known model of inflation. Also it is important to note that in this context, ∆N is defined as: Here ∆N is defined as: From different models of inflation one can explicitly compute e-foldings at horizon exit, which is given by the following expression: Consequently, the value of ∆N from observation can be estimated as: Now, to give a numerical estimate of the reheating temperature let us consider the following simplest monomial model: where V 0 fix the overall scale of the potential and p is the degree of the monomial which depends on the characteristic of the model. For this model the field value during reheating can be expressed as: The reheating scale is quantified in terms of the number of e-foldings as: Consequently, for the monomial model the reheating temperature can be quantified as: Reheating bound from model : Here V inf is the scale of inflation which is quantified by the following expression: Upper bound on inflationary scale : As a result, we get the following bound on the reheating temperature: Upper − bound on reheating temperature from inflation : which is true for any models of inflation. From the Planck 2018+BICEP2/Keck Array BK14 data the upper bound on the tensor-to-scalar ratio (primordial gravitational waves) is restricted to: where k * ∼ 0.05Mpc −1 is the pivot scale of momentum. This implies that the upper bound of reheating temperature from the Planck 2018+BICEP2/Keck Array BK14 data is given by: Here to writing down this expression for reheating temperature it is important to consider the following assumption: (a) Contribution from the kinetic term of the field which is mainly responsible for reheating is neglected.
(b) We also assume that reheating is described by scalar field.
This further implies that depending on the background particle physics model reheating temperature actually varies in a wide range and one cannot able to determine exactly its value as there is no such universal bound available earlier in this context. This is the main shortcoming of the phenomenological prediction of reheating temperature in the context of early universe cosmology.
On the other hand, just only considering the dynamical details of quantum chaos one can express the reheating temperature in terms of the Lyapunov exponent: Universal lower − bound'on reheating temperature : which is an universal lower bound on reheating temperature in the present context of discussion as it is not involve any model dependence from the background theory. This implies that the universal bound on quantum chaos in OEQFT restrict us to fix an universal model independent lower bound on reheating temperature. Combing the obtained bound in this paper and the upper bound obtained from inflation one can restrict the reheating temperature within a specified range. Additionally, the present analysis helps us put an unique upper bound on the Lyapunov exponent in terms of the scale of inflation (or tensor-to-scalar ratio) as:

What is OTOC?
Now it is important to note that the universal bound on quantum chaos can be achieved by computing the out of time ordered correlators (OTOC), which in general can be expressed in terms of commutators. In the study of quantum chaos, specifically in the context of Butterfly effect one can introduce two time dependent operators W (τ ) and V (τ ) from which one can define a commutator, [W (τ ), V (0)], where the operators are in general Hermitian in nature and they have introduced with time separation ∆τ = τ − τ = τ with τ = 0. This commutator actually captures the effect of perturbation by the operator V (0) on the later time measurement on the operator W (τ ) and the converse statement is also true. The time dependence of the operator W (τ ) in this context of discussion can be expressed in the Heisenberg representation as: The strength of such chaotic effect is characterised by the following measure: Four point quantum operator where the expectation value is in general the thermal averaged 17 , which is defined as: Here Z is the partition function which is defined as: and H is the Hamiltonian of the chaotic system under consideration. Here it is important to note that to construct the chaotic OTOC measure instead of using the two point operator, [W (τ ), V (0)] (commutator), here we have actually used the four point quantum operator, [W (τ ), V (0)] 2 (square of the commutator). The specific reason for such choice is following. To describe this let us first assume that we replace the commutator bracket by the Poisson bracket by considering the semi-classical limiting situation. In such a case the Poisson bracket shows typically an exponential growth, exp[λτ ], where λ is the Lyapunov exponent. But the signature of its coefficient can be anything, either positive or negative. Now further if we take the thermal averaging over this two point operator then both the contributions are cancelled each other in the semi-classical limit and will not contribute to describe chaos. From the quantum mechanical perspective, the two point thermal averaged operator, [W (τ ), V (0)] actually captures the description of correlation between the quantum Hermitian operators W (τ ) and V (0)., which decays in the large time limit (τ → ∞) and cannot describe the chaotic behaviour. On the other hand, the four point quantum operator after transforming it to the Poisson bracket in the semi-classical picture don't show any ambiguity in the signature of the co-efficient as it takes only positive value. After taking thermal average we get non vanishing result using which one can describe quantum chaos. In the quantum mechanical picture the four point thermal averaged operator, [W (τ ), V (0)] 2 not decays exponentially at the leading order in the large time limit (τ → ∞). Now, in the quantum mechanical description of the Butterfly effect predicts the following result: for any mathematical structure of the operators V (0) and W (t). Here it is important to note that, V (0)W (τ )W (τ )V (0) contribution is not directly effected by the quantum chaos. Also it is important to note that, in the present context for the sake of simplicity we additionally assume that: i.e. both the one point function or the thermal averaged expectation values of these operators vanishes.

Estimation of scrambling and dissipation time scales from OTOC
In the context of quantum chaos two important time scales are associated:

Scrambing time:
Here the associated time scale where the operator C(τ ) is relevant is known as the scrambling time scale τ * . Sometimes in literature this is known as the Ehrenfest time scale. A possible distinction between the classical and quantum description of chaos can be described by the Ehrenfest time scale in which the previously mentioned OTOC don't grow with respect to the associated time scale and saturates at the same scale.
In the next section we have provided a alternative chaos bound on OTOC (i.e. SFF in our case) from which we have further give an estimate of the bound on the Ehrenfest time scale.

Dissipation time:
Another time scale for chaos is the exponential decay time scale τ d in which the two point thermal correlation function behaves like V (0)V (τ ) . Sometimes in this literature it is known as the dissipation time scale or the collision time scale. In the context of strongly coupled quantum field theories at finite temperature it is expected that the dissipation time scale τ d ∼ β. It is also expected that for large time limit the more general form of the OTOC during this time scaled as: where · · · represent higher order terms which are more suppressed by the dissipation time scale τ d .
In the present context, additionally one can predict the connection between the quantum mechanical operator C(τ ) and quantum chaos by considering the semi-classical limit of a chaotic system which involves a single particle. To demonstrate this argument one can consider semi-classical billiards as a toy example. In the semi-classical limit one can take, V (0) = p(0), W (τ ) = q(τ ), where p and q is the generalized momentum and coordinate respectively. As a result in the semi-classical limiting approximation one can map the previously defined commutator bracket to the Poisson bracket, as given by: which can be treated as the classical analogous version of the quantum mechanical Butterfly effect. It is also expected that for such a system the nearby dynamical trajectories scale as, q(τ ) ∼ q(0) exp [λτ ] , where λ is the Lyapunov exponent. It is in principle divergent in nature for large time limiting situation. Now at the dissipation time scale τ d it is also expected that, τ d ∼ 1/λ, for which the nearby trajectory is convergent and is of the order of e. On the other hand, the prescribed OTOC can approximately expressed in semi-classical limit as 18 : is the area of the stadium and v is the velocity of the particle. Then the classical OTOC can be expressed as: where Z cl is the classical partition function defined as: . (4.44) Further taking A = 1 for simplicity we get: Further taking the limit t >> √ β we get the following simplified answer for classical OTOC for billiards:  Now at the scrambling time scale, τ * and dissipation time scale, τ d the OTOC approximately in the semi-classical limit scaled as: from which the scrambling time scale, τ * can be estimated as: This further implies that, in the semi-classical limit the scrambling time scale, τ * and dissipation time scale, τ d are related by the following expression: which explicitly shows that both the time scales for quantum chaos is different from each other and the fractional difference is given by the following expression: which is actually a large amount of hierarchy at the semi-classical limit as → 0. Now, the OTOC in the present context actually quantify the temporal growth of the Hermitian quantum mechanical operator W (τ ) is is introduced earlier in this section. In the general prescriptions of quantum field theory (QFT) such OTOC can be expressed in terms of the addition of simple type of operators, which span the quantum basis. Now, if the OTOC is large 19 then in such a situation with non-local interactions the scrambling time scale, τ * can be estimated as: where N bit is the number of qubits. Similarly for local interactions he scrambling time scale, τ * can be estimated by computing the separation between the quantum operators W and V . Additionally it is important to note that, quantization of a classical chaotic system may accommodate positive Lyapunov exponent from the OTOC mentioned earlier.
To quantify quantum chaos also the nearest neighbour distribution (NDD) for the spectrum of the energy is alternatively used 20 . Except Lyapunov exponent, in the present context of discussion OTOC (in our discussion it is SFF) play crucial role to quantify quantum chaos to explain dynamical features in the early epoch of universe.

Quantum chaos from RMT: An alternative treatment in cosmology
In this section, we will try to generalize spectral form factor (SFF) for any order of even polynomial potential. To serve this purpose, one can create such ensemble, such that all possible interaction between energy levels of many-body Hamiltonian would be accounted for by various matrices in the ensemble. If the Hamiltonian is time-reversal symmetric the required distribution will be invariant under orthogonal transformation. Else, it is invariant under unitary transformation.
In the thermodynamic limit (N → ∞) eigen value of density of random matrices showed a universal behaviour characterised by Wigner's Semicircle law. The results seemed to be applicable to a varied class of quantum system displaying chaotic behaviour. Chaos was also a hallmark of a few-body Hamiltonian (N finite), but better diagnostic for quantum systems was devised in which nearest neighbour spacing distribution (NNSD) of eigenvalues of the system will be chaotic if distribution is Wigner Dyson type: Here it is important to note that, here β is fixed at, β = 1 for Gaussian orthogonal ensemble and β = 2 for Gaussian unitary ensemble. In the present context of discussion Spectral Form Factor (SFF) is a tool for characterising spectrum ( i.e. discreteness of energy spectrum) of quantum system under consideration and defined by the following expression: Here Z(β) is the partition function of the quantum system and β = 1/T . For β = 0, the expression pick out contribution only form the difference between nearest neighbour energy eigenvalues at very late times. SFF when averaged over Gaussian random matrices, has very particular behaviour at large N with initial decay followed by a linear rise and then after a critical point saturation.This approach can relate a saturation limit for large N which can be treated as bound on chaos. Additionally, it is important to note that quantifying chaos through finding SFF is very useful when one cannot have a specific time dependent mass profile during cosmological particle production. In terms of scattering problem in the conduction wire if we don't know precisely the structure of interaction potential, then one can quantify chaos in terms of SFF rather than using Lyapunov exponent, as we have used in the previous section. Here we will discuss general approach to find SFF to quantify chaos for various even polynomial potential.

Quantifying chaos using RMT
Gaussian matrix ensemble is a collection of large number of matrices which are filled with random numbers picked arbitrarily from a Gaussian probability distribution. See refs. [93,94] for more details.  In table (2), we have explicitly mentioned the properties of the each elements of the Gaussian matrix ensemble in Random Matrix Theory (RMT).
Further, the joint probability distribution of such random matrix, which is characterized by the Gaussian potential is given by the following expression: where N represents the rank of the matrix M . If we consider any ensemble of matrices to keep this measure invariant under similarity transformation: such that it satisfies the following constraint: Here U being an orthogonal or unitary matrix. Then for most generalized ensemble one can implement the concepts of time independent Random Matrix Theory [95] in the present context of discussion. Now here integrating over the random matrix measure one can construct the following expression for the partition function for the Gaussian matrix ensemble, as given by: Further, using similarity transformation one can diagonalize the random matrix M as: On the other hand, ensemble in basis of eigenvalues of the matrix the partition function can be written as: where the action S(λ i ) is defined as 21 : Here we fix β = 1 for GOE and β = 2 for GUE. The overall 1/N come from scaling of eigenvalues by factor √ N . To find a solution we need to extremize the action w.r.t λ i , such that we get: Now we need the method of resolvents to derive the expression for the partition function (Z(β)) in the present context. In continum limit of eigenvalues we can use density of states (eigen values) ρ(λ), which gives the number of eigen values lying in between λ and λ + dλ. Therefore, saddle point of V (λ i ) is given by the following expression: This formalism is very useful when we can't able specify the particle interaction in the effective action. More precisely, in this situation when we really don't have any information about the particle interaction one can't able to define the action in terms of the usual language. Additionally it is important to note that, in our computation we consider that gravitational background is classical and non dynamical. However it will not explicitly appearing in the action for the distribution of eigen values of random matrices. Also during reheating since one can neglect the contribution from the expansion of our universe, then considering only the representative action for random distribution is sufficient enough for our discussion when we don't have any knowledge about the particle interactions at the level of action. In such a situation gravitational background is treated to be not evolving with time during reheating.
Here Pr represents the principal part of the integral. Solution of the principal part of the integral Eq (5.11) gives the eigen value density ρ(u) at large N limit. Now, we can define resolvent as given by: further, using Eq (5.12) we compute the following function: Next, we use the following resolvent identities for our computation performed in this paper: Here R denotes the resolvent and A, B both defined over same linear space. Consequently, Eq (5.13) can be recast into the following simplified form: Here we define: which implies that, here ω can be expressed as: Finally, in terms of newly defined function Ψ as stated in Eq (5.17), one can further recast Eq (5.16) as: Further, comparing the two equivalent definition of ω(x) we get the following differential equation for Ψ in terms of the eigen values of the random matrix, as given by: Therefore , the solution for Ψ(x) is given by the following characterestic polynomial : Here it is important to note that, the solution obtained in large N limit can be compared with the solution obtained using WKB approximation in Schrödinger equation. Then we can neglect the term 1 N ω (x) in Eq 5.16 and write down the following approximated algebraic equation of ω(x), given by: where we have introduced two new quantities ω(x) and ρ(x), which are defined as: Then solution of ω(x) is given by the following expression: Here for our discussion ω + (x) is redundant and only acceptable solution for our purpose is given by the following expression: Additionally, it important to mention that in large N limit we can write, is the density of eigen values from Wigner's semi-circle law. Consequently, the solution obtained in Eq (5.26) can be recast in the following simplified form in the large N limit as:ω This implies that, just by knowing the even polynomial structure of the potential V (x) one can able to find out the solution for the distribution of ω(x) in terms of the random variable x. In this context, which further implies that the Wigner's semicircle law is defined as the probability density function of eigen values of many random matrices is a semi-circle as N → ∞. On the other hand, for finite N , Schrödinger equation gives the corrections comparing with calculated result obtained in Eq (5.27), which is given by the following expression: scaling factor 1 2a .The semicircle nature predicted from Eq. (5.146). Consequently, one can write: where, L is the Lagrange multiplier and 1 denotes the total density. Now, we can generalize it to normal matrix model whose eigen value belongs to V i (union of contours). To characterize this here we introduce filling functions, which are described by the symbol i and consider the contours as: where d =dimension and n i eigen values are integrated over γ i . Further, we define Consequently, from Eq (5.30) one can write: which will be helpful for further computation. Now, for a contour, which is represented by: Consequently, Eq (5.29) can be recast into the following simplified form: Now the Fourier transform of the density function ρ(x) can be written as: using which the second term of Eq (5.37) can be written in Fourier space as: Now we know that the saddle points can be computed by imposing the following condition: During this computation one can further define the effective random potential, which is given by the following expression: Then one can recast Eq (5.37) in terms of the effective potential as: Further imposing the saddle point condition we get: which can be further written in terms of the eigen values of the random matrices as: Therefore within supp of ρρ one can write: On the other hand outside the supp of ρ since ρ(x) → 0, then in the large N limit one can write:ω Therefore jump (discontinuity) on real line along the support ρ(x) is given by the following expression: Then using Eq (5.48), we get the following simplified expression for the jump (discontinuity) (5.49) Now one can introduce a new function P (x) of random variable x as: which is analytic on C as it gives zero value of the jump. This is explicitly shown in the following: Additionally, it is important to note that, using the previous results we get: Here the most general solution for the density function is given by the following expression 22 : where both M (λ) and σ(λ) are polynomial in λ are defined as: Here we consider n number of intervals on which ρ(λ) is supported and a 2i−1 and a 2i are the end point. Further we consider a general case where instead of the specific form of the mass profile we only know the polynomial structure of interaction random potential V (M ) which is characterized in terms of the random matrix M . For our purpose we take it to be even polynomial potential written in the following general form: For n = 1 we also get the following simplified expressions for the polynomial M (λ) and σ(λ): Consequently, for n = 1 we get the following expression for the density function on semicircle: Now we use this ρ(λ) in ω(λ + i0) and Taylor expand in the limit λ → ∞ we get: which implies that all coefficients of λ r for r > 0 is zero and this gives n number of equations. This finally gives the full equation of M (λ) in terms of the coefficients C 2i . Solving these equations we get: Further equating the coefficients on both sides of the Eq (5.65) we get: and it will continue upto term by term giving all a n and we get the unique polynomial M (λ). We will verify this generalization for n = 1, 2, 3, 4, 5 and check their SFF in this work accordingly. For more general discussions see ref. [96,97] also.

OTOC in Random Matrix Theory (RMT)
In earlier section we have introduced OTOC and its application to cosmology. In this subsection, we will discuss about OTOC appearing in the context of RMT.

Two point OTOC
For this purpose, we start with two point correlation functions for the GUE which is described by the following equation: where the operator O 2 (τ ) in Heisenberg picture can be expressed as: Here it is important to note that the GUE measure dH is represented by the Hamiltonian H., which is invariant under the following unitary conjugation operation, which is described by Here U is the unitary matrix. Consequently, the GUE two point correlation function can be further expressed as: where dU is the Haar measure appearing in this context. After integrating over the Haar measure we get the following expression for the GUE two point correlation function: where the connected two point correlation function O 1 O 2 is defined as: Now we consider a special case where O 1 and O 2 are described Pauli operators. In such a situation, the GUE two point correlation function can be expressed as: Here I represents the 2 n dimensional Hilbert space in the present computation. To derive this above mentioned expression we have not assumed any additional assumption expect the fact that the Haar measure of GUE dH is invariant. This is a very useful information to study the physical characteristics of chaotic Hamiltonian at macroscopic scales.

Four point OTOC
Now we discuss about the four point OTOC for the GUE prescription. Here the fourth point OTOC can be expressed in terms of fourth Haar moment: where we consider (4!) 2 = 576 terms in this expression for four point OTOC. Now we consider a special situation, where all these operators appearing in the expression for the four point OTOC for GUE are described by Pauli operators. In such a case, the four point OTOC for GUE can be simplified as: where SFF 4 (τ ) is the four point SFF for GUE, which is defined by the following expression: and this is derived only by considering the leading order behaviour of four point SFF. Here additionally it is important to note that if we fix: This will give rise to non-zero expression for the four point OTOC for GUE. For other situations, where we get zero contribution to the four point OTOC for GUE. One can further generalise this statement for any arbitrary 2p point OTOC for GUE, which is given by the following expression: Generalizing the previous argument one can conclude that the final result for the 2p point OTOC for GUE is non zero when we have the following constraint: Here one can further show that for the GUE we get: 84) which indirectly implies that GUE is not sensitive to the fact that the operators as appearing in this context are out-of-time ordered or something else. Additionally, it is important to note that, if we compute the expression for the OTOC correlation function for a specified class of Hamiltonian operators, which are in general invariant under the operation of conjugation on the unitary matrix U . In such a situation from OTOC one can further express the OTOC in terms of SFF. This is a very well known technique in the study of many -body QFT systems, where particularly to study the underlying physics of thermalization and quantum quench []. In the next subsection we will provide an analytical proof of the equivalence of the two point SFF and the two point OTOC, which can be further generalized to any arbitrary 2p point correlation functions.

Spectral Form Factor (SFF) from OTOC
From the traditional perspective the idea of quantum chaos is used in the context of study of spectral aspects of statistical field theory. Recent developments are made in the context of black hole theory and quantum information theory where using OTOC one can quantify quantum chaos. However in this paper our one of the prime objective to apply the concept of quantum chaos to study early universe cosmology, which is obviously another new direction of future research area. In this subsection, our air is to give a formal proof which establish the connection between Spectral Form Factor (SFF) and OTOC in OEQFT. First of all we consider a limit where β = 1/T = 0 in which distribution of quantum operator insertions around a thermal circular path is very straightforward.
Let us consider a quantum mechanical Hamiltonian operator H operating on an I = 2 n dimensional Hilbert space and consists of n number of quantum bits (qbits). Next, we consider the two point correlation function O(0)O † (τ ) using which one define the following averaged two point correlation function: Here we assume that O is the Unitary operator which is integrated over a Haar measure on U(2 n ). Also it is important to note that the integral over the Haar measure can be translated in terms of the Pauli operators O k and I 2 = 2 2n = 4 n represents the total number of Pauli operators for this quantum n qubit system. Further, it is important to note that, to derive the expression for SFF from the present context additionally we need the first moment of the Haar ensemble, which is defined as: which can be be equivalently expressed in terms of the language of Pauli operator as: Next using Eq (5.87) in Eq (5.85), we get the following simplified result: where the two point SFF at infinite temperature is defined in terms of the quantum Hamiltonian H as: Here the result obtained in Eq (5.88) implies that the quantum averaged OTOC is proportional to the two point SFF at infinite temperature of the present context. This prescription can be further generalised to make the connection between any arbitrary 2p point quantum OTOC and 2p point SFF in this context. To establish this connection let us consider a 2p point quantum OTOC, which is described by: Now taking the average over such 2p point OTOC we get: where Q p is defined as: Here one can consider a special case where Consequently, the average over such 2p point OTOC can be further simplified to the following form: (5.94) Here it is important to note that the terms appearing in the are not symmetric because is an non-Hermitian quantum operator. This result establishes a direct connection between the spectral physics in statistical field theory and other physical observables. Apart from theoretical perspective one can use two point SFF as a good estimator for experimental measure. For this purpose one can consider the following standard deviation (or experimental estimation error) of the unitary operator O given by: By choosing the Haar unitary operator O as a random Clifford operator one can find a good estimator of two point SFF.
To give the similar proof at finite temperature let us consider the energy eigenvalue representation of OTOC, which is given by the following expression: where the time dependent expansion coefficient can be expressed as: Here we have used the fact that, H|n = E n |n . Consequently we get: This establishes the connection between OTOC and two point SFF at finite temperature

Two point SFF and thermal Green's function in RMT
In this subsection our prime objective is to explicitly compute the expression for SFF for different even polynomial potential of random matrices. This is very useful to quantify chaos when we have no information about the interaction or time dependent effective mass profile which will finally give rise to scattering in conduction wire in presence of impurity or cosmological particle creation during reheating.
Let us now consider a Thermofield Double State (TDS) associated with canonical quantum mechanical state at finite temperature. The time evolution of the TDS can be expressed as: Using this information one can define Spectral Form Factor (SFF) as: Here E n and E m correspond to the n -th and m -th level of the quantum system under consideration. Here the Boltzmann factor β = 1/T , where T is the temperature associated to the system. Apart from temperature dependent Boltzmann factor the definition of SFF also involves involves conformal time τ , which we have define in earlier section of this paper and during reheating τ ∝ t. Here t is the physical time scale and the proportionality factor is constant in space time. Now at very high temperature (β = 1/T → 0) and low temperature (β = 1/T → ∞) we get the following limiting behaviour of SFF, as given by: It is also observed that in τ → ∞ limiting situation the nearest neighbour energy spacings contribute only to the quantification of SFF. This implies that the concept of SFF also helps in understanding the time dynamics of the quantum system under consideration and also very useful tool to analyze the discreteness in energy spectrum. Chaotic system satisfy Wigner's formula which makes SFF a good observable for quantifying chaos.
In usual prescriptions, SFF is averaged over an statistical ensemble of random matrix. This is a very particular feature of SFF and can be directly linked to the quantification of quantum chaos. Before going to discuss further here it is important to note that, all distribution for eigenvalues are different from each other but quite similar at small scales. This is a very crucial information for the computation of SFF to quantify chaos. Now in the present context we define a new function G(β, τ ), which is represented by the following expression: Here, D(λ) = ρ(λ) =eigen value density. In the present context, G(β, τ ) characterize the two point correlation function which measures SFF. Now, one can divide the total Green's function G in two parts (connected and disconnected part of the Green's function), as given by: where disconnected part of the Green's function G dc and connected part of the Green's function G c can be expressed as: Now, for further analysis we consider the high temperature limit (β = 1/T → 0) and also can divide the total Green's function G in two parts (connected and disconnected part of the Green's function), as given by: where disconnected part of the Green's function G dc and connected part of the Green's function G c can be expressed as: Here we define the connected two-point correlation function, which is given by the following expression: To quantify this explicitly one can define the eigen value density function D(λ) in the neighbourhood of extremum of level density (ρ(λ)) as: where D(λ) is the average of the eigen value density function over the statistical ensemble of eigen values of the random matrices and δD(λ) represents the quantum fluctuation on D(λ).
Consequently, using this fact the two point correlation function reduced to the following form: and using this connected part of the Green's function G c can be further simplified as: Additionally, it is important to note that, the mean level density can be normalised in a semi circle using the following two conditions: Here D(λ) actually represents the number of eigen values lying between the small interval (λ, λ + dλ) and in the present context it is proportional to O( √ N ). On the other hand, ρ(λ) is the density which we get by extremising the action and treated to be free from all factor of N and all eigen values which are just O(1). In this context, the two variables λ and σ are related by the follwing expression: To compute the G dc and G c part of SFF explicitly let us first start with the one point function on the semi-circle as given by: On the other hand at very low temperature limit (β = 1/T → ∞) we get: simplified as:: Here it is important to note that, for different polynomial random potential we will get different expressions for the integral measure. Now we need to find the specific point after which properties of SFF drastically changes. We define this points as critical points. For general even order polynomial potential one can write down the following expression for the density function of the eigenvalues of the random matrices: n k=1 a n−k λ 2(n−k) ∀ even n (5.119) Further substituting Eq (5.119) in Eq (5.116) we get the following simplified expression for the one point function on the semi-circle: n k=1 a n−k λ 2(n−k) ∀ even n = n k=1 a n−k a 2 −a 2 −2k 4 n−k e 2iπk + e 2iπn a 2(k+n) Γ −k + n + 1 2 ; a 2 (β ± iτ ) 2 ∀ even n. Repeating the same calculation in high temperature (β = 1/T → 0) limit we get: ; −a 2 τ 2 ∀ even n. (5.121) For different polynomial potentials we can actually calculate the expansion coefficients a n−k and get the exact form of Z(β ± iτ ).
At finite temperature the disconnected part of the Green's function (G dc (β, τ )) can be expressed as: ; a 2 (β + iτ ) 2 × n m=1 a n−m −a 2 −2m 4 −m e 2iπm + e 2iπn a 2(m+n) Γ −m + n + 1 2 Further taking high temperature limit we get the following simplified expression for SFF as given by: ; −a 2 τ 2 (5.123) Next, we will consider late time limiting behaviour of the one point function, which can be expressed as: lim and at the high temperature (β = 1/T → 0) limit we get: Now, it is important to note from the previous discussion on SFF that, the connected part of the Green's function G c part of SFF depends on the two point correlation function δD(λ)δD(µ) and from RMT the exact from of this two-point function near the centre of spectrum (mean) of the eigen values is known and can be expressed in the following form: which can be derived using the method of orthogonal polynomials for Gaussian ensembles. This is true for any polynomial potential measure whose matrix (operator) is of single trace. Various polynomial potentials change only the eigen value distribution near edges of the distribution. There are two parts and they give different measures: 1. 1/N 2 part with sine squared function gives the ramp and have subdominant contribution.
2. 1/N part with Delta function gives the plateau and dominant.
Next, using Eq (5.126) in Eq (5.112) we get the following simplified expression for the connected part of the Green's function G c as given by: where we have used the fact that: To perform the integral present in the expression for G c we further substitute, λ + µ = E, λ − µ = ω. Consequently, the measure can be expressed as, dλ dµ = dE dω. Then at high temperature using this fact Eq (5.127) can be recast as: Then, at finite temperature the connected part of the Green's function can be written as: where δ(β) is the Dirac Delta Function, which is defined as: Since the integral over E gives trivial Diarc Delta function we choose our working region for which E = 0 (at hight temperature limit). Then the remaining integrand is only over ω and it finally gives: which gives us finally the following simplified expression: From the obtained result it is clearly observed that we get the linear growth in the region τ < 2πN and the constant plateau type behaviour in the region τ > 2πN . Also it is important to note that change in behaviour from region τ < 2πN to region τ > 2πN is abrupt. To show the behaviour of SFF explicitly we define argument of sin function as: as we choose N → ∞ and ω → 0. In this limiting situation we get the following results: 1. For large x(>> 1), sin x x → 0 and only the Dirac Delta function remains intact. So in this specific limit the vanishing of sin term implies that the oscillatory fluctuations don't contribute in the final expression for SFF. This limiting situation is called spectral rigidity.
2. For small x(<< 1), sin x x → 1. In this limiting situation the integral gets maximum contribution from the ω = 0 region. And this part contributes in ramp region.
We can also measure dip-time and it will give the change of decay behaviour exactly at the critical point. A direct relation between fall-off behaviour of the SFF and the edge behaviour of level density, at critical points can be established using Paley-Wiener Theorem [59]. Now we consider a function g(ζ) which is defined on a compact spatial support and its Fourier transform F (η) has a lower bound on the rate of decay is given by the following expression: Here N is a rational number and γ N is a real constant. A direct relation between the decay of the SFF and the edge effect of mean level density is given by the following expression: For the proof of this statement see ref. [59]. For decay behaviour of SFF at late time we use asymptotic behaivour of the solution appearing in the ref. [98]. Now to compute SFF we need to add both connected and disconnected part of the Green's function G(= G c + G dc ). Therefore, for different even polynomial potential we get finally the following expression for SFF at finite temperature: 137) where SFF(τ ) is defined with proper normalization.
After substituting the expression for G dc (β, τ ) we get the following expression for the SFF at finite temperature: ; a 2 (β + iτ ) 2 × n m=1 a n−m −a 2 −2m 4 −m e 2iπm + e 2iπn a 2(m+n) Γ −m + n + 1 2 Further taking the high temperature limit we get the following simplified expression for SFF -68 -as given by: Let us start our discussion with Gaussian random potential given by: Now for a single interval (n = 1) with end points −2a and 2a (semi-circle) we get: and we get the following expression for density function for eigen value of the random matrix M as given by: Further, Taylor expanding ω(λ + i0) we get the following expression: Then comparing the both the sides of above expression we get: Then the density function in terms of the eigen value of random matrix M is given by the following expression: and one point function of the partition function in presence of the Gaussian random potential can be expressed as: where 0F1 (A; B) is the regularized Hypergeometric function. Now here substituting τ = 0 we get: Further taking high temperature limit we get the following simplified expression for the one point function: which can be further simplified by taking the limit T = √ N τ → ∞ as: Now for the quadratic random potential disconnected part of the Green's function can be computed at finite temperature as: Further taking the limit T = √ N τ → ∞ we get the following simplified result: Now to compute SFF we need to add both connected and disconnected part of the Green's function G(= G c + G dc ). Therefore, for quadratic polynomial potential we get finally the following expression for SFF at finite temp: where SFF(β, τ ) is defined with proper normalization and in our prescription it gives the total Green's function as mentioned above. Further simplifying the result for high temperature limit we get the following expression for SFF, as given by: , (5.155) Further taking the limit T = √ N τ → ∞ we get the following simplified result for SFF:  From fig. 23(a), fig. 23(b), fig. 23(c) and fig. 23(d) we see that SFF at finite temperature decays with increasing τ and reach zero. But with changing β, SFF values remains almost same initially (for higher β or lower temperature).
For both the plots we have shown that SFF decays to zero for finite temperature. In

For Quartic random potential
Here we consider quartic random potential which can be written as: For a single interval (n = 1) with end points -2a and 2a (semi-circle) we get the following expression for density function for eigen value of the random matrix M as given by: Now, for the quartic random potential ω(λ + i0) can be expressed as: Now Taylor expanding ω(λ + i0) near λ → ∞ gives the following expression: Therefore equating both sides of the above equation gives: Then the density function in terms of the eigen value of random matrix M is given by the following expression: Further solving the constraint we get: and Here a 2 has imaginary value for g − 1 48 and the critical value is given by:  In fig. 25(a) and fig. 25(b) density function ρ(λ) for quartic potential is plotted with a = 1. The curve follows from Eq. (5.164). When g = 0 it matches with Wigner law. For g > 0 the curve shows a plateau region whereas for g < 0 it preserve the semicircular nature with minor deviation.The plateau region denotes the deviation from Wigner law even at very less effect of quartic term (as g is chosen to be small). The plateau region though converge with semicircle at end point. At g c = − 1 48 the curve deviates but converge to semicircle at end points where as for g < g c the curve never converge to semicircle one supporting its non-existence. (See eq. (5.166) for details.). Now we will calculate the one point function of the partition function for quartic random potential, which is given by the following expression: where I n (x) is the modified Bessel function of first kind with order n.
Further taking the high temperature limit we get the following simplified expression for the one point function: Therefore the first term vanishes exactly at the critical point g c = − 1 48 which gives: Now taking the limit T = √ N t → ∞ we get finally the following simplified result for the one point function: Now for the quartic random potential disconnected part of the Green's function can be computed at finite temperature as: which can be further simplified in the high temperature limiting situation as: Further taking the limit T = √ N τ → ∞ we get the following simplified result: Now to compute SFF we need to add both connected and disconnected part of the Green's function G(= G c + G dc ). Therefore, for quartic polynomial potential we get finally the following expression for SFF at finite temp: 174) where SFF(β, τ ) is defined with proper normalization and in our prescription it gives the total Green's function as mentioned above. Further simplifying the result for high temperature limit we get the following expression for SFF, as given by: Further taking the limit T = √ N τ → ∞ we get the following simplified result for SFF: Further simplifying the result for high temperature limit we get the following expression for SFF, as given by: , (5.177) From fig. 26(a) and fig. 26(b), we see that SFF at finite temperature decays with increasing τ and reach zero. But with changing β SFF values remains almost same initially (for higher β or lower value of temperature).
From both the figures we have shown that SFF decays to zero for finite temperature. In

For Sextic random potential
In this subsection we consider sextic random potential, as given by the following expression: For a single interval (n = 1) with end points -2a and 2a (semi-circle) we get the following expression for the density function in terms of the eigen value of random matrix M : Also for sextic potential ω(λ + i0) can be expressed as: In fig. 28(a) and fig. 28(b) for sextic potential behaviour of density function ρ(λ) is shown. The curve follows from Eq. (5.185). Again choosing g = h = 0 will produce the Wigner law. Deviating g and h by small amount shows deviation from Wigner semicircle law. For g > 0, h > 0 the curve shows plateau region though merge with semicircle at end points.   Further expanding ω(λ + i0) near λ → ∞ we get: Therefore, equating both the sides of the above equation we get: along with we get one additional constraint condition, as given by: Then, for the sextic random potential we get the following simplified expression for the density function in terms of the eigen value of the random matrix M , as given by: Solving the constraint condition we get, a 2 in terms of g and h. The real root for a 2 is given by the following expression: where we define the function F(g, h) as: Here we can check that putting h = 0 the constraint condition reduces to the following simplified form: 12ga 4 + a 2 = 1 (5.188) and the solution of this equation is given by the following expression: Here the critical value with h = 0 is given by: which is exactly same result as obtained for quartic potential in the previous subsection. Now the expression for the one point function for partition function at finite temperature can be computed as: Further in the high temperature limit the one point function for partition function can be simplified as: Next, simplifying the result for one point function in the limit T = √ N τ → ∞ we get: Now for the quadratic random potential disconnected part of the Green's function can be computed at finite temperature as: which can be further simplified in the high temperature limiting situation as: Further taking the limit T = √ N τ → ∞ we get the following simplified result: Now to compute SFF we need to add both connected and disconnected part of the Green's function G(= G c + G dc ). Therefore, for sextic polynomial potential we get finally the following expression for SFF at finite temp: −2 [(β + iτ )I 1 (2a(β + iτ )) (360a 2 h + β 2 (180a 4 h + 24a 2 g + 1) + 2iβτ (180a 4 h + 24a 2 g + 1) −τ 2 (180a 4 h + 24a 2 g + 1)) − 24aI 2 (2a(β + iτ )) (30h + (β + iτ ) 2 (15a 2 h + g))] where SFF(β, τ ) is defined with proper normalization and in our prescription it gives the total Green's function as mentioned above.
Further simplifying the result for high temperature limit we get the following expression for SFF, as given by: Further taking the limit T = √ N τ → ∞ we get the following simplified result for SFF: fig. 29(a) and fig. 29(b), we see that SFF at finite temperature decays with increasing τ and reach zero. But with changing β SFF values remains almost same initially (for higher β). In fig. 30(a), fig. 30(b), fig. 30(c), it is observed that SFF with variation in N get saturated at different value of τ . But with increasing N the value of the saturation point, will decrease. Subtracting the change of axis[SF F | τ =0 ] we get the predicted bound of SFF.

For Octa random potential
Here we consider octa random potential, as given by the following expression:  Then the function ω(λ + i0) can be expressed as: Further Taylor expanding ω(λ + i0) near λ → ∞ we get: Therefore, equating both the sides of the above equation we get: along with an additional constraint condition: a 2 + 12a 4 g + 60a 6 h + 280a 8 k = 1 (5.208) Solution of this constraint equation gives a 2 in terms of g, h and k. Since the solutions for a 2 are very complicated, we have not explicitly mentioned them here. Instead of writing full solution here we can check that putting h = 0 the constraint condition reduces to the following simplified form: 12ga 4 + a 2 = 1 (5.209) and the solution of this equation is given by the following expression: Here the critical value with h = 0 and k = 0 is given by the following expression: which is exactly same result as obtained for quartic and sextic (with h = 0) potential in the previous subsections. Then, the final expression for the density function in terms of the eigen value of the random matrix M can be written as: In fig. 31(a) and fig. 31(b) for octic potential behaviour of density function ρ(λ) is shown. The curve follows from Eq. (5.212). Again choosing g = h = k = 0 will produce the Wigner law. Deviating g, h and k by small amount shows deviation from Wigner semicircle law. For g = 0, h > 0, k > 0 the curve shows plateau region though merge with semicircle at end points. But choosing g > 0, h < 0, k < 0 and g > 0, h > 0, k < 0 show deviation from semicircle and don't converge even at end points. On the other hand, if we choose g > 0, h > 0, k > 0 then we get a valley region lying between two peaks of the maxima of the density distribution of eigen values of the random matrices under consideration. The same behaviour can be obtained by fixing g > 0, h < 0, k > 0, g = h = k = 1 and g = 0, h >> 0, k >> 0. Only slight difference can be visualised in the peak heights of the maxima and also in the spread in the valley region. But in all such cases in between it will not at all match with the Wigner semicircle law, but converge to the end points of the Wigner semi-circle , which is obtained by setting g = h = k = 0. Next, we compute the expression for the one point function of the partition function at finite temperature, which can be expressed as: where I n (x) is the modified Bessel function of first kind with n th order. Further, considering the high temperature limiting situation we get the following simplified expression for the one point function of the partition function: dλ √ 4a 2 − λ 2 80a 6 k + 6a 4 3h + 4kλ 2 + a 2 4g + 6hλ 2 + 8kλ 4 + 2gλ 2 + 3hλ 4 + 4kλ 6 + 1 2 e ∓iτ λ = a τ 6 ±J 1 (±2aτ )(1120a 6 kτ 5 + 60a 4 τ 3 (3hτ 2 − 112k) +24a 2 (gτ 5 − 15hτ 3 + 840kτ ) + τ 5 ) −24J 2 (±2aτ )(−30τ 2 (28a 2 k + h) + τ 4 (140a 4 k + 15a 2 h + g) + 1680k) . (5.214) Next, simplifying the result for one point function in the limit T = √ N τ → ∞ we get: Now for the octic random potential disconnected part of the Green's function can be computed at finite temperature as: 15a 2 hτ 4 + β 4 140a 4 k + 15a 2 h + g + 4iβ 3 τ 140a 4 k + 15a 2 h + g − 6β 2 140a 4 kτ 2 + 5h 3a 2 τ 2 − 1 − 140a 2 k + gτ 2 − 4iβτ 140a 4 kτ 2 + 15h a 2 τ 2 − 1 − 420a 2 k + gτ 2 + 140k a 4 τ 4 − 6a 2 τ 2 + 12 + gτ 4 − 30hτ 2 + |a| (β + iτ ) 3 I 1 (2(β + iτ ) |a|) −1120a 6 kτ 2 + 60a 4 112k − 3hτ 2 + 24a 2 15h − gτ 2 + β 2 1120a 6 k + 180a 4 h + 24a 2 g + 1 + 2iβτ 1120a 6 k + 180a 4 h + 24a 2 g + 1 − τ 2 which can be further simplified in the high temperature limiting situation as: Further taking the limit T = √ N τ → ∞ we get the following simplified result: (5.218) Now to compute SFF we need to add both connected and disconnected part of the Green's function G(= G c +G dc ). Therefore, for octic polynomial potential we get finally the following expression for SFF at finite temp: SFF(β, τ ) ≡ β 12 (β 2 + τ 2 ) 6 −24a 2 I 2 (2β |a|) β 4 140a 4 k + 15a 2 h + g + 6β 2 5h + 140a 2 k + 1680k + |a| β 3 I 1 (2(β) |a|) (6720a 4 k + 360a 2 h + β 2 (1120a 6 k + 180a 4 h + 24a 2 g + 1)) where SFF(β, τ ) is defined with proper normalization and in our prescription it gives the total Green's function as mentioned above. Further simplifying the result for high temperature limit we get the following expression for SFF, as given by: (5.220) Further taking the limit T = √ N τ → ∞ we get the following simplified result for SFF:  From these plots we can say that time variation of SFF follow oscillatory pattern initially but after certain time it has linear decaying amplitude for dominance of linear part. Then after τ > 2πN region SFF abruptly saturated due to second part of the connected part of the total Green's function G c . On the other hand, for τ < 2πN region SFF is decaying in amplitude and increasing with time. After τ > 2πN region the function will be constant thereafter.
Here it is important to note that, depending on the specific structure of the even polynomial random potential the upper bound on chaos very slightly changes (i.e. the amplitude for saturation of SFF is almost at the same order of magnitude for different even polynomial random potentials). But the late time behaviour for different random potentials are almost same as it shows complete saturation with respect to time. The saturation depends only on value of N . Also it is import to note from the plots that, for each even polynomial potential sudden transition from the random oscillatory behaviour to the perfect saturation of SFF take place at the unique time, τ = 2πN .

Estimation of dip-time scale from SFF
Now we introduce the concept of dip-time which denotes the change in fall-off behaviour of SFF near the critical points. It is estimated by comparing the initial fall-off behaviour with late time behaviour of the curve from which it starts the linear increase (ramp part).
For different even polynomial random model (see Eq (5.150 ) , Eq (5.170) , Eq (5.193), Eq (5.215)), we see that the fall off behaviour varies with τ − N 2 . Consequently, the disconnected part of the Green's function (G dc ) fall off as τ −N and we tabulated different fall off behaviour with linear increase rate in table. 3. Physical time t, conformal time τ and newly define time scale are dined as 23 : which specifically depends on different values of N , where it represents order of the even polynomial used in our paper to compute SFF. 23 Here to define the new time scale T we assume that the reheating approximation is perfectly valid. This implies that we can really neglect the expansion of the universe. This further pointing towards the fact the conformal time (τ ) and the physical time (t) are related by the following expression: where during reheating we have assumed that the conformal and physical time are almost same and the proportionality constant is the inverse of the scale factor a −1 , which is independent of both the time scales discussed in this paper. .

Potential 1st critical point 2nd Critical point 3rd Critical point 4th critical point
Gaussian Table 3. Fall-off behaviour near critical points for different even polynomial random potential. In the next section we will discuss about quantum correction of Fokker-Planck equation in the context of cosmological particle production. Here it is important to note that, particle production in cosmology can be treated as a chaotic random event and through our calculation we get quantum corrections due to the non-Gauaaian contribution in the probability distribution function. In this context the system can be treated as semiclassical. As a result, SFF shows a saturating behaviour on large time limit which implies that randomness in eigen value density has a upper bound though it is chaotic [59]. This relate that particle production also can have an upper bound which also confirmed using the computation of Lyapunov exponent.

5.6
Universal bound on quantum chaos from SFF and its application to cosmology In the previous subsection we have explicitly computed the analytical expression for SFF for generalized even polynomial random potential at finite temperature (β = 1/T =finite) in Eq:. (5.138) and at very high temperature (β = 1/T → 0) in Eq. (5.139). Now in this subsection our prime objective is to compute the analytical bound on SFF at long time interval i.e. τ → ∞. To derive the bound on SFF we first use the asymptotic behaviour of Hypergeomteric PFQ regularized function, which is given below: , −m + n + 5 2 ; a 2 (β ± iτ ) 2 ) = 0 ∀k = 1, 2, · · · , n, (5.224) This asymptotic behaviour of the Hypergeomteric PFQ regularized function remains same in the high temperature limit (β = 1/T → 0) also. Consequently, the asymptotic behaviour of the disconnected part of the Green's function can be expressed at finite temperature as well as in the limit β → 0 with finite N as: Similarly, the asymptotic behaviour of the connected part of the Green's function can be expressed at finite temperature as well as in the limit β → 0 with finite N as: Finally, adding the contribution from the disconnected and connected part of the Green's function in the asymptotic limit (τ → ∞) we get the following simplified expression for SFF at finite N as given by: (5.230) Also in the high temperature limit with finite N we get: (5.231) Here we considered the part of SFF only after τ > 2πN with finite N as we are considering τ → ∞ asymptotic limit. The main obstruction of the taking τ → ∞ asymptotic limit in the τ < 2πN with finite N divergent contribution in the connected part of the total Green's function G c as given by: (5.232) On the other hand, for the disconnected part of the Green's function we get the same result as obtained for τ (→ ∞) > 2πN with finite N case.
As a result, it gives divergent contribution to SFF at finite N is given by: For this reason, we will concentrate on the finite contribution on SFF coming from τ (→ ∞) > 2πN with finite N region. Finally, adding both the contribution from connected and disconnected part of the total Green's function for the asymptotic region τ (→ ∞) > 2πN (F inite N ) with 0 ≤ β(= 1/T ) ≤ ∞ we get the following upper and lower bound on SFF, as given by: (5.234) On the other hand, with large N limit one can consider the τ < 2πN region for the computation of the bound on SFF. To justify this statement we take the τ → ∞ asymptotic limit in the τ < 2πN with large N gives finite contribution in the connected part of the total Green's function G c as given by: (5.235) Similarly, for the disconnected part of the Green's function we get the same result as obtained for τ (→ ∞) > 2πN in previous case.
Finally, adding both the contribution from connected and disconnected part of the total Green's function for the asymptotic region τ (→ ∞) < 2πN (Large N ) with 0 ≤ β(= 1/T ) ≤ ∞ we get the following upper and lower bound on SFF, as given by: Bound on SFF from theory : behaviour of SFF with temperature and time for large and small N. At higher τ (Fig. ??),?? SFF decays with τ and at last goes to 1 πN . Else (Fig:-35(a),35(b),35(c),35(d)) SFF increases with τ for τ < 2πN and at last saturate to 1 πN . Now from the analytical solution we know that for large τ or large β G dc doesn't contribute to the SFF as the G c has τ (2πN ) 2 term.But for higher N and low τ there should be a minima( 1 πN − 1 N ) for this function within the range τ < 2πN . As soon as the point τ = 2πN is crossed the function change its form and get saturated. This way of calculation of bound conclude that nature of SFF at late τ shows same nature independent of type of potential. For infinite temperature SFF saturate at same level and at same τ value-for different potential with same N . For finite temperature it decays to zero irrespective of nature of potential. SFF is a measure of quantum chaos in a dynamical system. Bound on SFF prove that whatever be the interaction, every system at infinite temperature with same N saturated at same value. But at finite temperature randomness in the system decays to zero at late time and the system equilibriate within itself.
Here we consider Gaussian (Eq. Here we have Bessel's Function of first kind (I k (2aτ ) in which taking asymptotic time limit we get: Here n is order of the polynomial random potential. As a result for finite N and large N we get the bound on SFF s mentioned earlier.
Here it is important to note that, our prescribed bound on SFF is also the same as the saturation value of SFF for different potential for finite N and large N .
6 Randomness from higher order Fokker-Planck equation: A probabilistic treatment in cosmology 6.

Cosmological scattering problem
Here we discuss about the cosmological scattering problem due to the particle creation in the context of early universe physics (mostly during reheating epoch). For detailed derivation of the results see refs. [42,100], which we have followed in this discussion mostly. As we have already discussed in the first half of the paper that, the Klein-Gordon equation, which is the dynamical master equation of the particles created during reheating can be solved in the same way as Schrödinger problem by formulating it as scattering problem in presence of an impurity potential inside a conduction wire [101] and can be related to the phenomena of chaos [43,44,102,103]. In this section our prime objective is to establish this connection including the possible quantum effects (corrections) and we will try to develop a formalism to explain the quantum analogue of the chaos during cosmological particle creation.
To serve this purpose let us start with the solution of the Fourier mode of the field (created particle during reheating) after j-th non-adiabatic event, which can be expressed as: Fourier mode solution after j − th event : where β j and α j are the Bogoliulov coefficients. For the vacuum solution we set the initial condition as: Vacuum initial condition : β 0 = 0, α 0 = e iδ at τ = 0, (6.2) where δ is a phase factor. Now the vacuum initial condition implies: In this context the Bogoliulov coefficients satisfy following normalization condition: (6.4) This is analogous with scattering in presence of an impurity inside the conduction wire. Here we can relate Bogoliubov coefficients before and after the non-adiabatic event by transfer matrix as: When the wavelength of incoming mode is much larger than coherence interval of the nonadiabatic event, then the time dependent mass profile evolution m 2 (τ ), can't be resolved in wave. For this purpose, we take the following Dirac Delta profile of time dependent mass function: which is localized at time τ = τ j . Here j represents the number of non-adiabatic events and the total number of the events can be expressed as: Now, it is important to note that, these solutions will satisfy the following two fold junction conditions: In the present context, the transfer matrix M j can be expressed as: where λ j is defined as: . (6.11) Therefore in this computation the transmission and reflection co-efficient can be expressed as: Further, the local change in occupation number n j can be written in terms of transmission co-efficient as: Further, assuming the local change of occupation number is large only for k << m j i.e. λ j >> 1/2 then the transfer matrix can be simplified to the following polar form as: where the phase factors θ j and φ j are defined as: Further, using the transfer matrix in polar form we define the transmission, reflection probability and the total occupation number as: Additionally, it is important to note that, the total occupation number before (n(j)) and after (n(j − 1)) the j-th scattering are related by the following expression: n(j) = n(j − 1) + λ 2 j 1 + 2n(j − 1) + 2 n(j − 1)[1 + n(j − 1)] cos ∆ j + 2λ j n(j − 1)(1 + n(j − 1)) sin ∆ j (6.21) where ∆ j is the phase factor which is defined in terms of the Bogoliubov coefficients of the (j − 1) th events as: Now the vacuum initial condition demands that, n(0) = 0, (6.23) which further implies the following equation for n(−1): where we define A, B and C as: Using Eq (6.24) the solution for n(−1) can be written as: Similarly using Eq (6.21) recursively one can find out the expressions for occupation number for many scattering processes. Alternatively, using the concept of transfer matrices one can also compute the occupation number in the present context. To serve this purpose one can first write down the total transfer matrices for N s number of scatterer as: Using this the total occupation number can be expressed as: To model a phenomenological situation where width is finite, the scattering event is relevant and consider "sech" scatterers: Now, if we take the limit w j → ∞ then we get back the Dirac Delta mass profile as given by: Further, using the results obtained for transmission co-efficient we finally get the following simplified expression for the total occupation number: n j = cos 2 ( 1 2 π 1 + 2m j w j ) sinh 2 (πkw j ) . (6.33)

Fokker Planck Equation
In this subsection our prime objective is to construct Fokker Planck equation from the basic principles. To serve this purpose we start with the concept of probability density, which can be expressed in terms of Smoluchowski equation: Smoluchowski Equation : It actually explain the probability density for particle position of Brownian motion in a random system. For a Markovian process one can further express this in terms of Chapman-Kolmogorov equation where the probability density is conditional. It is important to note that, Smoluchowski equation describes a two point conditional probability distribution satisfying the following criteria: where we have added a small interval δτ to an existing interval τ to construct the probability density function. in such a situation the transfer matrix for the elongate interval can be expressed as: Here it is important to note that, to write this expression we have used the following set of rules: 1. First of all, we identify, M 1 = M τ and M 2 = M δτ .

Then we apply the composition law
3. Finally, we write M 1 = M τ as: This implies: Then, the time evolution of the probability density function can be expressed as: This gives Fokker Planck equation upto different order for occupation number n after appropriate parameterization of the transfer matrices and a marginalization over certain parameters.
Additionally it is important to note that, in order to reduce the complexity of the computation we suppress the wave number (k) dependence in our obtained results. As a consequence, the Smoluchowski equation as stated in Eq (6.1) can be simplified to the following form: Now, we Taylor expand both side after writing [n 1 , θ 1 , φ 1 ] in terms of [n, θ.φ].
Here one can write the occupation umber n 1 as: where the perturbed part of the occupation number δn can be expressed as: In this context, the right hand side of the above equation represents the perturbed Hamiltonian. Here we use perturbation theory to find the eigenvalues of occupation number n 1 in terms of the eigenvalues of occupation number n and the matrix elements of the perturbed part δn in the preferred choice of basis which diagonalizes the matrix n. Additionally, it is important to note that, the explicit expression for the perturbed angular parameter δθ is not very significant for our discussion. Instead of this, the angular difference (φ 2 − θ) play crucial role to quantify the perturbed contribution to the occupation number. Now, the conditional probability of getting Y at time t + τ in terms of probability of getting nearby to Y − ξ at time t and then to Y in time τ is given by: On the other hand, using Taylor series expansion of P 2 (Y 0 |Y, t + τ ) around τ = 0 we get: Here it is important to note that, in Taylor series expansion of the probability density function P 2 (Y 0 |Y, t + τ ) we truncate the series by considering upto the second term in the series. Further, comparing the right hand sides of Eq (6.45) and Eq (6.46), we finally get the second term of the Taylor expansion as given by: Next using Eq (6.45) in Eq (6.47) we get the following simplified result for the second term of the Taylor expansion as given by: (6.48) Now in this context the normalization condition for the probability density function can be written as: From Eq (6.48) we observe that it has scattering out and scattering in contributions respectively. Now, we can determine P 2 ≡ P (n 2 , θ 2 , φ 2 ; δτ ) by the condition that it maximizes the Shannon entropy: using the principles of maximum entropy ansatz. In the above expression, g 1 , g 2 and g 3 are the Lagrange multipliers. Here it is important to note that, U (θ 2 ) is an arbitrary function of θ 2 which has an extremum at the location θ 2 = 0 and can be explicitly determined by imposing additional constraint conditions i.e. symmetry arguments, consistency requirements and available knowledge of the microscopic sector of the system under consideration.
To apply the concept of maximum entropy ansatz we choose the following set of constraints, which are helpful to minimize Shannon entropy in the present computation:

Constraint I:
First of all, we talk about the Constraint I, which will fix the normalization condition of the probability density distribution as given by: This is obtained by setting the co-efficient of the Lagrange multiplier g 1 to zero.

Constraint II:
Using the Constraint II, it is possible to fix the local mean particle production rate , which is quantitatively defined as: This is obtained by setting the co-efficient of the Lagrange multiplier g 2 to zero.

Constraint III:
Finally the Constraint III demands that: 53) which basically implies that the addition of infinitesimal interval can't correspond to a finite significant change in transfer matrix.
To establish this statement we start with the following transfer matrix written in the polar form for j = 2: Now in the limit δτ → 0 we have: Consequently the transfer matrix can be simplified as: To impose this specific non-trivial constraint we assume that the following condition is satisfied: where U (θ 2 ) is a real valued and positive definite arbitrary function. This is possible if the function U (θ 2 ) has an extremum at θ 2 = 0 where e iθ 2 = 1. One can choose various types of function which can satisfy these constraints. For an example, as a phenomenological choice one can consider the following functional form: As a result, the probability density function reaches its maximum at θ 2 = 0 when the time interval δτ → 0. Further, extremizing the expression for the Shannon entropy we get the following expression for the probability density distribution function P (n 2 , θ 2 , φ 2 ; δτ ) as given by: where we introduce two new functions K(g 2 ) and K(g 3 ) which are defined as: dn 2 e −g 2 n 2 = 1 g 2 , (6.60) Further using Eq (6.52) and Eq (6.57), we get the following simplified expression for the probability density function: P (n 2 , θ 2 , φ 2 ; δτ ) as given by: Maximum Entropy Ansatz : = P (n 2 ; δτ )P (θ 2 ; δτ ) = P (n 2 , θ 2 ; δτ ), (6.62) which implies that the probability density function is independent of φ 2 after applying the maximum entropy ansatz. For weak scattering, this corresponds to scattering events being uniformly distributed. Now if we consider large number of scatterings, then applying Central Limit Theorem one can show that the final result is not sensitive to the probability density function P 2 . In this discussion we have explicitly provided the mathematical form of the probability density function, which is not very important to derive the Fokker-Planck equation.
Now if we use the fact that the probability density distribution function P 2 is completely independent of φ 2 one can further express the Smoluchowski equation in the following simplified form: P (n, θ, φ; τ + δτ ) ≡ P (n, θ; τ + δτ ) = P (n, θ, τ )P (dn + dn , dθ + dθ ; δτ )dn dθ = P (n + δn, θ + δθ; τ ) δτ (6.63) Further, integrating both sides of the above equation with respect to the parameter θ we get the following simplified expression: P (n; τ + δτ ) = dθ P (n, θ; τ + δτ ) = dθ P (n + δn, θ + δθ; τ ) δτ = P (n + δn; τ ) δτ (6.64) where during performing the integration over θ we explicitly use the information that the infinitesimal change in θ i.e. δθ is not functionally dependent on θ. Now, using Taylor expansion of P (n + δn; τ ) δτ with respect to the infinitesimal occupation number δn we get: ∂n q (δn) q δτ = P (n; τ ) δτ + ∂P (n; τ ) ∂n δn δτ + 1 2! ∂ 2 P (n; τ ) ∂n 2 (δn) 2 δτ + · · · = P (n; τ ) δτ + ∂P (n; τ ) ∂n ∂ δn δτ ∂τ + 1 2! ∂ 2 P (n; τ ) ∂n 2 (δn) 2 δτ δτ δτ + · · · , (6. 65) where in this context we have: P (n; τ ) δτ = P (n, τ ). (6.66) On the other hand, taking Taylor expansion of the probability density function P (n; τ + δτ ) with respect to the infinitesimal time interval δτ we get: Further, substituting Eq (6.65) and Eq (6.67) in Eq (6.64) and equating both the sides we get: Consequently, using Eq (6.44) one can define the following statistical moments: 2n(n + 1) n 2 + (1 + 6n + 6n 2 ) n 2 2 = 2n(n + 1)µδτ + (1 + 6n + 6n 2 )(µδτ ) 2 . (6.70) Here it is important to note that, for proper truncation of the moments we assume that the particle production rate is small locally i.e. µδτ < 1. For this reason the second factor is ignored in (δn) 2 δτ and finally we get: Diffusion term , (6.71) where in the mean particle production rate (defined earlier) we have restored the Fourier mode dependence i.e. µ = µ k . On the other hand, in the occupation number we have ignored the Fourier mode dependence. For more details see ref. [104]. Additionally, it is important to note that in presence of diffusion one can derive the Fokker-Planck equation from Langevin equation of the following mathematical form: where a(n) is defined in terms of an external deterministic force f (n) and contribution from frictional damping (which is characterised by the damping coefficient γ) as: and b(τ ) is the Gaussian random function (also known as Gaussian white noise), which satisfies the following criteria: In our prescription, the diffusion coefficient D(n) is given by the following expression: Diffusion Coefficient (Einstein s Relation) : Here it is important to note that, in a most generalised situation b(τ ) has a finite microscopic autocorrelation time for which it is a coloured noise defined as: where g(τ − τ ) is an arbitrary function of the time interval τ − τ . In such a situation it describes a non-Markovian process. One can also recast the Langevin equation in the following alternative form:

Langevin Equation II
: where, the Gaussian white noise satisfies the following criteria: Now we define: using which we finally get the following results: which is true for Q = 1. Now using the Chapman-Kolmogorov equation in the present context we get: Here upto the order we use the following fact: Also we have used the following well known identity of Dirac Delta function, as given by: Now, we expand the Dirac Delta function in the powers of , as given by: where it is important to note that we have Taylor expanded the Dirac Delta function of the order of B 2 , this is because of the reason that B ∼ O( √ ). Hence we use this result in Chapman-Kolmogorov equation and we get the following simplified integral: Further, performing the integration and using Eq (6.83) and Eq (6.84) we finally get the following simplified result of this integral: P (n; τ + |n 0 ) = P (n; τ |n 0 ) + ∂P (n; τ |n 0 ) ∂τ = P (n; τ |n 0 ) + − ∂ ∂n (a(n)P (n, τ |n 0 )) + ∂ 2 ∂n 2 (D(n)P (n, τ |n 0 )) + O( 2 ). ∂P (n; τ ) ∂τ = − ∂ ∂n (a(n)P (n; τ )) + ∂ 2 ∂n 2 (D(n)P (n; τ )) , (6.92) where we have used the notation P (n; τ |n 0 ) ≡ P (n; τ ) for simplicity.

Generalized Itô prescription:
In this situation for general Q one can write: where in the present situation the Stratonovich prescription corresponds to Q = 1/2. However, for our problem, we consider the simplest situation , where a(n) = 0, Q = 0 and D(n) = n(1 + n).

Now integrating the Langevin equation II over a small time interval we get:
n(τ + ) − n(τ ) = a(n(τ )) + τ + τ dτ D(n(τ ) b(τ ). (6.96) Now to deal with the product D(n(τ ) b(τ ) one can use various prescriptions and that will finally lead to different form of Fokker Planck equations. One of the possibility is to apply Stratonovich prescription, using which we can compute the integral as 24 : This will correspond to the following form of Fokker Planck equation, as given by: ∂P (n; τ ) ∂τ = − ∂ ∂n (a(n)P (n; τ )) + ∂ ∂n D(n) ∂ ∂n D(n)P (n; τ ) , (6.99) 24 In case of Itô prescription one can recast the integral I(n; τ | ) in to the following form: Now if we consider the contribution from the impurity potential is non vanishing then one can write: Generalized Fokker Planck Equation : ∂ ∂n D(n) ∂P (n; τ ) ∂n + β P (n; τ ) ∂V (n) ∂n = ∂P (n; τ ) ∂τ = − ∂J(n; τ ) ∂n . (6.107) then for equilibrium we set the Fokker-Planck current is zero 25 , for which we get finally the following result: from which we get the following Boltzmann probability distribution function for equilibrium: where P 0 = P (n = 0) is the normalization constant for the probability distribution. Now in the situation where the Fokker Planck current is non-vanishing one can use the following solution ansatz to solve the most Generalized Fokker Planck Equation, as given by: Using this ansatz one can write the following expression: (6.111) Further, substituting this result in the Generalized Fokker Planck Equation we get the following partial differential equation for the unknown function W (n; τ ), as given by: which can be recast in to the following simplified form: ∂ ∂n D(n) ∂W (n; τ ) ∂n − U (n)W (n; τ ) = ∂W (n; τ ) ∂τ , (6.113) where the effective potential U (n) is defined as: Effective Potential : Now let us only consider the time independent part of the Generalized Fokker Planck equation for which the wave function for the equilibrium is described by the following expression: where the solution is given by: where N eq is the normalization constant which can be fixed by the following normalization condition of the equilibrium wave function: , (6.117) where ψ 0 (n) physically represents the ground state wave function with energy eigen value E 0 = 0. All the exicied state have energy eigen value E p > 0 (for p > 1). To get the time dependence of the evolution equation we use the initial condition at time τ = 0 in terms of complete set of eigenfunctions ψ p (n) = n|p , which satisfy the following eigen value equation: which implies the following result: Here the expansion coefficient of the basis is defined as: where we have used the following expression: W (n, 0|n 0 ) = exp β 2 V (n 0 ) δ(n − n 0 ). (6.121) Now for the time dependent part the solution of the Generalized Fokker Planck equation can be written as: Then the probability distribution function can be expressed as: (6.123) In this paper we investigate the physical outcomes of the simplest possibility where V (n) = V (n 0 ) = 0 and U (n) = 0 for which one can write the following simplified expression: Now from Eq (6.127), one can write down the Fokker Planck equation in terms of the following operator equation:D FP P (n; τ ) = 0, (6.125) whereD FP is the Fokker Planck operator represented by: Now, we consider a special case where n >> 1, which gives the most simplest outcome in the present context. In such a situation one can approximately write down the following simplified form of the Fokker Planck equation, as given by: To solve this partial differential equation in the n >> 1 limit we use method of separation of variable, using which we can write: P (n; τ ) = P 1 (n)P 2 (τ ). (6.128) Further, using Eq (6.128) we get the following two sets of independent differential equations, as given by: Solution of Eq (6.145) and Eq (6.146) is given by: where we define a new constant: Additionally, A, B and C are arbitrary constants, which can be determined after imposing appropriate boundary conditions.
Consequently, the most general total solution for the probability density function in the limit n >> 1 can be expressed as: where we define A 1 and B 1 as: Now, after imposing the boundary condition it can be shown that in the large n limit (n → ∞) the solution obtained in Eq (6.150) can be expressed in terms of the following log-normal distribution. To check that explicitly, let us write the probability distribution as the Fourier transformation with respect to the occupation number n, which is given by: P (n; τ |n ; τ ) = 1 2π dk e iknP (k; τ |n ; τ ). (6.136) Using this one can write down the Fokker Planck equation in the Fourier space as: Then solution of Eq (6.137) can be expressed after imposing the initial condition as: Hence substituting back the above mentioned result into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 and further considering n → ∞ limit we get the following result for the probability distribution function, as given by: which is precisely a log-normal distribution function and σ is standard deviation of the log-normal distribution, which can be expressed as: One can explain the physics of this obtained result in the large n limit . In the earlier section we have discussed that the averaging over the phase factor of ln n can be expressed in terms of the logarithms of the occupation number of particles produced in each scattering events. Further using Central Limit Theorem, one can further interpret that ln n follows Gaussian profile (on the other hand, one can also say that in such a case n follows a log-normal distribution). But this physical explanation is only valid for large n limiting approximation.
In figure (36), we have shown the Evolution of the log normal probability density function with respect to the logarithm of the occupation number per mode ln(1+n k ), for a fixed time (µ k τ =fixed). Additionally, it is important to note that from the plot that for very large value of the occupation number n (large n limit) the log normal profile shows Gaussian features perfectly, which indicates the initial assumption regarding large n was consistent. Now, we consider another special case where n << 1 in the present context. In such a situation one can approximately write down the following simplified form of the Fokker Planck equation, as given by: To solve this partial differential equation in the n << 1 limit we use method of separation of variable, using which we can write: P (n; τ ) = P 1 (n)P 2 (τ ). (6.144) Further, using Eq (6.144) we get the following two sets of independent differential equations, as given by: Solution of Eq (6.145) and Eq (6.146) is given by: where we define a new constant: Additionally, D, E and F are arbitrary constants, which can be determined after imposing appropriate boundary conditions. Consequently, the most general total solution for the probability density function in the limit n >> 1 can be expressed as: where we define D 1 and E 1 as: Further using the Fourier transformation with respect to the occupation number n as mentioned in Eq (6.136), we get the following simplified expression for the Fokker Planck equation in n << 1 limit: which is obviously a simplest version of the Fokker Planck equation as it contains a single derivative with respect to time τ . Further imposing the previously used boundary condition for the limit n >> 1 in the present context we get the following result for the probability distribution function in the Fourier transformed space, as given by: Hence substituting back the above mentioned result into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 and further considering n → 0 limit we get the following result for the probability distribution function, as given by: which is not a log normal distribution function in n << 1 limit. One can explain the physics of this obtained result in the small n limit . In this situation one can observe deviations in the profile function. The prime reason for such deviations in small n limit is appearing due to the fact that, the total transmission probability is bounded by unity. In other words, on can say that this is only possible when n is bounded by zero in this context of discussion.
In figure (37), we have shown the evolution of the probability density function with respect to the occupation number per mode n k , for a fixed time (µ k τ =fixed). Additionally, it is important to note that for very small values of the parameter n we have observed from the plot that the deviation from Gaussian feature is observed. In other words, one can interpret that the deviation from log normal probability distribution function corresponds to the significant non-Gaussian features at small values of n. Apart from this one can also comment on the quantum mechanical origin of higher order non-Gaussian contributions appearing in Fokker Planck equation which are more appropriate at small values of n. In the next subsection we will discuss about the physical impacts of this additional higher order contributions in detail.

Corrected probability distribution profiles: Quantum effects from non-Gaussianity
In this subsection we get different order correction to the Fokker Planck equation that we have derived by Taylor expansion. As we already know that the Taylor expansion of the probability density distribution function is taken with respect to time τ . On the other hand, using Maximum entropy ansatz we have considered the Taylor expansion of the ensemble average of the distribution function with respect to the occupation number n. After that we equate both the results and comparing the coefficient of δτ from the both sides of the expansion (see previous Eq (6.68) for more details). Now without truncating both the sides of this expression one can get additional contributions in δτ and in its higher order. If we do the comparison including such additional contributions then it will give rise to corrected version of the Fokker Planck equation valid upto higher orders. Further solving these sets of differential equation order by order one can explicitly justify the validity of all such corrections in the Fokker Planck equation. In this paper we have investigated this possibility by considering the contributions upto fourth order. All such higher order correction terms are very useful to describe the non-Gaussian effects appearing during the process of cosmological particle production during reheating phase of early universe. On top of that, one can explain the origin of such higher order contributions in the quantum mechanical ground as it produces non vanishing significant effects in the expression for the higher order statistical moments directly originating from the various quantum mechanical correlations (one-point, two-point, three-point etc.) computed during cosmological particle production at the epoch of reheating of early universe. More precisely, the deviation from Gaussianity (in other words the deviation from log-normal distribution) in the present context can be directly linked with the quantum mechanical effects appearing during reheating epoch of early universe and for this reason one can interpret the higher order corrected version of the Fokker Planck equation as a quantum corrected Fokker Planck equation. Since in this paper we have provided the analytical correction upto the fourth order, one can say that in this derivation we have actually provided the fourth order quantum corrected Fokker Planck equation. The details of this derivations are explicitly discussed in the following sub sections, where doing the analysis we justify order by order that how such specific corrections will modify the log-normal distribution and its impact in the quantum mechanical ground.

First order contribution
In this context, our prime objective is to find out the first order contribution to the Fokker Planck equation and to solve this equation analytically, which will help us to understand the background physics related to the present formalism. To serve this purpose we equate both the sides of Eq (6.68) after Taylor expansion and compare coefficient of δτ . Consequently, we get the following partial differential equation: First order Fokker Planck Equation : 1 µ k ∂P (n; τ ) δτ = (1 + 2n) ∂P (n; τ ) ∂n + n(1 + n) ∂ 2 P (n; τ ) ∂n 2 . (6.155) Now to solve this partial differential equation we apply method of separation of variable, using which we can write the total solution in the following form: P (n; τ ) = P 1 (n)P 2 (τ ). (6.156) Further, using the solution ansatz stated in Eq (6.156) we get the following sets of independent differential equations, as given by: n(n + 1) d 2 dn 2 + (2n + 1) d dn + m 1 P 1 (n) = 0, (6.157) dP 2 (τ ) dτ + m 1 P 2 (τ ) = 0. (6.158) Solution of Eq (6.157) and Eq (6.158) is given by: Here C 1 , C 2 and C 3 are arbitrary integration constants which can be obtained by imposing appropriate boundary conditions. Additionally, we introduce a constant m 1 which is defined as: which will follow certain constraints in the present context. It is important to note that, to get real valued solution the constant m 1 satisfy the following condition: as Legendre polynomial has general form P N (x) with condition that N should be an integer greater than zero. For different values of N we get different m 1 following Eq (6.162). Consequently, the most general solution of probability distribution function P (n; τ ) is given by the following expression: 163) where we define two new constants, D 1 and D 2 by the following expressions: Here it is important to note that, this solution on limit n → ∞ converge to log-normal distribution as we have discussed earlier. In further section we will compare this result with obtained higher order calculations. From Eq (6.162), we see that the quantization property of m 1 eventually help us to predict the quantum nature of the present set-up. In other words, one can interpret N as a quantum number given by 26 : Quantum number I : N = 0, 1, 2, · · · , ∞ ∈ Z. (6.165) Further using Eq (6.162), one can further introduce another quantum number m 1 given by the following expression: Quantum number II : m 1 = 0, −2, −6, · · · , ∞. (6.166) Further using the Fourier transformation with respect to the occupation number n as mentioned in Eq (6.136), we get the following simplified expression for the Fokker Planck equation in n << 1 limit: ∂P (k; τ ) ∂τ = µ k (2n + 1)ik − k 2 n(n + 1) P (k; τ ), (6.167) which is obviously a simplest version of the Fokker Planck equation as it contains a single derivative with respect to time τ . Further imposing the previously used boundary condition used for the limit n >> 1 and n << 1 in the present context we get the following result for the probability distribution function in the Fourier transformed space, as given by: Hence substituting back the above mentioned result into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 we get the following result for the probability distribution function, as given by: P (n; τ ) = 1 2 µ k n(n + 1)τ π exp −n µ k (n + 1)τ + 1 4µ k τ (n + 1) + 1 , (6.169) which is coming from the first order contribution in the Fokker Planck equation. This expression is actually equivalent to the result obtained in Eq (6.163). 1st order contribution to probability distribution for fixed μτ Figure 38. Evolution of the first order contribution to the probability density function with respect to the the occupation number per mode n k , for a fixed time.
In figure (38), we have shown the evolution of the probability density function with respect to the occupation number per mode, for a fixed time (µ k τ =fixed). For very small values of the parameter n we have observed from the plot that the deviation from Gaussian feature is observed.

Second order contribution
In this context, our objective is to find out the contributions coming from second order in the Fokker Planck equation and to solve this equation numerically 27 . To serve this purpose we equate both sides of Eq (6.68) after Taylor expansion and compare the coefficient of δτ 2 . Consequently, we get the following partial differential equation: Second order Fokker Planck Equation : (1 + n) 2 ∂ 4 P (n; τ ) ∂n 4 + 2n 1 + 3n + 2n 2 ∂ 3 P (n; τ ) ∂n 3 Now to solve this partial differential equation we apply method of separation of variable, using which we can write the total solution in the following form: P (n; τ ) = P 1 (n)P 2 (τ ). (6.171) Further, using the solution ansatz stated in Eq (6.171) we get the following sets of independent differential equations, as given by: It is important to note that, the analytical solution of P 1 (n) is not possible for any arbitrary values of the constant m 2 , except the special case m 2 = 0. For this reason we use numerical technique to solve Eq (6.172). On the other hand Eq (6.173) is exactly solvable in the present context and the solution can be written as: where C 3 and C 4 are two arbitrary constants which can be fixed by choosing proper boundary conditions. Now we solve Eq (6.172) numerically for different values of m 2 along with given initial condition and also we consider the special case m 2 = 0 where we solve this equation analytically. Here it important to mention that, since arbitrary values of m 2 is allowed, one can consider integer as well as non integer values at the level of solution of differential equation. However, the only physically acceptable solution restrict us to only consider the integer values of m 2 because such second order corrected solution of the Fokker Planck equation is directly related to the quantum effects as we have mentioned earlier. As a result such integer values of m 2 can be interpreted as the (principal) quantum number i. e.
Quantum Number III : m 2 = 0, ±1, ±2, · · · , ±∞ ∈ Z. (6.175) For numerical solution we take the following assumptions: Here we assume that the particle production rate at low n(= 0.001) has a constant value and its derivatives also have same constant value for a given m 2 . Getting the numerical solution we plot (P 2 (n, τ ) vs n) them for some particular range of n. From figure 39 we can say that for m 2 = ±2, ±3 second order corrected probability distribution function almost overlap at lower values of the occupation number n but deviate significantly as n increases to large number. As a consequence, for low value of n particle production rate is independent on m 2 but as n increases they significantly deviate. It also implies that for higher values of n the integer m 2 constrains particle production rate. For m 2 = ±1 we found that the second order corrected probability distribution function significantly deviates from log normal (Gaussian) distribution and both of them explicitly show the signature of non-Gaussianity is the second order corrected distribution function. Finally, we have shown that for m 2 = 0 the amount of deviation from log normal (Gaussian) distribution is small compared to results obtained from the other values of m 2 . Now we discuss about the analytical solution of Eq (6.172) for the special case when we fix m 2 = 0. In this situation one can recast the Eq (6.172) in to the following simplified form: (1 + n) 2 d 4 dn 4 + 2n 1 + 3n + 2n 2 d 3 dn 3 + 1 + 6n + 6n 2 d 2 dn 2 P 1 (n) = 0 (6.177) Then analytical solution of Eq (6.177) can be written as: where a, b, c, b 1 , b 2 and u are all functions of i which is summed over and we have introduced them to use shorthand notation. In this context the functional dependence of all of these i dependent parameters are given by the following expressions: The solution contains generalized Hypergeometric PFQ Regularized function. Also C 1 and C 2 are arbitrary constant of integration which can be evaluated by imposing appropriate initial conditions.
For m 2 = 0 the total probability distribution function can be expressed as: For the special case m 2 = 0 and considering the large n limit (n → ∞) the Eq (6.177) reduces to the following form: After solving Eq (6.64) we get the following solution in the large n limit: where C 5 , C 6 , C 7 and C 8 are the arbitrary constants of integration which can be evaluated by imposing appropriate initial condition.
In the large n limit with m 2 = 0 the total probability distribution function can be expressed as: P (n; τ ) = 1 6 C 5 n 2 + 3C 6 n + C 7 n + C 8 C 3 e τ m 2 µ k + C 4 e −τ m 2 µ k , (6.183) which shows huge deviation from log normal (Gaussian) distribution. Further using the Fourier transformation with respect to the occupation number n as mentioned in Eq (6.136), we get the following simplified expression for the Fokker Planck equation at the second order: (1 + n) 2 k 4 − 2ink 3 (1 + 3n + 2n 2 ) − k 2 (1 + 6n + 6n 2 ) P (k; τ ), (6.184) which is obviously a simplest version of the Fokker Planck equation as it contains only two derivative with respect to time τ . In the present context we get the following result for the probability distribution function in the Fourier transformed space, as given by: where G(k; n ) is defined as: Additionally, C 1 and C 2 are arbitrary constants which is fixed by the following two fold boundary conditions, as given by: P (n; τ |n = 0; τ = τ ) = δ(n), (6.187) ∂P (n; τ |n ; τ ) ∂τ n =0,τ =τ = − δ(n) n , (6.188) which are necessary to solve the above mentioned second order differential equation.
As a result, we get the following set of constraints equations: Solving these equations we get: . (6.191) Using this solution we get the following probability distribution function in Fourier space with n = 0 and τ = 0, as given by: (6.192) Hence substituting back the above mentioned result into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 we get the following result for the probability distribution function, as given by: However this integral is not convergent within −∞ < k < ∞. For this reason we introduce a momentum cut-off −Λ C < k < Λ C . Consequently, we get the following regularised expression for the probability distribution function: In figure [40], we have shown the probability density function for second order correction with respect to the occupation number per mode, for a fixed time (τ = 0.15 µ k ). From this plot we have observed irregular oscillations with deviation from Gaussian feature.We have selected a high value for momentum cutoff[λ c ] for this particular plot.

Third order correction
In this context, our objective is to find out the contributions coming from third order in the Fokker Planck equation and to solve this equation numerically 28 . To serve this purpose we equate both sides of Eq (6.68) after Taylor expansion and compare the coefficient of δτ 3 . Consequently, we get the following partial differential equation: (1 + n) 3 ∂ 6 P (n; τ ) ∂n 6 + 3n 2 2 (1 + n) 2 (1 + 2n) ∂ 5 P (n; τ ) ∂n 5 +3n(1 + n)(1 + 5n + 5n 2 ) ∂ 4 P (n; τ ) ∂n 4 which can not able to solve analytically with any integer values of m 3 . We solve this equation for different values of m 3 numerically with assumed initial condition. Only for the special case, m 3 = 0 with large n limit we can able to provide an analytical solution in the present context. Now to solve this partial differential equation we apply method of separation of variable, using which we can write the total solution in the following form: P (n; τ ) = P 1 (n)P 2 (τ ). (6.196) Further, using the solution ansatz stated in Eq (6.171) we get the following sets of independent differential equations, as given by: It is important to note that, the analytical solution of P 1 (n) is not possible for any arbitrary values of the constant m 3 , except the special case m 3 = 0. For this reason we use numerical technique to solve Eq (6.197). On the other hand Eq (6.198) is exactly solvable in the present context and the solution can be written as: From Fig 41 one can say that the third order corrected probability distribution function for different value of m 3 = 0 overlap at lower n limit (n → 0) though deviate significantly at large n limit. At lower values of n, particle production probability is independent of m 3 and almost flat. But as soon as n reaches values greater than unity the distribution increase exponentially and for different m 3 they differ from each other. Now we discuss about the analytical solution of Eq (6.197) for the special case when we fix m 3 = 0. In this situation one can recast the Eq (6.197) in to the following simplified form: To find the analytical solution of Eq (6.202) one can further use the following simplification: where we introduce a new function Q 1 (n) which can be expressed in terms of P 1 (n) by following identification: On large n limit (n → ∞) the Eq 6.203 reduces to following extremely simplified form: The solution of this equation at large n limit can be expressed as: √ 71 ln n + 240n 15/2 ln n , (6.206) where C 1 ,C 2 and C 3 are arbitrary constant of integration and can be evaluated by imposing appropriate initial condition.
Here, in the large n limit the solution for P 1 (n) can be written as: where C 4 , C 5 and C 6 are arbitrary constant of integration and can be evaluated by imposing appropriate initial condition. Finally, in the large n limit the total probability distribution function can be expressed as: P (n; τ ) = 1 674880n 9/2 3040n 15/2 (37C 3 − 60 ln n + 110) which shows large deviation from log normal (Gaussian) distribution. Further using the Fourier transformation with respect to the occupation number n as mentioned in Eq (6.136), we get the following simplified expression for the Fokker Planck equation at the third order: which is obviously a simplest version of the Fokker Planck equation as it contains only three derivative with respect to time τ . In the present context we get the following result for the probability distribution function in the Fourier transformed space, as given by: where O(k; n ) is defined as: Additionally, C 1 , C 2 and C 3 are arbitrary constants which is fixed by the following three fold boundary conditions, as given by: which are necessary to solve the above mentioned third order differential equation.
As a result, we get the following set of constraints equations: Solving these equations we get: (6.220) Using this solution get the following probability distribution function in Fourier space, as given by: Hence substituting back into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 we get the following result for the probability distribution function, as given by:  In Fig:-42, we have shown the third order correction of probability density function with respect to the occupation number per mode, for a fixed time (µ k τ =fixed). From this plot we have observed primary gaussian feature followed by an exponential type increase.We consider the perturbative expansion to be valid. So the higher order correction contribute less than the previous order.To normalize the analytical solution we have used this. In latter part we can see the exponential type increase contribute in longer tail effect [ higher kurtosis].

Fourth order correction
In this context, our objective is to find out the contributions coming from fourth order in the Fokker Planck equation and to solve this equation numerically 29 . To serve this purpose we equate both sides of Eq (6.68) after Taylor expansion and compare the coefficient of δτ 4 . Consequently, we get the following partial differential equation: +20n(1 + n)(1 + 2n)(1 + 7n + 7n 2 ) ∂ 5 P (n; τ ) ∂n 5 which can not able to solve analytically with any integer values of m 4 . We solve this equation for different values of m 4 numerically with assumed initial condition. Only for the special case, m 4 = 0 with large n limit we can able to provide an analytical solution in the present context. Now to solve this partial differential equation we apply method of separation of variable, using which we can write the total solution in the following form: P (n; τ ) = P 1 (n)P 2 (τ ). (6.224) Further, using the solution ansatz stated in Eq (6.223) we get the following sets of independent differential equations, as given by: It is important to note that, the analytical solution of P 1 (n) is not possible for any arbitrary values of the constant m 4 . For this reason we use numerical technique to solve Eq (6.225).
Also considering the large n limit we have checked that Eq (6.225) is not analytically solvable. On the other hand Eq (6.226) is exactly solvable in the present context and the solution can be written as: where C 9 , C 10 , C 11 and C 12 are three arbitrary constants which can be fixed by choosing proper boundary conditions. Now to solve Eq (6.225) numerically for different values of m 4 along with given initial condition. Here it important to mention that, since arbitrary values of m 4 is allowed, one can consider integer as well as non integer values at the level of solution of differential equation. For numerical solution we take the following assumptions: Accoding to our assumption particle production probability has constant value at some particular small n value (n = 0.0001) and all its derivative has constant values and all those values are same. During the analysis we assume that the particle production probability and all its derivative has a constant same value for n = 0.0001, which is very very helpful for us to deal with the initial conditions during performing numerical techniques to solve the Eq (6.225).
From fig. 43 we observe that the fourth order corrected probability distribution function for different m 4 is almost flat upto n = 1 and after that the distribution function suddenly increases. Additionally, we observe that the fourth order correction has deviation from Gaussianity at small values of the occupation number. On the other hand, for large values of the occupation number we get a Gaussian like feature and that is shown explicitly in the mentioned plot.
Further using the Fourier transformation with respect to the occupation number n as mentioned in Eq (6.136), we get the following simplified expression for the Fokker Planck equation at the fourth order: (1 + n)(1 + 2n)(1 + 7n + 7n 2 )k 5 +(1 + 20n + 90n 2 + 140n 3 + 70n 4 )k 4 P (k; τ ), (6.230) which is obviously a simplest version of the Fokker Planck equation as it contains only four derivative with respect to time τ . In the present context we get the following result for the probability distribution function in the Fourier transformed space, as given by: where J (k; n ) is defined as: (1 + n )(1 + 2n )(1 + 7n + 7n 2 )k 5 +(1 + 20n + 90n 2 + 140n 3 + 70n 4 )k 4 . (6.232) Additionally, C 1 , C 2 , C 3 and C 4 are arbitrary constants which is fixed by the following three fold boundary conditions, as given by: P (n; τ |n = 0; τ = τ ) = δ(n), (6.233) ∂P (n; τ |n ; τ ) ∂τ n =0,τ =τ = − δ(n) n , (6.234) ∂ 2 P (n; τ |n ; τ ) ∂τ 2 n =0,τ =τ = 2 δ(n) n 2 (6.235) which are necessary to solve the above mentioned fourth order differential equation. As a result, we get the following set of constraints equations: 238) 239) (6.240) Solving these equations we get: (6.242) Using this solution get the following probability distribution function in Fourier space, as given by: Hence substituting back into the definition of Fourier transformation and setting the initial condition n = 0 and τ = 0 we get the following result for the probability distribution function, as given by: which is divergent within the interval −∞ < k < ∞. After introducing an IR and UV regulators, Q < k < L we can get a finite result, which we have not presented in this paper for its huge length. The origin of such corrections are the fourth order contribution in the Fokker Planck equation. In Fig:-44, we have shown the fourth order correction to the probability density function with respect to the occupation number per mode, for a fixed time (µ k τ =fixed) flow the intial gaussian then exponential decay type distribution.This decreasing feature suggest very low effect from fourth order correction which is also supported by perturbative expansion assumption.

Total solution considering different order correction
Now we plot total solutions with different order of correction with main solution.Then we check validation of different order of correction altogether and how they merge with each other at what limit. A. Upto second order correction:-Further we consider upto second order corrected solution with first order contribution at m 1 = 2 for different m 2 .  From Fig:-45(a), 45(b) we observe that at low values of n, P 1 + P 2 and P 2 are significantly different, but as we increase n they overlapped and P 2 effect is more over P 1 + P 2 so the second order solution dominate over the first order solution.On the other hand Fig:-45(c) shows how the second order contribution effect the primary gaussian curve. Distinct oscillation effect on gaussian curve shows the effect of second order correction which is not available from the numerical solution.

B. Upto third order correction:-
Here we add all the three previously derived contributions to produce the total probability distribution corrected upto third order.   fig. 46(a),46(b) we observe that P 1 + P 2 + P 3 , P 2 + P 3 and P 3 overlap at higher n limit though separated. At low n limit and P 1 + P 2 and P 2 overlap with each other but remain separated from P 1 + P 2 + P 3 for the complete range. This implies that third order contribution is dominant over the other two due to the non-linearities in the differential equation and behaves like a non-perturbative quantum effect at the level of solution. From  fig 46(c) we show the third order correction affects the tail of the gaussian and shift it bit higher.It is clear from fig 46(d). C. Upto fourth order correction:-Here we add all the four previously derived contributions to produce the total probability distribution corrected upto fourth order.  From fig. 47(a) we observe the final curve represented by P 1 + P 2 + P 3 + P 4 shifted the mean of the gaussian from its value at P 1 .Initial gaussian feature for P 1 +P 2 +P 3 and P 1 +P 2 are deviated at high n limit . Here P 1 + P 2 + P 3 + P 4 and P 1 follow exact gaussian nature but they are mirror image of one another. P 1 + P 2 + P 3 and P 1 + P 2 don't have this gaussian nature due to their divergence property as n increases. In fig 47(b) and fig 47(c) the effect of all order correction can be observed but effect of different order correction become evident from fig 47(d).Second order correction introduce the oscillating feature whereas third order correction increase the tail and responsible for higher kurtosis. Fourth order correction add a small positive effect to the previous correction without any shape change.
Previously it is shown that in ref. [42] the distribution will be log-normal at large n, considering the lowest order contribution coming from the solution of Fokker-Planck equation. We extend this result upto fourth order and shown the effect of different order correction.The oscillating nature and long tail can be specific feature of stochastic particle production in inflation epoch.This may be the effect of background field or noise at that time.The numerical solutions give different quantum numbers which support the quantum nature. In the next subsection we will calculate various statistical moments and from that we can discuss about the role of quantum effects and non-Gaussinanity from the Probability distribution profile for particle production in the context of early universe cosmology (mostly during reheating).

Calculation of statistical moments (or quantum correlation functions) from corrected probability distribution function
Here our prime objective is to compute the different statistical moments from the quantum corrected probability distribution function as obtain by Taylor expanding in order by order from Eq (6.68). From is corrected probability distribution function we compute the expression for n , n 2 , n 3 and n 4 and then calculate standard deviation, skewness and kurtosis for a given time. We have explicitly shown that the non vanishing contributions of skewness and kurtosis carries the signature of significant effect of non-Gaussianity. In this analysis the values of these moments are compared with predicted results obtained from log-normal (Gaussian) distribution and non zero values of kurtosis and skewness define the deviation from that.
To compute the moments we start with the following sets of master evolution equations valid in different orders, as given by: First Order Master Evolution Equation : where the first moment or the expectation value of the observable F is define as: Here F (n) is the physical observable in which we are interested in and P (n; τ ) is the corrected probability distribution function which is not necessarily log-normal (Gaussian) in nature.
In the present context of discussion, Eq (6.249) plays the role of generating function, which is commonly used in calculating various mathematical special functions. In our discussion, Eq (6.249) represents the statistical moment generating function in presence of quantum corrected probability distribution function. In terms of quantum mechanical language, Eq (6.249) signify the one point quantum correlation function and it is exactly equal to the statistical first moment in this discussion. Now, we explicitly compute the expressions for one point function (or first moment) of the occupation number i.e. n , two point function (or second moment) of the occupation number i.e. n 2 , three point function (or third moment) of the occupation number i.e. n 3 and four point function (or fourth moment) of the occupation number i.e. n 2 using the previously mentioned first, second, third and fourth order master equations. The detailed steps of the computations are appended bellow:

1.
Step I: First of all, we use the first order master evolution equation. Then we replace the function F by the occupation number n. Consequently, we get the following time evolution equation of the first moment or one point function n , given by: Step II: Secondly, we want to compute the expression for n 2 . To compute this we consider here the first and second order master equations, as mentioned earlier. Considering only the first order master equation we get the following analytical expression: 1 µ k ∂ n 2 ∂τ = 2n(1 + 2n) + 2n(1 + n) = 4n + 6n 2 = 4 n + 6 n 2 . (6.251) On the other hand, using the second order master equation we get the following analytical expression for the time evolution of the second moment or two point correlation: Finally, using the third order master equation we get the following analytical expression for the time evolution of the fourth moment or fourth point correlation: 1 µ 4 k ∂ 4 n 4 ∂τ 4 = 24(1 + 20n + 90n 2 + 140n 3 + 70n 4 ) = 480 n + 2160 n 2 + 3360 n 3 + 1680 n 4 + 24 . (6.259)

6.
Step VI: Using the result obtained in Step I and using the previously mentioned boundary condition we get the following expression for the one point function 30 (or first moment) of occupation number: First Moment (First Order) : which is further used to compute all the higher order moments from master evolution equation considering higher order Taylor expansion. Additionally, it is important to note that if we use higher order equations for the first moment then after imposing the boundary conditions we get the following results: Consequently, the total first moment can be written as: Total First Moment : n = n I + n II + n III + n IV = n I = 1 2 (e 2τ µ k − 1) .

7.
Step VII: Using the results obtained in Step II and using the previously mentioned boundary condition we get the following expression for the two point function 31 (or second moment) of occupation number: Second Moment (First Order) : Second Moment (Second Order) : which is further used to compute all the higher order moments from master evolution equation considering higher order Taylor expansion. Additionally, it is important to 31 It is important to note that, as far as quantum mechanical computation is concerned, it produces not exactly same result for the one point function and first moment of the occupation number. For two point function we actually get the following result: n(τ )n(τ ) = A(τ )δ(τ + τ ), (6.265) where A(τ ) is the amplitude of the two point function as given by the following expression: A(τ ) = n 2 = n 2 I + n 2 II . Consequently, the total second moment can be written as: Total Second Moment : n 2 = n 2 I + n 2 II + n 2 III + n 2 IV = n 2 I + n 2 II . (6.271) In fig. (49), we have explicitly shown the time dependent behaviour of second moment or amplitude of the two point function of the occupation number n 2 . As there is no contributions are coming from the third and fourth order moment generating master evolution equation for n 2 , the only contribution is coming from the first and second order master evolution equation. From this plot we see that for a fixed value of the parameter µ k = 1, 10, 100, at the lower values of the time the second moment or the amplitude of the two point function of the occupation number initially increase with time very very slowly. Then after a certain time when τ >> 1 it shows suddenly huge increment in the behaviour.

8.
Step VIII: Using the results obtained in Step III and using the previously mentioned boundary condition we get the following expression for the three point function 32 (or third In fig. (50), we have explicitly shown the time dependent behaviour of third moment or three point function of the occupation number n 3 . As there is no contributions are coming from the fourth order moment generating master evolution equation for n 3 , the only contribution is coming from the first, second and third order master evolution equation. From this plot we see that for a fixed value of the parameter µ k = 1, 10, 100, at the lower values of the time the third moment or the equal time amplitude of the three point function of the occupation number initially increase with time very slowly. Then after a certain time it shows suddenly huge increment in the behaviour. Most importantly, this plot shows the third moment or equal time amplitude of the three point function of the occupation number is not zero. This shows the second signature of the non-Gaussianity as we know for Gaussian probability distribution profile this is exactly zero.

9.
Step IX: Using the results obtained in Step IV and using the previously mentioned boundary condition we get the following expression for the four point function 33 (or fourth 33 It is important to note that, as far as quantum mechanical computation is concerned, it produces not exactly same result for the three point function and third moment of the occupation number. For three point function we actually get the following result: n(τ )n(τ )n(τ )n(τ ) = C(τ, τ , τ " , τ )δ(τ + τ + τ " + τ " ), (6.279) where C(τ, τ , τ " , τ ) is the amplitude of the four point function. If we fix τ = τ = τ " = τ (equal time) then we get the following expression: C(τ, τ, τ, τ ) = n 4 = n 4 I + n 4 II + n 4 III + n 4 IV .   The third and the fourth order corrected version of the fourth order moment equations are not exactly solvable analytically. For this reason we have applied numerical techniques to solve these differential equations.
Consequently, the total fourth moment can be written as:  This implies that the equal time amplitude part is exactly matches with the third moment of the occupation number in this context.
In fig. (51), we have explicitly shown the time dependent behaviour of fourth moment or amplitude of the four point function of the occupation number n 4 . From this plot we see that for a fixed value of the parameter µ k = 1, 10, 100, at the lower values of the time the third moment or the equal time amplitude of the four point function of the occupation number initially increase with time very slowly. Then after a certain time it shows suddenly huge increment in the behaviour.

Standard Deviation
Further, using the results obtained in the context of second moment or two point correlation function, in this subsection our prime objective is compute the expression for the Standard Deviation from the corrected version of the probability distribution function. In the present context of discussion Standard Deviation actually gives the spread of the peak of the corrected probability distribution function. Therefore, Standard Deviation considering upto first order is given by the following expression: S.D. uc = n 2 I − ( n I ) 2 = √ 2e 6τ µ k − 3e 4τ µ k + 1 2 √ 3 , (6.284) where the subscript "uc" stands for uncorrected.
On the other hand, after including the result from the second order the corrected expression for the Standard Deviation can be expressed as: From the fig.-52, we can see that the uncorrected Standard Deviation (first order) and corrected Standard Deviation (second order) has significant difference in low µ k τ limit and second order overlapped as they approach higher µ k τ . So for lower limit this second order correction is significant and for this reason during the computation of Kurtosis and Skewness we use total solution of standard deviation over the uncorrected one alone.In fig 52(a) the variance is with in the value 1 but in fig ?? variance has a enormous value.It is mere effect of tuning.The plots or values of variance can always be tuned to be within 1 using proper prefactor.

Skewness
In this subsection, our prime objective is to computed the expression for the Skewness from the corrected probability distribution function. Skewness actually measure the asymmetry of the probability distribution function of a real-valued random variable about its mean value. This measure can be positive or negative, or undefined. From positive skewness (for unimodal distribution) we can say normal curve has longer right tail.
Therefore, skewness without correction can be expressed as: Now from fig. 53, we can say that the corrected Skewness deviate significantly from the uncorrected one at low µ k τ limit. But we can see that at higher limit they overlap. Also Skewness is positive for the whole range which implies that the normal distribution curve has longer right tail. Moreover, there is a discontinuity of third order corrected Skewness in between the range 0.1 < τ < 1 and for the rest of the whole range of time Skewness decreased upto unity and then it is increased.

Kurtosis
Kurtosis is a measure of the tailedness of the probability distribution of a real-valued random variable. This is actually a descriptor of the shape of a probability distribution function and there are specific ways of quantifying it for a theoretical probability distribution and corresponding ways of estimating it from a sample from a population. It is important to note that, the Kurtosis of any univariate normal distribution is 3. For practical purposes it is common practice to compare the expression for Kurtosis of a probability distribution function to 3. Probability distributions with Kurtosis less than the value 3 are identified as platykurtic, although this information does not imply the distribution is flat-topped in nature. Rather, it implies that the probability distribution produces fewer and less extreme outliers than does the normal probability distribution. Probability distributions with Kurtosis greater than the value 3 are said to be leptokurtic. It is also common practice to use, the excess Kurtosis, which is the Kurtosis minus 3, to provide the comparison to the normal probability distribution profile. Like Skewness here also we calculate kurtosis from different distribution and get it at different order of correction.
Therefore, Kurtosis without correction can be expressed as: = 6 µ (−3e 4µ k τ + 2e 6µ k τ + 1) × 6(µ k − 1)e 2µ k τ + 3(5µ k − 2)e 8µ k τ − 2e 6µ k τ (3µ(4µ k τ + 3) − 5) − 3µ k + 2 .   fig. (54), we can say that uncorrected kurtosis deviate from corrected one in lower τ regime though overlapped in higher order [τ > 1]. So the contribution from the correction factor is important in lower regime. Also it is important to note that, Kurtosis is greater than the value 3 for the whole time regime, so the distribution is Leptokurtic and have fatter tails. Here fig 54(a) has a huge value for corrected kurtosis which can be bounded within a particular value using tuning parameter.
From the calculation of the higher-order statistical moments (or equivalently the amplitude of the quantum mechanical correlation functions) we get the following overall features to analyze the nature and physical outcomes of the corrected probability density function derived in this paper.
1. The Standard Deviation is significantly large for higher µ k τ , though very small for lower regime.
2. Skewness is positive throughout the time regime, though becomes vanishingly small at a specific time interval (0.1 < τ < 0.7).
3. Kurtosis is greater than 3 for the whole time regime.
4. The predicted Log-Normal Gaussian Distribution shows deviations at significant levels. Effects of the non-Gaussianity in the distribution function is clearly visualized.
5. The probability distribution has longer trailing ends and the trails go broad higher.
6. The probability distribution has a very low spread at lower µ k τ limit though highly spread out in larger limit.
Conequently, from the previously predicted result of [42] we show that distribution deviates from a Log-Normal distribution and the predicted distribution function appears to be leptokurtic and has a broad right trailing end.All of this portion have a simmilar connection with ref [105] 7 Conclusion In this paper we have addressed the issues which are appended below: • In this paper, we have provided the analogy between particle creation in primordial cosmology and scattering problem inside a conduction wire in presence of impurities. Such impurities are characterized by effective potential in the context of quantum mechanical description. On the other hand, in the context of primordial cosmology time dependent mass profile of created particles (couplings) mimics the same role.
• Specific time dependence of mass profile actually restricts the structure of the scattering effective potential. To establish the analogy between two theoretical frameworks we have further computed various characteristic features of conduction wire i.e. resistance, conductance (electrical properties), Lyapunov exponent (dynamical property), reflection and transmission coefficients (optical properties), occupation number and energy density (energetics) using the expression for Bogoliubov coefficients for different mass profiles which connects the ingoing and outgoing solution of the mode functions obtained in the context of particle creation process in cosmology.
• We have solved this particle creation problem using the following crucial steps: 1. First of all assuming that the interactions are well known we have studied the one to one correspondence between the particle creation problem in early universe cosmology with the scattering problem inside a conduction wire. Here we have additionally neglected the effect of the expansion of our universe and this is perfectly justifiable during the epoch of reheating. For this reason we call it as Reheating Approximation.
2. Secondly we have studied the same problem where the particle interactions are not known at all at the level of action. In such a situation, assuming the gravitational background is classical in nature and also assuming the previously mentioned Reheating Approximation we have demonstrated the problem with the help of Random Matrix Theory.
3. Further we have solved the dynamics of the particle creation problem by studying the higher order corrections in the Fokker Planck equation for previously mentioned random system where the interactions are not easily quantifiable at the level of action. We have constructed the fourth order corrected Fokker Planck equation from which we have provided the solution of the random probability distribution function. Such distributions are very very useful to study the dynamical systems when particle interactions are not well known. In our analysis we have identified all of these modifications as the quantum correction to the Fokker-Planck equation, the physical implications of which we have studied in detail in this paper.
• In this work, we have shown that the Lyapunov exponent varies inversely with the number of scatterers. Therefore, with an increase in the number of scatterers the Lyapunov exponent also reduces thereby reducing the amount of randomness in the system. This may be a hint to the fact that the Lyapunov exponent has a dependence on the momenta values of the incoming wave-function of the scalar field. Additionally, it is important to note that the upper bound of Lyuapunov exponent is restricted by the constraint λ ≤ 2π/β (where β = 1/T ), which is a generic bound on chaos obtained in the context of quantum field theory. As a consequence, one can find restriction on the upper bound on the reheating temperature for the different quenched mass profiles for which the chaos bound saturates. This is obviously a remarkable result in the present context as it can able to provide the explicit expression for the reheating temperature for a specified momentum scale, which was not predicted earlier in the detailed study of reheating. Most importantly, the bound on quantum chaos in terms of Lyapunov exponent directly restrict the value of reheating temperature without explicitly knowing the details of the particle interactions as appearing in the action. Just the knowledge of time dependence of the quenched mass profiles (in other words the knowledge of effective impurity potential as appearing inside the conduction wire) is sufficient enough to restrict the upper bound of reheating temperature due to quantum chaos.
• In this context we have also provided the expression for the two point quantum correlation function, which is known as Spectral Form Factor (SFF) for both in finite and zero temperature. Spectral Form Factor is actually a more strong measure to find chaotic behaviour of a dynamical system compared to Lyapunov exponent. We get saturating behaviour of SFF at late time scale, which indicates that it has an upper-bound. We can relate SFF for any potential (Even Polynomial Potential in this case). In the calculation of the Lyapunov Exponent for the specific time dependent mass profiles, we choose three different quenched protocols for mass profiles. Potential functions which can be represented by polynomial potential (Even only in our case) can be used to get the SFF-saturation. In this connection ,we have provided a model independent upper and lower bound of SFF, which is treated as the significant bound of quantum chaos (−1/N (1 − 1/π) ≤ SFF ≤ 1/πN ) in the context of particle production event in cosmology. This is obviously a remarkable result which we have explicitly computed in this paper.In [106] this same has been calculated for general polynomial with GUE.
• We have also presented the computation of quantum corrected Fokker-Planck equation which corresponds to the delta-scatterers. From this computation we have derived the corrected statistical distribution of the particle production events in cosmology. The distribution which has been predicted in [42] to be Gaussian doesn't retain its form when more correction terms are taken into account. This may be treated as a signature of non-Gaussian in particle production events during reheating (in cosmology).
The future discussions of the present work are mentioned in the following: • In this paper for our study of quantum chaos in the context of cosmology we have used a closed quantum system. As we have mentioned that the present computation has been performed for a massless scalar field which interacts with the heavy fields (which acts like scatterers inside the conduction wire). The entire calculation is being done for the set up when there is only a single massless scalar field that interacts with the scatterer. One may repeat the calculation for a large number of these scalar fields interacting with the scatterers which needs the introduction of the random matrix approach in a more generalized fashion.
• The system we have studied in this paper have no interactions with the background as the definition of the background in this set-up is itself an ill-defined one during reheating. To treat the entire system having being interacted with a background one needs to have a detailed description of the nature of background in the cosmological scenario. Then it will be possible to introduce the other non-linear and dissipative effects into the system introduced by the background itself. Such a treatment will be studied within the framework of an open quantum system interacting with the defined background set-up. One then needs to consider the entire system having being interacted with the background under a weak coupling limit. One such model as has been studied in [107].
• We have calculated the Lyapunov exponent and Spectral Form Factor in this paper which is a measure of chaos or non-linearity into the system. With the system prescribed in this work being treated as an open quantum system one may study the effects of dissipation being introduced into such a system which renders the system to be a stochastic one. With this, one may be able to study the effects of the nonlinearity being introduced into the system which may well be a good study to look for the behaviour of Lyapunov exponent and Spectral Form Factor.
• During the study of quantum correction in the Fokker Planck equation and the deviation from log normal distribution we have followed a specific approach in which we have considered the following possibilities: 1. We have neglected the contribution from the damping term in the Fokker Planck equation. One can include such effect and study its role in the context of cosmology (specifically during reheating).
2. During the computation we have followed a specific approach in which we have also neglected the effect of impurity potential at very high temperature during reheating. This will give rise to a simplest form of the Fokker Planck equation where only diffusion and drift contributions are appearing explicitly. But if we include the effect of impurity potential in presence of finite temperature then it will surely effect the final solution of the probability distribution function. One can include such additional effects and study its impact during reheating epoch of the early universe.
3. Furthermore, during the construction of the Fokker Planck equation from the basic principles we have followed a special approach in which the effect of diffusion and drift is appearing in a very simplified manner. However, in the study of statistical field theory Itô and Stratnovitch or more generalized prescriptions are used commonly to construct the Fokker Planck equation. Here it is important to note that in each case it will give rise to different mathematical structure of Fokker Planck equations. In the present context of discussion, one can follow such well known prescriptions to see its physical outcomes to solve the probability distribution function for the particle production and compare the results to check the appropriateness of these approaches during reheating. ∂P (n; τ ) ∂τ = − ∂ ∂n (a(n)P (n; τ )) + ∂ ∂n D(n) ∂ ∂n D(n)P (n; τ ) , (B.1) Here we take a(n) = 0 and D(n) = n(n + 1). Using this we get the following solution of probability distribution: P (n, τ ) = 1 2 √ π n(n + 1)τ µ k exp − 9(2n + 1) 2 τ µ k 16n(n + 1) (B.2) In fig. (56) we have shown the probability distribution function obtained from the Stratonovitch solution of the Fokker Planck equation. This solution is similar to the log normal distribution obtained from the present computation. From the plot we observe that for large value of occupation number n (n >> 1) the distribution function decays to a finite saturation value. On the other hand for small occupation number n (n << 1) we get peak in the distribution function for different values of µ k τ .      In fig. (58) we have shown the probability distribution function obtained from the solution of of the Fokker Planck equation derived in presence of finite temperature effective potential solution. From the plot we observe that for large value of occupation number n (n >> 1) the distribution function decays to a finite saturation value. On the other hand for small occupation number n (n << 1) we get peak in the distribution function for different values of µ k τ .