Circuit Complexity through phase transitions: Consequences in quantum state preparation

In this paper, we analyze the circuit complexity for preparing ground states of quantum many-body systems. In particular, how this complexity grows as the ground state approaches a quantum phase transition. We discuss different definitions of complexity, namely the one following the Fubini-Study metric or the Nielsen complexity. We also explore different models: Ising, ZZXZ or Dicke. In addition, different forms of state preparation are investigated: analytic or exact diagonalization techniques, adiabatic algorithms (with and without shortcuts), and Quantum Variational Eigensolvers. We find that the divergence (or lack thereof) of the complexity near a phase transition depends on the non-local character of the operations used to reach the ground state. For Fubini-Study based complexity, we extract the universal properties and their critical exponents. In practical algorithms, we find that the complexity depends crucially on whether or not the system passes close to a quantum critical point when preparing the state. For both VQE and Adiabatic algorithms, we provide explicit expressions and bound the growth of complexity with respect to the system size and the execution time, respectively.

How much does it cost to generate a target quantum state from another reference state?This is a rather general question that has been discussed in quantum information for obvious reasons.In quantum computation it is desirable to obtain the result with the minimum set of gates.This number is, roughly speaking, the computational cost and it is called circuit complexity (C) [1][2][3].It is, let us say, the quantum analog of the concept of computational complexity in computer science.Importantly enough, this cost builds upon a concrete physical architecture, i.e the available set of gates.Therefore, C not only depends on the reference and target states but on the restrictions for reaching the latter.This is quite natural if one thinks of an actual quantum computer where the possible operations have restrictions.Note that, if any unitary were allowed, a simple rotation would achieve the goal and every quantum state would be easily prepared, so that (essentially) the complexity would be a trivial quantity.Therefore, also in analytic calculations, the path between the reference and target is restricted to a set, e.g.gaussian states [4][5][6][7][8] .Beyond quantum computation, circuit complexity is a relevant concept in quantum gravity.In particular, for its consequences in holography [9][10][11].For those who are not experts (like us), we can say that holography describes quantum gravity within a region of space by looking at the boundary of that region, that is described by a non gravitational theory.Then, any bulk quantity in the gravitational theory is equivalent or dual to another quantity in the boundary of the non gravitational theory.One of the main problems of this duality is that the volume behind the black hole horizon keeps growing for a very long time while the entanglement at the boundary saturates at much shorter times.One possible solution is to conjecture that the dual of volume is not entanglement but complexity, via the identification Complexity = Volume.This is because we expect that volume is an extensive quantity, while entanglement (typically) fulfils an area law.Therefore, the calculations of complexity are beyond the quantum information community and different calculations in field theories have been discussed in the literature [12,13].
The notion of complexity (C) is much related to the geometry of states (or operators).It is a measure of the distance between two of them.Therefore, one possible choice for C is finding the geodesics in the Fubini-Study metric in the projected Hilbert space for the case of pure states.For mixed states different measures have been introduced via state purification [14] or distance measures for mixed states as the Bures distance [15].This geometric background is a powerful way to understand complexity, since it allows us to know how much it will cost to prepare a state by solving a geodesic equation.It is true, however, that the metric, in principle, can only be obtained in some cases: surely in exactly solvable models.And there we know how to prepare states.Thus, it is interesting to be able to predict the typical behaviour in general models.Here, we move in this direction.
In this article we are interested in a quite generic situation, i.e. when a critical point is crossed to reach the target state.In particular, we investigate what general statements about the behaviour of the circuit complexity we can make.We are not the first to calculate C in a quantum phase transition (QPT) [16].Recent papers discuss exactly solvable models as the topological Kitaev, Bose-Hubbard Lipkin-Meshkov-Glick ones and conformal field theories [17][18][19][20][21][22][23][24][25].Importantly enough, complexity has been shown to be a useful probe of topological phase transitions.Complementary to these calculations, in this work, we use that close to a transition point, the concept of universality emerges naturally, so we expect these universal properties to be inherited by complexity.If so, we can argue for its scaling laws or how complexity behaves regardless of model details or even on the particular chosen definition of complexity.In addition, we apply our theory for state preparation in quantum computers.This is a key and hard task [26].It is within the QMA complexity class [27], roughly speaking the NP-complete analogue for quantum computers.Nevertheless, quantum computers are expected to be better than classical methods such as density functional theory [28], density normalization group [29], tensor networks [30], quantum Montecarlo [31] or even ML-inspired techniques [32], in some instances.For a recent discussion of these issues, see [33].Heuristic quantum algorithms as adiabatic [34] or varational ones [35,36] can outperform classical calculations and serve for the generation of quantum states as quantum data, e.g.phase classification [37].Motivated by all of this, we discuss how useful the concept of complexity is and how much one can anticipate the difficulty of state preparation in variational quantum eigensolvers (VQEs) or adiabatic quantum algorithms (with and without shortcuts to adiabaticity).To challenge our theory we tackle both integrable and non-integrable models using numerical simulations and computing C.

A. Complexity overview
We find it convenient to discuss first the different notions of circuit complexity that we will use in this paper and the relationship between them.

Complexity à la Nielsen
The original notion of complexity is due to the works of Nielsen and collaborators [1][2][3].See [13] for a recent review.Restricting ourselves to unitary operations, target and reference states are related as T stands for time ordering.Notice that, A Cost function is formally defined as: with F some functional fulfilling some basic properties as continuity, homogeneity (F [U, λ U ] = λF [U, U ] for λ ≥ 0), positivity and the triangular inequality1 .If, in addition to these, smoothness is assumed and the Hessian of F as a function of U is strictly positive, the functional is a Finsler metric.Thus, C N is nothing but the geodesics.The suffix N stands for Nielsen.
Being a little more explicit, we can write that the evolution is given by with O n some operators and Y (n) (τ ) parameters.A usual functional is then given by, If we restrict ourselves to two level systems (qubits), F k (τ ) is a natural distance in SU (2 n ), such that d = t 0 dτ F k (τ ), Cf. with Eq. (3).What has been explained so far is the continuous version of complexity, that provides a lower bound for the number of gates needed to approximate |ψ T ⟩ from |ψ R ⟩ [1].The discrete version of C N can be computed introducing the functional (using the same notation as in the original [1]): where the Hamiltonian in this case is In the first sum, σ ranges over all possible one-and two-body interactions, that is, over all products of either one or two qubit gates.In the second sum, instead, the sum is over other tensor products of Pauli matrices and the identity.The factor p > 0 penalizes three, four, ... -body interactions.All put together, finding the geodesics in the continuum version is a good estimate of the resources needed to prepare a state.
At this point, we think it is necessary to emphasise something.If any unitary is possible, the complexity has a narrow utility, since its value is given by C = arccos(|⟨ψ R |ψ t ⟩|), i.e. of the order of one (it doesn't matter which state reference and destination are chosen).This can be verified by noting that the target state can always be written as A rotation in the subspace generated by {|ψ R ⟩, |ψ ⊥ R ⟩} does the job.Therefore, some restrictions on the possible unitaries or Hamiltonian (4) will be imposed.We will discuss this point in some depth later.

Circuit Complexity from the Fubini-Study metric
The functionals F discussed so far, see Eqs. ( 4) and (5), are not unique.Others can be chosen satisfying continuity, homogeneity, positivity and the triangular property.We want to discuss next another possibility where the distance between the reference and target states is given by the Fubini-Study metric.Originally proposed for Quantum Field Theories in [5], we prefer to study it here from a quantum information perspective.Let us time-slice the unitary (1) such that We have assumed that the unitaries and so the wave functions depend on the parameters λ.Then, for sufficiently small time step δτ := t n − t n−1 , the fidelity between two contiguous states is where χ F is denoted the fidelity susceptibility [38][39][40][41][42][43][44], see Ref. [45] for a review.Interestingly enough, we can relate χ F with the geometric tensor, in fact [Cf.Eq. ( 11)] with [5,46] g µν = Re (T µν ) .
Here, T µν is the quantum geometric tensor which is nothing but the Fubini-Study metric (FSM) on the CP n manifold, namely: with Another useful way of writing the metric tensor is as follows.Using a formal Taylor expansion for the states, the metric tensor can be written as, Setting now in the Hamiltonian (4), λν ≡ Y (ν) then |∂ ν ⟩ = O ν |ψ⟩, it is straightforward to see that [6], i.e. the fluctuations of the Hamiltonian operators O ν .
Using the fact that 1 − F n,n+1 is a distance, thus satisfying the properties we imposed for the F -functional, we have that we can understand C as the distance defined through the Fubini-Study metric: The suffix stands for Fubini-Study metric and the notion of distance is quite explicit.This is an alternative definition to that given by Eq. ( 3) that has some remarkable properties.The first one is that knowing the metric tensor the geodesics can be found, at least in principle, by solving the differential equation: Here, Γ are the Christoffel symbols: The second property of C FS is that, from its relation to the fidelity between states, F, its properties close to a QPT can be used when discussing the complexity, C, see also Eq. ( 13).

Some remarks comparing both approaches
The Nielsen complexity estimates the minimum number of gates required to transform an initial reference state to a target state.It operates by defining a set of available gates and a metric in the space of quantum circuits or unitary transformations.The optimization process within this space aims to minimize the number of gates.The Nielsen approach is operational, providing an estimation of the minimum number of operations based on the reference and target states, along with the available gates.On the other hand, the Fubini-Study approach considers a manifold of quantum states and defines the cost as the distance between pure quantum states using fidelity, as given in equations ( 8) and (17).In this sense, the Fubini-Study metric serves as a geometric measure, while the Nielsen complexity directly relates to the practical cost in a quantum computer.
In practical terms, there are notable differences between the two approaches.The Fubini-Study metric assigns a variable cost to specific gates based on the states they act on, while Nielsen assigns a fixed cost to each gate.Additionally, the Nielsen complexity may involve degenerate operations that leave the state unchanged (e.g., adding a global phase), resulting in a higher dimensionality of the space of unitaries.Different results are obtained using these approaches, depending on the specific application.Each form has its preferred use case.The geodesic equations provided by the Fubini-Study metric make it suitable for analytical calculations, while the Nielsen complexity is more relevant for cost estimation in quantum computation and state preparation.However, in some cases, both methods yield identical results, such as in the preparation of Gaussian fundamental states [6].

B. Main results and manuscript organization
For the exactly solvable systems that we discuss in this work, we find that C FS ≥ C N when crossing a phase transition.We understand this inequality as a consequence of the fact that the unitary space is larger, see previous subsection I A 3. In any case, C does not diverge at the critical point, its derivative does.For C FS we can characterize this divergence and its critical exponents in rather general circumstances.Let us remark, again, that throughout the paper we focus on the extensive part of C. Two models are studied in detail, namely the one-dimensional quantum Ising and Dicke models.After this general discussion, we focus on calculating the complexity when preparing a fundamental state in a quantum computer.Here, obviously, we compute C N in its discrete version.We explore two algorithms in detail.First, we discuss the circuit complexity in adiabatic algorithms with and without shortcuts to adiabaticity.We focused our study on one-dimensional spin lattices of different sizes.In this investigation, we found that using shortcuts does not significantly alter the complexity C N .However, we demonstrated that C N ∼ √ L × T , where L represents the system size, and T is the total time required to achieve a fixed fidelity, F, with the exact ground state (in our case, F = 0.9 was chosen).Thus, the complexity inherits the behavior of T close to a Quantum Phase Transition (QPT).Specifically, T is bounded by ∆ −2 , where ∆ represents the minimum gap between the ground state and the first excited state in the adiabatic algorithm.Then, we discuss the circuit complexity using VQEs.These algorithms are variational and do not need to cross the critical point even if the reference and target are in different phases.In such a case, C N is not necessarily aware of the QPT.On the other hand, if the target state is close enough to a phase transition, also in VQEs, the complexity grows.Importantly, we provide an explicit formula for C N , and by combining it with the correlation length generated using local Variational Quantum Eigensolver (VQE) ansatzs, we can show that C N ≳ L 3/2 .Therefore, this scaling poses challenges for our numerical capabilities, explaining the difficulties in finding reliable solutions around Quantum Phase Transitions (QPTs) when simulating the action of a VQE.
The rest of the manuscript is organized as follows.In the next section, II, we discuss the relation between circuit complexity, in this case C FS from Eq. ( 14), and the geometry of quantum states that allows extracting the critical exponents for the derivative of C.This is our first result that emphasizes that through phase transitions C FS has universal properties.In section III we perform explicit calculations for C FS and C N in two solvable systems, namely the one dimensional XY-anisotropic and Dicke models.We extract the critical exponents.Then, in section IV, we perform numerical simulations where C N is computed in two types of algorithms: variational and adiabatic ones.Concretely we benchmark with exactly solvable models as the Ising model, and we complement our discussion with non-integrable Hamiltonians as the ZZXZ model.Lastly, we discuss these results and conclude the paper in V. Some technical issues are left for the Appendices.The code used to obtain the numerical results is available upon request.

II. COMPLEXITY AND THE GEOMETRY OF STATES CLOSE TO A QUANTUM PHASE
TRANSITION.
In this section, we discuss general aspects for the complexity close to a QPT.To be as general as possible, we find it convenient to focus on C FS , Eq. ( 14).Within this geometric formalism, we see that, in general, the complexity is finite, but not its derivative, which can diverge when crossing a QPT.We study its finite size scaling obtaining the corresponding critical exponents.
A. Complexity and its derivative when crossing a QPT We have already argued in section I A 1 that if we are allowed to use any unitary, C is of the order of one.In the literature, several unitary restrictions have been used: considering one and two qubit gates or considering gaussian states when moving from reference to target states.In this subsection, we consider another kind of restriction, quite natural when talking about a QPT.We will consider that one (and only one) parameter, say λ, of the Hamiltonian model is varied to pass through the QPT, keeping other variables or parameters fixed.Thus the metric is unidimensional.We know, that in this case, the geodesic is given by: Therefore, Below, we will work some examples and we will see that T does not diverge at the QPT.However, if we compute the derivative instead: It is known that some components of the metric tensor can diverge, thus diverging the derivative of C. Equation ( 19) has two consequences.The first one is that, under quite general circumstances, the derivative of C close to a QPT is related to the metric tensor and inherits its universal properties.Thus, we can utilize the theory of the metric tensor and its behavior near a transition to automatically obtain the scaling for complexity.The second one is that this derivative can be used to witness and characterize QPTs.

B. Finite size scaling
Here, we review the scaling for the metric tensor [47], from which the behaviour of C follows directly, Cf.Eq. ( 19).Close to a critical point correlation length diverges as, with a a critical exponent.Similar relations occur for other thermodynamic quantities.In particular, and for what interests us, the metric tensor can be written as [47], with ∆ µν the corresponding critical exponent.Notice that, for the reasons already explained in section I A 3, from now on we will be interested in the intensive part of the metric tensor g µν → g µν /L d , with d the spatial dimensions.Near a phase transition, finite-size scaling dictates how quantities behave after scale transformations.Very briefly, after a length scale transformation x ′ = αx, time scales as τ ′ = α z τ , with z its critical exponent.This fixes the energy fluctuations ∆E∆τ ∼ 1 → ∆E ′ = α −z ∆E.Putting it all together, it is interesting to extract the value of the critical exponent ∆ µν above, which controls how the metric tensor behaves, in terms of other critical exponents that dictate more fundamental quantities.Looking at equations ( 4) and ( 12) and ( 13) and writing the scaling for the derivatives of the Hamiltonian as ∂ µ ′ H ′ = α −∆µ ∂ µ H we arrive to [47], Finally, merging, ( 20) and ( 21), we find that close enough to the transition, where the relevant length is given by the system size, L, we arrive to As a consequence of all of this and using (19), when a single parameter is varied across the QPT we have the scaling: It is remarkable that the complexity derivative scaling is dictated by universal exponents, whenever one parameter is varied to cross a critical point.In particular, if ∆ λλ > −2 the derivative is sub-extensive.If ∆ λλ = −2 it is extensive and if ∆ λλ < −2 is superextensive.

III. SOLVABLE HAMILTONIANS
Let us test the above ideas on a couple of solvable models: the one dimensional quantum Ising model [48] and the Dicke [49][50][51] model.

A. Quantum Ising model
The transverse field Ising model (Periodic Boundary Conditions will be assumed) is Hamiltonian ( 25) can be solved via the Jordan-Wigner transformation [48].This Ising model has a second order phase transition occurring at J c = 1(−1) in the L → ∞ limit.For J c > 1(J c < −1) the Z 2 symmetry is spontaneously broken and the g.s. is ferromagnetically (antiferromagnetically) ordered.W.l.o.g.we fix our attention in the paramagneticferromagnetic transition occurring at J c = 1.On top of that, the ground state can be written in terms of fermionic excitations (after the Jordan-Wigner transformation) as, For the rest of the section the metric tensor ( 12) is needed.It has been computed several times already [52,53] In the thermodynamic limit, the k-sum is an integral k → L/π and it can be computed explicitly, yielding

Complexity through QPTs
From Eq. (29) we see that g JJ diverges at J = J c .This is the reason behind the divergence in the derivative of the complexity at the QPT, Cf.Eq. (19).In figure 1, we plot C FS both in the continuum and for L-finite using either (29) or the sum (28).In both cases, the integral (18) is computed.It is evident that the complexity does not diverge at the QPT, but its derivative does, inheriting this behaviour from the metric tensor, Cf.Figs.1a and b.For the Ising transition, the exponent a = 1, Cf. Eq. (20).We know that ∆ hh /a = 1, so the complexity derivative diverges as ∼ L 1/2 at the Ising transition.

Relation between CN and CFS
Formula ( 26) is formally equivalent to the ground state for the 1D-Kitaev model.For the latter, C N has been computed in [18].If the reference, target and intermediate states have the same form (26), C N reads: where are the angles (27) at the target (reference) states.Following the same procedure as in [18] we checked that ∂ J C N ∼ log L, i.e. it diverges logarithmically.This must be confronted with the divergence (with critical exponent 1/2) for the case of ∂ J C FS .This is an important difference.While using the FS distance the complexity is associated with the fluctuations, cf.Eq. ( 13), the C N is more related to the angles difference and its divergence is therefore smoothed.

B. The Dicke model
The Hamiltonian for the ground state sector of the L-spin Dicke model can be written in terms of total spin operators of spin S = L/2 as [54] where the spin and ladder operators obey the canonical commutation relations [S z , S ± ] = ±S ± , [S + , S − ] = 2S z .This model can be solved in the thermodynamic limit, S → ∞, with a Holstein-Primakoff transformation on the spins yielding In the normal phase of the Dicke model we can obtain an effective Hamiltonian for S → ∞ by neglecting terms with 2S in the denominator in the Hamiltonian of Eq. ( 36), resulting in a completely symmetric model of coupled harmonic oscillators, one corresponding to the physical oscillator and the other corresponding to the spins within the Holstein-Primakoff transformation In the superradiant phase, the bosonic modes must be displaced to accommodate the macroscopic occupations that the spins and field develop in this phase.Once the displacements are introduced, terms with powers of 2S in the denominator are again neglected in the thermodynamic limit, yielding where µ = ω z Ω/ 4λ 2 and ā, b are the displaced bosonic operators [55].We omit the expressions of the displacement as they are irrelevant in the following.Both the normal and superradiant effective Hamiltonians can be diagonalized in the real space basis, where they present a gaussian profile given by where R = (x, y), x and y are the real-space coordinates associated to the modes a(ā) and b( b), A = U −1 M U with U a unitary matrix, M = diag [ϵ − , ϵ + ] and ϵ ± are the eigenmodes of the system [56].The overlap of two different ground states is given by This allows us to compute the components of the quantum metric tensor for the Dicke model exactly in the thermodynamic limit.We combine this with finite size results from exact diagonalization of Hamiltonian (31).The results are shown in Fig. 2. Just like we showed for the case of the Ising model, there is no divergence in C FS , the only signature of the phase transition is a non-analiticity that is only noticeable in the L → ∞ case.This non-analiticity, or its precursor in the case of finite L is best revealed as a divergence in the derivative of the complexity, which is naturally the square root of the metric tensor.Here we are considering the complexity along a λ-path and the divergence is revealed in ∂C FS = √ g λλ .We perform a finite-size scaling analysis of the metric tensor by fitting the maximal values (∂C FS ) max (L) and critical parameters at said maxima λ max (L) to their respective scaling laws The resulting critical exponents ν = 0.655(22) ≊ 2/3 and δ = 0.6711(15) ≊ 2/3 are in agreement with values reported in the literature [57].

IV. COMPLEXITY IN A QUANTUM COMPUTER, THE CASE OF GROUND STATE PREPARATION
In this section, we compute C N when preparing ground states in a quantum computer.We study both adiabatic algorithms and variational quantum eigensolvers (VQEs).Two versions of the former algorithms are discussed: with and without shortcuts to adiabaticity.
In both cases, the initial state is the "trivial zero" |0⟩ ≡ |00 • • • 0⟩4 .Some gates are applied to prepare the ground state of a given Hamiltonian.Here, we are especially interested when this initial state (that can be understood as the ground state in the paramagnetic phase) is in a different phase than the final one.In addition, we discuss whether or not a QPT is crossed during the algorithm.Finally, notice that in quantum computing applications it seems natural to compute C N and, in particular, its discrete version (the number of gates needed), Cf.Sec.I A 1. Thus, through this section, we compute C N .

A. Adiabatic algorithms
A systematic way of finding the ground state of a given Hamiltonian is by adiabatic passage or annealing.Let us consider the time-dependent Hamiltonian: Here, H 0 has a trivial ground state (easy to prepare), and H T is the hamiltonian from which we want to obtain its ground state.Consider that the time-dependent function λ(t) runs from λ(t = 0) = 0 to λ(t = T ) = 1, where T is the final time of the algorithm.At t = 0 the state is prepared in the ground state of H 0 .If λ is sufficiently small compared to the instantaneous gap, by means of the adiabatic theorem the final state is the ground state of H T [34].On the other hand, the adiabatic condition alerts us that as the gap closes, for example in continuous phase transitions, the execution time, i.e. the circuit depth, scales with the inverse of this gap, thus also C.
Importantly enough the adiabatic condition can be relaxed by introducing counter-diabatic terms.Generally speaking, instead of H(τ ) (whose ground states are |ψ(t n )⟩) what is evolved is the "modified" Hamiltonian [58,59]: The last term ensures that the time evolution exactly matches the instantaneous ground state of H(τ ) no matter how fast the evolution is.This is known in the literature as shortcuts to adiabaticity and H CD is called counter-diabatic Hamiltonian.There are different ways of writing H CD .In its original form we can write: with ∂ µ H ≡ ∂H/∂λ µ .We emphasize that at times 0 and t, |ψ R ⟩ and |ψ T ⟩ are ground states of H(0) and H(t) respectively.Explicitly Here τ means time, cf. with Eq (1).To connect this evolution with the previous sections, we note that the fidelity susceptibility can be written in terms of the H CD (τ ) fluctuations [Cf.Eq. ( 9) and ( 13)]: In practice H CD is difficult to find.Therefore, a systematic although approximate way of writing is convenient.Following [60] it can be rewritten as, Here, A is the adiabatic gauge potential that can be approximated as: where (l) is the "degree of approximation".On top of that, the {α k } are found variationally by minimising the action [61]: In many cases of interest, in the adiabatic protocol, H(τ ) is expected to be a local Hamiltonian, in particular a two body one.Notice that due to nested commutators, the higher the order (l), the longer the range of interaction.Following the functional (6) three, four, or higher order body interactions will be highly penalised.Thus, in what follows, we will restrict ourselves to l = 1 that introduces two body interactions at most.This, in turn, provides a systematic way of preparing, via trotterization, quantum states.

Complexity in adiabatic algorithms
As has been previously discussed, in order to compute the complexity as defined by Nielsen [Cf.Sec.(IA 1)], we only need to express our unitary operation as the time evolution of some Hamiltonian.In the present case, it is straightforward, with and without shortcuts, as the Hamiltonian ( 44) is given explicitly.We study the Ising model in transverse field and the ZZXZ model.For both, we use the function λ(t) = sin 2 π 2 sin 2 πt

2T
to drive the evolution from the initial Hamiltonian, H 0 , to the target one, H T .
For the Ising model, we start from H 0 = h x i σ x i , leaving the transverse field fixed at value h x = 1 and switching on the spin-spin interaction until we reach H int = J i σ z i σ z i+1 .The counter-diabatic Hamiltonian follows from equations ( 45), ( 48) and ( 49) with l = 1 yielding H CD (t) = λ(t)α(t) i (σ y i σ z i+1 +σ z i σ y i+1 ), where α(t) is the variational parameter in Eqs.(48) and (49).The full-time-dependent Hamiltonian reads FIG. 3. Complexity study for the Transverse Field Ising model using the adiabatic algorithm.(a) Evolution of the fidelity obtained with shortcuts to adiabaticity (solid lines) and without them (dashed lines) for increasing time lengths of the full algorithm and different target Js for L = 12.At shorter times the shortcuts provide better results, being identical to the simple case (without shortcuts) for the longest times.(b) Evolution of the gap between the ground state and the first excited state during the algorithm for the same values of J as in (a).The gap closes with an increasing value of |J|, explaining why longer times are needed for the larger |J| to obtain the same fidelity.(c) Complexity computed for different sizes with shortcuts (solid lines) and without shortcuts (dashed).As the gap closes, more gates are needed to achieve the fidelity threshold (0.9 in this case) but we do not find relevant differences between applying shortcuts or not in the final result for the complexity.(d) Complexity dependence on the chain size (root squared) and the preparation time.Each point corresponds to the preparation time for each |J| in (c).
Before discussing the complexity in adiabatic algorithms, let us investigate how much time is needed to reach the ground state and its relation to the Hamiltonian gap.It is well known that T determines the practicality of the algorithm, and in this section, we will explore its relationship with C N .To implement the unitary evolution via Trotter decomposition, we need to split the total time T into T /δT steps, where δT is the time discretization employed.A smaller δT ensures a more precise implementation but requires more gates, increasing the computational cost of the implementation.In the simulations presented here, we use δT = min(0.1,T /30).
In Figure 3a, we plot the fidelity between the final state obtained adiabatically and the target state.As expected, longer times result in better fidelity.Additionally, we confirm that at lower times, higher fidelities are achieved thanks to the counter-diabatic term.Figure 3b shows the gap evolution within the adiabatic algorithm.As |J| increases and the gap between the first excited state and the ground state closes, more time is required for the preparation.
Following Eq. ( 6), the Nielsen complexity C N for the adiabatic algorithm is given by: Here, we set the final time T as the time required for the adiabatic algorithm to reach a certain fidelity threshold, which in our case will be F = 0.9.From the above formula, we can extract a √ L. Figure 3c shows the actual Nielsen complexity values.Reflecting the fidelity behavior, the complexity jumps around the transition as the gap is closing.From Eq. ( 51), C N is proportional to T , and T depends on the gap, so this jump in complexity is expected when approaching the critical region.
More interesting is what we show in panel 3d, where C N / √ L shows universal behavior independent of L when plotting it as a function of T .Additionally, we observe that C N ∼ T .This can be understood as follows.First, we can ignore the case with shortcuts, as it barely affects the complexity (see Figure 3b).By doing so, Eq. ( 51) is simplified as follows: Here, J 0 (•) is the Bessel function.In the J-range studied here, the first approximation for the integral is quite accurate.This introduces a dependence on J 2 , however, this dependence is still small, and the overall dependence is very well approximated by C N ∝ T √ L, as confirmed by our numerical results (dashed line in Figure 3 d).As a consequence, we can use some results from the adiabatic theorem relating the final T and the minimum gap in Eq. ( 55).It has been proven that, in the best case, the preparation time scales with the minimum energy gap as T ∼ ∆ −2 [34], with ∆ ≡ min j̸ =k |E j − E k |.This dependency is inherited by the complexity and is shown in Appendix C. There, we verify that ∆ −2 holds as long as the gap is not too small, where the dependency is lost.This may be attributed to several reasons, such as the precision of our numerics or the fact that our 0.9-fidelity threshold cannot discern in those cases.
In order to check if this holds in other models, we also study the so-called antiferromagnetic ZZXZ model: Due to the combination of longitudinal and transverse fields, this is a non-integrable model.It is ideal, then, to explore the phenomenology of complexity beyond the exactly solvable models considered so far.In Fig. 4 we draw the phase diagram of the model at zero temperature as a function of the fields applied to the spins and the exchange constant [62].The critical line separates paramagnetic and antiferromagnetic phases.For our particular purposes, keeping the same initial Hamiltonian, H 0 , we set the transverse field, h x = 1 and the target longitudinal field to h z = 0.75.We thus study the quantum phase transition appearing when moving to different target values of J.This path is shown as the red line in Fig. 4, where the final point marks the maximum value simulated for the target J.Therefore, the transverse field is going to be fixed while we turn on both the longitudinal field and the magnetic interaction.The counter-diabatic term can be computed in the same fashion as before, getting the same result as in [60].The time-dependent Hamiltonian reads and the complexity acquires the following expression In figure 5 we show the results obtained for the different values of J and the chain sizes, L. The behaviour is equivalent to the previous model except that for sufficiently large values of J, the gap decreases sharply, closing completely (see figure 5b), causing the counter-diabatic terms to cause more error than the simple evolution itself, as we can see in panel (a) of the same figure.This is a consequence of the fact that our expression for the counter-diabatic term is FIG. 5. Complexity study for the ZZXZ model using the adiabatic algorithm.The phenomenology is essentially the same as for the TFI model.(a) Evolution of the fidelity obtained with shortcuts to adiabaticity (solid lines) and without them (dashed lines) for increasing time lengths of the full algorithm and L = 12.In this case we see that, for sufficiently large values of J, no applying shortcuts works better than applying them.This is explained by the gap closing much more abruptly than in the TFI model, as can be seen in (b).(c) Normalised complexity computed for different sizes with shortcuts (solid lines) and without shortcuts (dashed).As the gap closes, more gates are needed to achieve the fidelity threshold (0.9 in this case).(d) Complexity dependence on the chain size (root squared) and the preparation time.Each point corresponds to the preparation time for each |J| in (c).
not exact, but a first-order approximation of a general expression [cf Eq. ( 48)].This scenario serves to illustrate that the design of shortcuts can be tricky.Solutions to this problem could be to go to higher orders in the development of the CD term or to explore other λ(t) functions as in [63,64], where geometric arguments are used to obtain the optimal way to vary the time dependent parameters in the adiabatic evolution.
We can repeat the study in Eq. ( 52) (for the dependencies of the complexity on the system size, preparation time and model parameters) for this model.As one can see in Eq. ( 56), the only change is in the integral.
Despite obtaining an analogous analytical expression, we do not recover a linear dependence of C N / √ L on T in practice, as shown in Fig. 5d.Two effects contribute to the deviation.First, since we are now reaching larger values of J, the integral becomes more relevant than in the previous model.Thus, the complexity depends on J in two ways: through the preparation time T and through the integral.As a result, the scaling of the complexity with T is now superlinear, as a change in T is influenced by an underlying change in J which also contributes to the complexity through the integral.The second effect is numerical.The preparation time for this model is longer than for the TFI model, as shown in the other panels of Fig. 5.As a consequence, our discretization in preparation time is now insufficient to resolve the true dependence of C/ √ L on T .This effect manifests as bunching: we observe several different values of complexity for the same value of T .This is because the complexity is changing with J through the integral but the preparation time is stuck at the closest larger value.This effect should disappear with finer numerics.
The inclusion of shortcuts does not provide any significant advantage in terms of complexity reduction.This is because we have constrained these shortcuts to be as local as possible, in our case l = 1 in (49), introducing two body interactions at much.It is expected that by introducing long-range terms in (45) the complexity decreases as the system approaches to the QPT.This can be compared to the previous section II, where there was no restriction to local operations, so both C F and C N remained finite despite crossing the QPT.Other paths investigated in this work are sent to App.B.

B. Circuit Complexity in VQEs
VQEs, introduced in [36], use the fact that any quantum state can be written in terms of a unitary operation as where U ( ⃗ θ) is a parameterized unitary that transforms the initial state into the desired wave function |ϕ( ⃗ θ)⟩.This unitary can be implemented in a quantum circuit as a set of quantum gates.The expectation value of the Hamiltonian where we encode our problem (H) results The optimization process consists on minimizing the average energy of the parameterized state: The algorithm can be divided into three different stages.First, we need to choose the trial wave function (see Eq.( 57)).Choosing the unitary U ( ⃗ θ) is equivalent to constructing the quantum circuit that transforms the initial state into the parameterized wave function.The circuit used to achieve |ϕ( ⃗ θ)⟩ is called the ansatz and can be represented as, Choosing an appropriate ansatz is crucial for the optimization process.This choice depends completely on the model we are simulating and the set of gates available.We will dig into our choice of unitary below.The next step is constructing the Hamiltonian of the problem.Since this Hamiltonian is going to be evaluated later, Eq. ( 58), it must be written in terms of Pauli strings {I, σ x , σ y , σ z } ⊗L .Pauli operators are related to spin observables, which are suitable for direct measurement in quantum devices [65].With the Hamiltonian and the wave function defined, we can measure the energy of the state, which is the cost function.To compute this cost function, the expectation values of the Pauli observables are measured determining the value of the energy.Since the technique uses quantum and classical processors, VQEs are cast as hybrid algorithms.Our results are numerical and our Python code simply computes the product of the matrices U ( ⃗ θ) † HU ( ⃗ θ) previously defined and then projects onto the zero state obtaining ⟨0|U ( ⃗ θ) † HU ( ⃗ θ)|0⟩.We will not discuss its measurement overhead.Here, we are interested in the circuit complexity for reaching the desired ground state.The final step is to minimize this cost function through the variation of the parameters θ in the wave function.At the end of each iteration we obtain the value of the energy (59).Then, a classical optimizer determines the best direction of variation of the parameter vector ⃗ θ to minimize this value.We use as many iterations as needed until we converge to a final solution for the coordinates of the parameter vector.Ideally, this solution is the absolute minimum in the space of parameters.Still, obtaining this minimum is not an easy task.The optimizer can get trapped in local minima which will imply serious limitations in the minimization process.This problem and others have been previously discussed in the literature [65,66] and are out of scope for this work.Summarizing, we assume a given ansatz, the set of available gates in U ( ⃗ θ) in ( 57) and the hybrid algorithm finds the optimal solution.C N counts the number of gates, and once the VQE circuit is chosen, it can be done systematically.

Local VQE ansatz
We focus on a fixed geometry that is suitable for one-dimensional systems with single and two-qubit gates, besides the two-qubit gates act only on contiguous qubits.This ansatz can be interpreted as a Trotter approximation of continuous evolution by a local 1D Hamiltonian [67].In this case, we can separate the terms of the Hamiltonian that act on even and odd links and obtain two sets, each made of mutually commuting gates.In particular, the circuit is given by i.e. it consists of fundamental blocks (or layers) (separated by dashed lines above).Each layer is made out of single-qubit rotations R y (θ) and control-Z gates (CZ).At the end of the circuit, we add a final column of rotations (R y ).
For computing C N we rewrite the CZ gates in terms of Pauli operators, count the gates and use equation ( 6).This is a routine process that we send to Appendix A. Here, we just give the final result:

VQE complexity through QPTs
As before, we focus on Ising and ZZXZ models, Eqs. ( 25) and (53).In Figure 6 we summarize our results for the Ising Hamiltonian.In panels a-c) we plot the complexity using the local VQE ansatz to obtain the ground state at a given J for different chain sizes.We see that C N grows when the ground state approaches the QPT, that in this case is given by J c ∼ = 1 5 .In fact, close enough to the transition, the VQE cannot reach an acceptable ground state for a maximum depth of 8 (in our simulations).This can be checked in panel d) where the fidelity between the state obtained within the VQE algorithm and the exact ground state falls below 0.9 in the gray region of panel a) for L = 12.
It is possible to understand why the VQE fails around the QPT (gray zones in Figure 6).The local ansatz in Eq. ( 61) is a trotter-like decomposition of a time-dependent two-body interaction Hamiltonian [see Appendix A].As a result, the local ansatz can generate states whose correlation length, ξ, grows linearly with the number of layers (circuit depth), as described in detail in [67], which is rooted in Lieb-Robinson bounds.Thus, we have: On the other hand, for finite systems close to the phase transition, the correlation length saturates, ξ → L. Therefore, the lower bound for the circuit depth is L. Additionally, from the calculated VQE-complexity, Eq. ( 62), it scales as, Consequently, close to the QPT, This bound applies to the local ansatz and is specific to 1D systems.Our numerical simulations become intractable for d ≥ 8, so for reasonable sizes the VQE ansatz cannot approximate the ground state, explaining the gray zones.The grey zone indicates that the VQE does not converge for points inside that region in a reasonable number of layers to the fidelity threshold (0.9).(d) Fidelity obtained for different numbers of layers for points inside the grey box in (a) and in its vicinity for L = 12.For those points whose fidelity is above the threshold (0.9) only the best result has been plotted, for clarity's sake.
One final comment.Non-local ansatzes are expected to reduce the circuit depth, as discussed in [68].Further research on the complexity in this context would be interesting.It is worth noting that in the original proposal, long-range interactions are penalized in the functional F (τ ), see Eq. ( 6).Now, let us explain what happens when the target state is far away from the quantum phase transition (QPT).Before delving into the specifics, it is important to note that in finite simulations, deep within the ferromagnetic phase, the Z 2 symmetry is not broken.Therefore, the ground state manifold found by exact diagonalization is spanned by the states 1 2 (|0, ..., 0⟩ ± |1, ..., 1⟩).However, through energy minimization, the VQE reaches one of the fully polarized states, either |0, ..., 0⟩ or |1, ..., 1⟩, given that they are degenerate with the symmetric ground state.Our convergence criterion is based on reaching a fidelity of 0.9 between the state generated by the VQE and the result of exact diagonalization.Due to the discrepancy in the ground states obtained by both methods in the ferromagnetic phase, the fidelity is capped at 0.5, and the convergence criterion is never satisfied.To address this discrepancy and align with the physics of actual QPTs in the thermodynamic limit, where the symmetry is (spontaneously) broken, we decide to introduce a small bias, ϵ σ z i , in (25).That being said, it is remarkable that when the target state is far from the critical point, the complexity decreases, even though the target point and the reference state may belong to different phases.This is because, unlike the adiabatic algorithm, the VQE does not necessarily need to traverse states in the transition region to transition from |ψ R ⟩ to |ψ T ⟩; it can bypass criticality.This phenomenon is easy to understand in the Ising model since, in the paramagnetic phase, the ground state is approximately given by |+, ..., +⟩ (|+⟩ = 1/ √ 2(|0⟩ + |1⟩)), as shown in Eq. ( 25).This state is straightforward to prepare since it can be obtained through single-qubit rotations from the reference state |ψ R ⟩ = |0, ..., 0⟩.On the other hand, in the ferromagnetic phase, with the added bias, the ground state is either |0, ..., 0⟩ or |1, ..., 1⟩, which can also be easily obtained within the VQE.
We now consider the ZZXZ model, Hamiltonian (53) 6 .Here, we are not going to explicitly break the symmetry in order to discuss the scenario in which the symmetric ground state is sought.In the ZZXZ model, the QPT separates paramagnetic (PM) and antiferromagnetic (AFM) phases.In the PM phase, the behavior is analogous to the Ising model, Cf.Figs. 6 and 7. Deep in the AFM phase, the ground state manifold is spanned by the states |ψ AFM ⟩ ∼ = 1 √ 2 (|1, 0, 1, 0, ...⟩ ± |0, 1, 0, 1, ...⟩).Following the previous discussion, the VQE does not reach the symmetric ground state.Therefore, we see that C N grows as it approaches the phase transition (with our parameters J c ≲ 1, see Fig. 4) but does not decrease afterwards.At some point near criticality, the VQE cannot produce a ground state with a fidelity larger than 0.9, see panel d) for the L = 12 case, similar to the TFI model scenario.Here, however, the state remains difficult for the VQE after the near-transition region is surpassed.This is further confirmed in figure 8. There, we can see that although the total magnetization is well reproduced by the VQE (also the energy, in panel c), once we enter the antiferromagnetic phase the VQE generates either |1, 0, 1, 0, ...⟩ or |0, 1, 0, 1, ...⟩, as can be seen The grey zone, as in the TFI model, indicates that the algorithm fails to achieve fidelity over 0.9 for points within that region.(d) The fidelity behaviour with the depth of the ansatz shows that, again, once the QPT is crossed the algorithm cannot reach fidelities over 0.9.In contrast to the TFI model, here we don't recover high fidelity once we are fully in the antiferromagnetic phase, reaching a maximum value of 0.5 for the highest values of J (L = 12).In blue it is computed the fidelity as the overlap between the state generated by the VQE and the exact ground state; in red it is computed as the projection onto the subspace generated by the ground state and the first excited state.(c) Energy accuracy obtained for the same configurations displayed in the other panels computed as 1 − E rel , being E rel the relative error between the energy obtained from VQE and the exact value.
by computing the magnetization per site, which should be close to 1/2 in the exact ground state.However, the VQE gives 0 (1) for the even (odd) sites.To conclude our characterization, we see that all this is consistent with obtaining a F = 0.5, as well as a F ∼ = 1 if we compare the state generated by the VQE with the projection onto the subspace generated by the ground state and the first excited state.

V. DISCUSSION
Knowing in advance how much a computation will cost, even if only approximately, is of great help.Unfortunately, this estimation can pose a great challenge.Computer science has traditionally categorized problems into different complexity classes, allowing one to know whether a given problem is tractable on a classical computer.For a quantum computer, we can ask a similar question to know if the task we want to tackle is going to be feasible with the architecture we have at hand.For this purpose, the concept of circuit complexity was invented.Again, knowing the complexity of each task in any architecture seems too general to be able to give a concrete answer.On the other hand, we can shed some light on generic situations where some kind of general statement can be made.This is the idea that motivated us to write this manuscript.We have studied the situation in which a critical region is crossed in the process of preparing a state.
Our work has shown that, regardless of the type of complexity one chooses, and for diverse models, it appears that complexity grows if the algorithm visits states near a phase transition.We have further proven that this is a characteristic trait of typical algorithms for state preparation such as VQE and adiabatic evolution.The degree of divergence does depend on the definition of complexity used and on the allowed gates.In the case of local ansätze or evolutions, C tends to diverge as the system size grows, specifically as ∼ L 3/2 .Importantly, we have shown that VQEs, to the extent that they can go "directly" from the reference to the target state, can potentially avoid the divergence in complexity even if the reference and target states lie in different classes.Whether this is possible depends on the model, as it is determined by the degree of entanglement of the target and reference states.In the case of adiabatic algorithms, we have shown that complexity is bounded by the T , wich represents the total time of the algorithm.Therefore, for keeping the complexity down seems to be a matter of allowing non-local gates in the evolution, to fully exploit shortcuts to adiabaticity.This is supported analytically in Sec.III.Here, the Ising critical point is traversed along a restricted path of states of the form (26). Despite this restriction, these states are sufficiently non-local for C N to remain finite.
The impact of our work on the preparation of states in a quantum machine seems straightforward.What our results mean in the field of holography is another matter.Unfortunately, we do not have the knowledge to anticipate anything, but it would be interesting to think in this direction.Other ideas not discussed here would be the use of other types of complexity such as Krylov [23,[69][70][71][72] or mixed states and their behavior in thermal phase transitions.We leave this for future work.
Note Added in Proof.-Whilewe were finishing writing this manuscript, the paper [73], which discusses the importance of local and non-local gates in the computation of complexity, appeared in the arXiv.12. Relation between the Nielsen complexity in the adiabatic preparation and the system gap for the ZZXZ model.

)
FIG. 1. Study of the complexity for the Transverse Field Ising model.(a) Complexity for different sizes of the chain computed using the Fubini-Study metric.The discretization in J used is δJ = 1e−3.(b) Derivative of the Fubini-Study complexity for different L, δJ = 2e−3.(c) Finite size scaling of the maximum in the derivative of the Fubini-Study complexity.See that this maximum diverges polynomially with the size of the chain.(d) Study of the Nielsen complexity, δJ = 3e−4.(e) Derivative of the Nielsen complexity for different L, δJ = 3e−4.(f) Finite size scaling of the maximum in the derivative of the Nielsen complexity.See that this maximum diverges logarithmically.

FIG. 2 .
FIG.2.Fubini-Study complexity (a) and its derivative with respect to λ (b) across the phase transition of the Dicke model as a function of the system size (numerical results) and in the thermodynamic limit (analytical results).Plots on the right showcase the fits of λmax(L) (c) and (∂CFS)max(L) (d) (extracted from center plot) to their respective finite size scaling laws.All results are at resonance ωc = ωs = 1 and with a discretization of dλ = 10 −3 .Numerical results were obtained with a cutoff for bosonic excitations of Nexc = 30 .

FIG. 4 .
FIG. 4. Phase diagram of the ZZXZ model for zero temperature.The black dotted line signals the critical region between phases for different ratios of the fields (hx, hz) to the magnitude of the exchange interaction (J).The coloured lines depict the path followed for the adiabatic algorithm (red) and the values computed in the VQE (blue) [Cf.Sec.IV B].

70 FIG. 6 .
FIG.6.Transverse Field Ising model with bias, ϵ = 0.001.(a-c) Complexity as a function of J for sizes L = 10, 12, 14.The grey zone indicates that the VQE does not converge for points inside that region in a reasonable number of layers to the fidelity threshold (0.9).(d) Fidelity obtained for different numbers of layers for points inside the grey box in (a) and in its vicinity for L = 12.For those points whose fidelity is above the threshold (0.9) only the best result has been plotted, for clarity's sake.

50 FIG. 7 .
FIG.7.ZZXZ Ising model.(a-c) Complexity as a function of J for sizes L = 10, 12, 14.The grey zone, as in the TFI model, indicates that the algorithm fails to achieve fidelity over 0.9 for points within that region.(d) The fidelity behaviour with the depth of the ansatz shows that, again, once the QPT is crossed the algorithm cannot reach fidelities over 0.9.In contrast to the TFI model, here we don't recover high fidelity once we are fully in the antiferromagnetic phase, reaching a maximum value of 0.5 for the highest values of J (L = 12).

FIG. 8 .
FIG.8.VQE state characterization in the ZZXZ model for L = 12.(a) Magnetization of the spin chain as a function of J obtained from the states generated by VQE.The solid black line represents the total magnetization per site that the spin chain should have (obtained via exact diagonalization) whereas the dashed black line sets the magnetization per site in even/odd sites.(b) Evolution of the best fidelity obtained as a function of J.In blue it is computed the fidelity as the overlap between the state generated by the VQE and the exact ground state; in red it is computed as the projection onto the subspace generated by the ground state and the first excited state.(c) Energy accuracy obtained for the same configurations displayed in the other panels computed as 1 − E rel , being E rel the relative error between the energy obtained from VQE and the exact value.

FIG. 9 . 50 FIG. 11 .
FIG.9.Adiabatic evolution for the TFI model by switching on the transverse field instead of the spin-spin interaction.(a) Evolution of the fidelity obtained with shortcuts to adiabaticity (solid lines) and without them (dashed lines) for increasing time lengths of the full algorithm and L = 12.We see a clear difference with the plot in the main text, where the field is fixed and we vary the interaction, J.The gap closes much earlier for small field values (b), making the algorithm need much longer times to achieve high fidelity.