Squeezing, Chaos and Thermalization in Periodically Driven Quantum Systems: The Case of Bosonic Preheating

The phenomena of Squeezing and chaos have recently been studied in the context of inflation. We apply this formalism in the post-inflationary preheating phase. During this phase, inflaton field undergoes quasi-periodic oscillation, which acts as a driving force for the resonant growth of quantum fluctuation or particle production. Furthermore, the quantum state of the fluctuations is known to have evolved into a squeezed state. In this submission, we explore the underlying connection between the resonant growth, squeezing, and chaos by computing the Out of Time Order Correlator (OTOC) of phase space variables and establishing a relation among the Lyapunov, Floquet exponents, and squeezing parameters. For our study, we consider observationally favored $\alpha$-attractor E-model of inflaton which is coupled with the bosonic field. After the production, the system of produced bosonic fluctuations/particles from the inflaton is supposed to thermalize, and that is believed to have an intriguing connection to the nature of chaos of the system under perturbation. %By using this we calculated approximate lower bound of temperature ${\bar T}_{\rm MSS}$. We conjecture a relation between the thermalization temperature $({\bar T}_{\rm SS})$ of the system and quantum squeezing, which is further shown to be consistent with the well-known Rayleigh-Jeans formula for the temperature symbolized as ${\bar T}_{\rm RJ}$, and that is ${\bar T}_{\rm SS} \simeq {\bar T}_{\rm RJ}$. Finally, we show that the system temperature is in accord with the well-known lower bound on the temperature of a chaotic system proposed by Maldacena-Shenker-Stanford (MSS).


INTRODUCTION
Post-inflationary reheating is an important event in the early universe that is responsible for the thermalized universe we see today.Therefore, a theoretical understanding of this phase would be of paramount interest.Inflationary quantum fluctuations are responsible for producing curvature fluctuations, and their subsequent generation of structures is inevitably related to the process of transferring energy from inflaton to matters.In spite of having significant efforts , a proper theoretical understanding of this phase is still lacking particularly due to its complex non-linear evolution, and subsequent thermalization.
After a quasi-exponential accelerated phase of expansion, known as inflation [1][2][3][4][5][6][7], the universe is left with a super cold state of vanishing matter, entropy density, and a homogeneous inflaton field oscillating around its potential minima.When coupled with such oscillating inflaton, any fundamental field will experience a periodic driving force and have quantum mechanical production.Such production, therefore, leads to the transfer of energy from inflaton to standard model fields successfully populating the universe.Typically such a transferring process involves multiple stages of evaluations.Depending upon the coupling parameter, the initial stage may be non-perturbative [8][9][10][11][12][13][14][15][16][17][18] followed by perturbative production, or the process could be entirely perturbative [19][20][21][22][23][24][25][26].However, apart from non-linear nature of production, it is the thermalization of those produced particles in the expanding background which makes the whole process extremely complex and theoretically challenging.Studies of perturbative production and its thermalization have been discussed in the literature [27][28][29][30] taking a phenomenological approach with very little foundational progress.Considering the lattice simulation, studies have also been performed for non-perturbative production and their subsequent thermalization [31][32][33][34][35].However, to the best of our knowledge, a proper underlying mathematical framework to understand this thermalization process in the context of reheating is far from complete.It is in this realm, that we initiate formulating a theoretical framework that we believe to play an instrumental role in understanding such a phase.Apart from its theoretical motivation, it is important to stress that the production and thermalization during reheating may leave an indelible imprint on different cosmological observable related to Cosmic Microwave Background (CMB), baryogenesis, dark matter, and even small-scale structures in the universe.
In this submission, we particularly focus on the stage called preheating when particle production occurs with parametric resonance as a primary mechanism, and hence would naturally last for a very brief initial period of the entire reheating process.The formalism we adopt is based on the recent development [36][37][38][39][40][41] on the application of various intriguing concepts of quantum information theory in the realm of cosmology particularly centered around one of the fundamental questions in theoretical cosmology as to how the quantum correlation of fluctuations generated during inflation is dynamically decohered into post-inflationary classicalized correlation, and what would be the appropriate quantifier of their quantumness [36,42].During inflation, the primordial perturbation generically evolves into a squeezed state [43,44].Utilizing this squeezed state and to decode the quantumness of the inflationary correlations, a number of proposals have been put forth based on quantum information tools in the recent past namely, observing the Bell violation [37,45,46], quantum discord [36,[47][48][49], quantum stokes parameter [50].In the quantum field theory framework, on the other hand, the connection between squeezing and chaos has been known for quite some time [51,52].The growing interest in quantum information theory further unveiled the idea of OTOC, complexity as an important diagnostic tool to characterize the chaos first observed in the context of information scrambling by Black holes [53][54][55].Utilizing those tools quantum chaotic growth [38,39] of primordial perturbation has also been quantified through diagnostics such as OTOC [56][57][58], complexity [59][60][61], and are shown to have signatures in the correlations [62].
In this work, as described above, we shall make use of those different tools and explore interconnection among resonance, squeezing, and chaos during the preheating phase.We particularly show that the production under the background of periodic inflaton field during preheating is an ideal quantum chaotic system, where the OTOC of phase space operators exhibits an exponential growth ∼ e 2λt in time.Where, λ, the Lyapunov exponent characterizes the chaotic nature of the system under consideration.A few earlier isolated attempts in this direction can be found in [63,64], where the chaotic nature of the preheating phase has been indicated without detailed characterization.We performed a detailed study on this by calculating OTOC as a diagnostic tool for a class of α−attractor E-model of inflaton and different coupling with daughter fields.We extend our exploration further and shed light on how to quantify the possible thermalization temperature of the system in the parlance of squeezing formalism.From the perspective of a particular observer, the squeezed state is found to be thermal.Utilizing this observation, we arrive at an interesting relation between the squeezing and temperature of the system under consideration.At this juncture, we must remind the reader of a deep connection between chaos and thermalization as conjectured by Maldacena-Shenker-Stanford (MSS) [55] and this states that for any thermal quantum chaotic system, the characteristic Lyapunov exponent (λ) is constrained by its thermalization temperature λ ≤ 2πk B T /ℏ.Therefore, for any thermal quantum chaotic system, the temperature is set to have an upper bound fixed by the characteristic exponent λ.The temperature defined from the squeezing is indeed found to satisfy such an upper bound.
The order of construction of this paper is as follows.In Section II, we first give the model of inflaton background specifying different potential model parameters.We consider α-attractor E-model which is one of the observationally favored models in the context of early inflation, and reheating.In Section III, we construct our model in a two-mode squeezed state language.We derive three dynamical equations of squeezing parameter(r k ), squeezing angle(φ ) , rotation angle(θ k ), and the dynamics of these parameters signify the evolution of the quantum state of the produced fluctuation.We also obtain an equivalent expression of occupation number density in terms of squeezing parameter.In Section IV, we introduce the concept of OTOC, one of the diagnostics of quantum chaos, and we construct the quantum counterpart of classical symplectic matrix taking different combinations of field and its conjugate momenta.Finally, we shall calculate the OTOC of the system from the dominant eigenvalue of the matrix.In Section V, we discuss different inflationary parameters along with the dynamics of inflaton as a periodic driving source for three different models.The very dynamics of this driving force make the underlying chaotic nature and thermalization process distinguishable from one background model to another.In Section VI, we introduce two different interactions between the produced fluctuation and background.Here we consider two well-known three-leg and four-leg type interactions in the QFT framework.In the post-inflationary reheating phase, the production of elementary particles is widely studied by taking two different effects, perturbative and non-perturbative effects into account.Here we estimate a minimum bound on coupling strength above which fluctuation will be resonantly produced in a dynamical background.From Section VII, we calculate OTOC in two-mode squeezed state formalism for four-leg type interaction, and emphasize distinct features of its dynamics for three different background driving sources.In Section VIII, we repeat exactly the same analysis of Section VII for three-leg type interaction.In Section IX, we derive a relation between resonant growth index(Floquet exponent) and chaotic growth index(Lyapunov exponent) and also study the behavior of the OTOC spectrum.In Section X, we wish to have a semi-classical visualization of the quantum chaotic system by generating Poincaré section.Furthermore, we give an approximate estimate of the thermalized temperature of the system which is calculated from the squeezed state and Rayleigh-Jeans spectrum.We also check the consistency of these outcomes with MSS conjecture.Finally in Section XI, we conclude the paper stating some possible future direction of this work.

A.
Description of background: Models and equations We will construct our work on the basis of the consideration of a homogeneous, isotropic universe being described by the (observationally favored, spatially flat) Freedmann-Lametre-Robertson-Walker (FLRW) metric, where a(t) is scale factor in cosmic time t.At the end of quasi-exponential accelerated era, called inflation, inflaton field starts oscillating around its potential minimum.Because of the expansion of the background, the oscillation amplitude falls with time, resulting in a gradual decrease of Hubble parameter H = ȧ(t)/a(t) .The dynamics governing the background evolution are well-known Friedmann equations, where GeV is the reduced P lanck mass, and V (ϕ) is the inflaton potential.Solutions of the Eq.( 2), subject to slow roll initial conditions constrained by CMB observation (n s , r), give us the evolution of the background and Hubble scale after inflation.In order to proceed we will consider a class of observationally favored α-attractor E-model with the potential [65][66][67][68], By varying the exponent "n" we can achieve different power law forms of the potential namely, quadratic model(for n = 1), quartic model (for n = 2), and so on.The amplitude of the potential "Λ", which measures the energy content in the inflaton, is constrained by the CMB measurement, and is related to the scalar spectral index n s , the amplitude of the inflaton fluctuation measured as CMB normalization A s = 3 × 10 −9 and tensor to scalar ratio r.The model is favored by Planck 2018 observation [69], where n s = 0.9649 ± 0.0042 at 68% CL and 95% CL upper limit on tensor-to-scalar ratio r 0.002 is obtained as r 0.002 < 0.056 combining with the BICEP2/Keck Array BK15 data.We will be using these observations throughout.The tensor-to-scalar ratio r can be analytically expressed as, derived in [68].Another parameter (α) determines the shape of the potential.The energy scale of inflation related to the parameter Λ can be analytically expressed in terms of CMB parameters as [68].
In our work, we will obtain different backgrounds specifying different values of the exponent "n" in the given potential model (3).

B. Preheating Stage After Inflation
Just after the end of inflation, inflaton field amplitude is expected to be high(ϕ ∼ M pl ) and consequently, the energy stored in this field is also very high.In this stage, instead of behaving like an individual inflaton particle(quanta), the entire field behaves like a coherently oscillating inflaton field.This behavior can also be thought of as a condensate formed by the superposition of all zero momentum inflaton particles and coherently oscillating in the expanding background.Therefore, any other fields directly/indirectly coupled with this coherently oscillating inflaton will experience a periodic driving force, and if the coupling strength is appropriate, the field amplitude will be resonantly enhanced.It is this coupling regime where we will do the detailed theoretical exploration of the phenomena of particle production, resonance, squeezing, and chaos of the system that emerge in the early universe cosmology.Such phenomena is known as preheating, which has been studied quite extensively in the literature (See earlier references [8,9]).However, in this paper, we will redo this analysis from the perspective of quantum squeezing phenomena and its connection with chaos and thermalization.The span of this era is generically very short proportional to inverse decay width and happens within a few e-folds.In this short span explosive particle production takes place through parametric resonance, hence the energy is transferred very rapidly from the inflaton field to the daughter fields.In this paper, we consider the produced daughter fields to be bosonic in nature.This idea of explosive production through parametric resonance was first put forward and thoroughly studied in the phenomenal work by Linde, Kofman and Starobinsky [9].
To proceed toward the quantitative discussion of this production process, we tart with a general Lagrangian containing different interaction terms corresponding to different fields.The form of the general Lagrangian for inflaton (ϕ), daughter scalar (χ), and fermion (ψ) is taken to be, Where, V (ϕ) is inflaton potential, m χ is the bare mass of the produced scalar particle.F (ϕ, R) is the generic coupling term involving the daughter field and Ricci scalar (R) and background inflaton..The last term is the coupling between fermionic and background inflaton with the dimensionless coupling constant "h".Considering different interaction terms, we can describe the explosive production of different elementary particles from a time-dependent background.In our present work, we will mainly concentrate on bosonic preheating.Fermionic one we left for our future work.Our aim is to study the phenomena of preheating in squeezed state formalism and its theoretical implications in understanding non-equilibrium phenomena such as chaos, thermalization.

III. SQUEEZED STATE FORMALISM
In this section, we will build up the analytical structure of the model in a two-mode squeezed state framework which will be the mainstay of our numerical studies later (see earlier attempt [70]).Here, we mainly concentrate on scalar modes which are parametrically excited due to inflaton considering only bosonic interaction Lagrangian , Where √ −g = a 3 (t), m χ is the bare mass of produced particles, and F (ϕ) is a generic coupling function of the background inflaton (ϕ) and specific form of this function will be discussed later in Section VI.We ignore the Ricci curvature coupling for simplicity.Expressing the scalar field (χ) in terms of Fourier modes , the Lagrangian becomes, Using Eq.( 8) we can reach the following dynamical equation of each Fourier mode (χ ⃗ k ) as, In the following dynamical Eq.( 9) there is a damping term, 3H χ⃗ k which is non-zero in expanding background.
After rescaling the field χ ⃗ k and defining , and we obtain, Where Ω k is identified as a time-dependent frequency.It is important to note that due to inflaton, the terms such as F (ϕ), H, are oscillatory in nature.Influenced by these periodic driving force, each mode (X ⃗ k ) will be resonantly produced as particles.In the expanding background, the co-moving number density of such χ particles n k (t) with co-moving momentum mode k can be expressed in a well-known formula [8,9] In order to implement the two-mode squeezed state formalism, we express the Lagrangian (8) in terms of the rescaled field with two independent modes (X ⃗ k , X − ⃗ k ), and in terms of their canonically conjugate momenta, Π X − ⃗ k = Ẋ⃗ k and Π X ⃗ k = Ẋ− ⃗ k the Hamiltonian will be In the squeezed state formalism, the canonically conjugate phase space operators ( X⃗ k (t), ΠX ⃗ k (t)) are expressed terms of time-dependent creation and annihilation operators as Where time-independent frequency is defined as ω k = k 2 /a 2 + m 2 χ .The time dependent creation and annihilation operators, â⃗ k (t) and â † − ⃗ k (t) satisfy the standard commutation relation, â⃗ p , â † −⃗ q = (2π) 3 δ 3 (⃗ p + ⃗ q).Finally using Eq.( 14) in (13) we are left with the desired form of the Hamiltonian of this system ĤX = 1 2 In the above Hamiltonian, the time-dependent frequency Ω k encodes the interaction between the daughter field and oscillatory inflaton.The first part of the Hamiltonian (15) signifies the collection of free oscillators with equal and opposite momenta, and the second part is their interaction which creates or destroys them leading to wellknown two-mode squeezing phenomena between them.Significantly, the classical time-dependent background is responsible for such phenomena.Once we obtained the required Hamiltonian, the Heisenberg equation is derived as, where the new symbols are The Conventional notion tells us that time-dependent Hamiltonian of a system makes the concept of vacuum ambiguous.We shall justify this statement soon.In the present context, time dependence in the Hamiltonian (15) comes through the interaction with oscillating inflaton background as well as expanding spacetime.However to solve the above Eq.( 16), we shall resort to the idea of Bogoliubov transformation, Where, α k (t) and β k (t) are Bogoliubov coefficients satisfying the relation, |α k (t)| 2 − |β k (t)| 2 = 1 and t 0 is an initial time.The role of these Bogoliubov coefficients is to connect the late-time creation and annihilation operators with their counterparts at the initial time t = t 0 .Using the Eq. ( 18) in ( 16) we obtain With a view to determining the expression of two Bogoliubov coefficients α k and β k , we shall introduce two operators following the reference [36,38].The first one is two-mode unitary squeezing operator Ŝk (r k , φ k ) defined as, Ŝk (r The unitarity of the squeezing operator implies the anti-hermitian nature of B † k = − Bk as can be seen from Eq. (21).Two parameters r k and φ k in the Eq. ( 21) are quantifying the amount of squeezing of the quantum state and its squeezing angle respectively.To have complete evolution of the system, the second operator is the rotation operator, Rk (θ k ), defined as, Rk (θ k ) = e Dk with where, D † k = − Dk as it was in Eq. (21).Here rotation operator is characterized by a parameter θ k , called rotation angle.Combining these two operators Ŝk (r k , φ k ) and Rk (θ k ) one can construct a unitary time-evolution operator for every individual scalar mode as, One can easily check that the application of this unitary evolution operator on the creation operator generates its time evolution as A similar equation can be obtained for the annihilation operator by taking the complex conjugate.A comparison between Eq. ( 18) and ( 24) provides one of the forms of two coefficients α k (t) and β k (t) in terms of r k , φ k and θ k in this two-mode squeezed state formalism.The forms are as follows: Such choice naturally satisfies the relation, In light of this squeezed state formalism, we can alternatively define an expression of occupation number density n k which is equivalent to the expression (11).It is worth mentioning that the purpose of defining such coefficients namely Bogoliubov coefficients is to connect the late-time raising and lowering operators with the operators of initial time.In Heisenberg picture, as operator evolution is studied rather than state evolution, a vacuum state corresponding to the operators at initial time will no longer be a vacuum corresponding to the creation and annihilation operators at late time.
One can construct the late-time number operator( N⃗ k (t)) using (24) and conjugate, N⃗ and the vacuum expectation value of this late-time number operator in the so-called Bunch-Davis vacuum defined at the initial time(t 0 ) can be expressed as, The same expression can also be obtained for − ⃗ k mode, ⟨ N− ⃗ k (t)⟩ = sinh 2 r k .Interesting to note that the evolution of squeezing parameter r k influences the time-evolution of occupation number density for various modes ⃗ k.A closer look at the equations ( 24) signifies the fact that the mixing of raising and lowering operators in the evolution equation is solely responsible for a vacuum state in the remote past becoming an excited state in the distant future, hence justifying our previous statement regarding the ambiguity of the concept of vacuum for time-dependent Hamiltonian.
A. Two-mode squeezed states for bosonic particles We have already introduced key concepts of squeezed state formalism, and we proceed to identify the squeezed quantum state.It is a well-known fact that upon quantization, the states of a harmonic oscillator are coherent in nature and the states of a parametric oscillator(an oscillator having time-dependent frequency) are squeezed in nature.If one measures the uncertainty of two-phase space variables in a squeezed state corresponding to a particular system, then he will encounter the terms whereupon one grows exponentially and another one will die out exponentially [43,71].This can be clearly observed in an inverted harmonic oscillator(IO) state which is a beautiful example of a squeezed state.The Hamiltonian of IO having potential unbounded from below is given below: By defining x and p in terms of raising and lowering operators of the non-inverted oscillator, the inverted oscillator Hamiltonian becomes It can be easily shown that the squeezing of the quantum state for an inverted harmonic oscillator is due to the (a † ) 2 term.We notice that the second term of the Hamiltonian (15) is the same as that of the inverted oscillator Hamiltonian in its field theoretic form.Hence, we can infer that the presence of background driving force due to inflaton, the Hamiltonian (15) induces squeezing in the system.Introducing the fundamental concept for the well-known inverted harmonic oscillator, we shall now proceed towards the determination of the squeezed state by the application of the unitary evolution operator (23) on a two-mode unsqueezed vacuum state.Borrowing the idea as outlined in [36], we can factorize the full Hilbert space of the system ε into independent inner products of Hilbert spaces for two opposite momenta modes ⃗ k and − ⃗ k, ε = k∈R 3+ ε ⃗ k ⊗ ε − ⃗ k .Using this idea, we can write two-mode unsqueezed vacuum state as, Finally the application of two-mode unitary evolution operator (23) on the unsqueezed vacuum leads us to the momentum-preserving two-mode squeezed state1 : where two-mode excited state is given by This |ψ sq ⟩ ⃗ k,− ⃗ k is a two-mode squeezed state, well known in the context of quantum optics as an entangled state.
B. Dynamical equations of r k , φ k and θ k : According to the equations ( 25), we have the expression of two coefficients α k (t) and β k (t) in terms of r k , φ k and θ k .Plugging these expressions into the equations (19), one obtains three dynamical equations of r k , φ k and θ k .The utility of the solution of these dynamical equations lies in the study of not only the time evolution of the system but also the squeezed state (31).Those equations are as follows: Having analyzed the above three equations, we come to know that none of the dynamical equations of r k , φ k , and θ k depend upon the rotation angle θ k .So while studying the dynamics, the system will not exhibit any such sensitive dependence on the rotation angle θ k .

IV. OTOC: A DIAGNOSTIC OF CHAOS
OTOC, a nice acronym of out-of-time-order correlator is one of the diagnostics of specifying the presence of quantum chaos in a system.Unlike a classical system, studying and characterizing the chaos in quantum many-body systems is pretty challenging.For a classical system, if a small perturbation to the initial state of a system causes exponentially diverging phase-space trajectories, then such hypersensitivity to the initial state or initial condition of that system is a clear indication of the presence of chaos.This hypersensitivity can be quantitatively characterized by the Poisson bracket between position (q) and momentum (p) at unequal time: Where, λ n , the Lyapunov exponents are known to quantify the measure the chaos.In a quantum mechanical system, an analogous quantity can be defined between two operators with the classical Poisson bracket replaced by a Commutation bracket.The quantum mechanical analog of ( 34) is generically called OTOC q(t), p(0) , which behaves as a Poisson bracket in the semiclassical limit.Likewise, the exponential diverging trajectories in the classical case, the exponential amplification of this commutator is a potential signature of quantum chaos.
To quantify OTOC, the common practice is the use of a double unequal-time commutator Where, angle brackets ⟨...⟩ β stands for thermal average and β being the inverse temperature, where k B is the Boltzmann constant.Though we are taking commutation between q(t) and p(0), generically one can consider commutation between any pair of hermitian operators Ŵ , V , Therefore, according to the above discussions, a quantum chaotic system is also expected to exhibit exponential growth in OTOC, C(t) ∼ e 2λt , in a similar fashion to the classical one and we can define the analogous quantum Lyapunov exponent(λ) measuring chaos in the quantum system under consideration.In this work, we will calculate OTOC in the two-mode squeezed state parametrized described before for different classes of periodically driven systems in the context of preheating.For this we express the field and conjugate momentum (see Eqs. ( 14)) in terms of three squeezing parameters (r k , φ k , θ k ) as, Using the above conjugate variables the commutation relation assumes the following simple form, where f k (t, t 0 ) is defined as with (r 0 , θ 0 , φ 0 ) being the values of the squeezing parameters at initial time t 0 .The advantage of using such a complex conjugate pair is its being a c-number operator (Eq.( 38)).Using this property for every canonically conjugate pair of field variables in momentum space, we can straightforwardly compute the OTOC as Clearly, the appearance of the delta function advocates the momentum conservation in the system which is reminiscent of a homogeneous FLRW background.Ignoring those delta function terms, paying attention only to the amplitude part, we have In our subsequent studies, we will mainly focus on this amplitude part for various models.We have already mentioned that not only the commutation between position and momentum operator but OTOC, in general, can be taken as a double-commutator between any two hermitian operators.Hence, all possible combinations between field and conjugate momentum can capture the behavior of OTOC.We, therefore, construct a 2 × 2 matrix, M kk ′ , which is the quantum counterpart of classical symplectic matrix having unit determinant consisting of all such combinations of field momenta [38,72], as Where assuming the initial squeezing r 0 → 0, all the components of the c-number matrix are calculated to be From the symplectic matrix, we shall then extract the behavior of OTOC from the dominant eigenvalue of this matrix, and extract the information of "Lyapunov exponent"(λ) as discussed before.

V. INFLIONARY PARAMETERS AND ANALYTICAL SOLUTIONS
To proceed further towards the OTOC calculation, first, we need to prepare the background which represents a coherently oscillating inflaton field, and the Hubble scale which plays the role of the periodic driving force for the daughter fields.A reminiscence of the Eq. ( 2) enables us to determine the field dynamics as well as the Hubble scale numerically.As we are interested in three different equations of state, we study the background dynamics for three different n = 1, 2, 3. To study background dynamics numerically, the specification of proper initial condition of field value just after inflation demands special importance.During inflation, the inflaton satisfies usual slow-roll conditions.The usual condition for the end of inflation is set by one of the slow roll parameters ϵ ∝ V ′ /V 2 to be unity.In the context of α-attractor E-type potential model, we can derive an expression of the field value for general n at which inflation ends is where "end" stands for the value at the end of inflation or at the beginning of reheating and this will be used as an initial condition in our full numerical solution of the inflaton in the post-inflationary phase.Using this field amplitude in (3), we have the potential at the inflation end as A. Analytics solutions during preheating In our work, we shall use the α-attractor E-type potential potential model described earlier.Expanding the potential (3) around its minima we get the approximated power-law type potential for general "n", as In an expanding cosmological background, as the amplitude of the inflaton(Φ(t)) falls with time and generically assumes ϕ/M pl << 1 for significant duration of the reheating phase.Hence, we can safely neglect the higherorder terms in the expansion in the potential.In the post-inflationary era, inflaton field oscillates around the minima of this approximated power-law potential (46) with a decaying amplitude Φ(t).Oscillatory behavior is strongly influenced by the nature of the potential near its minimum.At this stage, the generic solution of the first equation of ( 2) can be expressed as, ϕ(t) = Φ(t)P(t) = Φ(t) ν̸ =0 P ν cos(νωt), where Φ(t) is timedependent inflaton amplitude in the dynamical background and P(t) is a quasi-periodic oscillatory function with the fundamental frequency for the periodic inflaton potential in Eq.( 46) [21], Note that the oscillation frequency is an explicit function of decaying inflaton amplitude.The energy-density, ρ ϕ and pressure, p ϕ , of homogeneous scalar field inflaton can be written as Subject to the approximated power-law potential (46), ignoring background expansion and taking oscillation average over one complete oscillation, we have average energy density ⟨ρ ϕ ⟩ = V (Φ), and pressure as ⟨p ϕ ⟩ = w ϕ ⟨ρ ϕ ⟩.Where, the average post-inflationary EoS is expressed as, w ϕ = (n − 1)/(n + 1).In order to obtain decaying amplitude Φ(t) for general "n", we utilise oscillation averaged energy-density and pressure in (48) together with (46) and the form of average EoS w ϕ in the continuity equation, ⟨ ρϕ ⟩ + 3H(1 + w ϕ )⟨ρ ϕ ⟩ ≃ 0, this yields, Where ϕ end and a end are inflaton amplitude and scale factor at the end of inflation respectively.Once the timevariation of inflaton amplitude is found for general "n", using (49), we can express the leading order behavior of Hubble scale in terms of decaying amplitude and quasi-periodic oscillatory function for general EoS as follows. with being the inflationary energy scale at the end of inflation.We recover the well-known leading behavior of the Hubble parameter ∼ a −3n/(n+1) .As an example, for quadratic potential with n = 1, we recover H ∼ a −3/2 .Substituting ( 49) in H 2 ≃ V (Φ)/3M 2 pl , we compute the leading order time variation of scale factor after inflation for general EoS.

1.
For n = 1: This is a model whereupon expanding the model potential about the minimum (i.e. at ϕ = 0) leads to quadratic power law potential (V ∼ ϕ 2 ).From the expression (46) for n = 1, we can obtain the inflaton mass m

√
3αM pl according to (52) which can be constrained by the CMB parameters(n s , A s ) through "Λ".The form of that well-defined mass term is m ϕ = 2Λ 2 / √ 3αM pl and subjected to the the central value of n s = 0.9649 and A s = (2.100± 0.030) × 10 −9 we get m (1) eff ∼ 1.23 × 10 −5 M pl .It can be shown both analytically and numerically that at the end of inflation inflaton amplitude ϕ end ∼ M pl , the condition 2/3α(ϕ end /M pl ) < 1 gives α > 0.5.The latest observation suggests α ≲ 10.For our numerical purpose, we assume α = 1.The model corresponds to an average EoS w ϕ ≃ 0 and hence the background inflaton field behaves like non-relativistic matter.An approximate analytic solution can be recovered from Eq. ( 49) and (50).Here field ϕ is taken in a unit of M pl and we will be using these dimensionless time-variable z1, z2, z3 throughout this work for n = 1, 2, 3 respectively.
For n = 1, non-negligible contribution comes from only the first Fourier component P 1 , which is From the above, we recover the leading order solution a(t) ∝ t 2/3 implying a dust-dominated universe.The analytic solution can indeed be observed from the numerical solution presented in the left panel of Fig. 1.In our numerical study, we have got the value of inflaton field ϕ end = 0.94M pl and Hubble scale H end = 4 × 10 −6 M pl at the end of inflation.

For n = 2:
For this case we have quartic power law potential V ∼ λϕ 4 around the minimum, with λ = (16Λ 4 )/(9α 2 M 4 pl ) ≃ 1.99 × 10 −10 with α = 1.Approximate analytic solution for n = 2 obtained from Eq. ( 49) and (50) as We numerically obtain |P 1 | = 0.47 and |P 3 | = 0.02.In terms of these Fourier components, P(z 2 ) takes the form P(z 2 ) We indeed recover the leading order solution for the scale factor a(t) ∝ t 1/2 , which corresponds to average EoS w ϕ ≃ 1/3 implying radiation-like behavior.The numerical solution for inflaton is presented in the middle panel of Fig. 1.The value of inflaton field is obtained as ϕ end ∼ 1.47M pl and the Hubble scale H end ∼ 3.6 × 10 −6 M pl at the end of inflation.
3. For n = 3: For this case we have sextic type potential V ∼ ϕ 6 , and approximate analytic solution derived from Eq. ( 49) and ( 50), are We numerically obtain |P 1 | = 0.46 and |P 3 | = 0.03.In terms of these Fourier components, P(z 3 ) takes the form We again recover the leading order scale factor a(t) ∝ t 4/9 , implying a stiff background with w ϕ ≃ 1/2.The numerical solution for the inflaton is presented in the third panel of Fig. 1, and we obtained ϕ end ∼ 1.8M pl using (44) and Hubble scale H end ∼ 3.5 × 10 −6 M pl at the end of inflation.Our study will be based on two different interactions, four-leg and three-leg type, and for each of these two, we will use three different background equations of state corresponding to n = 1, 2, 3 in the potential (3).

VI. INTRODUCING INTERACTION
Now we are on the verge of applying the above formalism in the context of the early reheating era, wellknown as preheating.During this stage, the inflaton field sets the background as a periodic driving force as just discussed above, and acts as a source for the massless daughter field.Our main objective of this paper would be to understand in detail the dynamical structure of this massless production process in order to gain some insight into the quantum chaos under periodic driving force in the cosmological setting.
In the present work, we are going to deal two types of interactions: i) Three-leg interaction and ii) Four-leg interaction with the massless daughter field say χ, where (σ, g) are the coupling constant.In the present study, the generic coupling function of inflaton, F (ϕ), will have the forms F (ϕ) = σϕ, g 2 ϕ 2 corresponding to three-leg and four-leg type interaction respectively.In the usual framework of perturbative field theory such interaction leads to the decay of inflaton ϕ → χχ, ϕϕ → χχ.
We consider coupling parameter regime where the χ-particle production will be non-perturbative [9,12,16].However, such non-perturbative production stops as time progresses.In our following discussion, we discuss the important condition under which the system will be in the non-perturbative resonance regime with the aforementioned inflaton interaction.
A. Dynamical resonance condition: stability-instability chart While constructing parameter space to study resonance phenomenon during the early reheating phase, one usually considers static background (a(t) = 1) neglecting the red shifting of momentum mode because of background expansion.Depending upon a specific background we are to construct this parameter space.The generic mode equation assumes the following form, As stated earlier the oscillating inflaton field is expressed as ϕ(z n ) = ΦP(z n ) with amplitude, Φ.This is the well-known generalized Hill equation [73].The solution of such an equation is strongly dependent upon two parameters (k, q).Where the emergent q-parameter, identified with the term associated with the purely oscillatory part of the equation is called resonance strength parameter, and in a static background those are, Employing Floquet theorem [73,74], solution of such equation shows exponential instability X k ∝ exp(μ k z n ) over time for a specific range of parameter values (k, g) forming a band for which Re μk > 0, where μk is the rescaled Floquet exponent, μk = 2πµ k /m (n) eff .Generating a parameter space we can have those sets of different values of the two parameters(k, g) from the unstable region for which the field solution will exhibit resonant instability.Consideration of a static Minkowskian background shows a clear distinction between stability and instability region in the parameter space (see Fig. 2).Any parameter set chosen from the instability region will lead to an unstable solution interpreted as resonant particle production, and vice versa.In reality, though, the background is dynamic, and hence parameters particularly the k mode instead of staying in a specific stability/instability band may transit from one band to another with the expansion.Thus in the case of expanding background, the distinction between stability and instability region disappears which is a familiar characteristic of Stochastic resonance.Because of the redshifting of each k mode as well as small coupling strength g, if we choose any mode from a narrow resonance regime in parameter space, we will not be able to see the resonance phenomenon as it turns out to be almost instantaneous.Therefore, to observe the exponentially growing solution in time, we need to choose a parameter set from a broad resonance regime with sufficiently high coupling strength g.Employing the usual Floquet analysis for the periodically oscillating background, we have constructed these parameter spaces as shown in Fig. 2 for three different models corresponding to background n = 1, 2, 3 for two different interactions.The structural difference in parameter space for two different interactions is very evident from the given figure.Tachyonic instability [75,76], due to switching sign of P(z) for tri-linear coupling, causes a very efficient and strong resonance, and because of that, in the lower panel of Fig. 2, we notice a very broad instability region in the parameter space for the three given model in comparison with quartic or four-leg interaction.
In the discussion so far we have talked about distinctive properties of flat-space resonance and resonance in the non-static background.But one thing that is not specified yet, is the precise condition of resonance for any general n.In the literature, conditions of resonance are usually specified by [9][10][11], considering q > 1, and that leads to broad resonance.Otherwise, q < 1 gives production in a very narrow resonance regime.Strictly speaking, such a condition defined in flat space is not appropriate, particularly when the driving force namely the inflaton dilutes very fast, Φ ∝ ϕ end (a/a end ) −3/n+1 depending on n values.In the background of such a fast-falling driving force, the resonance parameter q also reduces very rapidly in time as follows For those cases, the naive Minkowski resonance condition indeed overestimates the lower limit of the inflaton coupling parameters above which broad-resonance could occur2 .For the present case with a rapidly falling q-parameter, therefore, we propose a dynamical condition of resonance as follows: Resonant particle production is intimately tied with the violation of the adiabaticity condition, and that is expected to occur while inflaton crosses zero during its course of oscillation.To achieve significant resonant production within a certain period one needs to satisfy two important conditions.The first one would be to have the driving force, the oscillating inflaton, executing a few oscillations within the period of interest.The second condition is that within that period the resonance q-parameter should remain greater than unity.Combining the preceding two conditions, we state that for broad resonance to take place while the resonance parameter q time-evolving from its higher initial value to unity, it must complete at least one oscillation.Bearing this dynamical condition of broad resonance in mind, we derive the lower bound of the inflaton couplings for general background EoS.
To compute the number of oscillations required for q-parameter to change from its large initial value at the end of inflation to unity, we measure the dimensionless time-period of the oscillating inflaton as z where ω is calculated at Φ = M pl .With respect to this, the number of oscillations say N osc becomes, This dimensionless time-period z (ω) n matches very well with the numerically calculated one.Hence, to achieve efficient resonant production in a broad resonance regime, the minimum criterion that must be fulfilled is N osc > 1.Having this criterion, from (60) we get the lower bound of the coupling parameter for two interactions as [26] g > m It is very interesting to note that the lower limit of the coupling strength g to achieve efficient resonance also depends upon the time period of background oscillation besides its dependence on ϕ end and m The above condition sets the minimum possible coupling parameter that one requires at the beginning of the reheating to have efficient broad parametric resonance.However, due to expanding background, any particular resonant mode k also becomes non-resonant around its characteristic time scale.For example, in the present context, minimum dimensionless coupling strengths which are necessary to initiate an efficient resonance process are found to be g = σ/m (n) eff = 6 × 10 −5 , 4 × 10 −5 , 3.77 × 10 −5 for n = 1, 2, 3 respectively.It is important to mention that these lower bounds of coupling parameters are obtained for parametric resonance cases.But for three-leg type interaction, efficient production happens for coupling strength lower than the given bound because of Tachyonic growth present in the system.In the next section, while studying OTOC, we shall choose coupling parameters obeying the bound (61).

VII. SQUEEZING AND OTOC: FOR
In this section, we will concentrate on the study of chaos by explicitly calculating OTOC for n = 1, 2, 3 for four-leg or quartic interaction (g 2 ϕ 2 χ 2 ).The daughter quantum field χ which experiences periodic driving force due to a coherently oscillating inflaton field, will be observed to undergo resonant amplification.Such resonant amplification will be shown to be chaotic in nature.As discussed before, we compute OTOC on the time-evolving squeezed quantum state and show that such quantity will exponentially grow with time.The inverse time scale of such exponential growth will be identified as the Lyanupov exponent.We shall discuss the chaotic scenario and its multifarious characteristics with minute details for three different background dynamics as constructed in the last section separately.To this end, we would like to point out that in our present analysis, we will confine ourselves strictly to the regime where back-reaction is negligible.

A. Initial condition specification:
To perform a detailed numerical analysis the key dynamical equations are given in (33).For the appropriate initial conditions, we first write down the field mode in terms of three squeezing parameters.Replacing Eq. ( 14) in the mode expansion form of rescaled field X(t, ⃗ x) like (7) and also using (24), we get the following form The general field mode function can further be expressed in terms of squeezing parameters r k , φ k , θ k as, The proper specification of the initial condition can be associated with the well-known Bunch-Davies for each individual mode fluctuation.We can assume the solution of Eq. ( 10) by WKB ansatz of the form Now, we can choose the form of α k (t) and β k (t) as with the condition to achieve Bunch-Davies vacuum in the field mode X k and it's derivative is Comparing Eq. ( 66) with (25) at the initial time t = 0, we get the initial conditions for r k (t = 0) = 0, and θ k (t = 0) = 2lπ, where l is any positive integer including zero.We have the freedom to take any value of φ k (t = 0) since all the physical quantities are independent of squeezing angle φ k .Such choices essentially make the initial condition for the initial quantum state to be an unsqueezed one, and that is equivalent to the choice of Bunch-Davies vacuum.Once initial conditions are properly specified we are now in a position to proceed with our numerical studies.

B. Computing OTOC for n = 1
We shall now present our numerical results for r k , φ k , θ k solving three first-order coupled differential equations (33).We will be interested in those parameter (k, g) that belong to the instability region of the stabilityinstability chart (see Fig. 2), and hence the resonant solution exists.It is a well-known fact that long-wavelength modes (IR modes) are efficiently amplified during the preheating era.For pictorial representation and comparison, we have chosen three different long wavelength modes and inflaton daughter field coupling constant g = 5 × 10 −4 .With those parameter choices, we solve for the squeezing parameters shown in Fig. 3.During the initial phase of the evolution, the squeezing parameter, r k grows linearly in time z 1 and subsequently saturates.As expected, the slope of the growth increases with decreasing momentum.Long wavelength mode takes longer time to saturate, that depends on the decreasing amplitude of the background driving force due to Hubble expansion.All these features can indeed be seen in Fig. 3. Subject to the variation of r k as shown in Fig. 3 and using Eq. ( 27), we obtain occupation number density, n k (See Fig. 4), grows exponentially in time for all the three modes under consideration.However, generic dynamical features of squeezing (r k ) and number density ln(n k ) are that they initially grow and reaches a peak value and finally the system relaxes to saturation with small periodic fluctuations.Furthermore, as expected the long wavelength modes are easily excited than those of the higher momentum mode.These features nicely resemble the thermalization process of any thermodynamic system under perturbation.For every individual mode, the thermalization associated with the saturated region can be said to be achieved at their respective time scale.Long wavelength modes thermalize faster in time than short wavelength ones.From Fig. 3 we obtain different values of saturation time scale for three different modes having different wavelengths.For k = m eff , we have z 1 ∼ 9.9, for k = 3.5m eff , we have z 1 ∼ 10.7 and for k = 5.5m (1) eff , we have z 1 ∼ 12.For the present system, the reason behind such variation of saturation time scale is hidden in the dynamics of the background periodic driving force due to the coherently oscillating inflaton field.The dimensionless resonance parameter which controls the dynamics is q 2 g = (g/m (1) eff ) 2 Φ(t) 2 .Where, Φ(t) is decaying as ∼ a −3/2 (see Eq.( 53)).The condition of stochastic resonance in expanding background must satisfy [9], q 4  g m eff ≳ H, we have Considering CMB observed value of scalar spectral index n s = 0.9649 ± 0.0042, we have H end ∼ 10 −6 M pl .Post-inflationary dynamics leads H/m eff )(z end 1 /z 1 ) very small throughout the entire phase with m (1) eff ∼ 10 −5 M pl .On the other hand, the amplitude of Φ also decreases, Φ ∝ z −1 1 with the background expansion.Therefore, we have competition between the aforesaid two effects.As the dependence of the above condition (69) on H is weak (H 1/4 ), It has been observed that after a certain number of inflaton oscillations the above resonance condition (69) fails to satisfy and eventually resonance process ceases.This causes the saturation in all the physical quantities namely occupation density (n k ), squeezing amplitude (r k ), etc as shown in Fig. 3 and  Fig 4.
So far we discussed in detail the dynamics of quantum states of the produced particles in terms of their squeezing parameters and phenomena of chaotic resonance in terms of the field modes.In the following, we will study in detail how such non-trivial behavior of resonance is encoded into OTOC considered as an interesting diagnostic tool of chaos.We first numerically construct the symplectic matrix in terms of commutators of field mode and its conjugate momentum.From the dominant eigenvalue of the calculated symplectic matrix (42) we identify the OTOC dynamics which is parametrized by the Lyapunov exponent(λ k ).
For a quantum chaotic system the OTOC measured in terms of the function C(t) is expected to behave as C(t) ∼ e 2λt with λ being the quantum Lyapunov exponent.Therefore, the slope of (1/2)ln(OTOC) in time will give us the information of Lyapunov exponent λ.In our context we have been using dimensionless time-variable eff t/2π, and the associated dimensionless Lyapunov exponent will be λ ≡ 2πλ/m (1) eff which implies C(z 1 ) ∼ e 2 λz1 .In Fig. 4, we plotted (1/2)ln(OTOC) vs z 1 , and linear growth in time can be clearly observed for a finite period of time.Therefore, the initial chaotic nature is manifested through the linear logarithmic growth of OTOC of a quantum field mode which is under the periodic driving force.Our chronological discussions so far seem to suggest that for a quantum chaotic system squeezing of quantum state parameter r k , logarithmic growth of occupation number density ln(n k ) and logarithm of highest OTOC amplitude have some intriguing relation, which we will discuss later.Furthermore, as expected chaos remains present in the system until the resonance ends or saturation begins which essentially hints probably the obvious fact that chaos and instability are intimately tied with each other [63].Now our task is to determine the LE corresponding to different long wavelength modes.As one observes from Fig. 4, growth of 1 2 ln(OTOC) can be fitted by a straight line.We have the expression C(z 1 ) = Ae 2 λz1 such that ln C(z 1 ) = lnA + 2 λz 1 ⇒ 1 2 where A is any proportionality constant.For a particular k mode the λk = 2πλ k m (1) ef f in above expression (70) has many different zeros (red colored points) as in FIG. 5.For a particular k mode growth of OTOC depends upon a particular value of the Lyapunov exponent.However, Fig. 5 shows a time-varying nature LE, and we need to determine an effective value of this LE.To do so we prescribe the following procedure: We numerically compute the oscillation average of the function between two alternative zeros as marked by point P 1 at z 1 = 5.5 and P 2 at z 1 = 6 in Fig. 5. Repeating this process from the initial time when the growth of a mode starts to the time when it saturates, we define the averaged effective LE ( λeff k ) as,

⟩| oscillations Number of oscillation
Using this procedure, we compute effective LE for various k modes fixing the coupling strength g and also for different coupling parameters g with fixed k mode (see Table I).Finally in Fig. 9 we have given pictorial representation of Lyapunov spectra, that is λeff eff and λeff k vs g.From the figure, it is evident that the variation of LE with k/m (1) eff is not monotonic but rather chaotic.In spite of having this chaotic behavior, effectively with increasing k resonance effect decreases as expected.On the other hand in coupling space, there seems to exist an effective critical coupling around which the system becomes maximally chaotic.

C. For n = 2 :
We have explicitly shown all the necessary steps of our study in the previous case for n = 1 model.We present the numerical solution of r k , φ k , θ k using the equations (33) in the Fig. 6.Unlike the model n = 1, here squeezing parameter shows continuous growth over time without any late time saturation.As the particle  number density is tied with the behavior of squeezing parameter r k , the very nature of r k will have an impact on n k and this essentially causes an uninterrupted growth of n k over time.This feature is special for n = 2, when the inflaton field behaves like a relativistic fluid w ϕ = 1/3 rendering vanishing Ricci scalar.Due to this simple fact, background inflaton field can be made to satisfy Ψ ′′ + λΨ 3 = 0 [10] in conformal coordinate, which is independent of expansion.The solution of this equation is an oscillatory elliptic cosine function with constant amplitude.Such a feature is responsible for the uninterrupted growth of particle number over time for n = 2 model as opposed to other n values.In reality, such an uninterrupted growth process can be terminated by taking into account the backreaction of the produced particle into the background (for detail see [10,35].In our present study, we are not incorporating those additional effects.For this case, therefore, in a similar manner we evaluate the Eq.( 43) and substituting the obtained results in (42) we present the behavior of OTOC amplitude in Fig. 8.We further computed λeff k with k and coupling strength g.Taking various k modes and coupling strength g given in the Table II, in Fig. 9, we show the variation of approximate effective LE with dimensionless momentum mode k ≡ k/m (2) eff and coupling strength g.Effective LE takes a prominent peak at some value of k for a given coupling strength, then it gradually decreases with the increase of k.This sharp peak indicates the presence of a particular mode, for which the system is maximally chaotic.For n = 1 model, this peak is less sharp compared to n = 2. Variation of effective LE with g follows the more or less same behavior as it is with k.For a given k, system peaks at a particular coupling, then it falls with the increase of g with some fluctuations.

D. For n = 3 :
This particular value of n induces a stiff fluid of inflaton with an equation of state w ϕ = 0.5.In this background, the fluctuations evolve in a somewhat different manner compared to the models n = 1, 2. The  behavior of squeezing parameters is depicted in Fig. 7.Here we observe that the growth time scale of squeezing parameter r k is very large and occurs over a very large number of background oscillations as opposed to the previous cases.In Fig. 1, it is clearly seen that with the increase of background EoS, the decay rate of inflaton amplitude or in other words, the decay rate of the strength of the driving source gradually falls causing a slow energy transfer from background to the produced fluctuation.As long as the decay of inflaton amplitude takes place only because of background expansion, this slow energy transfer will essentially cause r k amplitude as well as OTOC amplitude(see Fig. 8) to grow very slowly over a long time.
Likewise, n = 2 case, the absence of any saturation of OTOC amplitude within a finite time makes the determination of exact λeff k for various k and coupling g a bit difficult.In order to have a rough estimate of the growth index of chaos, λeff k , for different k and g, we take the first prominent peak just before a plateau type region.For example, in Fig. 8, we take first peak around z 3 ∼ 820 for k = 0.1( See vertical dashed line in the right figure) and around z 3 ∼ 380 for k = 0.4.Sticking to this idea of local saturation, we have given approximate numerical values of λeff k corresponding to a few resonant modes for two different couplings in Table III.We expect that the nature of the variation of λeff k with two parameters( k, g) will remain unaffected by the choice of the aforesaid time scale.Stochastic nature, being an inherent property of the system results in a non-monotonic variation of λeff k with k and g as seen in Fig. 9.

VIII. SQUEEZING AND OTOC: FOR THREE-LEG INTERACTION (σϕχ 2 ) :
In this section, we shall repeat the same discussion as given in Section VII for three different n = 1, 2, 3 taking three-leg interaction σϕχ 2 .Before presenting our numerical results for different models subject to this interaction, we need to give a short discussion on the parametric resonance process for this particular interaction.The structure of the resonance greatly relies on the nature of the interaction.To study the effect of tri-linear interaction on resonance process, the relevant part in the time-dependent frequency of massless case Ḣ (see Eq.10) would be the first two terms.From the evolution of background inflaton in Fig. 1, it is evident that when ϕ < 0 during one half of each oscillation, modes satisfying k 2 < σ|ϕ(t)|a 2 leads to negative Ω 2 k (t).In that case, the well-known T achyonic instability sets in, and the modes experience exponential amplification.Therefore, for linear periodic driving force, the system encounters parametric as well as tachyonic instability which is a distinctive feature of tri-linear interaction.This process is known as T achyonic resonance [75][76][77].Both the instability leads to a very efficient production of fluctuation which causes the preheating phase to end within a few background oscillation.In the later discussion of this section, we will come across this distinctive feature of this particular interaction very closely.
Nevertheless, we will follow the same analysis and the behavior of the three squeezing parameters is depicted in Fig. 10 for different values of n.Both qualitative and quantitative differences can be observed in the smooth nature of the growth of squeezing and OTOC along with its amplitude.For n = 1, saturation time scale turns out to be for k = 1, z 1 ∼ 5.2, for k = 2, z 1 ∼ 6.9 and for k = 3, z 1 ∼ 9.2.For n = 2, due to effective flatness perceived by both background and fluctuations, every mode is expected to have eternal growth in r k .For all the cases the growth seems to be less stochastic compared to the four-leg interaction.This property will make an impact on the chaotic nature of the system.Using Eq.( 42), we also plotted the nature of OTOC to capture the chaotic dynamics of the system for three models(See Fig. 11).In Fig. 11, for n = 1 model, we notice that for k = 1, OTOC amplitude peaks around z 1 ∼ 5.2 and z 1 ∼ 9.2 for k = 3 and then it fluctuates without further growth.Finally employing the same methodology that was described in the last Section VII, we estimate the effective LE λeff k for different values of k with fixed g, and for different g with fixed k (see Fig. 12).In the Table IV    eff ) and coupling strength g for three models.

IX. RELATION BETWEEN FLOQUET EXPONENT AND LYAPUNOV EXPONENT, AND OTOC SPECTRUM
So far we mainly discussed the nature of OTOC and calculated the associated Lyapunov exponents.However, it is understood that such growth is intimately connected with the resonant production of modes.We, therefore, expect to have a direct relation between the Lyapunov and Floquet exponent which are the measure of chaos and resonance respectively under the influence of periodic driving force.To establish an approximate analytic relation between those two quantities, we shall resort to the method of successive scattering on parabolic potential discussed in detail [8,9].Here we quote the main results.If we investigate the resonance process more closely we notice that whenever the background inflaton crosses zero at the minimum, the particle number density shows a sharp growth.Otherwise, particle number density behaves as an adiabatic invariant quantity between two successive zero crossings.Such behavior allows one to obtain an approximate analytical expression of field mode in terms of Floquet exponent around each zero crossing time, say j-th crossing time is defined to be t j .If one consider ϕ 2 χ 2 interaction, near t j , massless field mode Eq.( 10) for n = 1 assumes following form in the parabolic potential around t = t j , where (Φ j ) is the time-dependent amplitude of the inflaton calculated at the instant of j-th crossing t j .For |t − t j | > 0, the adiabaticity condition is satisfied.Using the Bogoliubov method [8,9],one can compute the Where, κ 2 j = k 2 /(a(t j ) 2 gm (1) eff Φ j ), and θ tot j the total random phase accumulated by the field wave from the initial time to the instant t j .In order to obtain the nature of OTOC spectra as well as the possible relation between Floquet and Lyapunov exponent, we need to evaluate first the commutation relation between field(X k ) and conjugate momenta(Π k ) as outlined in Section IV.We closely follow the reference [9] in deriving the amplitude of the commutation relations using the asymptotic form of the field mode X k as in Eq. (64).We write the expression of OTOC amplitude as follows: Where V is a constant phase, t j = πj m 1 + e −πκ 2 .m (1) We have used Eq.( 72) to reach the final form of the spectra in Eq.( 74).Now associated with the Lyapunov exponent as outlined in section IV and using (64) we write, Lyapunov exponent at the instant of j-th crossing  as, Here the out-of-time order commutation of one combination field and conjugate momenta has been defined between t 0 and the j-th crossing instant t j .The above expression is true at the point of a particular zero crossing t j .We, therefore, can have an effective relation where the average is taken over time from the initial to the saturation time scale.So practically we take average over the number of oscillations executed by the background within the time required for saturation of a particular mode.In the above expression, the phase V is random around π/2. So, after taking the average over the entire saturation time scale for a particular mode, the third term will have a vanishing contribution.For the second term, we have numerically found a smaller contribution compared to μeff k .Finally, we find λeff k and μeff k are more or less the same which is consistent with the numerical result.From Eq.( 74) it is clearly seen that the spectral nature of OTOC is governed by the nature of the Floquet exponent.The growth index of particle number in Eq.( 72) being a function of random phase in expanding background causes somewhat random growth of number density after every zero crossing or every half of a period.This non-monotonic behavior is also reflected in the spectral behavior of OTOC amplitude as shown in Fig. 13.The random oscillatory behavior of OTOC amplitude in momentum space is also caused by the randomness of the Floquet exponent in the context of stochastic resonance as discussed in detail in [9].Out of three different background driving sources, only for n = 1 the system reach a saturation state3 after going through chaotic dynamics within a finite time.While evaluating OTOC spectra numerically for n = 1 model we have chosen the saturation time corresponding to the maximum mode that will be produced in a broad resonance regime.Unlike n = 1 model, no such saturation is observed for n = 2 model as shown in Fig. 8 and Fig. 11(Reason is discussed in Section VII).For this case, we computed OTOC spectra by truncating the evolution at any arbitrary instant of time for each momentum mode.With the increase of potential exponent n, slow decay of the amplitude of inflaton causes slow growth of OTOC amplitude (See Fig. 8 for n = 3), and for this, we have considered a local saturation time corresponding to the maximum mode chosen in the spectra in Fig. 13, where the coupling parameters are chosen from the broad resonance regime (See the lower bound in Eq.( 61)).Despite several fluctuations, one common feature observed in the spectra for all the models is that the higher the momentum, the lower the OTOC amplitude which actually mimics the nature of resonant particle production.The existence of different momentum bands for different values of resonance strength parameter (q g , q σ ) can be observed to be nicely imprinted in the broad oscillatory feature in the OTOC spectrum, particularly for quadratic coupling.However, for tri-linear coupling due to tachyonic and broad resonance, such oscillatory nature smoothens out.

A. Poincaré section: Semi-classical visualization
At the outset of this section, we gain some insight into how the chaotic nature of the system is connected to the squeezing of states.It is known that positive LE signifies the chaotic nature of a system.It actually measures the sensitivity of the system to its initial condition.The well-known method to explore chaotic behavior is to generate and analyze the Poincaré sections, particularly for non-integrable systems.Each field mode of our present system is periodically driven one dimensional parametric oscillator which does not have any conserved quantity.For this system, we identified the semi-classical phase space variables by calculating their expectation values in the squeezed quantum state.We obtain the Poincaré section or a surface of section, by mapping those semi-classical phase space variables projecting onto two dimensions.We generate this Poincaré section by sampling the data points in a 2D plane at a regular interval set by the periodicity of the driving force namely the oscillatory inflaton field.We present Poincaré section for n = 1 model for both types of interactions(threeleg and four-leg), besides this, we also give the parametric plot of X k and X ′ k which carries the information of squeezing of the system for both models(See Fig. 14).So for a non-chaotic system, we expect that all the data points will lie in an orbit and will form closed orbits, but here we see that the phase space points are scattered in the plane and are forming a cloud of points, which signifies the presence of chaos in the system.This typical behavior also tells us about the squeezing of the system.Looking at Fig. 3 we readily understand that the squeezing parameter of this system corresponding to particular parameter values (k, g) grows over a certain duration and then gets saturated, that means after the point of saturation there is no further squeezing in the system.Furthermore, the growth of OTOC in Fig. 4 also exhibits the same behavior.Connecting these two outcomes together, we can state that as long as squeezing is there, the system will remain chaotic which vindicates the statement that highly squeezed states are prime for quantum chaos.In the following section, we would like to shed light on the possible thermalization temperature of such a chaotic system.
B. Defining thermalization temperature and Maldacena-Shenker-Stanford (MSS) bound MSS bound: Chaos and thermalization have been the subject of investigation over a long period of time in non-equilibrium field theory.In this regard an intriguing bound on chaos has been proposed in [55] which we call MSS bound in brief, and the Lyapunov exponent is universally conjectured to be bounded as λ ≤ 2π(k B /ℏ)T , where T is assumed to be the temperature of the chaotic system when it becomes thermalized after the relaxation from chaotic instability.What it suggests the growth of chaos of a thermal quantum system is bounded by the system temperature T .If we reverse this argument, one may arrive at the following universal lower bound on the temperature of a quantum chaotic system T ≥ (ℏ/(2πk B ))λ.In our present system, we apply the proposed bound and estimate the minimum possible temperature that the system can achieve.So far we have discussed the chaotic nature and calculated the value of the associated Lyapunov exponents of a quantum scalar field when coupled with the oscillating inflaton.Further, we have performed our analysis for individual momentum modes and computed the Lyapunov exponent for each mode.Therefore, instead of getting the information about the temperature of the entire system, we will have the information of each mode.Using the bound we may write, Where λ k and T k are LE and the corresponding minimum possible temperature perceived by the mode.A more rigorous estimation would be to do the analysis in real space by using lattice simulation technique and identifying MSS bound.We will defer the lattice studies for our future work.Using the approximate values of effective LE given in Table I and Table IV for ϕ 2 χ 2 and σϕχ 2 interaction respectively, we can identify the lower bound of the temperature of the system by averaging over the resonant momentum band as follows, When the chaotic system reaches thermal equilibrium, the minimum temperature achieved by the system will be given by the following bound in Eq. (78).By exploiting this we obtain the lower bound of temperature for different coupling strengths for two different interactions.Using Eq.( 78), for example, we found TMSS = 2.62 × 10 11 (4.76×10 11 ) GeV for g = 3×10 −4 (5×10 −4 ) in case of g 2 ϕ 2 χ 2 interaction and TMSS = 2.33×10 11 (4.59×10 11 ) GeV for g = 8 × 10 −5 (1.6 × 10 −4 ) in case of σϕχ 2 interaction.We find that the lower bound of temperature is indeed approximately the same for different coupling and interaction as represented in Fig. 15.However, given the intimate relation between chaos and squeezing, can the squeezing of the daughter particle states also encode the information about the temperature of the system under consideration?In the following we conjecture that squeezed state indeed can encode the temperature of the system under consideration as follows: Temperature from Squeezed State: We have already obtained the squeezed state by the application of unitary evolution operator Ûk on an unsqueezed two-mode vacuum in the Eq.(31).For that squeezed state we can define the density matrix for two modes as ρ Since one of the two modes is always inaccessible to the observer, we can calculate the reduced density matrix ρ( ⃗ k) from the full one by tracing out the mode − ⃗ k and get We assume the produced particles are instantaneously thermalized.The above expression of a quantum state can be identified as a thermal squeezed state with inverse temperature defined as [36] We can therefore identify this as a system temperature under consideration.As the squeezing parameter r k is a function of time, βk will also vary with time.The very nature of r k shows a prominent peak near the saturation time scale as seen for n = 1 model observed in Fig. 3 and Fig. 10.The temperature of the system being tied with r k will be maximum at an instant, and for different k, this temperature peak position will be different.
Collecting only those maximum values, and following the definition of effective MSS bound stated in Eq.78, we define the average system temperature defined for the thermal squeezed state (SS) over resonance band as We assumed Boltzmann constant k B = 1 in natural unit.Our motive is to check the consistency of the variation of this average temperature for different coupling with the MSS bound and we anticipate that this temperature variation will abide by the MSS bound for any coupling strength.In the absence of any kind of self-interaction or back-reaction, as late time saturation of r k is obvious in n = 1 model for both the interaction, we concentrate only on n = 1 case in order to define an effective temperature of the system.Using the relation (82) for two different interactions, we determine an approximate effective temperature that the system can achieve after saturation for different inflaton coupling.For example TSS = 6.75 × 10 14 GeV for g = 2 × 10 −4 in case of ϕ 2 χ 2 interaction and TSS = 4 × 10 13 GeV for g = 3 × 10 −5 in case of ϕχ 2 interaction, which are greater than the MSS bound we calculated before.We indeed found the MSS bound to be satisfied for different inflaton coupling and parameters as depicted in Fig. 15.To this end, we should remind the reader that the temperature defined above is not the reheating temperature that is defined at the end of reheating.
Rayleigh-Jeans Temperature: We would like to stress another consideration in the context of thermalization of a system having a very large occupation number.In this limit when the system is in thermal equilibrium the Rayleigh-Jeans spectrum can be defined as n k E k ≃ T k , for occupation number n k and associated energy E k corresponding to a mode k mode [31,35].In this equation, T is again the temperature of the system under consideration.In our present context energy per mode can be calculated to be using the Hamiltonian H X of the produced massless fluctuation in Eq. (15).As the occupation number is usually very high during the preheating phase, to have a better understanding of the thermalization across all the resonant modes for a particular coupling, we can use the combination n k E k .Following the Eq.( 82) here also we can define an average temperature using the Rayleigh-Jeans(RJ) formula as Most interestingly, we find that the equilibrium temperature predicted by the thermal squeezed state in Eq.( 82) matches with the thermalized temperature of the system predicted by the Rayleigh-Jeans formula given in Eq.( 83) pretty well as shown in Fig. 15 by red and green dots accordingly.This indeed establishes a beautiful connection between the non-equilibrium squeezed quantum system with the non-equilibrium preheating phase in the early universe.
Since our present analysis is far from complete, and hence the estimation of different temperatures may not strictly be taken as actual system temperature.At this juncture, therefore, it would be worth pointing out the standard expectation, and we would like to close this section by indicating a naive estimation of the maximum possible temperature the system may attain right after the end inflation, and that can be calculated by utilizing the following relation H 2 end = 1/(3M 2 pl )ρ = π 2 g/(90M 2 pl )T 4 with g ∼ 107 being the effective number of degrees of freedom.In the expression, we considered Stephan-Boltzmann law temperature of radiation energy density ρ = (gπ 2 /30)T 4 .Assuming H end ∼ 4 × 10 −6 M pl , one gets T ∼ 2.6 × 10 15 GeV which anyway falls within the estimation displayed in Fig. 15.The over and under estimation of the temperature with respect to the perturbative temperature just mentioned above can be evaded by the full lattice simulation which we consider for our future study.

XI. CONCLUSIONS
Understanding the non-perturbative dynamics of particle production is of immense importance, particularly in the quantum field theory framework.In the realm of cosmology, the early universe reheating phase is an interesting laboratory to understand such phenomena.During the reheating phase, inflaton plays the role of periodic driving force which may lead to resonant particle production if the appropriate conditions are satisfied.Extensive exploration has been performed to understand such non-perturbative phenomena with the motivation to successfully reheat our universe.In this submission, we re-explored such non-perturbative production in the preheating stage and analyzed its underlying physics in relation to quantum squeezing, chaos, and their intriguing connection with thermalization.To the best of our knowledge along that direction very little attention has been paid (see early references [27][28][29][30]) mainly because of a lack of understanding of the possible connection between the reheating and its subsequent evolution to standard cosmology.However, preheating can be assumed as an ideal phase of non-equilibrium processes which demand very special attention from the point of view of theoretical understanding of early universe cosmology.The mechanism of thermalization is of great theoretical importance on its own which is not very well understood.Particularly during reheating, it is believed to have played a significant role in some well-studied cosmological phenomena such as producing matter-anti-matter asymmetry (baryogenesis), and dark matter abundance that we observe today.In the context of perturbative reheating and Lattice studies, the issue of thermalization has been discussed [31,32,34,35].However, a theoretical understanding of this is still far from complete.In this paper for the first time we initiate a theoretical exploration of the phenomena of squeezing, chaos, and their deep connection with thermalization during preheating, and we advocate through this submission that we have some promising results that are worth exploring with greater detail in the future.For our present study, we have considered two well-studied inflaton interactions with the daughter field namely ϕχ 2 and ϕ 2 χ 2 .After the inflation, the inflaton field oscillates near the minimum of its potential ϕ n and acts as a quasi-periodic driving force leading to chaotic growth of quantum fluctuation which is observed to be imprinted in the OTOC of phase space variables.We have computed the associated effective Lyapunov exponent ( λeff k ) characterizing the chaos in the system.During this chaotic growth, the quantum state of the produced particles evolved into a squeezed state.We have further explored the underlying connection between the chaotic growth and the squeezing of the quantum states and showed the connection between Lyapunov exponent and squeezing parameters.Thermalization is believed to be deeply connected to the chaotic behavior of a system, which is beautifully conjectured [55] by proposing an inequality relating Lyapunov exponent and system temperature under consideration λ ≳ T .By using this we calculated the approximate lower bound of temperature TMSS .We further conjectured a relation between the system temperature and quantum squeezing averaged over phase space and calculated the temperature TSS (see Eq.82) which indeed satisfies the MSS bound TSS > TMSS .Finally to validate our aforesaid conjecture we resort to Rayleigh-Jeans definition of the temperature of a system being in thermal equilibrium.In the classical limit when the occupation number is very large, a thermalized system satisfies the well-known Rayleigh-Jeans formula nE ∼ T with (n, E) being the occupation number and particle energy accordingly.Surprisingly, the temperature TRJ derived from the Rayleigh-Jeans formula turned out to be TRJ ≈ TSS , which is quite consistent with our definition.To this end, we find that among three different background models, n = 1, 2, 3, only n = 1 appears to have a close resemblance to the thermalization process of a typical thermodynamic system under perturbation.The process of thermalization itself is a complex phenomenon.In particular, how a quantum system initially prepared in far-from-equilibrium states can evolve into a thermal equilibrium state is still not a completely understood subject.Eigenstate Thermalization Hypothesis [78][79][80][81] is a well-known proposal towards this direction.In this context, it would be interesting to investigate this ETH in the reheating era which we have left for our future endeavor.
t/(2π) by multiplying effective inflaton mass parameter m (n) eff at the end of inflation.

3 ϕFIG. 1 :
FIG.1: Figure represents the post-inflationary evolution of inflaton background.Here field ϕ is taken in a unit of M pl and we will be using these dimensionless time-variable z1, z2, z3 throughout this work for n = 1, 2, 3 respectively.
this limiting value varies from one EoS to another depending upon ϕ end , m (n) eff and time-period z (ω) n .

35 FIG. 2 :
FIG. 2: Upper Panel: Figure represents the flat space parameter space for the model n = 1, n = 2 and n = 3 respectively from the left corresponding to g 2 ϕ 2 χ 2 interaction.Lower Panel : Figure represents the flat space parameter space for the model n = 1, n = 2 and n = 3 respectively from the left corresponding to σϕχ 2 interaction.

FIG. 3 :
FIG.3: Figure represents the time evolution of squeezing parameter(r k ), squeezing angle(φ k ) and rotation angle(θ k ) for three different momentum modes(k) and the coupling strength is chosen to be g = 5×10 −4 .Here Squeezing angle(φ k ) and rotation angle (θ k ) are given in units of radian.Here we have taken three different modes for which squeezing parameters grow efficiently giving a considerable production of particles.

FIG. 4 :
FIG. 4: Left Panel:Figure represents the growth of occupation number density n k for three efficient k modes with g = 5 × 10 −4 .Right Panel: Figure represents the time variation of logarithm of OTOC amplitude for two different k modes corresponding to significant production with g = 5 × 10 −4 .Two black dashed lines for two modes indicate the average straight line-like growth of ln(OTOC) and from the slope of those lines we estimate the average growth indices of chaos of the system.

FIG. 6 :
FIG. 6: Figure represents the time evolution of squeezing parameter(r k ), squeezing angle(φ k ) and rotation angle(θ k ) for three efficient k modes, and the coupling strength is chosen to be g = 4 × 10 −5 .Here Squeezing angle(φ k ) and rotation angle (θ k ) are given in unit of radian.

FIG. 12 :
FIG. 12: Figure represents the Variation of λeffk with dimensionless modes k and dimensionless coupling g for three models in three-leg type interaction.

FIG. 14 :
FIG. 14: Upper Panel: Figure on the left represents the parametric plot of X k and X ′ k where k = 1 and figure on the right represents the Poincaré section for n = 1 model for ϕ 2 χ 2 interaction.Different colors correspond to different initial conditions.In the Poincaré section, Blue colored points are phase space points for k = 1 and Red colored points are phase space points for k = 3.Here dimensionless coupling strength g = 5 × 10 −4 .Lower Panel: Figure on the left represents the parametric plot of X k and X ′ k where k = 1 and figure on the right represents the Poincaré section for n = 1 model for σϕχ 2 interaction.In the Poincaré section, Blue colored points are phase space points for k = 1 and Magenta colored points are phase space points for k = 2.Here dimensionless coupling strength g = 8 × 10 −5 .

FIG. 15 :
FIG. 15: Left Panel: Figure represents the variation of effective temperature with coupling strength for four-leg type interaction.Right Panel:Figure represents the variation of effective temperature with coupling strength for three-leg type interaction.It is seen that effective temperatures for different coupling are consistently following the lower bound of temperature for both interactions.Here red dots indicate the temperature obtained through the relation (82) and green dots indicate the temperature obtained through Rayleigh-Jeans formula given in (83) for both the interactions.

TABLE II :
Variation of effective LE( λeff k ) with k and coupling strength "g" for n = 2(g 2 ϕ 2 χ 2 ) , Table V and Table VI, we have given some numerical values of λeff k for different momentum modes k and coupling constant g ≡ σ/m

TABLE III :
Variation of effective LE( λeff k ) with k and coupling strength "g" for n = 3(g 2 ϕ 2

TABLE IV :
Variation of effective LE( λeff k ) with k and dimensionless coupling strength "g" for n = 1(σϕχ 2 )

TABLE V :
Variation of effective LE( λeff k ) with k and dimensionless coupling strength "g" for n = 2(σϕχ 2 )

TABLE VI :
Variation of effective LE( λeff k ) with k and dimensionless coupling strength "g" for n = 3(σϕχ 2 )