Quantum aspects of chaos and complexity from bouncing cosmology: A study with two-mode single field squeezed state formalism

$Circuit~ Complexity$, a well known computational technique has recently become the backbone of the physics community to probe the chaotic behaviour and random quantum fluctuations of quantum fields. This paper is devoted to the study of out-of-equilibrium aspects and quantum chaos appearing in the universe from the paradigm of two well known bouncing cosmological solutions viz. $Cosine~ hyperbolic$ and $Exponential$ models of scale factors. Besides $circuit~ complexity$, we use the $Out-of-Time~ Ordered~ correlation~ (OTOC)$ functions for probing the random behaviour of the universe both at early and the late times. In particular, we use the techniques of well known two-mode squeezed state formalism in cosmological perturbation theory as a key ingredient for the purpose of our computation. To give an appropriate theoretical interpretation that is consistent with the observational perspective we use the scale factor and the number of e-foldings as a dynamical variable instead of conformal time for this computation. From this study, we found that the period of post bounce is the most interesting one. Though it may not be immediately visible, but an exponential rise can be seen in the $complexity$ once the post bounce feature is extrapolated to the present time scales. We also find within the very small acceptable error range a universal connecting relation between Complexity computed from two different kinds of cost functionals-$linearly~ weighted$ and $geodesic~ weighted$ with the OTOC. Furthermore, from the $complexity$ computation obtained from both the cosmological models and also using the well known MSS bound on quantum Lyapunov exponent, $\lambda\leq 2\pi/\beta$ for the saturation of chaos, we estimate the lower bound on the equilibrium temperature of our universe at late time scale. Finally, we provide a rough estimation of the scrambling time in terms of the conformal time.


Introduction
The idea of circuit complexity [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] has recently gained huge attraction of the theoretical physics community and is recently used as a diagnostic for Quantum chaos [19][20][21][22][23][24][25][26][27][28][29]. The absence of a proper tool to develop a wholesome understanding about the AdS/CFT correspondence [30] in certain black hole settings is what motivated the high energy theoretical physics community to apply this computational concept in the context of Quantum Field Theory (QFT). The information about the bulk geometry that can be extracted from the boundary Conformal Field Theory (CFT) remains very much incomplete and is one of the toughest challenges that one faces when probing black hole physics beyond the horizon. One of the main difficulties in boundary field theories is that it reaches thermal equilibrium very quickly while the Einstein-Rosen bridge continues to grow. These challenges motivated Leonard Susskind and collaborators to propose the Complexity=Volume and Compexity=Action conjectures to probe gravity beyond the horizon of black holes and have led to the development of enormous new ideas about the application of complexity and other information theoretic measures in the gravity sector [1,[31][32][33][34][35]. The CV conjecture suggests that the holographic complexity of the boundary field theory is equal to the volume of an extremal codimension one surface extending the boundary time slice into the bulk whereas the CA conjecture suggests that the complexity of a boundary state is dual to the gravitational action evaluated on the Wheeler-DeWitt patch. However the traditional way of computing complexity has certain shortcomings when applied to holography and QFT states. Generally in these contexts, one considers a continuum of states and a proper way to define complexity in this continuum of states faces several questions that need to be addressed. To name some of them, selecting the initial reference state, a set of infinitesimal unitary generators or quantum gates, a proper measure for understanding the role of these gates in minimizing the distance function and the procedure it follows. One of the proposals for facing these issues is to compute quantum complexity using the path length obtained by integrating the Fubini study line element joining the reference and the target state. The reference state is mainly chosen to be Gaussian because the ground states of free field theories are in general Gaussian. For Gaussian quantum states, a geometric way of computing the complexity was given in [36][37][38]. It includes two different methods commonly known as the wave-function approach [2] or the covariance matrix approach [39,40]. The wave function approach has been found to be the most insightful one to probe the underlying physics, especially in the context of time evolution.
Sharing an intimate relation with the Out-Of-Time-Ordered-Correlation functions [41][42][43], abbreviated as OTOC, these two measures has been the recent tools to probe quantum randomness and chaos in various quantum mechanical systems. OTOC's which first appeared in literature in the context of superconductivity [41] soon became popular as a theoretical probe to explore the out of equilibrium phenomenon in finite temperature field theories, bulk gravitational theories and many-body quantum systems. A lot of investigation has followed since then to conclude that whether OTOC's can be considered as a good measure to study stochastic randomness and chaos of quantum systems at out of equilibrium phase. Together with OTOCs, complexity is now considered to be an integral part of the machinery used in the diagnosis of quantum randomness and chaos. Both of these measures have been found to provide information like Lyapunov exponent, scrambling time etc., which are by far the most essential quantities required to comment on the chaoticity of any quantum mechanical system.
In this work, our attempt will be to apply this quantum information theoretic measure to the framework of bouncing cosmological paradigm. Bouncing cosmology is gaining traction to resolve the problem of Big Bang Singularity in recent years . A solid model in bouncing cosmology can resolve the Horizon problem, Flatness problem, the CMB Inhomgeneity and other problems that are prevalent in the current model of Big Bang and Inflationary cos-mology . One way of getting a non-singular ghost free bouncing models is through non-local infinite derivative gravity theories with an addition of appropriate non-local function in the Einstein-Hilbert action in the ultraviolet regime that captures all the derivative terms [111][112][113][114][115]. Moreover non singular bouncing solutions of a positive cosmological constant can make inflation geodesically complete [116]. The primary motivation to apply the formalism of cosmological complexity in bouncing background is that the study of complexity can give great insight about a given model in bouncing cosmology and the explicit calculation of the Lyapunov exponent and the corresponding lower bound on equilibrium temperature [19] during the bouncing period can be very useful in our understanding of primordial cosmology. In this paper, we intend to apply this concept of cosmological complexity under a squeezed state formalism with scalar cosmological perturbations to two well known bouncing solutions -the cosine hyperbolic bounce [116,117] and the exponential bounce [113], which we have derived from usual Einstein gravity with two different models of dynamical scalar matter field embedded in spatially flat (k = 0) Friedmann-Lemaitre-Robertson-Walker (FLRW) cosmological background in 3 + 1 dimensions. However, the exact same solutions can also be derived from higher derivative non-local gravity theory admitting isotropic and homogeneous bouncing universes in the absence of matter [116,117].
We have developed a framework for bouncing cosmology from potentials derived from String theory descriptions at very high energy scale, that can be treated with the squeezed state formalism [118][119][120][121][122][123][124], and using that result the cosmological complexity can be further analyzed. We write a generalized scalar perturbation in the framework of bouncing cosmology and expressed the action, and its parameters including the dispersion relation without truncating higher order terms initially and then give the limiting solutions in the sub-Hubble, Horizon crossing and the super-Hubble regions. The Hamiltonian is also written in its most general form, as compared to [118] before fixing the initial conditions at the horizon crossing scale at kτ = −1 and formulating the squeezed states with a next-to-leading order time dependent slowly varying term in the dispersion relation that we found after the quantization of the Hamiltonian to be more relevant in the context of bouncing cosmology. 1 Other works in cosmological complexity [11,12], have only considered the leading constants in the dispersion relation and squeezed state formalism under the assumptions of stationary background space time. We have then focused our further analysis with bounce in the sub-Hubble region (−kτ 1) to get a better analysis of the quantum fluctuations as compared to the super-Hubble region which falls under the classical domain. This is where the necessary approximations to the dispersion relation is made and the complexity cost functions based on an early general description of family of cost functions is derived. A universality relation between the OTOC and the complexity has also been given under certain conditions. We make certain key observations from our numerical analysis including: • Observation I: Behaviour of squeeze parameters in and around the bounce and at late times.
• Observation II: Initially fluctuating complexity that grows at later times and achieves a saturation at very large time.
1 Note: In refs. [11,12,118] the authors have not considered the slowly varying contribution in the evolution in the sub-Hubble region (−kτ 1) in their computation. During describing the inflationary paradigm all of them have considered the exact de Sitter solution, which is in realistic cosmological analysis is not very useful and also appropriate. The prime reason is using exact de Sitter solution one cannot able to stop the inflation at all in the evolutionary time scale or equivalently in the field space. To stop inflation in an appropriate field space one needs to include slow-roll parameters, which basically considering the small but significant deviation from exact de Sitter solution. When the slow roll parameter reaches the unity the end of inflation is ensured. Evolutionary scale ⌘ ⌧,a(⌧), N < l a t e x i t s h a 1 _ b a s e 6 4 = " d J I 8 1 7 n w O t y 6 7 / 9 2 Q N b B B f 9 k S V s = " > A A A C l n i c d V F b a 9 s w F J b d X b r s 0 m x 9 G e x F L C w k E I L t m a Z 5 2 C g N Z X s a G T R t I Q p B l o 8 T U f l S S Q 4 1 x v 1 H + z N 7 2 7 + Z H G d j 1 w + E P n 3 f O Z y j c 4 J M c K U d 5 5 t l 7 9 2 7 / + D h / q P W 4 y d P n x 2 0 n 7 + 4 U G k u G c x Y K l J 5 F V A F g i c w 0 1 w L u M o k 0 D g Q c B l c T 2 r / c g N S 8 T Q 5 1 0 U G i 5 i u E h 5 x R r W R l u 0 v 3 e 7 5 O x J J y k q 3 K k k A m l Z k B T c / J Y 9 k v G p e I R E J E R D p X k l k j C d p n A m 4 5 b q o i O S r t e 5 X Z d h Y Z 5 t U 5 H U B K o s 7 x a i A q l 8 N C G l 1 u 3 / 7 u A k g c J P z D S a a 5 g N M e / X d H + C S G B N / q u 6 2 W L Y 7 z n B 8 f O T 5 R 9 g Z O s 7 I 9 d y a e C P / r Y 9 d o 9 T o o B 2 m y / Z X E q Y s j y H R T F C l 5 q 6 T 6 U V J p e b M 1 G y R X E F G 2 T V d w d z Q h M a g F u V 2 r B V + Y 5 Q Q R 6 k 0 J 9 F 4 q / 6 a U d J Y q S I O T G R M 9 V r 9 6 d X i v 7 x 5 r q P j R c m T L N e Q s K Z Q l A u s U 1 z v C I d c A t O i M I Q y y U 2 v m K 2 p W Y E 2 m 2 y Z I f z 4 K f 4 / u f C G r j 8 c f / Y 7 J 6 e 7 c e y j V + g 1 6 i E X j d A J + o i m a I a Y d W i N r V N r Y r + 0 3 9 t n 9 o c m 1 L Z 2 O Y f o N 9 j T 7 z p J z K U = < / l a t e x i t > • Observation III: There exists a smooth transition between the non-equilibrium growing phase and the equilibrium saturating phase.
• Observation IV: The saturation at late times indicates a bound on chaos, which makes it possible to describe the Lyapunov exponent and the lower bound of the equilibrium temperature using the well known, Maldacena(M) Shenker(S) Stanford(S) saturation bound on quantum Lyapunov exponent [19], 2

MSS Bound
where T is the equilibrium temperature corresponding to saturation of quantum chaos at the late time scale.
• Observation V: The two different measures used for complexity point to Lyapunov exponent whose fractional deviation is under ten percent, and hence it is safe to assume that our universality relation holds perfectly, in the context of our study.
• Observation VI: We have also very roughly estimated the scrambling time [126,127] for both models and found them to be decent indicator of the time that the OTOC may take to attain equilibrium in both cases.
We expect the bound on quantum chaos and hence the resulting Lyapunov exponent from the two measures of complexity to be much more closer in value by doing the analysis with a full dispersion relation given in the paper. We had initially done the numerical analysis against scale factor for simplicity, but to connect with the observational constraints we have extended the analysis of the complexity in the bouncing background with respect to the co-moving Hubble radius as well, which can be further expressed in terms of the number of e-foldings. It is expected from the present study that this theoretical formulation and the corresponding analysis of cosmological complexity and, its connection with quantum chaos through OTOC could act as a very strong theoretical indicator for future observational probes for studying non-equilibrium physics within the framework of bouncing cosmology.

Organization of the Paper:
• In Sec. 2 a brief review of the concept of circuit complexity has been given and how it can be used to probe new areas of physics in the context of Cosmology.
• Sec. 3 introduces the reader to the framework of Bouncing cosmology and the models that we have considered for the computation of complexity.
• In Sec. 4 a detailed computation of the cosmological scalar perturbations in the bouncing cosmology framework has been provided along with the origin of the squeezed quantum states and its various solutions.
• In Sec. 5 a discussion on the complexity for the squeezed quantum states has been given.
• Finally in Sec. 6 the computational details of the considered models has been provided with all the relevant discussions. We conclude with all our major observations and future prospects in this direction.

Circuit Complexity for dummies
The concept of circuit complexity was primarily used in the field of Computer Science to know the depth of different circuits. It is basically defined as the effort required to carry out a given task or the difficulty in implementing a given task. The task at hand is essentially to prepare a desired quantum field theoretic target state from a reference state. It is generally an optimization technique. Technically, it refers to the minimum number of unitary operations required to implement a given task. The process of carrying out the task involves constructing a unitary transformation that takes a given reference state to the desired final state. The unitary operator being referred to here usually represents the sequences of quantum gates {g i 1 , g i 2 , ....g i n } required to achieve the desired the target state.
Of course, there exist infinite such sequences that produce the desired target state from the given reference state, but the complexity of a quantum circuit provides the sequence which requires the minimum number of gates to do so. This optimal number will depend on the choice of the reference state, |ψ R 〉 and the gate set {g 1 , g 2 , ....g n }. The construction of the unitary operator involves finding a time-dependent Hamiltonian that produces the desired U.
The unitary operator is then constructed from a continuous sequence of parametrized pathordered exponential of the chosen Hamiltonian, The variable s parametrizes a path in the space of unitaries. The Hamiltonian H(s) can be expanded in terms of generalized Pauli matrices i.e.
where M I are the basis in which the Hamiltonian is expanded and the coefficients Y I (s) are the control functions that decide the gate acting at certain values of the parameters. The control function basically represents a tangent in the space of unitaries and acts as the Hamiltonian in the Schrödinger equation satisfied by the unitarity operator U, The idea then is to define a cost for the various possible paths, minimizing which leads to the identification of the optimal circuit. The cost functional is defined as follows: where F is a local cost function depending on the position U(s) and the tangent vector Y I (s).
Once the concept of cost function is introduced, the problem is identical to finding the trajectory of a particle by minimizing the action from the Lagrangian F (U(s), Y I (s)). There are certain desirable features for F to be a cost functional [2] viz. smoothness, positivity, triangle inequality and positive homogenity. Some of the simplest cost functionals which satisfy the above properties and the ones which we have considered in this paper are the linear and the quadratic cost functionals defined as [2]: Linear cost functional : Quadratic cost functional : where the degree of homogeneity is 1 for both of them.
To be precise the cost function F 1 comes closest to counting the number of gates required to make the optimal circuit. The measure F 2 however brings in a notion of proper distance in Riemannian geometry and converts the problem of constructing the optimal circuit to finding the shortest curve connecting the initial and the final states in that geometry. Some other types of cost functionals are also discussed in [2,9].
On the other hand, a general class of inhomogeneous and homogeneous family of functionals are represented by the following expression [9]: F κ family cost functional : [F κ ] 1 κ family cost functional : where for all family members, the degree of homogeneity is represented by the superscript, κ > 1. Here the inhomogeneous family of functionals, F κ was introduced to match the results obtained from both the leading order UV divergences appearing from the well known, "com-plexity= action" [34,35] and "complexity=volume" [34,35] conjectures proposed within the framework of holography. Though the results agreed with with the holographic complexity results, these cost functions do not satisfy the homogeneity property i.e Eq. (6) is not invariant under the reparametrizations of s for these κ family cost functions. Apart from these previously mentioned measures, one can further introduce the following sets of basis independent and state independent cost functionals, which are given by [9]: Trace norm cost functional : Schatten norm cost functional : The Schatten norm cost functional helps to express the circuit complexity in a basis independent way, a problem which occurs with the general κ family of cost functions including the linear and the quadratic ones. Further, one can construct few more state dependent cost functionals which are given by the following expressions [9]: In the context of cosmology, using the quantum squeezed state formalism in the perturbation picture enables one to compute the expression for the cosmological complexity. Scalar perturbations on an expanding background can naturally be described with the formalism of squeezed quantum states. The ground state is chosen as the reference state while the mode is inside the horizon, and a target state consisting of the time-evolved cosmological perturbation on the expanding background. Thus squeezed state formalism gives an elegant way of defining the reference and the target state between which the circuit complexity can be computed. The squeezed state formalism also enables to translate the entire problem in terms of just two quantities known as the squeezed state parameter and the squeezed angle. The whole idea of squeezed state formalism can be easily understood using the well known model of inverted harmonic oscillator. In this formalism, the wave function is squeezed with a large uncertainty in one direction and with a small uncertainty in another direction. Similar observations can be found if one looks into the phase space trajectories of an inverted harmonic oscillator. The presence of one growing and one decaying solution produces a squeezing effect even in the classical level. The main idea behind the squeezed states is to re-parametrize the unitary operator as the product of a squeezed and a rotation operator. The squeezed and the rotation operator can further be expressed entirely in terms of the creation and annihilation operators.
The significance of the rotation operator is less as it mainly produces a phase factor. However, the squeezing operator is of prime significance as the entire problem and all the important observables can eventually be expressed in terms of two quantities, the squeezing parameter and the squeezing angle. Thus, the squeezed state formalism not only gives an elegant way of finding the target and the reference state but also helps to express all the important observables in terms of only two quantities as will be seen in the upcoming sections.
Here, the complexity can be defined in terms of all the previously mentioned different types of cost functionals and one can test as to which ones do give out the best features in terms of the study of quantum chaos for a given cosmological model of our universe. However, in this paper we have restricted our computation by considering only the cost functionals, F 1 and F 2 from which we compute the expression for cosmological complexity. Using the universality relation, we have further computed the expression for OTOC, Quantum Lyapunov exponent and the lower bound on the equilibrium temperature of a system within the framework of bouncing cosmological paradigm.

A simple framework for Bouncing Cosmology
In this section, our prime objective is to construct a bouncing cosmological framework that can further participate in the computation of cosmological complexity. In the present context, we start with the following representative action, given by: where we have fixed the reduced Planck mass M p = 1 for the simplification of the computation. We have introduced a single scalar field with a kinetic term that is minimally coupled with the classical gravitational background. Here V (φ) is the effective potential for the scalar field φ in 3 + 1 dimensions from which we will describe pre-bounce, bounce, and post bounce scenario. We consider here two models which can serve our purpose: Model I : Cosine Hyperbolic model Model II : Exponential model where r 1 is the dimensionless parameter in the Planckian units for both of the bouncing models and V 0 represents the overall energy scale of the potential which mimics the role of Cosmological Constant at a very high energy scale. In Fig. 2, the potentials of the two bouncing cosmology models considered in this paper have been studied with respect to the field variable φ. For both the models, the potential for the pre-bounce region decreases exponentially to negligible values as the value of the field variable increases. An exponential increase in the potential is seen in both the cases as φ is increased. However, for the bouncing region, the behaviour of the potentials is widely different. For the Cosine hyperbolic model, the potential of the bouncing region is negative and goes to large negative values for a slight change in the field variable φ. For the Exponential model, the potential of the bouncing region does also vary even though it is not apparent in Fig. 2.
The aforementioned potentials used to describe the pre-bounce region, bouncing region and post-bounce region can be derived from String Theory descriptions at a very high scale. On the other hand, one can think of another equivalent situation where without introducing a scalar field in the classical gravitational background, one can also study the cosmological bouncing framework. Originally, the concept of cosmological bounce was proposed to resolve the coordinate intrinsic singularity of space-time at the time scale of the Big Bang, which is t = 0. This is because the inflationary paradigm cannot resolve this issue. Not only that, the well known Swampland Criteria and Trans-Planckian Censorship Criteria  which are very useful to construct a physically consistent Effective Field Theory framework at a relatively lower scale than the very high UV cut-off scale of quantum gravity, commonly fixed at the Planck scale, can be described by bouncing paradigm more consistently than the inflation. Additionally, the bouncing cosmological paradigm can be done in presence of higher derivative quantum gravity corrections to the Einstein-Hilbert action. If such corrections are only a function of Ricci scalar then it is known as, f (R) gravity, and within this class R + αR 2 , which is known as the Starobinsky model is the most famous one. 3 One can show that using 3 In the Jordan frame one can actually compute the corresponding mathematical form of the f (R) gravity by making use of the following equations in M p = 1 unit: this model, along with infinite derivative non-local correction to the gravity sector of the form, R + RF (2)R , [116,[214][215][216] and a Cosmological Constant term Λ, can produce the same type of bouncing solution in the spatially flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric in 3 + 1 dimensional space-time, which is described by the following line element: where τ is the conformal time coordinate which is related to the physical time coordinate t through the following replacement relation in the line element: The prime objective to include such non-local correction was to produce a ghost-free renormalizable theory of gravity whose classical limit will be consistent with the local Einstein-Hilbert gravity contribution. Apart from this, the bouncing framework is very important in the context of primordial cosmology because the Big Bang singularity can be removed from the theory by imposing the bouncing condition on the related scale factors in the spatially flat FLRW background, which can be explicitly computed by making use of the Friedmann equation and the Klein-Gordon equation for the scalar field φ. At the cosmological bounce scale t = t B one has to satisfy the following constraint conditions to find out the appropriate dynamical solutions of the field equations: Bouncing condition I : Bouncing condition II : This same condition for the bounce at the conformal time scale τ = τ B can be further translated in the following simplified form: Bouncing condition I : Bouncing condition II : This implies that the mathematical structure of the bouncing conditions remains the same in physical time and the conformal time coordinates, though they are not exactly the same as we have pointed earlier. One can also write constraint conditions on the potential function at the point of bounce, which is given by the following expressions: For an example, for the potential V (φ) = 1 2 m 2 φ φ 2 , with φ 1 we get, f (R) = R 2 and with φ 1 we get, f (R) = R. So by considering both the limiting contribution one can construct a f (R) function which is basically made up of both R and R 2 contributions and they are appearing with appropriate coefficients i.e., f (R) = αR+βR 2 . For φ 1, we have α β and for φ 1 we have α β.
Consequently, around the point of bounce if we expand the potential function in Taylor series in the field space, we get: where the first three terms are the renormalizable contributions and other · · · represent nonrenormalizable terms. From the previously mentioned models the scale factors can be computed in terms of the physical time coordinate as: Model I : Model II : a Post exp 9 2   Model I :

Pre-Bounce
Model II : The scale factors have been plotted in the Logarithmic scale to show the rising values near the boundary of the bouncing region as can be seen in Fig. 3. We also expect the models to satisfy the bouncing conditions given before, in Fig. 5 and Fig. 6, and for both the models H B = 0 and H B > 0 can be verified. As we will see later, the behaviour of H and the difference in signs on either side of the point of bounce will require two different squeezing parameters that describe the bouncing region, one before the point of bounce and one after it, which will result in differing behaviour of the complexity.
In the next section using these solutions our prime objective is to perform the cosmological perturbation and find the explicit role of these class of solutions to construct the squeezed vacuum sates.

Perturbation with squeezed quantum states in Bouncing Cosmology 4.1 Scalar perturbation in Bouncing Cosmology
In this section we will study squeezed state formalism within the framework of cosmological perturbation theory [217][218][219][220][221][222] for FLRW spatially flat background specifically for post-bounce, bounce and pre bounce region. In this context one needs to consider the following perturbation in the scalar field in the De Sitter background: and to express the whole dynamics in terms of a gauge invariant description through a variable: At the level of first order perturbation theory in a spatially flat FLRW background metric, we fix the following gauge constraints: which fix the space-time re-parametrization. In this gauge, the spatial curvature of constant hyper-surface vanishes, which implies curvature perturbation variable is conserved outside the horizon.
Applying the ADM formalism one can further compute the second-order perturbed action for scalar modes. The action, after gauge fixing, can then be expressed by the following: Now, to re-parametrize the above mentioned second-order perturbed action expressed for primordial scalar perturbation, we introduce the following space-time dependent variable: which helps transform the perturbed action to that of the familiar mathematical form of canonical scalar field. In the cosmology literature, this is known as the Mukhanov variable, in terms of which we will perform the rest of the computation. Additionally, it is important to note that the newly defined quantity, ε(τ) is the conformal time dependent slowly varying parameter, which is defined as: Consequently, the new version of the second order perturbed action for the scalar perturbation after re-parametrization in terms of the Mukhanov variable can be written as: (39) Now, we explicitly compute the following crucial conformal time dependent contribution, which plays a significant role to explore various unknown physical facts of the primordial universe: Our job is now to further convert the second-order perturbed action for the scalar degrees of freedom in terms of the Fourier modes, by implementing the following ansatz for the Fourier transformation: using which one can compute the following contributions from the time and space derivative of the perturbed field variable appearing in the second-order action : After the substitution of all the aforementioned expressions, the simplified version of the second-order perturbation for the scalar modes in Fourier space can be further recast as: where it is important to note that: Now after varying the second-order perturbed action with respect to the perturbed field variable expressed in the Fourier space, we get the following equation of motion: This is commonly known as the Mukhanov-Sasaki equation and actually represents the classical equation of motion of a parametric oscillator where the frequency of the oscillator is conformal time dependent and in the present context of discussion, be explicitly given by : where we have introduced a conformal time dependent effective mass in the present computation, which is quantified by the following expression: where for the purpose of simplification of computation we have introduced a conformal time dependent mass parameter, ν B (τ), which is defined as: where · · · is the contribution which is varying very slowly in the context of our present discussion.

Scalar mode function
As a result, the Mukhanov-Sasaki equation can be translated into the following simplified form: The most general analytical solution of the above equation can be expressed as: where are Hankel functions of the first and second kind,respectively, with argument −kτ and order ν B . During this computation, we have also used the fact that the conformal time-dependent quantity ν B is varying very slowly with respect to the evolutionary time scale of our universe. Additionally it is important to note that, the two integration constants, C 1 and C 2 can be fixed by the choice of the initial quantum vacuum state in the present context. In this work, we choose the most popular and the simplest initial quantum vacuum state, which is known as Bunch Davies vacuum or Hartle Hawking vacuum or Chernkov vacuum, and can be fixed by choosing C 1 = 1 and C 2 = 0.
Consequently, we get the following solution: Upon further considering −kτ → 0 and −kτ → ∞ asymptotic limits, one can write the following simplified form of the Hankel functions of the first kind: Using these asymptotic results of the Hankel functions of the first kind the most general solution for the perturbed field can be expressed as: In the present solution, the slowly varying time-dependent mass parameter ν B (τ) is a completely model-dependent one. For this reason, to fix the value and the behaviour of the slowlyvarying function with respect to the underlying conformal time scale we need to explicitly compute this expression for different models which are describing the pre-bounce, bounce, post-bounce, and the away from the bounce region. 4 One can further consider two asymptotic cases, super-Hubble and the sub-Hubble which might be extremely useful to study the physical impact of the mode function obtained for the scalar fluctuations in the two different physical regions as mentioned before. In terms of the representative dynamical scale, the super-Hubble and the sub-Hubble limit is described by −kτ 1 and −kτ 1, respectively. Additionally, it is important to note that in this context of the discussion, the cosmological horizon crossing is described by −kτ = 1. Now we shall implement all the discussed limits to get simplified results from the scalar mode function obtained previously within the framework of bouncing cosmological paradigm. These limiting results are appended below: Horizon crossing solution : 4 Note: Here it is important to note that during inflation the mass parameter ν B = 3 2 , if we exactly follow the De Sitter expansion in the spatially flat FLRW background. But in order to stop inflation, one needs to consider a slight deviation from exact De Sitter expansion during inflation, and technically this slight amount of deviation has been taken by considering the slowly varying time dependent slow-roll parameters. So it is expected that for exact De Sitter expansion, the factor ν B − 3 2 will exactly vanish, and for the quasi-De Sitter expansion, this difference will be proportional to the amount of deviation from the exact De Sitter expansion. But in the present context we are interested in the pre-bounce, bounce, post-bounce, and away from bounce, where it appears to us that the analytical solution of the scalar mode function appearing from the cosmological perturbation in the spatially flat FLRW background is identical to the structure that one may compute by solving the equation of motion of the scalar mode fluctuation, which is the Mukhanov-Sasaki equation in the context of inflation. The significant difference can be observed clearly if we look into the mathematical structure and the leading , sub-leading order contribution appearing in the expression of the mass parameter in both of the cases separately. For inflation, this value is slightly larger than 3 2 , which as we told can demonstrate the quasi-De Sitter expansion. On the other hand, for the alternative to the inflationary paradigm -which is described by pre-bounce, bounce, post-bounce, etc., it is expected that the value of the mass parameter will be completely different from 3 2 and the amount of deviation from the exact De Sitter is very large. This is because the slowly varying parameter ε and its derivatives are significantly large compared to the value obtained for this parameter, which is smaller than unity during inflation and approximately unity at the end of inflation. Apart from this underlying significant difference, for the sake of consistency with the previous works and their findings, we have expressed the solution of the scalar mode function for the pre-bounce, bounce, post-bounce, and away from the bounce phases like the result obtained from inflation.

Quantization of Hamiltonian for scalar modes
Using these solutions, one can further compute the expression for the derivatives of these field variables with respect to the conformal time scale, which will be helpful for the further computation in the present context: As mentioned in the previous subsection, one needs to further consider two asymptotic cases, the super-Hubble and the sub-Hubble limiting situation which might be extremely useful to study the physical impact of the obtained mode function for the scalar fluctuations in the present context. In terms of the representative dynamical scale, the super-Hubble and the sub-Hubble limit is described by −kτ 1 and −kτ 1, respectively. Additionally, it is important to note that in this context of the discussion, the cosmological horizon crossing is described by −kτ = 1. By following the same logical reasoning one can write down the following expressions for the conformal time derivative of the mode functions from scalar fluctuations which will explicitly contribute further in the expression for the canonically conjugate momenta associated with these scalar modes: Horizon crossing solution : Super − Hubble limiting solution : Now, our next objective is to construct the classical Hamiltonian function studied for the present parametric oscillator problem. For this purpose, we need to find out the expression for the canonically conjugate momentum for the classical cosmologically perturbed scalar field variable appearing previously in the second order action perturbed action of the system that  Figure 7: The real part of v k for the Cosine hyperbolic and the exponential bounce has been plotted. We can see the behaviour of v k in the sub-horizon region |kτ| 1. Around the Horizon crossing, the Re(v k ) is slightly negative for the Cosine hyperbolic bounce whereas it is almost zero for the exponential bounce. Far from the horizon crossing, the behaviour of both the models is almost identical. It slowly starts rising before becoming highly oscillatory with increasing amplitude and frequency as they approach the point of bounce kτ b ≈ 30. The field variables in all regions behave almost similarly. However, a noticeable difference between the two models lies in the fact that the amplitude of the field variable for a particular region is in contrast for both the models. If the amplitude is maximum for a particular region for model-I then it is minimum for model-II.
we have mentioned earlier in this section, and it is given by the following expression: Further, using the above mentioned results one can construct the expression for the classical Hamiltonian function from the present problem set up, which is given by: where the time dependent mass µ 2 (k, τ) of the parametric oscillator is given by the following expression: Next, using the previously mentioned solution of classical mode function we can further construct the quantum mechanical operators in the Heisenberg picture:  Figure 8: The imaginary part of v k for the Cosine hyperbolic and the exponential bounce has been plotted. We can see the behaviour of v k in the sub-horizon region |kτ| 1. Around the Horizon crossing, the I m(v k ) starts off positive for the Cosine hyperbolic bounce, whereas it starts off negative for the exponential bounce. Far from the horizon crossing, the behaviour of both the models are almost identical. It slowly starts rising before becoming highly oscillatory with increasing amplitude and frequency as they approach the point of bounce kτ b ≈ 30. The field variables in all regions behave almost similarly. However a noticeable difference between the two models lies in the fact that the amplitude of the field variable for a particular region is in contrast for both the models. If the amplitude is maximum for a particular region for model-I, then it is minimum for model-II.
Now using the above mentioned quantum operator one can finally express the canonical Hamiltonian for the parametric oscillator in the following quantized form (See Appendix A): where we define Ω k (τ) and λ k (τ) by the following expressions: Here Ω k (τ) represents the conformal time dependent dispersion relation in the present bouncing cosmological set-up, and λ k (τ) basically captures the slowly conformal time varying function ln z(τ), where z(τ) = a 2ε, is the Mukhanov variable, which appears during the computation of cosmological perturbation for scalar modes in the bouncing set-up. For the details of the computation, please refer to Appendix B.

Fixing the initial condition at horizon crossing
Here it is important to note that one can fix the initial condition in such a way that, at the time scale τ = τ 0 , we get the following normalization: provided we have imposed a constraint that, kτ 0 = −1, which basically represents the horizon crossing scale. Following this fact it is further expected that at any arbitrary later time scale τ in the Heisenberg picture one can write the associated quantum operators for the present problem as:v where both the creation and the annihilation operators at time τ can be expressed in terms of the results obtained from the initial time scale τ = τ 0 using the following unitary similarity transformation in the Heisenberg picture: Our next job is to determine the expression for the above mentioned unitary operator in the context of cosmological primordial perturbations of the scalar modes and to determine this expression, the well known squeezed state formalism used in the context of quantum mechanics will play a significant role.

Squeezed state formalism in Cosmology
The unitary evolution operator U, produced by the previously mentioned full quadratic quantized Hamiltonian function, can be factorized by following the proposal given in refs. [118,119] and can be written as: where R is the two mode rotation operator, which is defined as: andŜ is the two-mode squeezing operator, defined as: Here the squeezing amplitude is represented by the time-dependent parameter, r k (τ) ,and the squeezing angle or the phase is represented by the time-dependent parameter φ k (τ). Additionally, it is important to note that, the two-mode rotation operator,R produces an irrelevant phase contribution exp(iθ k (τ)) while acted upon the initial quantum vacuum state and can be ignored from our current analysis to avoid the appearance of unnecessary junks. By recognizing that the interaction of the cosmological perturbation with the conformal time-dependent scale factor in the spatially flat FLRW background leads to a conformal time-dependent frequency for the canonically normalized parametric oscillator, the appearance of a squeezed quantum mechanical state for cosmological primordial perturbations is quite natural. The quantization of the conformal time dependent parametric oscillator is then described in terms of two-mode squeezed state formalism as introduced in ref. [118].
For our further computation we choose the ground state of the free Hamiltonian as the initial quantum mechanical state:â which is basically a Poincare invariant vacuum state in the present context of discussion. Now we are going to use the squeezed quantum operatorŜ which acts on the above mentioned initial vacuum state and produce a two-mode squeezed quantum vacuum state, as: with the following two-mode excited or usually known as the occupation number state given by the following expression: Consequently, in the present context of discussion the full quantum wave function can be expressed in terms of the product of the wave function for each two-mode pair as k, −k given by the following expression:

Time evolution in squeezed state formalism
Now we go back to the previous discussion where we have written the creation and the annihilation operators of the conformal time dependent parametric oscillator in the cosmological perturbation theory at any arbitrary time using the Heisenberg picture. This will help us to explicitly identify the time evolution of the perturbation field variable operator corresponding to the scalar modes and its associated canonically conjugate momentum operator. In terms of the above mentioned squeezed quantum state description one can further express the creation and annihilation operators in the present context as the unitary operator for the time evolution in the Heisenberg picture. The unitary operator can in turn be factorized in terms of the two-mode rotation operator and two-mode squeezed quantum state operator as we have discussed earlier. After performing the unitary similarity transformation in terms of the twomode rotation and squeezed operator, one can write down the following expressions for the creation and the annihilation quantum operators at any arbitrary time scale τ as: (83) Consequently, the quantum operator associated with the cosmological perturbation field variable for the scalar fluctuation and the its canonically conjugate momenta can be expressed as: Here we identify the classical mode function and the associated canonically conjugate momentum in terms of the squeezed parameters as: π k (τ) = π k (τ 0 ) cosh r k (τ) exp(iθ k (τ)) + sinh r k (τ) exp(i(θ k (τ) + 2φ k (τ))) . (87) Further, the time evolution of the conformal time dependent quantum operatorsR andŜ are described by the Schrödinger equation, which gives the following set of differential equations for the squeezing parameters in the present context: where the time dependent factors, λ k (τ) and Ω k (τ) in the squeezed state picture in the sub-Hubble region (−kτ 1) can be recast as: Here it is important to note that in the sub-Hubble region the factor Ω k (τ) is mainly controlled by the momentum scale of the scalar mode of the perturbation, k, and the slowly varying time dependence is taken care of by the conformal time dependent mass parameter ν B , which can be approximately written by considering the contribution upto the next-to-leading order as: where we have neglected the contributions of all higher order small correction terms for the computational simplicity. Now after substituting the above mentioned expression for the mass parameter ν B one can further write the following simplified form of the factor, Ω k (τ) in the sub-Hubble region, as: where γ E is the Euler-Mascheroni constant, which is γ E = 0.577. For a more detailed discussion on dispersion relation please refer to Appendix C.

Quantum complexity from squeezed quantum states in Bouncing cosmology
In this section, we compute the complexity from the squeezed cosmological perturbations studied in the previous section for the bouncing framework. We use the wave function formalism of computing circuit complexity developed by [2,3] and used extensively in [10][11][12]. Computing the circuit complexity involves choosing a certain reference state and a target state. In the case of cosmological perturbations, a commonly chosen reference state is the two-mode quantum initial vacuum state |0〉 k,−k , as mentioned in the previous section. The target quantum state is the squeezed two-mode vacuum state |Ψ sq 〉 k,−k . In ref. [2,3] the authors expressed the reference and the target states as Gaussian wave-functions. We follow an identical approach in this paper for further computation. We use the following field operator and its associated canonically conjugate momentum operator as: where v k (τ 0 ) and π k (τ 0 ) fix the initial condition on the classical scalar mode and its associated canonically conjugate momentum at the horizon crossing scale, −kτ 0 = 1. We have computed their explicit expressions in the previous section. Additionally, we have also computed the expressions for the associated quantum operators at any arbitrary time scale τ in terms of the squeezed conformal time dependent parameters r k (τ) and θ k (τ) in the Heisenberg picture of quantum mechanics. At any arbitrary time scale τ, these cosmological quantum operators satisfy the following well known equal time commutation relation (ETCR), given by: The two-mode vacuum state wave function, which we choose as our reference state is defined as:â which has the following usual Gaussian structure: where we have used the expression for Ω k in the sub-Hubble region, the approximated analytical expression of which we have already derived explicitly in the previous section. The wave function of the target or the squeezed quantum state for the cosmological perturbation can be calculated by noting that a particular combination of the squeezing parameters along with the creation and annihilation operator annihilates the two mode squeezed vacuum state, constructed in the previous section. That particular combination is written as: The cosmological perturbed field space representation of the wave function is given by the following expression: where the coefficients A(τ) and B(τ) are the functions of the squeezing parameter r k (τ) and the squeezing angle φ k (τ), and are explicitly given by the following expression: The vacuum reference and the target squeezed state written in 98 and 100 is eventually used to calculate the complexity from two types of cost functions namely the "linear weighting" (C 1 ) and the "geodesic weighting" (C 2 ) respectively within the framework of Cosmology and represented by the following expressions: A trivial generalisation of the complexity measure of the homogeneous and inhomogeneous family of cost functionals can be also be done in the present context. The expression of the complexity for the homogeneous family is given by: Similarly the expression of the complexity for the inhomogeneous family can be written as: where we define the following functions: It might happen that in some particular context, the measures C 1 and C 2 are not good enough to probe the underlying chaos and randomness of the system. Complexity calculated from the homogeneous and the non homogeneous family might come in handy in that scenario and may bring out some essential features which remains unidentified by C 1 and C 2 . In this paper, we have mainly focused on the complexity measure calculated from the linear and the geodesic weighted cost functionals to comment on the chaoticity of the universe from the bouncing cosmological framework. We can express the complexity measures in terms of the squeezed state parameters. Substituting 101 in 103, 104, 105 and 106 the complexity measures for the bouncing set up for two mode squeezed vacuum state can be written as: One can further derive some approximate analytical expressions for the cosmological complexity in different limiting situations, which are discussed below: For small r k (τ) and φ k (τ) one can write: In this limit, we have the following simplified formulae of cosmological complexity for the bouncing set up for two mode squeezed vacuum state: 2. Large r k (τ) & Large φ k (τ): For large r k (τ) and φ k (τ) one can write: Consequently, the cosmological complexity for the bouncing set up for two mode squeezed vacuum state reduces to the following simplified expressions: which will finally lead to the following approximated connecting relationship between the two cosmological complexities computed from different cost functions: 3. Small r k (τ) & Large φ k (τ): For small r k (τ) and large φ k (τ) one can write: Consequently, we have the following simplified formulae of cosmological complexity for the bouncing set up for two mode squeezed vacuum state: which will finally lead to the following approximated relationship between the two cosmological complexities computed from different cost functions: For large r k (τ) and small φ k (τ) one can write: Consequently, we have the following simplified formulae of cosmological complexity for the bouncing set up for two mode squeezed vacuum state: In the next section, we have done a detailed numerical analysis with the already introduced models of bounce to study their physical impacts on cosmological complexity from two types of cost functions and interpret the physical outcomes from those models.

Numerical results and interpretation: Connection with quantum chaos
In this section our prime objective is to numerically solve the time evolution equations of the conformal time dependent squeezed state parameter r k (τ) and squeezed angle θ k (τ), given in Eq. (88) and Eq. (89). However, instead of using the conformal time τ as the dynamical variable, we have chosen the scale factor a(τ) to make the computation simpler and physically justifiable. To perform the change in variable from τ to a(τ) we have to replace the following differential operator in the above mentioned evolution equations using the chain rule, as: In general quantum field theory literature we usually identify such type of variable transformation as field redefinition. One can treat the scale factor a(τ) as a classical field and the same interpretation is valid in this context. Consequently, the evolution of the squeezed state parameter r k (a) and squeezed angle θ k (a), can be recast in terms of the newly defined dynamical variable a(τ) as: In the above set of evolution equations, since we do not need to care about the explicit conformal time dependence we have written the scale factor a(τ) as a, where a itself is treated as a new dynamical variable. Once we numerically solve the evolution of the squeezed state parameter r k (a) and squeezed angle θ k (a) in terms of the scale factor a, we can construct the target quantum state out of a Gaussian initial state. This will further help us to numerically compute and understand the quantum complexities in Eq (103) and Eq (104) within the framework of primordial cosmological perturbation theory, where the effects of the quantum fluctuations is treated in terms of the squeezed state parameter r k (a) and squeezed angle θ k (a) in the squeezed state formalism. For the explicit computational details, we suggest the readers to look into the previous two sections very carefully where we have explicitly shown why and how these interesting connections can be established. Now since we have a good understanding of both the complexities, C 1 (a) and C 2 (a), we will compute them from the two previously mentioned cost functions and analyze the behaviour, from C 1 (a) and C 2 (a), vs scale factor a plots, specifically in the exponentially rising region. Now, by studying the exponential rise in the complexities, C 1 (a) and C 2 (a), one can write the following approximated expression for the complexities: which are valid only in the domain of exponential rising with respect to the scale factor a. Additionally, it is important to note that, though the exponential growth feature is common in both the complexities, we have written the expressions for the two complexities separately because the overall amplitudes, which are represented by c 1 and c 2 , and the slope of the previously mentioned plots, quantified by two factors, λ 1 and λ 2 , are different which can be confirmed by comparing the features of both the plots. This can be demonstrated as: where a grow is the specified value of the scale factor from the region where exponential growth feature can be explicitly visible from the complexities vs scale factor plots. Most importantly, Eq (139) is a conjectured relationship which we have proposed by seeing and comparing the numerical behaviour of the obtained plots from this analysis. For this reason, we have written ≈ symbol instead of using =. To know the complete evolution one needs to solve the system numerically which will give us an exact result, valid in all evolutionary regions of the scale factor a, and not only in the exponentially rising region. On the other hand, by doing the explicit computation of out-of time-ordered correlation (OTOC) functions obtained from the classical field a and its canonically conjugate momenta π a one can find the following relationship: which is again valid in the specific region of interest. Here λ is identified to be the Quantum Lyapunov Exponent which captures the effect of chaos in the quantum regime, and in ref. [19], the authors, Juan Maldacena (M), Stephen Shenker (S) and Douglus Stanford (S) have found that for a generic quantum chaotic system Quantum Lyapunov Exponent has to be bounded by the following saturation upper bound, as given by: where β is the inverse equilibrium temperature of the chaotic system during saturation of the OTOC at large evolutionary scale. The equality symbol in the MSS bound represents the maximal saturation of chaotic OTOC. Now further using Eq (139) and Eq (141) together we get the following simplified results: Now after studying the above mentioned equations we can arrive at the following conclusion: which implies that the connection between OTOC and the two different measure of complexities are not strictly same. Additionally, in the present context we have the following restriction: However, to have an universal feature it is expected that the following fact is also true in the present context: which is only true in the limit, λ i → λ, c i → c ∀ i = 1, 2. In this limit, precisely we have: Negligibly small sub−leading effects (148) This further implies that if we neglect all the extremely small sub-leading contributions, and restrict our attention to only the leading order term then it is possible to write down the following universal highlighting relationship between all possible measures of complexities and the OTOC, as: Here it is important to note that, the above mentioned universal relation is perfectly consistent with the ref. [11]. The only difference is that, here we have achieved the universality using the dynamical variable, scale factor a and in ref. [11], the authors have pointed such universality using the physical time variable t. Though, both the discussions hold good in their preferred choice of dynamical variables, ultimately both of them support the same chaotic behaviour during the exponential rise. Also it is observed that, when the universality is achieved we expect to get a saturation in the behaviour of complexities as well as in the OTOC with respect to the dynamical scale a. Now to have a precise agreement with consistency condition, which is described by the well known MSS bound, one needs to satisfy the following constraint, which will provide a cost function model dependent lower bound on the Lyapunov exponent appearing from the definition of the complexities: If the maximal saturation is achieved, then from this relation one can further get a lower bound on the equilibrium temperature of the quantum system of our universe under study during bouncing scenario, and this is given by: Finally, when the universality as well the maximal saturation both have been achieved simultaneously in the above mentioned expression, the equality gives the exact estimation of the equilibrium temperature of the quantum system of the universe studied during bounce, which is valid at very large values of the evolutionary scale represented by a. Now from the above bound since λ i ∼ λ and that λ i = λ, it is also expected that the lower bound on the equilibrium temperature can have two predictions in terms of the two possibilities of the complexities originated from two possible cost functions in the present context. However, the numerical order of both of the predictions computed from the plots will be same and somewhat in a broader sense support the universality criteria, which tells us that both of the predicted temperature will not be much different. From the above obtained lower bound on the equilibrium temperature one more important aspect we want to point here is that, this result does not depend on any particular particle content or a specific model available during bounce and gives us a generic estimation of the equilibrium temperature. Now if we are thinking about the more realistic cosmological observation then it is not very good to study the evolution with respect to the scale factor, because in the context of realistic cosmology the scale factor is not the direct physical observable which one can probe in the observation for various cosmological missions running (or supposed to run in the near future) to test the signatures of the primordial cosmological paradigm. In that case instead of using the scale factor a one can consider a more physically realistic variable, which is the rescaled number of e-foldings, N , which one can use as a direct probe in various cosmological observations. In this specific situation one needs to use the following transformation for which the linear differential operator appearing in the evolutionary equations of the squeezed parameter and the squeezed angle will be modified as: where we have used the following couple of facts for the above mentioned transformation: Here, N is the actual number of e-foldings, N is the number of e-foldings in terms of the re-defined variables, and ε(τ) is the slowly varying conformal time dependent parameter. Consequently, the evolution of the squeezed state parameter r k (N ) and squeezed angle θ k (N ), can be recast in terms of the newly defined dynamical preferred choice of suitable variable N as: In this context, r co = (aH) −1 or r co = H −1 represents the co-moving Hubble radius, which is extremely important quantity in terms of which the newly re-defined number of e-foldings have been expressed in terms of the good old definition of the number of e-foldings. So instead of solving these sets of first order coupled differential equations in terms of the dynamical variable a here our further objective is to study the evolution numerically with respect to the re-defined dynamical variable, N . Here, we additionally want to point out that by replacing the dynamical variable a in terms of the re-defined expression for the number of e-foldings N , we can write down similar type of conclusion which we have written earlier to interpret the exponential growth and then the saturation in the large scale. Here one can write: Similarly one can derive the universality relation which will be same as the previous one. In the next two subsections, we will explicitly numerically solve the previously mentioned dynamical equations of the squeezed parameter and squeezing angle with respect both the dynamical variables, scale factor a and the re-defined number of e-foldings N for cosine hyperbolic and exponential bouncing models that we have introduced in the first section of the paper. The explicit details of the analysis and the corresponding physical interpretation of the numerical results and the plots are discussed in the following two subsections. Discussion of the differential equations with respect to different dynamical variables is given in Appendix D.
Another important aspect that one can estimate numerically from our present set up, is the well known scrambling time scale. Within the framework of quantum chaos this time scale plays very significant role to understand the underlying behaviour of the physical systems. There are several definitions associated with this quantity, in the theoretical physics community in different contexts for physical interpretation of various unknown phenomena. We will now quote the most frequently used definitions, and follow one of them in the present context to numerically estimate the order of scrambling time scale from the bouncing cosmological scenario:

Definition I:
According to this definition this is the time which the OTOC takes to equilibriate. This is a very modern definition and directly associated with the phenomena of quantum mechanical chaos. 5

Definition II:
According to this definition this is the time which a system takes starting in an arbitrary tensor product state to become nearly maximally entangled.
Now according to Leonard Susskind [126] and later pointed in many other refs. [223] for the first scrambler the scrambling time scale can be computed as: where β is the inverse of the equilibrium temperature of the physical system which corresponds to the saturation of quantum chaos and N represents the very large number of configurations. Now making use of the MSS bound one can further simplify the above mentioned expression and obtain a lower bound on the scrambling time scale in terms of the quantum Lyapunov exponent: Here the equality holds good for the maximal saturation of chaos. Now, within the present framework we have used the conformal time dependent scale factor a and/or the number of e-foldings N as dynamical variable using which we have studied all the evolution of cosmological complexity and the OTOC in this paper (for the details see the next two subsections.). One can then ask a very justifiable question in this case that how we then define the scrambling time scale within the framework of cosmology? Following the previous logical discussions and interpretations of the universality relation between the cosmological complexity and cosmological OTOC by replacing the time with the scale factor one can define the scale factor at scrambling time scale or scarmbling scale factor, which is given by: Here the index i = 1, 2 is used to differentiate between the value of the scale factors obtained from the two definitions of complexities used in this paper. To hold the universality between the cosmological complexities and the OTOC we have previously shown the deviation from the results obtained from both of the definition has to lie within a very small numerical error range. It is expected that the same argument also holds here perfectly and in the next two subsections we are going to investigate this very carefully from the numerical plots to ensure the justifiability of this statement. Now, we have already computed the expression for the scale factor in terms of the conformal time for both the models and also most of the quantum chaotic predictions are appearing ( for the details see the next two subsections.) from the bouncing solutions. For this reason using those definitions one can extract the information of the associated scrambling time scale in the conformal coordinates within the framework of bouncing cosmological paradigm. Additionally, since we also know the connecting relationship between the physical time scale and the conformal time scale, then using this it is further possible to determine the scrambling time scale in terms of the physical time coordinate in cosmology. In the next two subsections, for two different bouncing models we are going to estimate this time scale from the numerical plots. Finally, there is a confusion regarding the fact that in the present cosmological set up how can one give a numerical estimation of the factor N which represents number of physical configurations. We are now going to give an estimate of this factor in the present context in terms of the known parameters. To obtain this estimate we start with the following relationship: Using this relation and truncating the expression for OTOC in the second term we get: where in the usual quantum chaos literature one can identify: The one can further write the expression for the scrambling scale factor in terms of the known parameters as: From the numerical plots given in the next two subsections one can estimate both λ i and c i (for i = 1, 2) from both the bouncing models and from this relation it is possible to give a numerical estimation of the scrambling time scale from the models of bouncing cosmology discussed in this paper. Additionally it is important to note that, in this connection the equivalent result can be obtained by considering the number of e-foldings as the dynamical variable instead of the scale factor within the framework of cosmology.

Cosine Hyperbolic bounce
We have numerically plotted the squeezing parameters and the derived complexity measures for cosine hyperbolic in four different regions -pre bounce boundary, pre point of bounce, post point of bounce and post bounce boundary against the scale factor 6 . From Sec. 3 we can see that at present time and at a time much before the boundary (τ → −∞) the value of scale factor a = 1. We have taken the value of pre-boundary and post-boundary parameters r k (a = 1) = 1, φ k (a = 1) = 1 to set our initial conditions, and ensured continuity at a boundar y as initial conditions for the bouncing region parameters for numerically solving differential equations with respect to scale factor (Eqs. D. 15 and D.16). For the analysis of Cosine hyperbolic bounce we have taken −kτ b = 30 and the range of −kτ goes from 0 to 60. The parameter r 1 , appearing in the expression of the scale factor is related to the cosmological constant by the relation r 1 = 2Λ/9. We have fixed the value of Λ to 10 −4 for our numerical analysis.
For the squeezing paramater plotted in Fig. 9 • the pre-boundary and the post-boundary behaviour is oscillatory with decreasing amplitude as it approaches a = 1(very early times in case of pre bounce boundary and present time in case of post bounce boundary), • while inside the bouncing region we see highly oscillatory behaviour near the point of bounce(with very high amplitudes) that saturates into a given value as it nears the boundary. This saturation behaviour of squeezing parameter near the boundary might point to saturation behaviour of the Complexities as we will see in further analysis.
The squeezing angle and the sine of twice its value are also important to understand the Squeezing operator. See Fig. 10 and Fig. 11.
• The has an exponential increase even against a logarithmic scale, with the rate of increase falling down while approaching the bounce boundary from earlier times in the pre-boundary region. The frequency of the sin 2φ k corresponds to the rapid rate at which it increases, initially oscillating really fast to slow spaced oscillations at the boundary. • Upon entering the bouncing region the angle just has a sturdy exponential rise till the point of bounce after which it exponentially increases with a slow rate till the boundary after crossing the point of bounce. The sine again behaves similarly with slowed down oscillation at the boundary, where we see saturated rate of change in the angle.
• Outside the boundary the angle exponentially decreases at a rapid rate and the sine value correspondingly increases in oscillations as we approach present time.   good measures of the circuit complexity, the linearly weighted complexity C 1 shows similarity to the calculations from holographic side. Both the complexity measures have very similar behaviour Fig. 12 and Fig. 13, • The value of complexity outside the bouncing boundary based on the respective squeezing parameters defined there for very early times and nearing present times is oscillatory with smaller frequency at the boundary.
• Inside the bouncing region prior to the bounce both the complexities cross the boundary with a sturdy rise and then go on to become highly oscillatory and spike up. The post point of bounce values are of greater interest as they show a sturdy rise and a nice saturation when extrapolated to present time.
• Even though the post bounce boundary behaviour looks oscillatory it is important to note that the growing behaviour of complexity at post point of bounce. We see a sudden exponential rise near the boundary. The analysis of growing complexity observed by extrapolating the post point of bounce values at late times shows saturation after an initial rise across the boundary. We have written down the values in Table 1.
We can see extremely high complexity values at the point of bounce. This points to highly complex transformations taking place between the reference and target quantum state during the bounce.
The slope of logarithm of complexity at the point of rise directly corresponds to the value of the quantum Lyapunov exponent as mentioned in Eq (140). To predict the slope of the logarithmic value of complexities we consider the change of y-axis value over the range of the x-axis value i.e. between point of rise and point of saturation. Mathematically it is represented by For this we have plotted the logarithm of complexity values in Fig. 14 and Fig. 15. We observe the qualitative features to be same as that of the complexity graphs, showing corresponding oscillatory, rising and saturation at the respective regions. We calculate the Lyapunov exponent from the post point of bounce case as it shows exponential and saturation at late times and this gives an estimation on the lower bound of temperature.  The Lyapunov exponent can be calculated from these values given in Table 2: The estimated lower bound on the temperature from the calculated values of the Lyapunov exponents are Using Eq(166), we have numerically calculated the lower bound of scrambling time scale in terms of scale factor and conformal time. We have considered the region of saturation and taken the values of complexity at what we have perceived as the starting point and the ending point of the region of saturation. We have then calculated ∆a i (τ sc ), which will give us the lower bound of the scrambling interval in terms of the scale factor. We have converted this in terms of conformal time for easy physical interpretation. In our numerical analysis we have extensively used the conformal time and we have normalized all other numerical measures with respect to conformal time in both models whereas the physical time is not normalized with respect to our numerical analysis. Hence calculating scrambling time scale in terms of physical time will not make much sense quantitatively in our case without appropriate normalization (and redoing complete analysis). In our cosine hyperbolic case we have normalized conformal time in such a way bounce occurs at τ b = −3000, and present day time is 0, and hence we can interpret the values given in Table 3, qualitatively in terms of physical time too. We get conformal scrambling time scales around one-tenth of the time period since bounce till present. This roughly points to the time taken for OTOC to attain equilibrium as can be seen from the graph. A quicker scrambling time scale points to smoother saturation of complexity. A sense of time period in terms of physical time can then be qualitatively understood using this argument. The OTOC values plots calculated from the universality relation mentioned in Eq (149). The behaviour observed is very similar to the complexity behaviour at different regions -from being random oscillations outside the bouncing region to settling at the boundary to again oscillating and spiking at the point of bounce. The OTOC values at different points have been   written in Table 4. One noticeable observation is the really small value of the OTOC at the point of bounce from both the complexity measures. We know that more than the scale factor the Number of e-foldings(N = log aH) is a measurable and interesting value. Using the simple relation we have also plotted Complexity and the OTOC against N . One can see from the behaviour of a, H from Sec. 3 and Sec. 3, that the direction of log aH will be different inside the bouncing region and outside the bouncing region. For this reason the plots have been made separately to ensure readability. 7 7 Reading graphs vs N : Inside the bouncing region the value of N at bounce and boundary can be obtained and it is seen that N bounce < N boundar y and extrapolating the same one can get value at present time which is greater For both linearly weighted complexity C 1 and geodesically weighted complexity C 2 : • the outside bouncing region behaviour is oscillatory and random. The oscillations decrease near the boundary. We have extrapolated their graphs inside the boundary to show that their oscillatory behaviour drops to a certain/very low complexity. But this is not what is actually expected from complexity measures previously done inside the bouncing region. See Fig. 19 and Fig. 21.
• Inside the bouncing region we see that the pre point of bounce graph starts at a high value at boundary and oscillates randomly till bounce. From Fig. 18 we observe continuity at point of bounce(at a higher value than whatever the extrapolated outside region lines pointed at) and see the post point of bounce graph follow random oscillations till nearing the boundary from which it starts rising and slowly goes on to saturates at late times when extrapolated.
• From Fig. 20 the C 2 behaviour is seen to be similar although the oscillations have a single defined smoother trough and peak. We do not see continuity at point of bounce, but we see the similar rising behaviour of post point of bounce line as it approaches boundary and saturates upon extrapolation. We can also observe the behaviour of the predicted OTOC values from the complexities. The behaviour is very similar to that of the one we saw with scale factor.
• outside the bouncing region the values are highly oscillating and produce random fluctuations upon extrapolation. The oscillations are smoother(smooth peaks and troughs) in the case of OTOC from C 2 . See Fig. 23 and Fig. 25.
• Inside the bouncing region as before we have extrapolated the post point of bounce line that decreases steadily from near the boundary to the extrapolated value of N pr esent where it saturates as seen in Fig. 22 and Fig. 24 This is the signature of a chaotic system. The pre point of bounce line shows random fluctuations and an increasing value of OTOC if extrapolated.

Exponential bounce
We have numerically plotted the squeezing parameters and the derived Complexity measures for Exponential bounce model in four different regions -pre bounce boundary, pre point of  bounce, post point of bounce and post bounce boundary against the scale factor. In Sec. 3, the scale factor of the exponential model has been plotted with respect to the conformal time. It can be seen that at present time and at a time much before the boundary (τ → −∞) the value of scale factor a = 1. We have taken the value of pre-boundary and post-boundary parameters r k (a = 1) = 1, φ k (a = 1) = 1 to set our initial conditions, and ensured continuity at a boundar y as initial conditions for the bouncing region parameters for numerically solving differential equations with respect to scale factor (Eqs. D. 15 and D.16). For the analysis of exponential bounce we have taken −kτ b = 30 and the range of −kτ goes from 0 to 60. We have fixed the value of Λ to 10 −4 for our numerical analysis in this case as well. • In Fig. 26 the behaviour of the squeezing parameter r k has been plotted with respect to the scale factor of the model for four different regions of interest. We have taken the value of pre-boundary and post-boundary parameters r k (a = 1) = 1 to set our initial conditions, and ensured continuity at a boundar y as initial conditions for the bouncing region parameters for numerically solving Eq(D.15), and Eq(D.16).
• The pre-boundary and the post-boundary behaviour of r k is oscillatory with decreasing amplitude as it approaches a = 1.
• The behaviour of r k for the post bounce region can be seen to be highly oscillatory near the point of bounce and the amplitude of oscillations reduces significantly near the  boundary. However the behaviour of the squeezed parameter r k is almost constant with minor fluctuations in the region between pre bounce boundary and the point of bounce.
In Fig. 27 and Fig. 28, the squeezing angle and the sine of twice its value has been plotted with respect to the scale factor.
• The squeezed angle parameter φ k slows an exponential increase starting from almost zero, for the pre-bounce boundary region, with the rate of rise decreasing as the parameter approaches the boundary of the bouncing region from very early times. The sine of twice the angle of the squeezed parameter in this region is a periodic function with the frequency of oscillation decreasing as it approaches the boundary of the bouncing region.
• φ k shows a asymptotic rise for the pre point of bounce region and the sine of twice the angle of φ k shows wild oscillation in this region.
• The squeezed angle parameter φ k slows an exponential increase for the post point of bounce region and the sine of twice the angle of φ k shows oscillatory behaviour with the frequency of oscillation larger near the point of bounce than near the boundary.
• An exponential decay of the squeezed angle parameter can be seen for the post bounce boundary region with the value approaching zero for the present day. The sine of twice the angle shows regular oscillatory behaviour with the frequency of oscillations increasing as one approaches the present day. In Fig. 29 and Fig. 30 the complexity measures C 1 and C 2 have been plotted with respect to the scale factor. The behaviour of both the complexity measures are almost identical.
• The value of complexity outside the bouncing boundary based on the respective squeezing parameters defined there for very early times and nearing present times is oscillatory with smaller frequency at the boundary. The amplitude of oscillation however remains almost identical for the early and nearing present times to that at the boundary.
• Inside the bouncing region the complexity before the pre-bounce shows wild oscillations.
• The behaviour of the post point of bounce shows interesting features on extrapolation to the present day. Although within its domain i.e inside the bouncing boundary region the behaviour is oscillatory as the pre-point of bounce behaviour but extrapolation to the present day shows exponential increase in the complexity after a certain time. This is different to what we observed in the cosine hyperbolic model where the rising behaviour of complexities was seen inside the bouncing region before the post bounce boundary itself.
In Table 5 we have presented the values of complexity at various time scales as observed in this exponential bouncing cosmology model.  From the complexity values shown in Table 5, one can interpret that the system tends towards a highly chaotic behaviour near the point of bounce, which is understood from the maximum complexity value at that point. In the language of squeezed quantum states, the rapid oscillation of the squeezing parameters near the point of bounce may be an indirect way of signifying chaos. Since we observe most interesting features from the post point of bounce plots(on extrapolating to the present time), the prediction of Lyapunov exponent from that case is extremely useful as it gives an estimation on the lower bound of temperature. The slope of logarithm of complexity at the point of rise directly corresponds to the value of the quantum Lyapunov exponent as mentioned in Eq (140). To predict the slope of the logarithmic value of complexities we consider the change of y-axis value over the range of the x-axis value i.e. between point of rise and point of saturation. Mathematically it is represented by For this we have plotted the logarithm of complexity values in Fig. 31 and Fig. 32. We observe the qualitative features to be same as that of the complexity graphs, showing corresponding oscillatory, rising and saturation at the respective regions. We calculate the Lyapunov exponent from the post point of bounce case as it shows exponential feature and this gives an estimation on the lower bound of temperature. In Table 6 we have written the numerical values of the logarithm of the complexity values at the region of a where complexity shows an exponential rise. The Lyapunov exponent calculated from these values are: The estimated lower bound on the temperature from the calculated values of the Lyapunov exponents are Using Eq (166), similar to the cosine hyperbolic case we have also computed lower bound of scrambling time intervals for the exponential case. The main difference is that as we have seen for the exponential case the complexity values do not actually saturate even at much later times. Hence we have calculated the scrambling time in the region of rise (the same region that we have numerically considered for calculating the Lyapunov exponent in exponential case). It is unclear whether the physical interpretation of the scrambling time will remain same as there is no given region of saturation in the exponential case as we had for the cosine hyperbolic case. Nevertheless an estimated value for the same in the region of rise is given in Table 7. Since our normalization for conformal time at the exponential case is different with bounce at τ b = −150 and present time at 0, we can see that the scrambling period is more than one-half of time from bounce to present day. Such a high scrambling time (more time for OTOC to attain equilibrium) can be due to the lack of saturation and late and perpetually rising complexity values, and a never hence a never saturating OTOC. Hence the interpretation of scrambling time in rising region might point to fact that it takes really long (almost never) for OTOC to attain equilibrium hence hinting at the lack of saturation in the chaotic behaviour that we have seen in complexity in the case of cosine hyperbolic model.  Fig. 33 the predicted OTOC from complexity measure C 1 has been plotted with respect to the scale factor. We observe that the outer bouncing boundary curves shows similar features and is oscillatory with the amplitude of post bounce boundary region always lying above the pre-bounce boundary region. Inside the boundary region OTOC shows rapid oscillations for both the pre and the post point of bounce. However interesting features is observed from the post point of bounce curve. The frequency of oscillations starts decreasing near the boundary of the bouncing region. On extrapolation to the present time the OTOC predicted from C 1 actually shows an exponential decay.
• In Fig. 34 the predicted OTOC from complexity measure C 2 has been plotted with respect to the scale factor. We observe that the outer bouncing boundary curves shows similar features and is oscillatory but unlike the OTOC predicted from C 1 . The amplitude of post bounce boundary region is almost identical to the pre-bounce boundary region.
Inside the boundary region OTOC shows rapid oscillations for both the pre and the post point of bounce. The frequency of oscillations for the post point of bounce curve starts decreasing near the boundary of the bouncing region. On extrapolation it to the present time the OTOC predicted from C 2 actually shows a similar exponential decay as OTOC predicted from C 1 . In Table 8 the numerical values of the predicted OTOC's from both the complexity measures from our present analysis is written. Again at the point of bounce the OTOC shows a drastic reduction in the values.
As discussed in the previous section it is better to relate the complexity with some quantity which is observable. The number of e-foldings is one such observable quantity. Again we plot the complexity corresponding to within the bouncing region and the outside boundary region separately as done for Model-I.
• In Fig. 35 the linearly weighted complexity measure C 1 has been plotted with respect to the number of e-foldings inside the bouncing region. We observe that the behaviour   near the boundary is some random fluctuations of negligible amplitudes, however as it approaches the point of bounce it fluctuates wildly. Similarly the behaviour post point of bounce is arbitrary and random near the point of bounce whereas it takes a regular periodic shape on approaching the boundary. However the interesting part can be realised on extrapolating the post bounce behaviour to the present times. We see a sudden exponential rise as it approaches the present times.
• In Fig. 36 we have plotted the complexity measure C 1 as a function of the number of efoldings for outside the bouncing region. For both the pre and the post boundary region we see a smooth, regular and periodic behaviour of the complexity even on extrapolating it inside the boundary region. However an important feature to notice is that the frequency of oscillation for the pre boundary region decreases when it approaches the boundary whereas the post boundary behaviour shows a contrasting behaviour as the complexity approaches the present time. Also the value of the complexity for the prebounce boundary region is always greater than the post bounce boundary region. • In Fig. 37 we have plotted the geodesically weighted measure of complexity with respect to the number of e-folds for inside the bouncing region. We observe that the behaviour of C 2 in this region is almost identical to the behaviour of the linearly weighted measure of complexity. Even the extrapolated behaviour of the post point of bounce is identical to that of C 1 .
• We observe similar periodic behaviour in the complexity measure C 2 as C 1 . However an important observation is that unlike C 1 , C 2 has equal values for all regions even when extrapolated inside the bouncing region. Table 9 contains all the key features of the complexity measures C 1 and C 2 and the Out of time ordered correlation functions predicted from them for both of the Bouncing cosmological models.

Conclusions
From our study of the complexity measures computed from the linearly weighted and geodesically weighted cost functionals within the framework of bouncing cosmology we have the following final remarks: • Remark I: The complexity measure calculated from two different types of cost functionals has an overall identical behaviour for both the models of bouncing cosmology under consideration with some noticeable differences. Though their behaviour is identical, it is evident  from the plots that for a particular value of the scale factor, the linearly weighted complexity measure (C 1 ), is always greater than the geodesically weighted one (C 2 ). We see this feature in both the models.
• Remark II: We observe that the complexity measure calculated for the post-point of bounce for both the models is the most interesting one. Though we observe random quantum fluctuations inside the bouncing region, once extrapolated to the present day we observe an exponential rise in the complexity measures followed by its saturation. An important point worth noting is that for the Cosine Hyperbolic bounce model, the starting point of the exponential rise in the two complexity measures is observed inside the bouncing region itself and the saturation is well observed on extrapolation. However, the rise in the two complexities for the Exponential bounce model is observed only after extrapolating it outside the bouncing region. We do not observe saturation as such in the complexity measures even on extrapolating it to the present day for the Exponential bounce model. We observe similar behaviour of the complexities when the analysis is done concerning the observationally important quantity known as the number of e-folds.
• Remark III: The behaviour of the complexity measure outside the bouncing region is not of prime significance as we observe smooth, well behaved and periodic nature in those regions. Though the periodicity may not be equal near the boundary of the bouncing region and the present time, it does not give us any random or chaotic behaviour in those regions. However, the behaviour of complexity measure in the pre-point of bounce region is of some significance because we observe random oscillations of the complexity in this region. For the exponential case, these oscillations are extremely wild and it can be attributed to quantum mechanical fluctuations.
• Remark IV: We observe an exponential decay in the predicted OTOC's computed from the complexity measures C 1 and C 2 , which is in accordance with the recently established predictions of OTOC's in the context of Cosmology [43]. This behaviour of the OTOC's is observed not only with respect to the theoretical measure of scale factor but also with the respect to the observational measure of the number of e-folds.
• Remark V: The Lyapunov exponent calculated from both complexities at the region of rising only differ in the second decimal place for both models. The fractional variation between the Lyapunov exponents is observed to be less than ten percent. We expect the variation to be much lesser if we consider higher-order terms in the dispersion relation for the numerical analysis.
• Remark VI: Since we have solved the same dynamical equations with respect to the Number of efoldings, another calculation of Lyapunov exponent with respect to the number of efoldings, is expected to have the same order and hence a separate calculation is redundant.
• Remark VII: We get a theoretical prediction for the lower bound of temperature for both models in the region of rising which falls from before the boundary till late times for the hyperbolic cosine case and completely outside the boundary and at much late times for the exponential case.
• Remark VIII: The choice of initial conditions at the horizon crossing that we have chosen for evaluating the perturbed action can have significant changes to the obtained results in the form of complexity, OTOC, and the Lyapunov exponent.
• Remark IX: The scrambling time scale for both cases has been estimated. For the cosine hyperbolic case (estimated in the region of saturation) it turns out to be one-tenth of the conformal time period from bounce to present day, whereas for the exponential case the estimated scrambling time (in the region of rise), turns out to be more than one-half of the conformal bounce-to-present day time period, which might signify the lesser time required for OTOC to attain equilibrium in the cosine hyperbolic case (also pointing towards saturation of complexity), whereas a long time required for exponential OTOC to attain equilibrium hence hinting at a never saturating complexity value at late times.
The future prospects of the present work are appended below: • Prospect I: The framework for bouncing cosmology that we have mentioned along with the generalized perturbed action, dispersion relations, and the Hamiltonian can be used in a more general way without any truncation of higher-order terms for investigating cosmological complexity in any models. Though we have focused on the sub-Hubble region due to its quantum fluctuations, numerically one can use the full solution, rather approximating it in the sub-Hubble region. One can also apply the same analysis for the generic inflationary paradigm, which has not been considered in any work yet in an appropriate way.
• Prospect II: We are working currently on applying the same framework to certain models in Island Cosmology as discussed in [224,225]. We will be planning to work out on this framework with different quantum initial conditions which will appear very soon in an upcoming paper.
• Prospect III: An interesting study would be to see how this complexity can be used to study the non equilibrium phenomenon and chaoticity in various entangled systems [226][227][228][229][230][231][232]. This measure might be used to see if the long-range correlation between systems induces chaoticity and quantum randomness or not.
• Prospect IV: It is naturally expected that chaos and randomness might be an inherent property of Open quantum systems [233][234][235] depending on the properties of the quantum dissipation and its impact. Complexity finds another use in this direction of research where one might be inclined to study chaos in nature from a more realistic point of view.
• Prospect V: It is also possible to find a general representation of squeezed state formalism for the multi-field interacting scenario in the context of cosmological perturbation theory written in a spatially flat FLRW background. But till now there are no concrete results available on this in the cosmology literature because the result is completely dependent on the type and strength of the interaction and extremely model dependent. For this reason, it is extremely difficult to deal with such types of computations within the framework of the quantum field theory of cosmological perturbation theory. But apart from having these mentioned difficulties, if one can write a general structure of the squeezed state formalism at least for two interacting fields in the spatially flat FLRW cosmological background by considering all possible general renormalizable coupling and interactions in the perturbative regime of the quantum field theory then it is possible to compute many physical observables out of those results. Using the interacting two field squeezed state formalism it is also good to understand the role of quantum mechanical chaos and complexity within the framework of cosmological perturbation theory. The future aim should be to carry forward such computations and explore some of the important unknown important underlying physical features of cosmological perturbation theory in presence of interacting quantum mechanical fluctuations.
• Prospect VI: In this paper, we have restricted our analysis only for scalar mode quantum fluctuations generated from a cosmological perturbation in the spatially flat FLRW cosmological background. However, a similar analysis can be extended for the primordial gravitational waves appearing from the tensor mode fluctuations in the same cosmological background set up. It would be really nice and also important to check how the primordial gravitation waves and related tensor mode fluctuations get affected by two modes squeezed state formalism and also important to study how that will further put a stringent constraint on the phenomena of quantum mechanical chaos and complexity.
• Prospect VII: Generally, people try to explain this concept of chaos and complexity through various models. However a completely model-independent notion of complexity can be given from the perspective of effective field theory [236,237]. It is in general possible to start from the single EFT action and derive all the models under various constraints satisfied by the parameters of the action. Squeezed state formalism for such a universal action can be developed to generalize an give a model-independent prescription of complexity.

Note:
The background Gaussian image in the table of contents is inspired from [238].

A Quantization of Hamiltonian for scalar modes in terms of squeezed parameters in cosmological perturbation theory
In this appendix, we present the details of the quantization of the Hamiltonian for scalar modes obtained from the cosmological perturbation theory within the framework of bouncing cosmology. To serve this purpose, let us start with the expression for the classical Hamiltonian function already derived earlier in the paper and it is represented by: where the conformal time dependent mass µ 2 (k, τ) of the parametric oscillator is given by the following expression: Also one can express the field velocity with respect to the canonically conjugate momentum density in the Fourier space as: Next, using the classical mode function we can further construct the quantum mechanical operators:v Now we evaluate the following quantities: Now using the above mentioned quantum operator one can finally express the canonical Hamiltonian for the parametric oscillator in the following quantized form: Contribution from the free term Contribution from the Interaction term where we define Ω k (τ) and φ k (τ) by the following expressions: Here Ω k (τ) represents the conformal time dependent dispersion relation in the present bouncing cosmological set-up and φ k (τ) is squeezing angle appearing in the squeezed state formalism as discussed earlier in the text.

B Hamilton's equations in the Heisenberg picture in cosmological perturbation for scalar modes
Next, using the previously mentioned solution of classical mode function we can further construct the quantum mechanical operators in the Heisenberg picture: Now, our objective is to find out the fact that whether the mode functions for both the field variables and their associated momenta in Fourier space satisfy the well known, Hamilton equations or not. In the Heisenberg picture one can write down the following where the field operator and the corresponding canonically conjugate momentum operator can be expressed in terms of the creation and annihilation operators of conformal time dependent parametric oscillator as:v 8 Here we have explicitly used the following operator identity which is valid in the Heisenberg quantum mechanical picture: As a result we get the following simplified form of the equation of motion: which is the generalized version of the well known Mukhanov Sasaki equation. In the sub-Hubble region (−kτ 1) one can simplify the expression for the dispersion relation, Ω k (τ), which is explicitly discussed in the next section. Now considering the leading order contribution we get the following expression for the conformal time dependent frequency parameter: which is exactly appearing in the Mukhanov Sasaki equation.

(C.4)
Here we have chosen the initial condition at the time scale τ = τ 0 by considering the horizon crossing scale, −kτ 0 = 1. We impose this condition on the perturbation field variable and on the canonically conjugate momentum obtained for scalar fluctuation. We finally get: Neglecting the phase factors in the above equation and also noting that ν B ≈ 1 2 + · · · , we get a pretty simplified expression for Ω k (τ), i.e.

C.1 Sub-Hubble limiting result
In the sub-Hubble limit, −kτ 1, it is expected to have very small contribution from the squeezed parameter, r k (τ) for which one can use the following approximations: cosh r k (τ) ≈ 1, sinh r k (τ) ≈ r k (τ). (C.8) Consequently, in the limit r k (τ) → 0, we get the following result for the dispersion relation in the sub-Hubble region: which is basically dependent on the co-moving wave number and a very slowly varying time dependent quantity ν B at the sub-Hubble scale. In the previous ref. [11,12] the authors have not considered this additional slowly varying time dependence appearing from the parameter ν B , which is not appropriate if we want to extract the information regarding quantum correlation in the out-of-equilibrium phase where random fluctuations play significant role. Now we will explicitly show how the slow time dependence is appearing in the parameter ν B . In the sub-Hubble region the conformal time dependent mass parameter ν B can be approximately written by considering the contribution upto the next-to-leading order as: where we have neglected the contributions of all higher order small correction terms appearing as · · · for the computational simplicity. But out of all the terms in the correction part, H /H 2 term gives the most significant contribution. Due to slowly varying feature with respect to the in the super-Hubble region: = Ω Sub k (τ) 1 + 2 sinh r k (τ) sinh r k (τ) − 1 6 , (C. 16) where Ω Sub k (τ) is the dispersion relation derived in the previous sub section in the sub-Hubble region. In the super-Hubble region one needs to consider the contributions appearing in the bracketed terms in the above derived expression during the study of the evolution with respect to any dynamical parameters involved in the system.
One can further consider a more simpler situation in the super-Hubble region, which is described by very small value of the squeezed parameter, r k (τ) which is fixed by the following contribution: sinh r k (τ) ≈ r k (τ). (C.17) Here we consider r k (τ) to be small but not approaching zero and also we neglect the quadratic contribution in r k (τ) due to smallness approximation. As a result, finally we get the following simplified answer for the dispersion relation in this specific situation: where one need to consider the contribution from the second term in the evolution equations and this ensures the fact that, Ω Sup k (τ) = Ω Sub k (τ) in this limiting situation.
(C. 22) In the above mentioned discussion it is clearly evident that to match the dispersion relation obtained from the sub-Hubble and super-Hubble region at the horizon crossing the squeezed parameter has to be either, r k (τ 0 ) = nπ, ∀n ∈ , or r k (τ 0 ) = sinh −1 1 6 .

D Equivalent representations of the evolution equations in twomode squeezed state formalism
In this section, we will discuss about three equivalent representations of the evolution equation of the squeezed parameter and squeezed angle using which one can study the impact of the two mode squeezed state formalism in the present bouncing cosmological set up which is described in the spatially flat cosmological FLRW background. The details of each of the three representation has been discussed in the following three consecutive subsections respectively.

D.1 Representation I: In terms of conformal time
The time evolution equations of the conformal time dependent squeezed state parameter r k (τ) and squeezed angle θ k (τ) are given by: The above set of evolution equations, are coupled differential equations of squeezed state parameter r k (τ) and squeezed angle θ k (τ) where in both conformal time derivatives are involved. We choose the initial condition is at the horizon crossing scale, −kτ 0 = 1 at τ = τ 0 and also consider the sub-Hubble (−kτ 1) region for the computational purpose, where the scalar modes for two momenta k and −k having all possible values becomes quantum in nature. Using these information one can numerically solve these equations to construct the target quantum state out of a Gaussian initial state. This will further help us to numerically compute and understand the quantum complexities in Eq (103) and Eq (104) within the framework of primordial cosmological perturbation theory, where the effects of the quantum fluctuations is treated in terms of the squeezed state parameter r k (τ) and squeezed angle θ k (τ) in the squeezed state formalism. Now we will discuss about the strong and the weak coupling region and the behaviour and the physical outcome of these evolution equation:

Strong coupling region and freeze-out phenomena :
In the strong coupling region the effect of squeezing phenomena become maximum because:

Weak coupling region and oscillation phenomena :
In the weak coupling region the effect of oscillation phenomena become maximum because: As a result we approximate: where β k = sin −1 λ k Ω k 1. (D.9) As a result we get following simplified form of the evolution equations: Consequently, we get the following analytical solution: For the cosmological models when the modes appearing from the cosmological perturbation lies within the horizon, the above mentioned solutions works perfectly well. On average, the squeezing parameter r k (τ) during this time is almost constant and the perturbation do not grow at all.

D.2 Representation II: In terms of scale factor
In this section instead of using the conformal time τ as the dynamical variable, we have chosen the scale factor a(τ) to make the computation simpler and physically justifiable. To perform the change in variable from τ to a(τ) we have to replace the following differential operator in the above mentioned evolution equations using the chain rule, as: . (D.14) Consequently, the evolution of the squeezed state parameter r k (a) and squeezed angle θ k (a), can be recast in terms of the newly defined dynamical variable a(τ) as: Once we numerically solve the evolution of the squeezed state parameter r k (a) and squeezed angle θ k (a) in terms of the scale factor a, we can construct the target quantum state out of a Gaussian initial state. This will further help us to numerically compute and understand the quantum complexities in Eq (103) and Eq (104) within the framework of primordial cosmological perturbation theory.

D.3 Representation III: In terms of co-moving Hubble radius/ number of efoldings
Now if we think about the more realistic cosmological observation then it is not very good to study the evolution with respect to the scale factor, because in the context of realistic cosmology the scale factor is not the direct physical observable that can be probed in the observation for various cosmological missions running (or supposed to run in the near future) to test the signatures of the primordial cosmological paradigm. In this specific situation one needs to use the following transformation for which the linear differential operator appearing in the evolutionary equations of the squeezed parameter and the squeezed angle will be modified as: where we have used the following couple of facts for the above mentioned transformation: Here, N is the actual number of e-foldings, N is the number of e-foldings in terms of the re-defined variables, and ε(τ) is the slowly varying conformal time dependent parameter. Consequently, the evolution of the squeezed state parameter r k (N ) and squeezed angle θ k (N ), can be recast in terms of the newly defined dynamical preferred choice of suitable variable N as: In this context, r co = (aH) −1 or r co = H −1 represents the co-moving Hubble radius, which is extremely important quantity in terms of which the newly re-defined number of e-foldings have been expressed in terms of the good old definition of the number of e-foldings.