Quantum Circuit Complexity of Primordial Perturbations

We study the quantum circuit complexity of cosmological perturbations in different models of the early universe. A natural measure for the complexity of cosmological perturbations is based on the symplectic group, allowing us to identify complexity with geodesics in the hyperbolic plane. We investigate the complexity of both the mode functions and the physical perturbations, arguing that the latter often provides a more insightful description of the physics involved. In all models the total complexity reached is rather large. Inflationary perturbations may be represented by a comparatively simple quantum circuit, while the perturbations during a matter-dominated contracting phase present the most rapid growth in complexity. Ekpyrotic perturbations reside in the middle and are distinguished by the smallest growth of complexity before horizon exit. Our analysis serves to highlight how different cosmological models achieve the same end result for the perturbations via different routes and how all models show a pronounced sensitivity to initial conditions.


Introduction
The oldest optical record that we have of the universe is the cosmic microwave background radiation (CMB), which shows us the state of the universe about 380,000 years after the start of the hot big bang expansion. The CMB contains temperature fluctuations that, using known plasma physics, can be extrapolated back in time to primordial density fluctuations with an almost scale-invariant spectrum over the range of wavelengths that can be observed. In these calculations, it is an excellent approximation to treat the primordial perturbations as classical density perturbations.
A striking idea is that the ultimate origin of the primordial perturbations lies in quantum fluctuations that were amplified and rendered effectively classical in the early universe. The best known scenario of this type is inflation [1][2][3][4], which can simultaneously amplify quantum perturbations, render them classical and explain their seemingly acausal correlations by tracing them to earlier causal processes. It came as a surprise to many that there exist alternative scenarios, in particular ekpyrotic cosmology [5,6] and a contracting matter phase [7], that can achieve the same goal with entirely different physics. In all cases, however, the underlying idea is that quantum fluctuations are turned into effectively classical density perturbations. Thus in all these models the universe acts as a quantum computer, processing an initial quantum state (usually taken to be the vacuum state) into a state that can explain what we observe in the CMB. In this paper we want to explore this alternative description in terms of quantum computation, in particular by calculating the complexity of the involved computation. As the name suggests, the complexity of a quantum computation may be thought of as the difficulty in building a quantum computer performing the same task. More precisely, the complexity is taken to be the minimum number of quantum gates (from a specified set) that the task requires. In other words, we are asking how complicated a quantum computer would have to be in order to simulate the perturbations in various early universe models.
As we will see, for the states relevant to early universe cosmology, we will quite naturally be led to depict their evolution in hyperbolic geometry, see Fig. 1. The cosmological phases mentioned above turn quantum fluctuations into effectively classical fluctuations by combining two effects: they amplify the perturbations and also turn the quantum state into a highly squeezed state where the uncertainty in momentum is vastly smaller than that in amplitude. In this ampli cation squeezing Figure 1. We may usefully display the evolution of the quantum state of cosmological perturbations on a Poincaré disk, which forms a representation of hyperbolic geometry. Starting from the Bunch-Davies vacuum state in the centre of the disk, the arrows indicate the directions of evolution in which the perturbations are amplified or squeezed. A final state which is both amplified and squeezed, as required to match the observations of the CMB, corresponds to a state very close to the boundary directly right from the centre, as indicated by the star. manner the quantum states become equivalent to a statistical mixture of classical perturbations [8][9][10][11][12][13][14][15][16]. On the Poincaré disk in Fig. 1 we have indicated by arrows the directions of evolution corresponding to amplification and squeezing, starting from the vacuum state in the centre. Even though all viable early universe models must end up with the same final state, they employ rather different routes to get there. Fig. 2 illustrates this by showing the evolution in the different early universe models that we investigate in this paper (the actual numerical calculations performed to produce this figure will be explained in due course). As is immediately apparent from the figure, in some models amplification and squeezing occur separately (in particular in inflation), while in other models they occur simultaneously (e.g. for isocurvature perturbations during ekpyrosis). Thus it may already be guessed that the complexities of the corresponding computations will also turn out to be different, and this is indeed what we find.
In fact we find that the complexity depends quite strongly on the cosmological model. In all cases the final complexity is very high, due to the vast range of scales that need to be processed. We actually find that it is useful to characterise the evolution of complexity in terms of the number of e-folds, rather than in terms of physical time. Broadly speaking, inflation turns out to be a simpler quantum computer than contracting cosmologies, with a contracting matter phase being the most complex. Ekpyrotic perturbations behave in an intermediate regime, though they have the feature that complexity grows the slowest while perturbations are still on subhorizon scales. An interesting aspect is that the complexity is sensitive to the total duration, and thus also to the beginning, of the cosmological phases in question. This offers the hope that complexity may help in further elucidating the conditions required at the beginning of the early universe models considered, in order to see how they may eventually form a part of a complete cosmology.
We will start the paper with two review sections, one on circuit complexity (section 2) and one on the quantisation of cosmological perturbations (section 3). Readers may skip them if they feel comfortable with the subjects, though we review them with a particular application in mind, and hence the approach may be interesting even for the expert reader. We will then apply these methods to investigate circuit complexity for different cosmological models. Inflation is considered in section 4, while we divide the analysis for ekpyrosis up into single-field models (section 5) and two-field models (section 6), given that the involved physics is significantly different. We then compare these results with matter contraction and adiabatic ekpyrosis in section 7, where we will see that these models provide the extremes on the spectrum of complexity.  Fig. 1. The colour code depicts the time evolution as a function of e-folding number N for a given mode from its subhorizon, ultraviolet regime (in dark blue), through horizon exit at N = N exit (in chartreuse colour), all the way to its super-horizon, infrared regime (in red). The differences between inflation (where amplification and squeezing occur in succession) with for instance isocurvature perturbations during ekpyrosis (where amplification and squeezing occur simultaneously) becomes manifest. Not all perturbations achieve a highly amplified and squeezed end state: in particular curvature perturbations during ekpyrosis do not, as the evolution does not get close to the boundary. Meanwhile, a contracting matter phase reaches the boundary extremely fast, with much of the distance covered already sub-horizon.
Our results are discussed in section 8, and we include a brief comparison with earlier proposals in appendix A. We use natural units with = 1 and 8πG N = 1 throughout.

Brief review of quantum circuit complexity
In quantum computation, complexity describes how difficult it is to build a circuit that transforms a given reference state |Ψ R into a target state |Ψ T . A conceptually straightforward measure for the difficulty of performing a computational task is simply to look at how many quantum gates one needs in order to perform the required transformation. Complexity thus provides a quantitative way of evaluating how much a wave function has changed. The fact that we will only be interested in Gaussian perturbations greatly simplifies the analysis (a useful exposition, which we will partly follow here, can be found in [17]). Let us assume in the present section that we are working with a one-dimensional harmonic oscillator in position space x. We will further assume that the reference and target states are given by the wave functions with frequencies ω and Ω, which at first we take to be positive real numbers. The evolution between the reference and target states will be unitary, The question is then how many gates one needs in order to implement the unitary operator U or, rather, what the minimum number of such gates might be. We can proceed by first discretising the evolution by considering small steps of size . Our gates will be elementary unitary operators. Useful examples are given by the following operators, displayed here along with the effect that they have on a wave function |Ψ(x) : Here p x denotes the momentum operator conjugate to x, that is to say p x = −i∂ x . In the examples above, H effects a phase change in the wave function, while J shifts the position. Most useful to us is Q, which leads to a scaling in the position (and also includes a normalisation factor). These are just a few examples -one could in principle consider many more operators. A general circuit is built from a number of such gates performed in succession, i.e., where the αs are natural numbers. General circuits would contain operations that are later undone again by other gates. We are interested in the minimal number of gates required, however. For our specific example in Eq. (2.1) it is clear that we will only need the scaling (Q) gate, so that for us a useful circuit will be of the form Thus we can obtain the required transformation from |Ψ R to |Ψ T as long as 2α = ln(Ω/ω).
In general one should now assign a measure to the 'depth' of the circuit. A popular method is to assign a metric to the space of unitary operations, after taking the continuum limit → 0. This is useful, as the shortest circuit will then correspond to a geodesic in this space, and this geometric method allows one to use the power of differential geometry (see [18]). We will turn to such an example shortly. Our present example is however so simple that it is sufficient to directly equate the complexity C to the number of gates, where we have to re-scale the measure by in order to obtain a well-defined limit as → 0: Note that the complexity evolves logarithmically, hence even small numerical changes in the complexity correspond to significant evolution of the wave function. As we will see, the early universe provides a laboratory that is surprisingly efficient at producing complex quantum circuits. In cosmology, the frequencies ω, Ω are in general complex valued since we are dealing with generic Gaussian states. We therefore have to find an appropriate generalisation of Eq. (2.6). Perhaps the most straightforward generalisation is simply to let α become a complex number too, which one should then think of as encoding the effect of two quantum gates, one for the real and one for the imaginary part of α. This has been formalised in [19] (and used in the cosmological context, see e.g. [20,21]), with the complexity now defined via (2.7) We will refer to this definition as the analytically continued (a.c.) complexity. There exists another generalisation, however, which leads to an appealing geometric interpretation, both of the states and of their complexity. This generalisation was developed in [22,23] and takes as its starting point the covariance matrices associated to Gaussian states. Let us first rewrite the target state as with a, b being real valued. Moreover, we group the position and momentum into a combined coordinate ξ m = (x, p) with m = 1, 2. Since we are working with a Gaussian state, all information is contained in the quadratic combinations where the antisymmetric matrix encodes the canonical commutation relations, while the symmetric covariance matrix G T has entries given by the expectation values (2.11) The evolution equation (2.2) now becomes The transformationsξ m =M m n ξ n that preserve the canonical commutation relations are elements of the symplectic group,M ∈ Sp(2, R). This suggests the use of gates that belong to the algebra sp (2). In fact, it is precisely quadratic combinations of x and p that form the corresponding generators, As noticed in [22], it is sufficient to consider the sub-algebra formed by V and W -moreover, this restriction will lead to a useful connection with hyperbolic geometry, as we will now review. A matrix representation of V, W , along with the associated group elements/gates is provided by These gates are sufficient to take a reference matrix G R = diag(1/ω, ω) to the target (2.11). We may slightly simplify the task by using a squeezing operation S = diag( √ ω, 1/ √ ω) on the reference and target states, with the effect that Thus the reference has become the identity.
A general circuit will again consist of a sequence of gates. It is here that we will make the transition from a discrete to a continuous description, which will allow us to make a connection with differential geometry. In the continuous description a circuit is represented by a path-ordered exponential 16) with M I = (V 2×2 , W 2×2 ) and Y I (s) describing switches that turn the respective gates on or off. The full circuit then runs from the identity at s = 0 to the targetG T at s = 1.
It is straightforward to verify that a general element of such a circuit can be parameterised by Note that the transformation law (2.12) involves the transpose, thus also generating the required first row/second column entry. The connection to geometry is found by inverting (2.16) according to where the factor 1/2 comes from the normalisation of the gates M I . From this we can immediately obtain a circuit geometry specified by the line element Here, following [22], we chose g IJ = diag(1, 1/2), but other choices are equally simple to implement. The line element above may be recognised as the metric on the hyperbolic upper half plane H 2 . Optimal circuits then correspond to geodesics, which on the hyperbolic plane are given by arcs of circles that are perpendicular to the boundary z = 0. Our reference and target states correspond to the coordinate locations (2.20) Then the complexity may be defined as the hyperbolic distance between these points, We will refer to this definition simply as the "complexity", though for the sake of comparison with the analytically continued complexity (2.7) (such as in appendix A for instance) we shall sometimes specifically refer to it as the "hyperbolic complexity". Finally, we mention that it can be useful to map the hyperbolic plane to a finite representation, in particular to the Poincaré disk. This mapping is most easily expressed in terms of the complex coordinate Z ≡ y + iz. Then a point Z on the hyperbolic half-plane is mapped to the point Z−i Z+i on the disk. In particular, our reference state (0, 1) corresponds to Z = i and gets mapped to the origin of the disk. Target states that are far from the reference state, and thus obtain a large hyperbolic complexity, can then be found very close to the edge of the Poincaré disk. This is the representation that is used in Figs. 1 and 2. It thus becomes clear how amplification and squeezing of the perturbations (ω/a and b/a growing, respectively) moves one around the Poincaré disk, getting closer to the edge and the point (1, 0) on the disk as |Z| grows.

Quantisation of cosmological perturbations and complexity thereof
The easiest way to model both inflation and ekpyrosis is by considering the dynamics of scalar fields coupled to gravity and evolving in an appropriate potential. We will start our analysis with models involving only a single scalar field σ with an exponential potential V (σ). The action is given by Note that is always taken to be positive, but V 0 can be positive or negative. In a flat Robertson-Walker spacetime, the equations of motion reduce to When is constant (i.e., when the scalar field equation of state (EoS) is constant), there exists an exact scaling solution to the equations of motion, given by For inflation, one would take V 0 > 0 and assume the universe to be expanding (i.e., t > 0), with EoS < 1. This then corresponds to accelerated expansion of the scale factor and slow evolution of the scalar field. For ekpyrosis, one would take V 0 < 0 and assume the universe to be contracting (i.e., t < 0) with > 3. This would correspond to slow contraction and fast evolution of the scalar field. Note that the EoS (slow-roll/fast-roll) parameter can be re-expressed in various ways, Below, we will work in conformal time τ , defined by dτ ≡ dt/a, in terms of which the scaling solution becomes where a prime denotes a derivative with respect to conformal time. Conformal time naturally runs over negative values for both inflation and ekpyrosis. We will also often use the e-folding number to characterise the time dependence later. We define the number of e-folds of evolution via dN ≡ d ln(a|H|) = d ln |H| since it is the appropriate definition with regard to the flatness problem. For inflation H is roughly constant so that instead one often uses dN = d ln a, but this is not useful for ekpyrosis where a is approximately constant. Thus, we will use dN = −d ln(−τ ) = (1 − )d ln a, which is valid generally. We will be primarily interested in perturbations of these models. For this, we will work in comoving gauge, in which the scalar field fluctuation vanishes, δσ = 0, and the scalar degree of freedom is represented by the comoving curvature perturbation R. The metric on spatial hypersurfaces is then given by h ij = a 2 (1 + 2R)δ ij . We will only consider scalar perturbations -an analogous calculation could be performed for gravitational waves. At quadratic order, the action for the comoving curvature perturbation is remarkably simple and given by The corresponding equation of motion is In order to quantise the perturbations, it is useful to define the Mukhanov-Sasaki variable v ≡ zR , Switching to conformal time, the action (3.7) now becomes canonically normalised: The action is quadratic and will thus lead to a linear equation of motion. This implies that it will be useful to expand the perturbations into Fourier modes: The equation of motion for each Fourier mode of wavenumber k is then The linearity of the equation implies that each mode evolves independently, and there is no mode mixing.
We can quantise the system in the Heisenberg picture by promoting the mode functions to operators and writing these new operators as a linear combination of annihilation and creation operators,v Here the f k (τ ) are time-dependent (complex) solutions of the equations of motion (3.12), which because of the spatial isotropy of the background depend only on the modulus k ≡ |k|. Note that the definition above implies the relationv −k =v † k , which ensures that the comoving curvature perturbation is real valued, as it should. We then require the annihilation/creation operators to satisfy the canonical quantisation condition, This condition implies that the field operatorv and its conjugate momentumπ ≡v − (z /z)v satisfy the canonical equal-time commutation relation, as well as the trivial commutators, The Wronskian is a constant of motion, where we have fixed the right-hand side in such a way as to ensure the canonical normalisation of the mode functions.
The above quantisation procedure in the Heisenberg picture is standard in early universe cosmology. But in order to investigate the analogy with circuit complexity, it is more useful to work directly with the wave function, i.e., in the Schrödinger picture. The wave function will be of Gaussian form since at our level of approximation the perturbations are governed by a quadratic action. It will in fact be a product of Gaussians for each wavenumber k, and thus we may focus on a single wavenumber and write where the proportionality constant is determined upon normalisation. We are interested in the vacuum stateâ|Ψ = 0. If we rewrite the annihilation operator as iâ = (f * − z z f * )v − f * π , then withπ → −i∂ v we can deduce an expression for the correlator A σσ in terms of the mode functions as We also added the correlator of the comoving curvature perturbation, whose expression follows immediately upon using v = zR. The equation of motion for the correlator follows from the Schrödinger equation i|Ψ =Ĥ|Ψ , whereĤ is the Hamiltonian operator, or equivalently from the Heisenberg equation of motion, We can now solve for the mode functions defined above. On a scaling solution, the mode equation (3.12) becomes It can be solved exactly by rewriting it in the form The solution approaching the Minkowski vacuum e −ikτ / √ 2k in the far past is given in terms of a Hankel function of the first kind, The phase is immaterial in what follows, and hence we will drop it. The explicit solution to the mode functions allows us to evaluate the correlator. Using the formula d dx H (1) α (x) we obtain the compact expression Note that the term involving z /z has disappeared from A σσ . This cancellation follows directly from the definition α 2 σ = 1/4 + τ 2 z /z, which implies that z behaves as a power law z ∝ (−τ ) n for some real n, in turn implying that α 2 σ = (τ z /z − 1/2) 2 . Thus for any single field model, where the fluctuations are adiabatic, i.e., those of the field that drives the background, the term inversely proportional to the comoving horizon, z /z ∝ (−τ ) −1 ∝ |H|, disappears from the correlator. This is important for models where the comoving horizon shrinks, such as inflation and ekpyrosis, since a potential source of amplification of perturbations is removed. As we will see, the consequences of this fact for inflation and ekpyrosis differ drastically.
At early times (τ → −∞), the correlator approximates its Minkoswki vacuum value A σσ k. In order to obtain the late-time limit, τ → 0 − , one has to use the asymptotic expansion of the Hankel function, The leading real and imaginary parts of the correlator depend on the value of α σ and are given (in the late-time limit) by Given that z 2 = 2 ā 2 0 |τ | 1−2ασ when 0 ≤ < 1 or ≥ 3, we can also find the corresponding expressions for the correlator for the comoving curvature perturbation R, We will analyse the implications of these expressions for relevant single field models of the early universe in the next sections. Given that during inflation and ekpyrosis the wave function evolves a lot, we can expect the complexity to grow significantly too. In previous works, the complexity 1 was analysed with regard to the canonically normalised variable v [20,21]. In that case it is natural to use the early Minkowski correlator A σσ k as reference state. The target state is taken to be the late-time super-horizon state in which the perturbations find themselves at the end of inflation or at the end of the ekpyrotic scenario, just before reheating occurs in every case. Based on the definition (2.21) -so using the hyperbolic measure rather than the analytically continued one -we may thus define the complexity of the canonical variable v as 2 As we will see below, it may make more sense however to work in terms of the physical variable, which is the comoving curvature perturbation. Then the complexity must be defined with respect to the state at some early time, which one may think of as the start of inflation or ekpyrosis. We will take the initial time to be τ i and the final time at reheating τ f . Implicit in this prescription is the assumption that inflation and ekpyrosis had a start and did not reach back arbitrarily far into the past. This assumption seems well justified in light of recent works analysing the beginning stages of inflation [24][25][26][27][28][29][30][31][32][33][34], and it is well founded for ekpyrosis, in particular in the context of its cyclic realisations [35][36][37][38][39][40][41]. Starting again from the definition in Eq. (2.21), we may thus write the complexity of the curvature perturbation as , (3.29) in the time interval τ i ≤ τ ≤ τ f . This provides the basic ingredients to analyse the evolution of complexity in inflation and ekpyrosis in the following sections. A few observations can immediately be made: as long as one remains deeply sub-horizon (−kτ 1), we notice that A σσ k implies X v 1 and C v 0, so complexity in terms of the canonical variable does not grow. In terms of the curvature perturbation, the corresponding early-time sub-horizon limit is A RR (τ ) 2 ka(τ ) 2 , and so one has where a i ≡ a(τ i ). As the universe expands (a a i ) or contracts (a a i ), either one of the two terms will dominate. Thus, as X R 1, the complexity may be approximated as C R 2 −1/2 ln(2X R ), and so one finds where the + sign holds for expansion (when 0 ≤ < 1) and the − sign for contraction (when > 1). Therefore on sub-horizon scales, as we will confirm in the next sections, inflation with 1 has C R √ 2N and ekpyrosis with 3 has C R √ 2N / .

Inflation
Inflation is characterised by accelerated expansion, hence 0 ≤ < 1 and consequently α σ > 3/2. In the slow-roll limit, where 1, we may approximate α σ ≈ 3/2 + ; we will illustrate our results in this limit. One can see from Eq. (3.26b) that the dispersion, of the canonically normalised modes v is growing as inflation proceeds, i.e., as τ → 0 − . This means that these modes are strongly amplified. Meanwhile, since the ratio is also growing fast, we can see that the wave function evolves into a highly squeezed state. Another perspective is offered by the correlator for the physical perturbation, namely the comoving curvature perturbation, in Eq. (3.27b) -see Fig. 3. From its late-time asymptotic expression, we can see that the real part evolves to a constant. This is a reflection of the fact that on large (super-horizon) scales, the comoving curvature perturbation becomes constant. The spectrum of the perturbations can be read off from the k dependence, which implies a spectral index n s = 4 − 2α σ ≈ 1 − 2 with a small red tilt. Since the state is highly squeezed, | Im(A σσ )/ Re(A σσ )| 1, the fact that the curvature perturbation is approximately conserved on large scales implies that on these scales the momentum is known to be small with high precision and consequently the bulk of the uncertainty resides in the amplitude.
Using the full expression for the correlator [Eq. (3.24)], one can compute the complexity given by either (3.28) or (3.29). A numerical example of the evolution of the two complexities,    C v for the mode functions and C R for the curvature perturbations, is shown in Fig. 4. It is striking how differently they evolve: the most obvious difference is that C v only grows on superhorizon scales, while C R grows throughout. In fact, for C v we have the following leading-order approximate expressions for modes that have not yet exited the horizon and for those that already have (using Eq. (3.26b) and Taylor expanding Eq. (3.28) as τ → 0 − ), where τ indicates the conformal time at horizon exit, while the complexity of the curvature perturbation is well approximated by The above expression on sub-horizon scales was already derived in Eq. (3.31), while for superhorizon scales we used Eq. (3.27b) and expanded Eq. (3.29) as τ → 0 − . These expressions can be recast as the growth of complexity (∆C) in the different regimes as a function of the duration of that regime (∆N ) as These approximate expressions are very accurate, as may be seen by the closeness of the numerical curve to the dashed lines in Fig. 4. The 'canonical' complexity grows roughly as √ 2 times the number of e-folds that a given mode spends outside of the horizon. For modes that exit the horizon say 50 e-folds before the end of inflation, the complexity will thus end up being approximately 50 √ 2. Meanwhile, modes that leave just before the end of inflation would have essentially vanishing complexity by the time of reheating. By contrast, the 'curvature complexity' already grows (again roughly as √ 2∆N ) while a mode is still sub-horizon. This is because the curvature perturbation is significantly amplified as it is drawn out of the vacuum. Once the perturbation exits the horizon, the complexity keeps growing (roughly at the same rate) due to the imaginary part of the correlator, which is still growing. Thus, even though the real part of the correlator becomes constant, signalling that the curvature perturbation is conserved outside of the horizon, the wave function has a rapidly growing imaginary part. In other words, the wave function for the curvature perturbation is increasingly of Wentzel-Kramers-Brillouin (WKB) form, with a rapidly changing phase and a slowly changing amplitude. The phase drives the complexity to even higher values. Overall, the curvature complexity is thus much higher than the canonical complexity. If an inflationary phase lasts say 80 e-folds then the same mode considered above would undergo a growth of complexity of about 30 √ 2 ≈ 42 before horizon exit, followed by an additional 50 √ 2 ≈ 71 units of growth after horizon exit. More generally, if the inflationary phase lasts N tot e-folds in total, then the complexity of a mode exiting N e-folds before the end is given by The canonical complexity is thus ignorant about the beginning and duration of inflation. By contrast, the curvature complexity is sensitive to the total duration of inflation and hence also to the pre-horizon-exit evolution. The above result for the curvature complexity in inflation admits an additional interesting interpretation. Approximating inflation by de Sitter spacetime and thinking of de Sitter as a thermal system with temperature given by the inverse of the horizon radius (i.e., T ∼ H), one can ask the question of how fast the system can 'scramble' [42,43] perturbations over the Hubble horizon. From the point of view of a dynamical system, trajectories moving apart exponentially fast as a function of time (∼ e λ L t ) can be diagnosed as chaotic, and the strength of the sensitivity to the initial conditions is characterised by the Lyapunov exponent λ L . The time at which trajectories have moved apart by an O(1) factor is then representative of the scrambling time (t * ), hence t * ∼ 1/λ L . For a thermal quantum system, it has been conjectured that there exists an upper bound on the growth of chaos given by the temperature, namely λ L T [44], and black holes as well as de Sitter are potentially among the fastest scrambling systems in the universe, given that they saturate the conjectured bound [42,43]. For generic quantum systems, there exist various approaches and techniques to quantify the chaoticity and correspondingly compute the Lyapunov exponent and scrambling time (for de Sitter we certainly expect λ L ∼ T ∼ H, see e.g. [45][46][47]). Complexity might be a promising quantity in that respect. Indeed, as the evolution of complexity depends on the logarithm of the correlator, linear growth in complexity is actually indicative of exponential separation of trajectories in 'field space' (the space of quantum states) and thus of chaos [48][49][50][51][52][53]. One could therefore attempt to identify λ L ≡ dC/dt, which would yield a scrambling time of order the Hubble parameter, λ L √ 2H (on super-horizon scales for C v , but on all scales for C R ). There might however be different roles for the canonical and curvature complexities: the canonical complexity may be a good measure of chaos in the sense that it is the effective mass squared of the mode functions that transitions from positive to negative near horizon exit, implying an enhanced sensitivity to initial conditions as the horizon is crossed. On the other hand, the physical complexity is sensitive to the entire inflationary evolution, and thus to the vast separation of scales achieved over the course of an entire inflationary phase. This may also lead to connections with the trans-Planckian censorship conjecture [30][31][32]54]. It will certainly be of interest to explore the connections between complexity and chaos in much more detail in the future.

Single-field ekpyrosis
Ekpyrosis is a phase of slow contraction with ultra-stiff EoS > 3. This EoS is required in order to make the ekpyrotic phase an attractor, such that anisotropies do not grow despite the fact that the universe is contracting. During this phase, the scale factor is almost constant, but the Hubble rate grows quickly in magnitude. Since the scale factor evolves little, the wavelength of perturbations changes equally little. However, the rapidly growing Hubble rate implies that the comoving horizon 1/|aH| ∝ (−τ ) shrinks rapidly as the universe contracts, τ → 0 − . Consequently perturbation modes with ever shorter wavelengths successively exit the horizon, just as in inflation.
A crucial difference with inflation is that the adiabatic modes associated with the ekpyrotic scalar field do not get amplified much and do not develop a nearly scale-invariant spectrum. One may see this by inspection of the correlator for the curvature perturbation. Since the EoS is ultra-stiff, we have that α σ < 1/2 (for large we have α σ ≈ 1/2 − 1/ ) and thus Eq. (3.27a) implies that at late times both real and imaginary parts of the correlator A RR tend to constants. A numerical example is shown in Fig. 5. Thus the ratio of the imaginary part to the real part does not grow, implying that no squeezing occurs and that these perturbations cannot be given a stochastic classical interpretation. Moreover, the spectral index n s = 4 − 2α σ > 3 remains very blue.
Given these observations, we may conclude that the wave function of adiabatic modes does not evolve much during an ekpyrotic phase, and thus we expect the complexity to remain rather small. A numerical example is shown in Fig. 6 Fig. 5. The canonical complexity keeps growing even after horizon exit, while the curvature complexity reaches a constant value. As discussed in the main text, the curvature complexity seems to give a more faithful representation of the salient physical effects.
the left panel) keeps growing even after horizon exit, due to the fact that the correlator for the Mukhanov-Sasaki variable keeps growing. Equation (3.26a) shows that it grows in magnitude even at late times, with real and imaginary parts growing at the same rate as a result of the absence of squeezing (specifically A σσ ∼ (−τ ) −2/ , and so C v ∼ √ 2N / , which is depicted by the dashed orange line in the left panel of the figure). Our discussion in the previous paragraph indicates that the curvature complexity represents a more sensible reflection of the physical processes during the ekpyrotic phase, and thus constitutes the preferred measure of complexity (as was the case for inflation).

Two-field ekpyrotic models
In the previous section we saw that the adiabatic modes are not significantly amplified during an ekpyrotic phase, and moreover they retain a blue spectrum. Thus, put simply, they cannot explain the primordial temperature fluctuations seen in the cosmic background radiation and are utterly unimportant on large scales. However, the ekpyrotic background solution does possess a scaling symmetry, in close analogy to inflationary spacetimes. The problem, most clearly discussed in [55], is that during ekpyrosis the scaling symmetry projects out of the curvature perturbation. A resolution of this issue is to consider two-field models. If there is a second ekpyrotic field, then the difference of the two fields (which inherits the scaling symmetry from both fields) is a gauge-invariant variable and may serve as the seed for primordial (nearly) scaleinvariant perturbations. This has been discussed at length in several papers, e.g., [56][57][58][59][60][61], hence we will only provide a brief review here.
In addition to the adiabatic scalar σ driving the background we will consider the transverse field s, whose fluctuations represent entropy (or isocurvature) perturbations. By definition we will assume that the ekpyrotic background evolution remains unchanged, so that at the background level s = 0. (In the case where one starts with two ekpyrotic fields, this corresponds to rotating the field basis such that σ points along the background trajectory.) The perturbations δs will however be important. They are sensitive to the shape of the potential, which we take to be given by where κ 2 is a positive parameter. In a theory of two scalar fields with two exact exponential potentials, the above potential arises after the above-mentioned field rotation and leads to κ 2 = 1.
Thus κ 2 may be thought of as parameterising the deviation from exact exponential potentials. Note that the potential is unstable in the s direction -this crucial feature will be responsible for the amplification of δs perturbations. It also leads to a sharpening of the issue of initial conditions, with possible implications investigated in [39,40,62]. Non-Gaussian corrections are encoded in higher-order terms in the potential -these lead to interesting observational consequences (see e.g. [63][64][65][66][67][68][69]), but will not be important for our present work. There exists a second class of models, again with two scalar fields, but where instead of having a potential that is unstable in the transverse direction one considers a non-minimal coupling between the two scalar fields. A judicious choice of coupling can allow the adiabatic field to transfer its scaling symmetry to the second scalar field, without introducing instabilities in the background dynamics. These models, and their observational consequences, were explored in [70][71][72]. They have the advantage that they lead to non-Gaussian signatures that are significantly smaller, yet within reach of near-future observations [73,74].
In all of these models, the ekpyrotic phase flattens the universe, suppresses anisotropies and amplifies entropy perturbations. In the adiabatic direction, the potential is approximately a negative exponential. This cannot grow indefinitely towards ever larger magnitudes, and one thus expects the potential to turn off at some point. This will mark the end of the ekpyrotic phase, and with the influence of a potential having disappeared, the evolution afterwards is that of a phase dominated by the kinetic energy of the adiabatic field. We will sometimes refer to this kinetic phase as 'kination'. Kination is usually envisaged to be a rather short phase, immediately followed by a bounce (and reheating) into the standard hot big bang expanding phase of the universe. The details of the bounce and reheating are highly model dependent (see e.g. [75][76][77][78][79]), just as for reheating after inflation. Nevertheless, the long wavelength modes we are interested in remain essentially unchanged during a non-singular bounce (see e.g. [80][81][82]). Thus we will end our analysis with the kinetic phase, the same way we ended the analysis before inflationary reheating in the previous section.
The analysis is greatly simplified by the realisation that the adiabatic/curvature and entropic/isocurvature perturbations remain decoupled during ekpyrosis and most (or all) of kination, even for models where the background scalar fields are coupled non-minimally [72]. The Lagrangian is thus a sum of terms involving solely the adiabatic mode v σ ≡ zR and terms involving only the entropic mode v s ≡ a δs (where from here on we will consider a single Fourier mode and drop the subscript indicating the wavenumber to lighten the notation): The equations of motion are thus and the canonical momenta are given by The fields can be quantised by following the usual procedure of promoting fields to operators, with f σ,s (τ ), g σ,s (τ ) being complex, linearly independent solutions of (6.3). The conserved Wronskian combinations are now slightly more involved and read All of this is simply the two-field generalisation of the discussion in section 3. One can now use the Wronskians to express the annihilation operators in terms of the field operators and momenta, This allows us to transition to the Schrödinger picture, where the vacuum wave function is defined viaâ|Ψ =b|Ψ = 0. Thus, with the commutation relations realised viaπ → −i∂ v , and up to an unimportant normalisation proportionality factor, the wave function is given by with the correlators During the ekpyrotic phase, the scale factor evolves as a function of conformal time τ the same way as in Eq. (3.5) and consequently The Fourier mode functions for the adiabatic (v σ ) and entropic (v s ) fields then have the solution where the integration constants have been fixed such that at early times the mode functions approach their expressions in the Minkowski vacuum (up to an unimportant phase). The order of the Hankel functions is given by There is no mixing between the modes. Hence for the adiabatic modes we will recover exactly what we had before in section 5. Note that for κ 2 ≈ 1 we have that α s ≈ 3/2 and thus we can expect that by contrast the entropic modes will behave much more like inflationary perturbations. The correlators may be written as .
At early times, τ → −∞, the correlators approximate their Minkowski vacuum values A σσ k, A ss k. At late times, we may again use the asymptotic form of the Hankel functions in Eq. (3.25) to obtain the same A σσ as in Eq. (3.26a) and where we have kept the leading real and imaginary parts. Given that the physical perturbation is δs = v s /a, it makes more sense to consider its correlator A δsδs = a 2 A ss , which during the ekpyrotic phase evolves to In the last line we made the approximation α s ≈ 3/2. Several features can immediately be read off: since > 3, the real part of the correlator shrinks, indicating that these modes will be amplified. Meanwhile, the imaginary part grows, the ratio between imaginary and real parts growing ever more rapidly as |τ | −3 . Thus the state becomes highly squeezed and the perturbations become equivalent to a stochastic mixture of classical perturbations. In contrast to inflation, the dispersion of the entropy perturbations (∼ 1/ Re[A δsδs ]) does not reach a constant value, but rather keeps growing as the ekpyrotic phase proceeds. This is because the ekpyrotic potential is steep and keeps evolving to larger magnitudes. In order for the entropy perturbations to be able to act as realistic seeds for primordial density perturbations, they must reach a magnitude of around 10 −4 , implying that the potential must reach the grand unified scale [59]. As discussed in section 5, the potential turns off at this point and a kinetic phase ensues. For the purposes of our present study, it is a good approximation to assume that the potential becomes zero abruptly. We just have to make sure that we match the scale factor and the mode functions at the ekpyrosis-kination transition (call it τ e-k ). This may be done as follows: during the kinetic phase, we may write the scale factor as a kin (τ ) = a kin,0 (τ c −τ ) 1/2 , allowing for the fact that the would-be crunch is shifted by an amount τ c compared to the ekpyrotic phase. Moreover, a general solution of the mode equation (6.3) during the kinetic phase (with V = 0) is given by a linear combination of Hankel functions of the first and second kind as Then we can match the field values and momenta at the matching time τ e-k by imposing a(τ e-k ) = a kin (τ e-k ) , a (τ e-k ) = a kin (τ e-k ) , g(τ e-k ) = g kin (τ e-k ) , g (τ e-k ) = g kin (τ e-k ) , (6.17) and solving for the constants a kin,0 , τ c , c 1 , c 2 . The evolution of the correlator A δsδs is shown in Fig. 7 as a function of conformal time. At the ekpyrotic-kinetic transition, the correlator is continuous, but its derivative is not. This is because the derivative of the correlator involves second derivatives of a and g, and these are sensitive to the abrupt turning off of the potential. Figure 7 shows that during the kinetic phase, the entropy perturbations are further amplified since the real part of A δsδs keeps decreasing. This is however just a logarithmic growth, which does not alter the ekpyrotic amplification significantly. As for the imaginary part, it is similarly reduced during the kinetic phase, implying that the amount of squeezing remains approximately constant. We can verify this using a simple analytic approximation. In the large-scale limit, −kτ 1, the mode function (6.16) during the kinetic phase can be approximated as g kin (τ ) √ −τ [C 1 + C 2 ln(−kτ )], where the complex constants C 1,2 are related to c 1,2 . At fixed k and in the limit τ → 0 − , it is thus dominated by g kin (τ ) ∼ √ −τ ln(−kτ ). Then, the correlator for the entropy perturbation scales as (to leading order for the real and imaginary part) , (6.18) whereC is a real constant related to C 1,2 (and c 1,2 for that matter), and so for the physical perturbations, we have .
The overall logarithmic correction is easily visible in Fig. 7.
We are finally in a position to investigate the circuit complexity of the entropy perturbations. Here we also have the choice of looking at the complexity of the mode functions v s or that of the physical perturbation δs. Their respective expressions are given by where τ i marks the start of the ekpyrotic phase. A numerical example showing the evolution of complexity is shown in Fig. 8 for = 12 and k = 0.01. The figure shows both adiabatic and entropy perturbations, during ekpyrosis and kination, for both definitions of complexity. The adiabatic perturbations were discussed in section 5, here the only novelty is that the evolution is followed into the kinetic phase. From the figure we can see that the canonical complexity grows significantly during the kinetic phase, even though the adiabatic modes evolve little, are not amplified and remain unsqueezed. This is again a reason to prefer the definition of complexity in terms of the physical variable, which remains essentially constant for the adiabatic modes throughout ekpyrosis and kination. Our real interest lies with the entropy perturbations. For these, the physical complexity grows significantly during the ekpyrotic phase, while being constant during the kinetic phase. With the aforementioned analytic series expansion for the correlators, we may find the following approximations for the complexity, , (sub-horizon) (6.21a) (super-horizon, ekpyrotic) (6.21b) C δs (τ ) C δs (τ e-k ) = const. , (super-horizon, kination) (6.21c) and hence the complexity growth can be written as , with respect to f σ and g s for adiabatic and entropy modes, respectively), which we argue in the text to be somewhat misleading. A better measure of complexity is that of the physical perturbations, shown in the right panel (so here C R denotes both C R and C δs , i.e., the complexity when taking the correlator to be A RR = z 2 A σσ and A δsδs = a 2 A ss for adiabatic and entropy modes, respectively). The complexity of the adiabatic curvature perturbation remains very small, but that of the entropy perturbation is strongly enhanced. At a later stage, the entropy perturbation is envisaged to act as a source for the curvature perturbation, so that in the end the curvature perturbations will inherit the complexity of the entropy perturbations.
where the second approximations assume a large EoS . Note that, sub-horizon, the physical complexity first increases a little simply due to the overall re-scaling caused by the background, but as the numerical example in Fig. 8 shows, this increase is eventually overwhelmed by the growth of the imaginary part of A δsδs on large scales, cf. Eq. (6.15). We may thus approximate the total growth of complexity during the ekpyrotic phase by 2 √ 2N in the large limit. We immediately notice that this is twice as large as the super-horizon growth of complexity in inflation [recall Eq. (4.5d)]. Of course, the total growth depends on the time spent by a mode on sub-and super-horizon scales. In particular, an entropy fluctuation that exits the horizon very late during ekpyrosis (so that its super-horizon evolution is very short) acquires very little complexity, especially compared to a similar inflationary curvature fluctuation. This is perhaps the most important difference between ekpyrosis and inflation: the sub-horizon growth of physical complexity is mitigated in ekpyrosis in the large limit, while it is significant for inflation, recall Eq. (4.5c). Another difference is that the physical complexity is further restrained due to the kinetic phase in the ekpyrotic scenario. Indeed, instead of continuing its growth on super-horizon scales, the complexity saturates to a constant value 3 . It is important to point out, though, that in realistic models the kinetic phase should last only a few e-folds before a bounce and reheating occur.
The authors of [21] conjectured the existence of an upper bound on the growth of complexity in contracting universes. This bound is clearly violated by the entropy perturbation studied here, Re(A δsδs ) Figure 9. Graphs of the real and imaginary parts of the correlator A δsδs = a 2 A ss for the EoS parameter decreasing in unit intervals from = 11.1 (lightest curve) to = 3.1 (darkest curve) as a function of conformal time. Here we take k = 0.001, and α s is taken to be 3/2 regardless of the value for since the formula for α s given in the text is not accurate enough for small values of . The vertical, dotted line is the time of horizon exit. As gets smaller, the imaginary part of the correlator grows by less and less (always in absolute value). The numerical evolution has been normalised here such that the correlator takes the same value at ekpyrotic-kinetic matching. In the plot there are always 20 e-folds of ekpyrosis followed by 5 e-folds of kinetic evolution.
because their growth depends just as much on the transverse potential as on the background evolution. Moreover, we would like to point out that analysing the growth of complexity as a function of physical time may be misleading. We saw that for de Sitter dC/dt √ 2H, but a similar derivative for super-horizon entropy modes in ekpyrosis yields dC/dt 2 √ 2 |H|, assuming is so large that the scale factor is approximately a constant and hence t ≈ τ . Thus, one could claim that the growth rate of the complexity of entropy modes in ekpyrosis can be made as large as wanted (it is unbounded by ), and furthermore it is not constant (it keeps growing as the universe contracts since |H(t)| grows). However, as we saw, the growth rate really is bounded for a fixed number of e-folds of ekpyrosis, hence we believe the more appropriate time variable is N . In that sense, if we think of the growth rate dC/dN as characterising the chaotic nature of the perturbations, then we obtain a hierarchy for the models on large scales as dC assuming the entropy fluctuations δs in ekpyrosis are later converted into curvature perturbations R. We will comment on the interpretation of this hierarchy in the discussion section.
It is interesting at this point to explore a little more the dependence of the complexity in ekpyrosis on the EoS parameter, in particular when it is only marginally in the ekpyrotic domain -a numerical example of this situation with decreasing in unit intervals from = 11.1 to = 3.1 is shown in Figs. 9 and 10. In all cases, the real part of the correlator becomes small, hence the entropy modes are amplified. However, for small , the imaginary part grows significantly less. This results in a slightly smaller growth of the complexity of entropy modes. When the EoS is only marginally larger than 3, e.g. = 3.1 for the darkest curve in Figs. 9 and 10, the evolution and behaviour of complexity starts deviating: the growth rate is still proportional to √ 2/ on sub-horizon scales, but this is now O(1) rather than being suppressed by a large ; and on super-horizon scales, while the growth rate is now smaller the smaller is, the behaviour changes upon the transition to the kinetic phase. While the complexity saturates (or grows extremely slowly) during kination for away from 3, when gets very close to 3 the complexity continues to grow as ln [ln(−kτ )] 2 . Indeed, as we can see from Fig. 9, | Im A δsδs (τ e-k )| is actually smaller than A δsδs (τ i ) for the darkest curve, so the constant term (Im A δsδs ) 2 /(A δsδs (τ i ) Re A δsδs ) is suppressed compared to the growing term A δsδs (τ i )/ Re A δsδs in X δs (and correspondingly in C δs ). Let us point out that if the imaginary part of the correlator changes too little, it will be difficult to achieve a high degree of classicality of the perturbations. In such extreme cases, one would have to perform a more rigorous calculation of decoherence and the quantum-to-classical transition, along the lines of [16]. Models with small background are nevertheless easier to construct in supergravity (see for instance [83] and the discussion in [84]). One should point out however that such models require a long ekpyrotic phase, with a very large field displacement, which may be difficult to implement in a reliable effective theory [84] (for a more general discussion, see e.g. [85][86][87] and references therein).
The evolution of complexity will change as the universe enters the bounce phase. Just before, during or just after the bounce, two-field ekpyrotic models envisage a process which uses the entropy perturbations as a source for the adiabatic perturbations. A simple incarnation of this idea, motivated by the original colliding brane ekpyrotic scenario [88], is that in the effective theory a bending of the trajectory on scalar field space will occur. Another possibility is that the timing of the bounce itself is modulated by the entropy perturbation, and thus the timing of reheating is modulated [77]. Whatever the details of the process may be, during this conversion process the adiabatic perturbations (which due to their blue spectrum were essentially absent on large scales) inherit the large-scale properties of the entropic perturbations, in particular the nearly scale-invariant spectrum, the large amplitude and also the complexity. In this way, largescale density perturbations are generated at the start of the hot big bang phase. If the universe reaches thermal equilibrium, then the entropy perturbations will be unobservable later on. This in itself constitutes an efficient process of decoherence, rendering the curvature perturbations classical [16]. Thus via the so-called entropic mechanism, two-field ekpyrotic models have the potential to explain the observed primordial perturbations -their complexity being at the root of the later complexities seen in the globally expanding and locally gravitationally collapsing universe.

Matter domination and adiabatic ekpyrosis as other alternatives
In our review of adiabatic perturbations in section 3, which was then applied to inflation ( 1) and ekpyrosis ( > 3), we always had the underlying assumption that the EoS is actually constant (or close enough to a constant as a leading-order approximation). In more generality, this might not hold, and one may as well construct different early universe scenarios in which the EoS has an important time dependence. The more general equation of motion for the canonical variable of a single adiabatic mode, v ≡ zR, is still given by Eq. (3.12), except with z ≡ a 2 /c s , dy ≡ c s dτ , and where a prime in the equation of motion now denotes a derivative with respect to the rescaled conformal time y. This way, one allows for a sound speed c s possibly different from unity and time dependent, as may arise with a scalar field having a non-canonical kinetic structure. If z ∝ |y| n , then z /z = n(n − 1)τ 2 , and the requirement for scale invariance, z /z = 2/τ 2 , is achievable if n = −1 or n = 2. There are thus many possible ways in which a(τ ), (τ ), and c s (τ ) may yield a scale-invariant power spectrum. If we consider c s to be constant for simplicity, this means y = τ , and so one needs a √ ∝ τ 2 or a √ ∝ |τ | −1 . If is also constant, then a ∝ τ 2 corresponds to the matter-dominated contracting scenario [7,89,90] (see also [91] for a recent exposition of the so-called matter bounce scenario and its issues), while a ∝ |τ | −1 is slow-roll inflation (de Sitter to this level of approximation), which correspond to the two 'standard' adiabatic structure formation scenarios. Conversely, we could explore cases where a is essentially constant (as in the original ekpyrotic scenario), but where is (strongly) time dependent. In this case, we see that √ ∝ |τ | −1 constitutes an interesting scenario in which the EoS starts small (inflation like, though the universe is very slowly contracting) at early times (τ → −∞) and increases to large ekpyrotic-like values at late times (τ → 0 − ). This is the basis of the adiabatic ekpyrosis 4 model [98,99] (note that in this model adiabatic perturbations behave differently from the adiabatic perturbations in ordinary ekpyrosis studied in sections 5 and 6). Matter-dominated contraction is an adiabatic scenario with constant EoS = 3/2 (so a(τ ) ∝ τ 2 ). Thus most expressions of section 3 are immediately applicable, in particular Eq. (3.26b) for the correlator A σσ on large scales with α σ = 3/2. Multiplying by z 2 = 3ā 2 0 τ 4 yields Re A RR ∼ τ 6 and Im A RR ∼ |τ | 5 , hence the amplification is very efficient as 1/ Re A RR ∼ τ −6 and so is squeezing with Im A RR / Re A RR ∼ |τ | −1 . Regarding complexity, from Eq. (3.31) one has C R 2 √ 2N on sub-horizon scales, while the evolution on super-horizon scales can be approximated as C R (τ ) C R (τ ) + 3 √ 2 ln(τ /τ ), hence the super-horizon complexity growth is ∆C R 3 √ 2∆N . We note that those are the largest growth rates of complexity encountered so far in this work. We present a full numerical calculation of the complexity during matter contraction in Fig. 11 where this is explicit.
Let us now turn our attention to adiabatic ekpyrosis. In a realistic adiabatic ekpyrosis scenario, the phase during which the EoS rapidly evolves can only last for a few e-folds at most. It is then followed by a standard ekpyrotic phase where the EoS settles to a (large) constant value. This phase can last longer, but perturbations exiting the horizon during this phase will not acquire a scale-invariant power spectrum. Thus, appropriate model building gives one a scaleinvariant power spectrum over the range of modes that are of observational interest, but which becomes blue outside of this range. Models that can satisfy all the constraints (observational and theoretical [99]) can be engineered with some level of tuning. For the purpose of the present analysis, we will rather consider a toy model from [98], which may not meet all constraints, but which will capture the essential features of adiabatic ekpyrosis with regard to complexity. Specifically, let us parameterise the evolution of the function z as where c and H 0 are constants. It is assumed that the scale factor is almost a constant throughout (e.g., a(t) ∝ (−t) 2/c 2 ≈ 1 with c 1), so that physical time and conformal time are approximately equal, t ≈ τ . Then at early times (τ → −∞), the EoS rapidly evolves as desired, ∼ 1/τ 2 , while at late times (τ → 0 − ) the EoS tends to a large constant, c 2 /2. From this, the function z behaves as in inflation (more precisely the de Sitter limit → 0) at early times, z ∼ 1/|τ |, which does not only imply the same scale-invariant power spectrum on large scales but also the same complexity growth for that time period on sub-and super-horizon scales, specifically C R √ 2N . Modes of observational interest have thus evolved to large scales when the EoS moves closer to its constant, late-time value. At that point, we expect the perturbations to be matched with the super-horizon mode solutions of standard single-field ekpyrotic cosmology and acquire the corresponding complexity evolution. In that limit, z is approximately constant, so from the mode equation in the far infrared (k → 0), R k + 2(z /z)R k 0, we can see that v k ∼ R k ∼ c 1 + c 2 |τ | for some integration constants c 1,2 that are obtainable upon matching. Therefore, to leading order the correlators A σσ and A RR and the corresponding complexities, C v and C R , all tend to constants at late times. This can be verified by taking Eq. (7.1) and solving the mode equation (3.12) by means of numerical methods. Subsequently, using the numerical solution for the mode function, one can compute the correlators and then the complexity. The result is shown in Fig. 11 in dark blue. In the rapidly evolving EoS phase, the complexity behaves as in inflation (light pink). On sub-and super-horizon scales (for N 11.5 and 11.5 N 17), the curvature complexity C R grows as √ 2N . Then, as the EoS tends toward its constant value (i.e., z ∼ √ ∼ const.), happening when N 17 (still on super-horizon scales), the complexity tends to a constant, i.e., it saturates. This is in agreement with the analytical approximation derived in the previous paragraph.

Discussion and conclusions
All viable early universe models must in some way be able to explain the fluctuations observed in the CMB. In the present work, we have shown how quantum circuit complexity provides a useful characterisation of the different ways in which cosmological models achieve this goal. The models we have analysed (inflation, ekpyrosis and a contracting matter phase) all rely on the amplification and squeezing of quantum perturbations. But the details of how these phases proceed differ markedly. A useful summary of our results is provided by Fig. 11, and let us also recall the following super-horizon evolutions: (8.1) Two main features are immediately obvious: the complexity that is achieved depends primarily (essentially linearly) on the number of e-folds of evolution. And the coefficient of proportionality depends on the cosmological model, i.e., it serves to distinguish the different models.
The growth of complexity is smallest for the most popular early universe model, namely inflation. (In adiabatic ekpyrosis the growth rate is the same, before it caps off.) Thus inflation acts as a comparatively 'simple' quantum computer in drawing quantum perturbations out of the vacuum and turning them into effectively classical density perturbations. Contracting models have a higher growth of complexity, hence they are more 'efficient' at quickly producing a complex system, and in particular a contracting matter phase leads to the largest complexity. Ekpyrotic models reside in between inflation and matter contraction, and they possess the distinguishing feature that on sub-horizon scales the growth of complexity is very small, so that essentially the entire complexity comes from super-horizon evolution, cf. again Fig. 11. It is interesting that the models come out as being so clearly distinguished. Note in particular that within each class of models the specific dependence on the equation of state is rather modest. Also, a dependence on the wavenumber comes about only through its influence on the time of horizon exit, and even this only when there is a significant difference in the growth of complexity before and after horizon crossing. Moreover, there is no explicit dependence at all on the energy scales involved (e.g. of the potential). Thus complexity provides a truly complementary characterisation of cosmological models, more attuned to their quantum properties.
In order to define complexity we have used a measure that was developed in particular for Gaussian states and that is related to the Sp(2, R) symmetry of the associated quantum mechanics [22,23]. This measure is conceptually appealing, as it provides a link with hyperbolic geometry -see Fig. 2 for a useful visual illustration of the evolution of cosmological correlations (using the same numerical models as in Fig. 11). Moreover, the hyperbolic measure has a structure that is sensitive to both amplification and squeezing, which are precisely the features that are important for early universe models. A comparison with another popular measure is provided in appendix A. It will be important to see how the present study can be generalised to non-Gaussian corrections, which are bound to play a significant role in future observations.
We have focused on the complexity of the physical perturbations and contrasted it to that of the re-scaled mode functions. This distinction ends up being rather crucial. The canonically normalised mode functions are not directly sensitive to the overall expansion or contraction of spacetime, and thus they miss the sometimes vast changes of physical wavelength that various cosmological models cause (however they are highly sensitive to the changes in physics occurring near horizon exit, or at junctions with different phases of evolution). Physical perturbations are not only the ones that are directly related to observations, but they depend much more crucially on the entire history of a cosmological phase, and in particular they are sensitive to the conditions at the beginning. Thus the complexity of physical perturbations offers the prospect of better characterising the initial conditions for the cosmological models that are studied, which will be important in terms of incorporating such phases into a complete cosmology. In this respect we suspect that useful links with the puzzles of trans-Planckian perturbations may also be developed in future work.
An additional theoretical avenue that deserves further exploration is the relation of complexity to chaos. For inflation, these issues are already understood to some extent, principally because inflation may be regarded as a thermal system. But for alternative cosmological models, in particular contracting models, such an identification is not available. Yet some chaotic features are certainly present in the dynamics of such models, and it would be interesting to see to what extent complexity may provide a useful diagnostic of these. As is often the case, the confrontation of ideas from different parts of physics is likely to lead to fruitful new insights.
A Comparison with the analytically continued complexity formula As mentioned in section 2, a straightforward generalisation of the simplest measure of complexity is obtained by analytically continuing that formula to the case where the frequencies involved may be complex. This approach, and its implications for some cosmological models, has been pursued in [19][20][21]. In many circumstances, the two approaches yield qualitatively similar results, but there are exceptions.
To illustrate this, it is useful to compare Fig. 10, showing the evolution of hyperbolic complexity C for ekpyrotic models with relatively small EoS, with Fig. 12, which shows the analytically continued complexity C (a.c.) for the same models. As can be seen from the figures, the evolution of hyperbolic complexity follows a natural progression as the EoS is lowered, while the analytically continued complexity starts showing bizarre features when the EoS approaches the lower bound = 3. This may be understood from the following heuristic rewriting of the two definitions of complexity. The inverse of the real part of the correlator determines the amplification A, while the ratio of the imaginary to the real part of the correlator is a measure of the squeezing S. Then the two definitions of complexities may be heuristically written as C (a.c.) ∼ 1 √ 2 ln 1 A 2 1 + S 2 . (A.1b) Thus we see that the hyperbolic complexity is separately dependent on the amplification and the squeezing and can grow when either of these properties evolves. By contrast, for the analytically continued complexity, amplification and squeezing may counteract each other to some extent. This is exactly what happens for ekpyrotic perturbations when is small, since one can see from Fig. 9 that in such a case the imaginary part of the correlator grows only slowly (though it is quite large), while the amplification proceeds without much change compared to the cases with a larger EoS. This has as a consequence that the analytically continued complexity is reduced again, even though the perturbations are both amplified and squeezed. This then leads to the misleading perception that the perturbations evolve just as much during a short phase of kination as during the preceding ekpyrotic phase. The hyperbolic definition avoids this pitfall and seems better suited to us in order to characterise cosmological perturbations.