Circuit Complexity as a novel probe of Quantum Entanglement: A study with Black Hole Gas in arbitrary dimensions

In this article, we investigate the quantum circuit complexity and entanglement entropy in the recently studied black hole gas framework using the two-mode squeezed states formalism written in arbitrary dimensional spatially flat cosmological Friedmann-Lema$\hat{i}$tre-Robertson-Walker (FLRW) background space-time. We compute the various complexity measures and study the evolution of these complexities by following two different prescriptions viz. Covariant matrix method and Nielsen's method. Independently, using the two-mode squeezed states formalism we also compute the R\'enyi and Von-Neumann entanglement entropy, which show an inherent connection between the entanglement entropy and quantum circuit complexity. We study the behaviour of the complexity measures and entanglement entropy separately for three different spatial dimensions and observe various significant different features in three spatial dimensions on the evolution of these quantities with respect to the scale factor. Furthermore, we also study the underlying behaviour of the equilibrium temperature with two of the most essential quantities i.e. rate of change of complexity with scale factor and the entanglement entropy. We observe that irrespective of the spatial dimension, the equilibrium temperature depends quartically on entanglement entropy.


Introduction
Circuit Complexity has become a helping hand to not only the high energy physics community but also to the people from other branches as well . This quantum information theory technique has been significantly used recently to probe many features which were previously difficult. Though this concept is a computation tool, its contribution in the field of physics of late is massive. It provides a way to probe physics behind the horizon of black holes through the use of the "Com-plexity=Volume" and "Complexity=Action" conjectures [25][26][27][28][29]. Since then it has been extensively used in quantum field theory and in studies involving AdS/CFT correspondence [30]. These holographic approaches connect a probe on the gravity side with a concept of quantum information theory.
In the recent past, along with the out-of-time ordered correlation functions [31][32][33][34][35], it has formed the web of quantum chaos. It has been found to reveal essential information like Lyapunov exponent, scrambling time etc. required to diagnose chaos in a system. Many interesting works have been done using this tool in wide areas of physics. It was studied for cosmological islands in [36], where the authors showed that entanglement entropy from circuit complexity via a famous relation proposed in [25] resembles the page curve in some particular regime where a universal relation between circuit complexity, OTOC and entanglement entropy can be written. It was used to study early universe chaos within the framework of bouncing cosmology [37]. People have also computed circuit complexity in the context of supersymmetric quantum field theory [38]. The connection between entanglement and emergence of space-time has been an active area of study, where the entanglement entropy is the minimum cross-sectional area of an Einstein-Rosen Bridge(ERB). However, classically the ERB continues to grow for a very long time, whereas the dual thermodynamic system comes to a thermal equilibrium quickly. This led Susskind to introduce a new variable namely Complexity which could be responsible for the ERB growth [25]. Complexity of a quantum system is a real quantity and its growth rate has been conjectured to be proportional to the entropy of the black hole based on these observations. It was very recently shown in [39] that there exists some relationship between entangling power and circuit complexity. Most importantly, if the entanglement entropy grows linearly with time, the geometric circuit complexity also grows linearly.
In this paper, we will study the evolution of complexity w.r.t. the entanglement entropy of the blackhole gas model [40]. An important feature of this model is that the total entropy of the blackhole gas is directly proportional to the volume of the system instead of the area and the system behaves like a thermodynamic gas. We will compute the most common measures for determining the entanglement between the squeezed states of blackhole gas, namely the von Neumann and Rényi entropy, which quantifies the amount of uncertainty linked to the density matrix. We will measure the entanglement entropy by constructing an effective thermal representation for the reduced single-mode state from the two-mode squeezed state which varies linearly w.r.t. to squeezing parameter r. Then by using the scale factor predicted by the blackhole gas model in the flat space-time metric as a dynamical variable we will study the evolution of complexity in three spatial dimensions and compare it with the entanglement entropy in terms of squeezed state parameters. One of the reasons we are interested in the blackhole gas is that its equation of state describes a universe right after the big bang and before the start of inflation, if one wants to avoid such an equation of state governing radiations of very high densities then one needs to have inflation right from the Planck scale.

Key Highlights
• The behaviour of circuit complexity calculated from two different approaches viz. Covariance matrix method and Nielsen's wave function method has been studied w.r.t. scale factor for the blackhole gas model. We observe interesting behaviours for different spatial dimensions.
• The behaviour of the Von Neumann entropy and the Rényi entropy has been studied w.r.t. the scale factor for different spatial dimensions. We observe similar features for d = 1, 2 whereas for d = 3 we observe slightly different behaviour.
• The behaviour of dC/da w.r.t. von Neumann and Rényi entropy has been studied for different spatial dimensions. It is shown that by no means is it a linear function.
• The behaviour of the equilibrium temperature for the blackhole gas model is identical compared to entropy for different spatial dimensions whereas it seems dependent on spatial dimension when compared to dC/da.
Organization of the paper The organization of the paper is as follows: • We begin by providing a review of the Black Hole gas given by Samir Mathur in [40] in the section A Short Note on Black Hole Gas. Solving the Friedmann equation for blackhole gas, we relate the scale factor a(t) with spatial dimension.
• In sec:2 Black hole gas perturbation theory in d+1 dimensions, we investigate the squeezed state formalism by perturbing the black hole gas geometry in FLRW spatially flat spacetime..
• In sec:3 A Short Note on Circuit Complexity, We review of the Circuit Complexity in the section. We discuss the geometric framework of circuit complexity developed by Nielsen and collaborators.
• In sec:4 Circuit Complexity of two-mode squeezed states, after introducing the notion of squeezed states, we calculate it's circuit complexity using two approaches: Complexity using the covariance matrix and Nielsen's method of wave functions.
• In sec:5 Entanglement entropy of two-mode squeezed states, we compute the entanglement entropy of two-mode squeezed states. We compute Rényi-entropy, Von-Neumann entropy and Rényi-2 entropy. We find entanglement entropy grows linearly with increasing squeezing parameters for the short time period. Then, we compare entanglement entropy with circuit complexity obtained in sec:1.
• In sec:6 Numerical Results, we numerically study the behaviour of the complexity measures and entanglement entropy separately for three spatial dimensions (d= 1,2 and 3).
• In sec:7 Conclusions, we conclude with some discussions.

A Short Note on Black Hole Gas
In this section, we review a model proposed in ref. [40], where the author have studied the state of a system moving towards maximal entropy S. Note that this system is unlike the inflationary model where we have a low entropy state after the inflation because the positive energy of the matter content is compensated by the gravitational potential. Here, we provide a quick derivation of the equation of state that describes the pre-inflationary universe. We will consider a configuration where we find entropy S(E, V ) of a system in a toridal box of volume V in the limit E → ∞. At low enough energies the matter phase corresponds to radiation whose entropy as a function of E and V in dimension d can be given by S ≈ V ρ d d+1 . By putting more energy into the box one can look for a configuration where the system turns into a blackhole of radius R whose entropy is given by One might suspect that for a given box of radius R, eqn.
(1) describes the state with maximum entropy for energy E = E bh , since throwing more energy into the black hole will only result in increasing the hubble expansion. However if we let go the constraint that the energy inside the box shouldn't be greater than the mass of the blackhole inside the box then one could arrive at a configuration where the entropy of the system is greater than (1). This could be achieved by putting N number of black hole each of radius R in a lattice instead of a single box of volume V . The number of black holes in such a configuration is given as Therefore the total entropy of the system is We notice that (3) is in contrast with (1) where the entropy is directly proportional to the volume rather than the area of the horizon. The energy leading to such a state is undoubtedly greater than E bh and could be expressed as follows substituting the value of R from (4) to (3) we get for ρ = E/V we get, This contrast in the definition of entropy is subjected to the constraint that microstates cannot expand freely to larger size unlike in asymptotically flat space where the entropy is given by the area law. Also the resulting lattice configuration with E > E bh , having N number of black holes wouldn't collapse to form one large blackhole as the entropy corresponding to the lattice configuration is larger than a single black hole state in a box of volume V . Now from the first law of thermodynamics we can show Now from (8) we see that equation of state takes the form ρ = wp with w = 1. The solution for the black hole gas model can be obtained by starting with the FLRW flat metric in 1+d dimensions which is given by the following line element Solving the Friedmann equation for black hole gas in the d + 1 dimensional flat metric with scale factor a(t) we get The above FLRW flat metric can be written in terms of the conformal time scale by using the following conversion relation Further integrating the both sides of the above equation we get the following relationship between the physical time t and the conformal time τ in the arbitrary d + 1 dimension black hole gas: This relationship is extremely useful for the present computation which helps us the directly translate the information in terms of desired conformal time from physical time.
In this conformal time coordinates, the flat FLRW line element gets transformed as Hence, the above solution of the black hole gas scale factor can be written in terms of conformal time as follows The corresponding Hubble parameter with respect to the conformal time scale can also be written as: Now we will describe an equivalent scenario inside a black hole gas where one can obtain the above mentioned scale factor and Hubble parameter in any arbitrary d + 1 dimensions. In this scenario we embed a scalar field within the framework of Einstein gravity having the previously mentioned spatially flat FLRW space-time, where the solution of the scale factor is exactly same as mentioned earlier. In this framework the representative action of the scenario is described as: where we have fixed the Planck mass M p = 1 for the present computation. First term in the action represent the usual Einstein Hilbert term, second term represents the kinetic term of the embedded scalar filed and the last term represents the potential function in a black hole gas V (φ) in any arbitrary d + 1 dimensions. We have found the the following two possibilities of the potential functions are allowed in the present context which can finally give rise to same same scale factor which we have mentioned earlier: Here it is important to note that, since we have embedded the scalar field in the homogeneous and isotropic spatially flat FLRW background it turns out to be the field is only function of the time coordinate. Further solving the Klein Gordon field equation in d+1 dimension spatially flat FLRW background the dynamical solution of the homogeneous and isotropic background scalar field φ in terms of the conformal time coordinate can be expressed as: This solution actually representing the dynamics of the field inside the black hole gas in d + 1 dimension spatially flat FLRW background.  the context of black hole gas: In Fig. 1 and 2, we have plotted the two choices of fields with respect to the conformal time. In choice I, we see that for d = 1 the field keeps on increasing monotonously and for higher dimensions i.e. d = 2, 3 the field value decreases slowly initially w.r.t. conformal time and then rapidly as we increase the value of conformal time. The exact opposite behaviour can be seen for choice II. For d = 1 the field value keeps on decreasing monotonously, where as for d = 2, 3 field value grows steadily and then at a faster rate with conformal time.
In Fig. 3 and 4, we have plotted potential for the two choices against the field. As can be seen from the Eq. (17), the potential function takes the form of exponentially decreasing and increasing function for first and second choices respectively.
In Fig. 5, we have drawn potential function with respect to conformal time. We have taken the logarithm for the potential. It can be seen that for d = 1, the potential goes on decreasing continuously in a straight line where as for d = 2, 3 the potential rises quickly after going through a slow rise period in the beginning.
In fig: 6, we have plotted the scale factor of the black hole gas model w.r.t the conformal time scale. The plots have been done by fixing the value of the constant a 0 to 1. We observe a significant difference in the behavior of the scale factor for the spatial dimension d=1 and the higher spatial dimensions. The scale factor corresponding to the spatial dimension d=1 shows an increasing behavior in the late time scales. The scale factor for higher spatial dimension shows a decreasing behavior. The decrease is linear for the spatial dimension d=2 whereas it is non-   infinity at late conformal time scales (near τ =0) for the other spatial dimensions. Moreover, the value taken by the Hubble parameter for the higher spatial dimensions is always negative.
One could also arrive to (6) by using the T and S duality symmetries of the string theory, this seems interesting because it doesn't require the state of blackhole gas to be state near the big bang since microstates of stringy theory could also describe the microstates of blackhole.

Black hole gas perturbation theory in d+1 dimensions
In this section we will study squeezed state formalism within the framework of black hole gas theory for FLRW spatially flat background. In this context one needs to consider the following perturbation in the scalar field: 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: δφ( x, τ ) = 0, 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 introducing the Mukhanov variable, defined as v( x, τ ) = z(τ ) ζ( x, τ ) where, z(τ ) = a d−1 2 (τ ) 2 (τ ), the second order perturbed action can be rewritten as: where the quantity (τ ) is known as the conformal time dependent slow roll parameter and is defined as For the black hole gas model, it is very easy to verify that the slow roll parameter is equal to the dimension in which the black hole model is being considered i.e.
which is finally independent of conformal time coordinate τ .
By implementing the following ansatz for the Fourier transformation: v( x, τ ) : the second-order perturbation for the scalar modes in Fourier space can be further recast as: Now, varying the above second order action, we get the following equation of motion The above equation is known as the Mukhanov-Sassaki equation with the frequency of the oscillator given by The conformal time dependent effective mass in the present computation is given by with the conformal time dependent mass parameter for the black hole gas in arbitrary dimension upto leading order terms is given by where we can clearly observe that the slowly varying conformal time dependence is appearing from the third term where we have truncated the expansion. During computing this mass parameter have used the following facts: where the explicit expression for the Hubble parameter in conformal time coordinate is quoted in the previous section. Finally using this result we get the following simplified expression for the mass parameter for black hole gas: A general solution to the equation of motion is written as νBHG (−kτ ) and H (2) νBHG (−kτ ) are Hankel functions of the first and second kind,respectively, with argument −kτ and order ν BHG . The two integration constants can be fixed by the choice of various initial conditions. In this paper, we restrict our choice to only the Bunch-Davies vacuum case in which one fixes the initial conditions, C 1 as 1 and C 2 as 0. However, it is generally difficult to work with these full solutions and one takes the asymptotic limits of the solutions which is given by where the function f d (k, τ ) takes the following form The normalization with the factor f d (k, τ ) has been done in such a way that for spatial dimension d = 3, it becomes functions with respect to the conformal time. We observe an identical behavior for both the real and imaginary part of the mode functions. It shows an oscillatory behavior with respect to the conformal time and the amplitude of oscillation decreases with the conformal time evolution. Also, it can be observed that the amplitude of oscillation decreases with the increase in the number of spatial dimension. From the mode functions, one calculates the conjugate momentum to the mode functions and thus constructs the classical Hamiltonian function. By promotion the mode function and the conjugate momentum to quantum mechanical operators in the Heisenberg picture, one quantizes the Hamiltonian which is written as where, the origin of the creation and the annihilation operators can be understood, when one promotes the mode functions and its conjugate momentum to quantum mechanical operators. where the symbols Ω k (τ ) and λ k (τ ) are defined by the following expressions: where the quantity µ 2 (k, τ ) is given by Imposing the initial condition at the horizon crossing (k/H = 1), given by time scale (τ = τ 0 ) which is given by the following conditions one can calculate the quantum operators at any arbitrary scale in the Heisenberg picture. Our next job is to determine the expression of the unitary operator in the context of cosmological primordial perturbations of the scalar modes where the concept of squeezed state formalism plays a significant role. In this approach, the unitary operator is factorized as follows : 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 timedependent parameter, r k (τ ) ,and the squeezing angle or the phase is represented by the time-dependent parameter φ k (τ ). The two-mode rotation operator,R produces an irrelevant phase contribution and is ignored. The ground state of the free part of the above Hamiltonian is taken as the initial quantum state, whereas the two mode squeezed quantum vacuum state obtained by acting the squeezed operator on the initial vacuum is taken as the final target state. The time evolution of the conformal time dependent quantum operatorsR andŜ, described by the Schrödinger equation, gives the following set of differential equations for the squeezing parameters:

A Short Note on Circuit Complexity
One of the challenges in Quantum Information processing is it to find out the efficient circuit for implementing a unitary operation U which can be used to solve a computational problem like Search Algorithm or Shor's factoring [41][42][43]. In computer science, similar term called complexity [44,45] exists, which can be defined as the minimum number of computational gates required to implement a certain algorithm. If we extend this definition to quantum version as minimum number of quantum gates out of basic unitary gates [46] in order to implement a unitary operation U , we get quantum complexity [47,48]. Therefore, understanding the difficulty of implementing such unitary operation U , as a sequence of logical gates, is very helpful and at the same time challenging.
In refs [49][50][51], the authors introduced a geometric approach to compute quantum circuit complexity based on the idea that finding optimal quantum circuit is equivalent to problems of computing geodesics in Riemannian geometry. Here, we define a Riemannian metric on the space of n-qubit operations, and distance d(I, U ) between identity and target unitary operation U is equivalent to the number of quantum gates which is then identified as the circuit complexity. It was shown that minimizing this distance d(I, U ) i.e. finding geodesic length, gives a good measure of complexity. Thus, one can employ well developed tools of Riemannian geometry such as Levi-Civita connection, geodesic, curvature, etc. to analyze the quantum circuit complexity. It is important to note that, even though in refs. [49][50][51] geometric techniques were introduced to study complexity, in ref. [52], the authors have previously used techniques from theory of symmetric spaces to study time-optimal control of quantum evolution.
Let U be a transformation which transforms reference state |ψ R to the target state |ψ T via: The unitary transform U , in the language of quantum computation has a order of unitary gates Q i such that U = Q 1 Q 2 ...Q d where d is the depth of the circuit. We can also introduce the tolerance which tells us whether the transformation is successful: This makes sense because in any practical implications, it is difficult to represent the unitary transformation U exactly as a combination of discrete unitary Q i operations.
Obviously, there exists infinite number of ways to achieve this target state |ψ T from the reference state. The circuit complexity is then the depth of the optimal circuit out of this infinite possibilities.
Motivated from the theory of Hamiltonian control problem, authors in refs. [49][50][51] introduced a geometric approach to compute this circuit complexity which was later used to compute complexity in various quantum mechanical and quantum field theoretic models. Instead of directly counting discrete set of gates required for constructing U , Nielsen's approach is geometric. In this method, with a time-dependent Hamiltonian H(t) one constructs unitary U as: Here, the hermitian operators O I forms the basis for time dependent Hamiltonian H(τ ). Path ordering operator P is another version of time ordering operator which indicates that the circuit, made out of non-commuting operators, is built from right to left. The right to left application of operators is a choice of convention. The control functions Y I (τ ) can be thought as particular gates added at a time s represented by Eq. (51).
One can then construct paths in the space of unitaries as: The most interesting case is when trajectory satisfy the boundary conditions U (τ = 0) = 1 and U (τ = 1) = U . These O I and Y I (τ ) satisfy: Eq. (53) is actually just a time-dependent Schrödinger equation in the form when one solves it via time ordered exponentials.
As we discussed before, there are infinite number of ways of achieving unitary transformation U . However, not all processes are optimal. In order to find out the optimal transformation, a cost function F (U, Y (τ )) is defined. This cost function F (U, Y (τ )) is a local functional along the trajectory of the U (τ ) and tangent vectors Y (τ )). Now, for each path the cost is defined as: Nielsen showed that, the variational geometric approach of minimizing this functional is basically finding the optimized quantum circuit. Standing on the physical grounds, the cost function F should satisfy certain properties. Those are appended below.
• Continuity: F ∈ C 0 i.e. F should be continuous. It is reasonable to assume continuity on physical grounds.
• Positivity: Based on the definition of cost function F , it is reasonable to expect where equality holds if and only if v = 0. The equality condition also implies that the reference and target is basically same.
• Positive homogeneity: For any positive real number α and any vector v, we get F (αv) = αF (v).
• Triangle Inequality: F satisfies triangle inequality, for all tangent vectors v and v . Special case is satisfied if and only if v and v are along the same ray coming out from the origin.
If one extends the continuity condition F ∈ C 0 with F ∈ C ∞ i.e F is smooth, then the manifold is known as Finsler Manifold. Nielsen's geometric approach of determining complexity is computing the geodesic in Finsler geometry, and the length of this geodesic gives the complexity. This definition of complexity is also called the geometric circuit complexity.
In literature, there are different choices of these cost functions F (U, v). These choices depends on how one defines the complexity, and the elementary gate sets for their setup. Some of the simple examples are: We can now give comments on various choices of these cost functions. F 1 , linear cost functional, measure is the nearest concept that is close to counting individual gates in the circuit. F 2 , quadratic cost functional, can be considered as the proper distance in the manifold. F 1p can be thought of as a cost function where penalty factors p I are used to favor certain directions over others. This becomes reasonable when one consider elementary gates coupled only to the neighboring qubits and discard those qubits which are non-local. Depending on the system, one can choose different cost functions to study the circuit complexity.
One can also introduce a general class of inhomogeneous and homogeneous family of functionals represented by: where k ≥ 1 represents the degree of homogeneity. F k was introduced in the context of hologrpahy to match the results obtained from "Complexity = Action" and "Complexity = Volume" conjectures.

Circuit Complexity of two mode squeezed states
A simple but a very rich example of entangled multimode field states is two-mode squeezed vacuum state. More details about two-mode squeezed states can be found in [53]. As already defined in the previous section, the two mode squeezing operator is given by: where, ξ = r k e iφ k , r k and φ k is known as squeezing parameters and 0 ≤ r k < ∞ and 0 ≤ φ k ≤ 2π. The two mode squeezed vacuum state, which will act as our target state is given by the action of two mode squeezing operatorŜ k (ξ) on two-mode vacuum state (initial state), In terms of number states, one can show that the state of two-mode squeezed states is given by: We are now in the position to compute circuit complexity of the reference and target state of two mode squeezed states and compare it to the entanglement entropy. For this, we need to write our reference and target quantum state as gaussian wave functions. The auxiliary position and momentum variables are: The reference state i.e the two mode vacuum state, in the position space can be expressed as a gaussian wave function as follows: The target state, two mode squeezed state, in the position space has the wavefunction: where, A and B are the coefficients and are functions of squeezing parameters r k and φ k : It is helpful to define three other terms for simplifying complexity calculation: Three methods of computing complexity was discussed in [54]. For our case, two methods i.e. computing complexity via covariance matrix method and Nielsen's method are relevant. As we discussed before in Eq. (57) and (58), complexity depends on the choice of cost functions. Let C 1 be the circuit complexity corresponding to linear cost functional F 1 , C 2 to quadratic cost functional F 2 and C k to k family of functionals F k .

Complexity via Covariance matrix method
This method is interesting because complexity from the covariance matrix method is independent of squeezing angle φ k . We will see later that entanglement entropy obtained is also independent of the squeezing angle φ k , so the comparison between circuit complexity and entanglement entropy is more visible in this approach. Since our reference and target states (64) and (65) are in Gaussian form, we can express it as covariance matrix. The covariance matrix for the reference state is given as: The covariance matrix for the target state is given as: where Σ k and Σ − k are defined in (67). The covariance matrix basically carries the same information as the wave function. In the context of covariance matrix approach, circuit complexity quantifies the number of gates to take covariance matrix of reference state to the covariance matrix of target state. We will factorize the covariance matrix G into two 2 × 2 matrices. This gives us the benefit that we can compute complexity for each block and sum over all Ω k to give total complexity. The two symmetric blocks for the reference states are: While, the two symmetric blocks for the target states are: The basis for each block is changed to make the calculation easier as follows: where, S is a specifically choosen matrix such that G s=0 = 1. In our case, matrix S is given by: This impliesG s=0 = 1 and We can assume that k is real. In the language of covariance matrix, the unitary transformation of wave functions can be expressed as: The unitary transformations is then parametrized with gates satisfying SL(2, R) algebra as: cos(µ(τ ))cosh(ρ(τ )) − sin(θ(τ ))sinh(ρ(τ )) − sin(µ(τ ))cosh(ρ(τ )) + cos(θ(τ ))sinh(ρ(τ )) sin(µ(τ ))cosh(ρ(τ )) + cos(θ(τ ))sinh(ρ(τ )) cos(µ(τ ))cosh(ρ(τ )) + sin(θ(τ ))sinh(ρ(τ ))   (77) where, µ, ρ, θ are the coordinates on the SL(2, R) group. Now, we will set following boundary conditions: This boundary conditions applied with the parametrized unitary transformations gives: , (79) In order to make the calculation simpler, we choose: Given this conditions the metric forŨ becomes: The simple geodesic is a straight line on this geometry which is given by ρ(τ ) = ρ(1)τ . From, the boundary conditions (79), we get: In order to get total circuit complexity, we have to sum over both values of k i.e. k and −k. Therefore, for linear and quadratic cost functions C 1 and C 2 , we will get: Using the explicit form of Σ k , Σ − k , the circuit complexity reduces to simple form which is independent of squeezing angle φ k : These two cost functions are then related by: Naturally, for small squeezing parameter r k → 0, C 1 ≈ 0 and C 2 ≈ 0. This makes sense as for small squeezing parameter r k , the reference and target states are basically same.

Complexity via Nielsen's wave-function method
Unlike covariance matrix method, Nielsen's approach using wavefunctions gives circuit complexity of two mode squeezed states that is sensitive to both squeezing parameters: r k and φ k . The general philosophy of computing circuit complexity is basically same as in covariance approach. However, instead of representing the wave function as a covariance matrix, we will directly compute the complexity using reference and target two-mode squeezed states, i.e. eq: (64) and eq: (65) respectively. Then, we will be able to write the circuit complexity in terms of squeezing parameters r k and φ k .
The exponent of the target state i.e. two-mode squeezed states eq: (65) can be diagonalized as: where, N is the normalization constant i.e. denominator in 65 and,M The unsqueezed state, reference state is also a guassian wave function represented by: Our two gaussian wave functions has a form: where, v = (q k , q −k ) and A τ is an 2 × 2 diagonal matrix. For the target state eq: (86), while for our reference state eq: (88), matrix A is A τ =0 . So, The unitary transformation eq: (52) acts like, The boundary conditions is given by: U can be parametrized as in eq: (51) such that at τ = 1, the required target state is achieved. Since, A τ =1 and A τ =0 can have complex elements, elementary gates are restricted to GL(2, C) unitaries. Tangent vector components Y I in eq: (53) are complex parameters while O I are the generators. Eq: (53) can also be expressed as: where, we have note that: and I, J = 0, 1, 2, 3. The metric is then given by: (96) For simplicity, we will choose penalty factors G IJ = δ IJ where we fix it to unity. The off-diagonal elements in GL(2, C) can be set to zero as they increase the distance between states. The U (τ ) will become: where, α i (τ ) are complex parameters and O diagonal i are generators with identity at i diagonal elements. The metric takes a simple form: where Re and Im indicates real and imaginary part of α k respectively. The geodesic is again a straight line in the manifold given by: for each (i ∈ k, −k) and (p = Re and Im). Given the boundary conditions, we will get, for each (i ∈ k, −k). Now, the circuit complexity for linear C 1 (Ω k ) and quadratic cost C 2 (Ω k ) functions can be derived as follows: C 1 (Ω k ) = α k,Re (τ = 1) + α −k,Re (τ = 1) + α k,Im (τ = 1) + α −k,Im (τ = 1) C 2 (Ω k ) = (α k,Re (τ = 1)) 2 + (α −k,Re (τ = 1)) 2 + (α k,Im (τ = 1)) 2 + (α −k,Im (τ = 1)) 2 Using explicit values of Σ k , Σ − k , ω k and ω − k from Eq.

Entanglement entropy of two mode squeezed states
In this section, we will compute the entanglement entropy for the two-mode squeezed states and compare it to the Circuit complexity. Not only our states are entangled, there is also a strong correlation between the two modes. |ψ sq k, −k is also an eigenstate of the number difference number operatorn k −n −k with eigenvalue 0, wheren k =ĉ † kĉ −k andn −k =ĉ † −kĉ k . Due to this strong correlation and symmetry between two modes, average photon number in each mode is same: The reduced density operators for the individual modes are given by: The probability of having n photons in a single mode k or −k is: Commonly used entanglement entropies are von-Neumann and Rényi Entanglment entropies. For a density operatorρ, von-Neumann entropy is given by: If the density operatorρ is pure, then S(ρ pure ) = 0, while for mixed states S(ρ mixed ) > 0. It is usually not trivial to calculate the entropy. However, for the basis in which density operator is diagonal such as in Schmidt basis, entropy can be calculated simply from the diagonal elements as: Since our two mode squeezed state 61 is already in the form of Schmidt decomposition, and we also have the form of reduced density operators of individual modes a and b, we can calculate the Von-Neumann entanglement entropy by realizing that the diagonal elements ρ kk is P (i) n . Then, the Von-Neumann entropy, a measure of degree of entanglement is: tanh 2n r k cosh 2 r k ln(tanh 2n r k ) − ln(cosh 2 r k ) = ln(cosh 2 r k )cosh 2 r k − ln(sinh 2 r k )sinh 2 r k (118) We have plotted the Von-Neumann entanglement entropy in fig: 10. It can be seen that entanglement entropy increases with increasing squeezing parameter r k . Note that we didn't calculate the entropy corresponding to the squeezed state Eq. 61 because naturally this entropy is going to be zero as it is a pure state. Instead, we have calculated entropy for the reduced density matrix.
We can now generalize Von-Neumann entropy to get Rényi entropy for the reduced density operator: where µ ≥ 0 is the Rényi Parameter and d is the Schmidt rank of the squeezed state Eq. 61 which is infinity. Again, we can see that Rényi entropy increases with increasing squeezing parameter r k . In fig: 11, we have plotted Rényi-entanglement entropy for various Rényi parameters µ. For very large squeezing parameter, we get: If we take the limit µ → 1, we get the Von-Neumann entropy 118. Meanwhile, Rényi-2 entropy is given by S 2 (r k ) = ln cosh2r k .
One can also calculate effective temperature of the source by computing the thermal distribution with an average photon number n i = sinh 2 r k . The average photon number of the thermal field is given by: Then, one can compute the effective temperature as: where, ω i = i/c is the frequence of the mode and i ∈ (k, −k).

Quantum Circuit Complexity Vs Entanglement
Now, the comment on comparison of entanglement entropy with circuit complexity is in order. It was very recently shown in [39] that there exists some relationship between entangling power and circuit complexity. Most importantly, if the entanglement entropy grows linearly with time, the geometric circuit complexity also grows linearly.
Generally, quantum circuit complexity and entanglement are different quantities. However, for small values of circuit cost and entanglement, one can use the entanglement entropy to bound the circuit complexity. The argument presented in [39] is that quantum gates that are close to the identity performs little entanglement from product or entangled states. One of the interesting corollary presented in [39] is that whenever entanglement entropy grows linearly in time, the circuit complexity also  grows linearly. Linear growth of entanglement entropy is a generic feature of several quenched many-body systems. Our analysis of complexity and entanglement entropy of two-mode squeezed states is in agreement with the result in [39]. Since both entanglement entropy and circuit complexity computed with Covariance matrix method are independent of squeezing angle φ k , the comparison is clear than with Nielsen's method of wavefunctions. The explicit form of circuit complexity with covariance matrix method is obtained in eq 84: While comparing the form of entanglement entropy eq. 118 and 119 with this circuit complexity, we get: In fig: 12, we have plotted the comparison between Von-Neumann entropy and Circuit complexity (computed using covariance matrix method). The circuit complexity C 1 and C 2 grows linearly just like entanglement entropy. Up to these distances, circuit complexity is indeed lower bounded by entanglement entropy. This result can have physical interpretation. Since entanglement entropy for the two mode squeezed states increases with increasing r k , entanglement entropy from the vacuum to distant states is large. Therefore, in the context of two-mode squeezed states, with proper circuit complexity cost, entanglement entropy could be used as a measure of complexity. So far, we have only compared circuit complexity obtained via Covariance approach. A more detailed numerical comparison of circuit complexity via Nielsen's approach with entanglement entropy and temperature will be discussed in the numerical analysis section.

Numerical Results
In this section, we do the numerical analysis of the circuit complexity calculated for the model of "Black Hole gas". To provide a wholesome and physically relevant discussion, we do the analysis in terms of the scale factor. We begin by solving the evolution equations of the squeezed state parameters given in eqn: 47.   To recast the above differential equations and study the time evolution in terms of scale factor, a simple change of variable is implemented, which transforms the above equation. This change of variable is sometimes called as field redefinition.
In fig: 13, we have plotted the evolution of the squeezed state parameter with respect to the scale factor. The behaviour of the squeezed state parameter r k is crucial   for understanding the behaviour of the circuit complexity and its evolution with the scale factor. From the behaviour of the squeezed state parameters, we see a widely different behaviour of the model in (1+1) dimensions i.e d=1 in the plots. The behaviour for the higher dimensions however looks to be pretty similar. The squeezing is large and growing at early times, however after a certain scale, the squeezing freezes and saturates at a con-stant value of squeezing. The increase in the squeezing grows up to a very large scale for spatial dimension 1 (d=1 in the plots) and the freezing of the squeezing effect is not observed even for extremely high scales. This makes the spatial dimension 1 markedly different from the higher spatial dimensions where the freezing effect in the squeezing is explicitly observed.
In fig: 14, we have plotted the squeezing angle φ k with    respect to the scale factor. For the model considered in this paper, we observe that for the spatial dimension 1, the squeezed angle rises for initial scales and is frozen and saturated at intermediate and late scales. However for higher dimension, the squeezed angle increases at the initial scales but shows a fall after a certain characteristic scale.
In fig: 15 and fig: 16, we have plotted the circuit com-plexity with respect to the scale factor calculated from the two different cost functionals using both Nielsen's and covariance approach. Let us make a comparative analysis of complexity obtained from Nielsen's and covariance approach. The structure of circuit complexity in the covariance approach has the similar pattern as the squeezing parameter r k in fig: 13 and has almost no feature coming out of the squeezing angle φ k of fig: 14.

T e m p e r a t u r e v s d C i / d a c u r v e T e m p e r a t u r e v s e n t a n g l e m e n t e n t r o p y c u r v e T e m p e r a t u r e d C i / d a a n d e n t r o p y V a r i a t i o n o f t e m p e r a t u r e w i t h e n t r o p y a n d d C i / d a a n d e n t r o p y f o r d = 1
(a) Using Nielsen's method.

T e m p e r a t u r e v s d C i / d a c u r v e T e m p e r a t u r e v s e n t a n g l e m e n t e n t r o p y c u r v e T e m p e r a t u r e d C i / d a a n d e n t r o p y V a r i a t i o n o f t e m p e r a t u r e w i t h e n t r o p y a n d d C i / d a a n d e n t r o p y f o r d = 1
(b) Using Covariance matrix method.

C i / d a a n d e n t r o p y V a r i a t i o n o f t e m p e r a t u r e w i t h e n t r o p y a n d d C i / d a a n d e n t r o p y f o r d = 2
(b) Using Covariance matrix method. This makes sense as the circuit complexity obtained from the covariance approach is independent of the squeezing angle. Irrespective of the spatial dimension, the circuit complexity C 1 amd C 2 gradually increases and saturates after some values of a.
In contrast to the covariance approach, Nielsen's ap-proach gives a different story of circuit complexity. This is mainly due to the reasoning that the circuit complexity in Nielsen's approach is dependent on both squeezing parameters: r k and φ k from fig: 13 and 14. This gives one to look at the detail of evolution of wave-function uniquely.
As already pointed out in the previous discussion, the

T e m p e r a t u r e v s d C i / d a c u r v e T e m p e r a t u r e v s e n t a n g l e m e n t e n t r o p y c u r v e T e m p e r a t u r e d C i / d a a n d e n t r o p y V a r i a t i o n o f t e m p e r a t u r e w i t h e n t r o p y a n d d C i / d a a n d e n t r o p y f o r d = 3
(a) Using Nielsen's method.

. 4 T e m p e r a t u r e v s d C i / d a c u r v e T e m p e r a t u r e v s e n t a n g l e m e n t e n t r o p y c u r v e T e m p e r a t u r e d C i / d a a n d e n t r o p y V a r i a t i o n o f t e m p e r a t u r e w i t h e n t r o p y a n d d C i / d a a n d e n t r o p y f o r d = 3
(b) Using Covariance matrix method. speciality of the spatial dimension 1 can be clearly understood from the complexity plots as well. The initial rise in the complexity measures is observed irrespective of the spatial dimension though the scale factor upto which the rise is observed is influenced by the spatial dimension. With the increase in spatial dimension, the rise in the complexities is observed till the lower scale factors. After a critical value of the scale factor the complexity measures show a gradual fall in the values. This rate of fall is found to be extremely less for the spatial dimension 1 where even at large value of the scale factor, only a small fall in the value of the complexity is observed. It can be noted that for higher dimension, complexity measure falls off quickly. The faster a complexity measure falls to a certain minima, it starts to oscillate soon after as seen in the graph. For d = 3 the oscillation starts early compared to d = 2. Also, such oscillatory behaviour is saturated at higher value of scale factor. The oscillations in fig 15 and 15 for higher spatial dimensions could be a hint of the quantum gravity corrections in the very early universe in terms of vaccum fluctuations of "virtual black holes" of radii R. Such fluctuations could effectively resolve the cosmological constant puzzle. The 'vecro component' which describes the part of wavefunctional associated to virtual black hole fluctuations could alter the overall vaccum energy giving us an effective value of cosmological constant Λ = (GR 2 ) −1 and hence resolving the issue. The other way to look at the oscillation of complexity is that at the minimum complexity regions, the distance in initial and evolved states are low as one state can be evolved to next with low number of quantum gates while it is opposite in the maximum complexity regions. So, the structure of the wavefunction in lower complexity regions are more close to the intial states than the one in high complexity region.
In fig: 17, fig: 18 and fig: 19, we have plotted the behavior of entanglement entropy with respect to the scale factor. The two curves in the plots correspond to the two types of entanglement entropy we have considered in this paper viz. von-Neumann entanglement entropy and Rényi entropy. Even though the overall behavior of both forms of entanglement entropy are identical, we still a minute difference. It can be seen that the von-Neumann entropy rises faster to a higher value compared to Rényi entropy. This feature is observed for all spatial dimensions. For spatial dimension d = 1, we observe an increasing behavior of the entropy through the entire range of the scale factor. But for the spatial dimension d = 2, 3, we observe an initial increase in the entropy which then starts to oscillate with its amplitude decaying for higher value of scale factor. It can be noted that with rise in the number of spatial dimension, the rise in entropy decreases and hence saturates to a lower value. We would like to relate the entropy calculated from the squeezed state formalism with the entropy of the black hole gas. One can comment about the entropy of the black hole gas from the entropy calculated using the squeezed state formalism because the information about the black hole gas is itself encoded in the squeezed parameter r k . To be more precise, the evolution equations for the squeezed state parameters written in eqn (47) has been solved using the solution of the scale factor of the black hole gas model as the dynamical variable. Hence, the information about the black hole gas model propagates through the squeezed state parameters to any quantity we calculate. Thus, the entnaglement entropy calculated from the squeezed state parameter is intimately related with the entropy of the black hole gas model.
In fig: 20, we have plotted dC i /da computed with both Nielsen's and Covariance approach w.r.t the von Neumann entanglement entropy to inspect the validity of the conjectured relation proposed by Susskind between complexity and entanglement entropy. We observe that for the spatial dimension d=1, in the initial values of entanglement entropy, the behavior dC i /da shows an increasing behaviour. However, at the intermediate scales, dC i /da shows a sharp fall followed by a rise and saturation at large values of entanglement entropy. Thus we observe a non-linear relation between dCi da and entropy. For low values of entanglement entropy, the difference in amplitude of dC 1 /da and dC 2 /da is higher in Nielsen's approach than in the covariance-approach. This can be because in Covariance approach C 1 and C 2 are related by C 1 = √ 2C 2 = 4 √ 2r k while in Nielsen's approach C 1 and C 2 has complicated relation.
In fig: 21 and fig: 22, we study the behaviour of dC i /da with Von Neumann entanglement entropy for the spatial dimension d=2 and d=3. We observe an almost identical behavior for the higher spatial dimensions with the behavior shown in spatial dimension 1.
In fig: 23, fig: 24, fig: 25, we have plotted dC i /da vs Rényi entropy (µ = 2) for different spatial dimensions. It is observed that the overall behaviour of dC i /da with respect to the Rényi entropy is identical to what we observe in the Von-Neumann entanglement entropy case. This identical nature in the behaviour of dCi da is observed for all spatial dimensions.
In fig: 26, 27, 28 we have plotted the behaviour of the equilibrium temperature of the black hole gas with respect to dC i /da and entanglement entropy for the spatial dimension d=1,2,3. The reason we have plotted the behaviour of the equilibrium temperature with respect to dC i /da and the entanglement entropy on the same plot was to get an idea of how it behaves with two most important quantity in our analysis i.e dC i /da and entanglement entropy. The motivation came from Susskind's conjectured relation where he connected the rate of change of complexity with the entanglement entropy and the equilibrium temperature. However, instead of using dC i /dt, we have used dC i /da as we have used the scale factor, which is cosmologically much more relevant quantity, as the dynamical variable of our analysis. The red curve in the plot shows the behaviour of the equilibrium temperature with respect to the entanglement entropy, whereas the black curve shows the behaviour with respect to dC i /da. It is clearly evident that irrespective of the spatial dimension, the behaviour of the equilibrium temperature shows an increasing behaviour. One can approximate the behaviour as follows: T ∝ S 4 . However, it can be seen that the behaviour of the equilibrium temperature is overall not identical in nature with dC i /da, for different spatial dimension and the measure to compute complexity, though some of the features do match. In Nielsen's approach, it can be observed that for negative values of dC i /da, for two values of dC i /da the black hole gas model attains same value of the equilibrium temperature. For spatial dimension d=1, in the intermediate and positive values of dC i /da, the equilibrium temperature is almost constant, but for higher spatial dimension the multivalue nature of dC i /da with the equilibrium temperature returns. However, in the covariance approach in all three spatial dimensions, the behavior is identical.
Thus we see that irrespective of the spatial dimension and the approach of computing complexity in which the black hole model is considered, neither dC i /da nor entanglement entropy has a linear relationship with the equilibrium temperature.

Quantum extremal islands vs Black hole gas
In this portion, we are going to give a comparative analysis of the quantum extremal islands with the black hole gas model from the perspective of circuit complexity.
• Circuit complexity calculated from the solution of Cosmological islands resembled the page curve in a specific parameter space [36] but for the black hole gas model we observe different behavior of the circuit complexity for different spatial dimension.
• In another parameter space the behaviour of the circuit complexity for the island model showed only a rising behavior which is also different from the one we observe for the black hole gas model.
• The entanglement entropy predicted from the circuit complexity in the cosmological island model again resembled page curve in a particular parameter space and showed a decreasing behavior in another parameter space whereas for the black hole gas model, the entropy showed a increasing behavior for the spatial dimension 1 and 2 and an increasing behavior followed by an oscillation for the spatial dimension 3.
• The oscillatory behavior of the circuit complexity at large values of scale factor, which is observed for the higher spatial dimensions for the black hole gas model is absent in the cosmological island model, even when probed to very high scales. The difference in amplitude of dC 1 /da and dC 2 /da is lower.
The difference in amplitude of dC 1 /da and dC 2 /da is higher.
dC i /da, Temperature and entropy(S) in Black hole Gas Model The behavior of Temperature with dC i /da and entropy is different in three different spatial dimensions.
The behavior of temperature with dC i /da and entropy is same in three different spatial dimensions.

Conclusions
Through analysis of the black hole gas model from the perspective of circuit complexity and entanglement entropy using the language of squeezed state formalism we arrive at the following conclusions: • The circuit complexity computed from Nielsen's wave function approach provides a much better understanding than that computed from the co-variance matrix method as it depends on both the squeezing angle and the squeezing parameter and hence can be related to the entanglement entropy.
• The behaviour of the circuit complexity for the spatial dimension d=1 is significantly different from higher spatial dimensions. Whereas complexity saturates or changes significantly slowly at large scale factors for d=1, it falls off rapidly and has an oscillatory behaviour for higher spatial dimensions.
• The behaviour of the entanglement entropy w.r.t the scale factor for different spatial dimension shows different features. For d=1, it is just an increasing function whereas for d=2 and 3, we observe an oscillatory behavior with the frequency of oscillation increasing with the increase in spatial dimension.
• We observe that for no spatial dimensions the quantity dC i /da varies linearly with the Von-Neumann entanglement entropy or Rényi entropy.
• For different spatial dimensions, the behaviour of the equilibrium temperature with dC i /da is peculiar and it is not possible to predict an approximate relation and one has to study different ranges of dC i /da separately to understand the behavior of the equilibrium temperature.
• For different spatial dimensions, from the behavior of equilibrium temperature with entanglement entropy, it can be understood that the relation between entanglement entropy and equilibrium temperature is not linear but goes as: • From the comparative analysis of the black hole gas model with that of the cosmological islands from the perspective of circuit complexity, we can conclude that Circuit complexity can be used as a useful tool to discover the underlying features of a model which are otherwise difficult to analyse.
The future prospects of the work can be written: • Circuit complexity has been studied for thermofield double states [20]. The process of thermalization can also be realized by a process known as quantum quench, where the states are expressed as the generalised Calabrese Cardy form. Hence one can explore the thermalization phenomenon using circuit complexity.
• People have studied circuit complexity as a deformation in the euclidean path integral for CFT's. This is mainly known as path integral optimization [19]. However, these deformations appear in the context of cosmological perturbation theory as well and one can try to extend this circuit complexity using path integral optimization in de Sitter space.