Circuit complexity in interacting QFTs and RG flows

We consider circuit complexity in certain interacting scalar quantum field theories, mainly focusing on the $\phi^4$ theory. We work out the circuit complexity for evolving from a nearly Gaussian unentangled reference state to the entangled ground state of the theory. Our approach uses Nielsen's geometric method, which translates into working out the geodesic equation arising from a certain cost functional. We present a general method, making use of integral transforms, to do the required lattice sums analytically and give explicit expressions for the $d=2,3$ cases. Our method enables a study of circuit complexity in the epsilon expansion for the Wilson-Fisher fixed point. We find that with increasing dimensionality the circuit depth increases in the presence of the $\phi^4$ interaction eventually causing the perturbative calculation to breakdown. We discuss how circuit complexity relates with the renormalization group.


Introduction
In the context of quantum information theory, circuit complexity is the number of unitary operations needed to perform a desired task in a quantum circuit [1][2][3][4][5][6][7][8][9][10]. Of course it is desirable that the number of steps is minimum to have the most efficient implementation of a quantum algorithm. In particular it is important to quantify this so that one can have meaningful comparisons with classical algorithms. The question of circuity complexity in the context of quantum field theories is still relatively novel with very few results. In [11][12][13], it was shown that the non-perturbative calculation of n-particle scattering amplitude in a scalar φ 4 theory on a quantum computer would have an exponential advantage over known algorithms which can be implemented on a classical computer which uses perturbative Feynman diagram techniques to perform such a calculation. This is quite remarkable and the question naturally arises how a quantum computer would compute other interesting quantities that are calculated by conventional means. Especially, is there a connection between renormalization group flows and circuit complexity?
Recently, observations due to Susskind and collaborators based on observations related to thermalization during black hole formation in holography have spurred activity in computing circuit complexity in quantum field theories [14][15][16][17][18][19][20][21][22][23][24][25][26]. In holography, it was observed that while the entanglement entropy asymptotes to a constant with time as the black hole thermalizes, the size of the Einstein-Rosen bridge in the context of the eternal AdS black hole keeps increasing. It was proposed that the analogous quantity that keeps increasing after thermalization is complexity. Two interesting proposals were given in the context of AdS/CFT [27][28][29][30][31][32]. The first one is the volume of a maximal codimension-one bulk surface extending to the boundary of AdS space time (complexity = volume). This particular slice is chosen in a way such that it asymptotes to a specific time slice on which the boundary state resides. The second one is to consider the so called Wheeler-DeWitt (WDW) patch which is basically the domain of dependence of a bulk Cauchy surface anchored at a specific time. Both these two objects have the potential to probe physics behind the horizon and both of them evolve with time even after the thermal equilibrium has been reached. These two proposals have been subjected to a host of interesting checks .
In [66][67][68] a geometric approach for circuit complexity was put forward which was studied in great detail in the context of free scalar field theories in [14]. It was proposed that to find the minimum complexity, one writes down a suitable cost function in the space of unitaries, which works out to be in general a Finsler space, and then minimizes it. The distance from the reference state (|ψ R ) to the target state (|ψ T ) is then the geodesic distance. To elaborate, we want where U (s) is a unitary operator. Here s parametrizes a path in the Hilbert space with the boundary conditions such that the reference state is located at s = 0 and the target state is located at s = 1. This is just a matter of convenience . We can redefine this parameter. Now every unitary operator can be written as follows, The basic notion of complexity (or more suitably in this case "circuit depth") is to provide a measure to count the number of "elementary" gates which can be combined to form the required U (s). In general, this is not a unique procedure and is difficult to accomplish. Following [14,[66][67][68] we try to find the shortest path or the geodesic connecting the reference and target state. To do so, we first define a suitable "cost function" F (U,U ) and the complexity is then defined as, We then minimize this cost function which then gives the geodesic connecting the two states. Then evaluating D(U ) on this geodesic, we obtain a measure for the complexity. Now there are various possible choices for the "cost function". But there are some desirable properties [66][67][68] that these cost functions should satisfy. They should be continuous, positive definite, homogeneous and satisfy the triangle inequality. These properties help us to identify these functions as legitimate functions measuring the distance between two points on the underlying manifold. In addition to these, if these functions are infinitely differentiable, then (1.3) gives the distance between two points on a Finsler manifold. Keeping these properties in mind we can choose various possible functions. We quote the ones advocated in [14,19,24]  (1.4) Here p I are some weights which at this moment are arbitrary. We observe that F κ=1 is directly related to the number of gates that has been used to achieve the target state (at least upto a certain tolerance |||ψ T − U |ψ R || < , where is a small and adjustable parameter). On the other hand F 2 with p I = 1 for all I is basically a distance function on a given manifold. These calculations have been generalized for free fermions in [17,19]. Interesting connections were found with holographic proposals in spite of the fact that these were free theories. It was suggested in [69] that tensor network approach maybe useful to derive qualitative features of holography, especially to understand the notion of emergent geometry from the field theory. One such useful tool is cMERA ("Continuous Multiscale Entanglement Renormalization Ansatz") [70] which provide us with several interesting features of holography [71][72][73]. In [74], the authors have computed complexity for the cMERA circuit for ground state of free scalar theory. Basically for the free scalar theory the ground state wavefunction furnished by cMERA can be parametrized as an SU (1, 1) coherent state. Then complexity was computed by first defining a Fubini-Study metric for SU (1, 1) manifold and computing the length of the geodesic with a suitable boundary condition. Further motivated by the tensor network representation of the partition function an alternative method of computing complexity for conformal field theory has been proposed in [75,76] based on an "optimization" procedure which basically determines how to represent most efficiently a partition function for conformal field theories. Then the complexity can be computed by evaluating simply the Liouville action [75][76][77][78]. This method has been generalized to include perturbations for 2 dimensional spacetimes [79]. These methods give the same leading divergent term for the complexity [14].
Our goal in this paper is to consider interacting scalar field theories. A primary motivation is to try to establish a connection with the renormalization group perspective. Our progress in this front will be modest. We will compute circuit complexity in a variety of interesting interacting scalar QFTs to leading order in perturbation. We will generalize the approach of [14] to the interacting case. In the process of doing so, we will encounter several subtleties. For starters, since our approach will be based on the group GL(N, R) we will not be able to use a purely Gaussian reference state for reasons we will explain. Rather we will forced to start with a slightly non-Gaussian reference state. This will lead to interesting technical complications. Namely we will find that we will be forced to make the cost functional dependent on the perturbative coupling so that we can smoothly interpolate to the free theory.
Taking into account these complications, we will then turn to evaluating circuit complexity in various dimensions including fractional dimensions to make a connection with the Wilson-Fisher fixed point in the epsilon expansion. We find that while the free theory depends linearly on the spatial volume, the interacting part shows a fractional volume dependence. Next we find that as dimensionality increases, the circuit depth also increases in the presence of the interaction (for positive coupling). In the RG paradigm, we know that the Gaussian fixed point is stable for d > 4 while the Wilson-Fisher fixed point is stable for d < 4. From the perspective of a potential quantum computer, we then find that since the circuit depth corresponding to turning on a coupling increases (and in fact diverges worse than the free theory eventually) with increasing dimensionality, it will be harder to perform the corresponding computation. Therefore, it appears to us that circuit complexity can be used as a diagnostic to analyse RG flows. Eventually, one can then hope that there could exist a monotonicity property much like the c-theorems in quantum field theory [80].
The fact that there must be an interesting connection between the renormalization group and complexity is not unexpected. In the context of the kind of calculations that were initiated in [14], one can easily see this as follows. If we turn on a perturbative coupling, then clearly the complexity answer will get modified by this coupling. In the context of renormalized perturbation theory, the coupling is a function of the RG scale. As a result we can consider writing down a differential equation for the circuit complexity in terms of this scale. Since we are considering first order perturbation theory, the differential equation will relate circuit complexity to the beta function of the theory as well as the flow equation for the mass parameter. While this is unsurprising, what is important to know is the nature of this relation-in particular, is the perturbative approach to circuit complexity well defined in any dimension? Can we identify fixed points and establish if a fixed point is stable or unstable by considering circuit complexity? For instance intuitively we may expect that if a fixed point was stable, then moving away from this fixed point would increase the complexity, while if it was unstable then the reverse would happen. We will attempt to take some modest steps in these directions. A potentially useful spin-off in our investigation is that we will come up with a general analytic method to perform the required lattice sums. This will enable us to consider even fractional dimensions. This method is outlined in Appendix (C) and will be heavily used in the paper.
The plan of the paper is as follows: In Section (2) we discuss the λ φ 4 on lattice as coupled Harmonic oscillators and solve its ground state wavefunction. Then we detail the complexity calculation for the the two oscillator case by generalizing the arguments of [14]. We discuss in detail all the subtleties regarding our approach and the construction of the circuit. In Section (3) we generalize this for the arbitrary N oscillator case for arbitrary spacetime dimensions d and take the continuous limit of the expression for the complexity. Then we derive flow equations for the complexity and study its implications. In Section (4) we generalize all these for theories with arbitrary number (N ) of scalar fields. In Section (5) we discuss briefly computation of complexity of φ 4 theory using a different set of gates. Then we end with a summary and a list of some of interesting future problems. All the supplementary materials which we have deemed useful for the reader have been placed in the appendices.

2
Circuit complexity with φ 4 interaction-the 2-oscillator case In this paper we will consider a massive scalar field theory with aλ φ 4 interaction term. We will follow the notation in [14] to facilitate an easy comparison. The Hamiltonian for the theory is, where d is the spacetime dimensions. We will assume the couplingλ 1 so that we can work in a perturbative framework.λ is dimensionful having mass dimension 4 − d. Next we discretize this theory on a d − 1 dimensional lattice. After discretization the Hamiltonian takes the following form, n denotes the location of the points on the lattice. Then using, we arrive at the following 5 , We will focus on evaluating complexity for the ground state of this Hamiltonian. It is evident that this system is nothing but coupled anharmonic oscillators. For simplicity, first we focus on two coupled oscillators. The Hamiltonian is: We have also set M = 1 here to make the analysis of the harmonic oscillator case more convenientthis will not affect the final answers which have to be dimensionally correct. Eigenstates of this Hamiltonian can be easily solved in normal mode coordinates. , (2.6) We define where, H n (x) s are the Hermite polynomials with H 0 (x) = 1. The expression for the ground state eigenfunction to first order in λ can be written as Note that [λ] = 4 and [λ] = 4 − d in our notation. While introducing Ω seems redundant, it will facilitate a comparision with the coupled harmonic oscillator case as in [14] and we will continue using it. and (2.10) For later convenience we will use where, where it is understood that the expression in (2.11) can be trusted only upto linear order in λ.
We will also use an approximate wavefunction 6 Defining the fidelity as we find that F (ψ, ψ) = 3 λ 32

Circuit complexity
We would like to evaluate the circuit complexity for this ground state wavefunction starting from a reference state. This needs several assumptions on our part. We will need to specify the reference state and available gates. We begin by writing the wavefunction in the following form: N s is a normalization factor. s is a running parameter and parametrizes the space of circuits. For s = 1 this coincides with the target state (2.11) with N s=1 = (ω 0ω1 ) 1/4 √ π exp(a 0 ), while s = 0 will be the reference state. Here the idea is to write the exponent of the wavefunction as a matrix A(s) conjugated by a basis vector v. If the state is just a Gaussian state then v = {x 0 ,x 1 } only. But for our case the states are not Gaussian so we have to extend the definition of this basis vector v. Below we discuss this explicitly for our setup.
A desirable property of the reference state is that it should not contain any entanglement in the original coordinates. Also we have to keep in mind that in order to represent it in the form (2.16), we have to allow non-linear terms (x 2 1 and x 2 2 ) in the basis vector. Keeping in mind these two facts, we choose the following, Here λ 0 is some parameter that we will fix later on. It parametrizes a non-Gaussianity in the reference wavefunction. Now going to the normal coordinates we get, Now we rewrite this state in the form (2.16). We also want to make sure that the matrix A is non singular. We have to choose the following basis (this is one of the choices and the minimal one as far as we can see) In this basis, b is arbitrary. Now given the basis vector v in (2.19) we can easily verify that both for the reference and the target state we can write the exponents in the form v a A a b v b as advocated in (2.16). Also we can see that to write the exponents in this form it is absolutely necessary to include quadratic terms in the basis (2.19).
Next we want to make the determinant of this matrix A positive which needs 2 < b < 4. Now further to make our analysis simple we choose b = 3 to kill the off-diagonal components in the reference state. Now the target state in this basis will look like, The coefficients are given in (2.12). Also we restrictb such that the determinant of this matrix is positive definite. Asω 1 >ω 0 , we have a 3 > a 4 > 0. Also a 5 > 0. These conditions together with fact that the eigenvalues are positive lead to, The upper limit is always positive but depending on the values ofω 0 ,ω 0 , 1 − 1 6 (ω 0 +ω 0 ) 2 ω 0ω1 can be both positive or negative. But we have to make sure thatb a 5 is also positive to make the determinant of A(s = 1) positive. So we restrict ourselves to, Now starting from this reference state the target state can be achieved via unitary evolution The unitary takes the following form, As explained before we have to act the reference state by the set of the operators O I s in a particular sequence. Note that, the Y I (s) s depend on the path, i.e., on the particular sequence in which O I 's are acting on the reference state. Our target is to find the shortest possible path such that we will achieve minimum complexity. To do so, we try to attach a geometrical interpretation to this process [14,[66][67][68]. Now looking at the structure of the matrix A in (2.21) we can consider U (s) to be an element of GL(5, R) with positive determinant. Then we can write U (s) in the following way From this we get, Hence, Next the infinitesimal distance i.e., metric in the parameter space defined by Y I 's can be written as, (2.31) At this stage we have various choices for G IJ . First, for simplicity we set G IJ as 25 by 25 identity matrix. Now U (s) is an element of GL(5, R). To do further calculations we have to chose a suitable parametrization for U (s). In general any element (g) of GL(5, R) can be written in terms of product of an orthogonal matrix, a diagonal matrix and an upper triangular matrix (G = KAN ), which is known as Iwasawa decomposition. But after closely inspecting our target state (2.21), we can further simplify our choice of parametrization. We can easily infer that our target state is of a block diagonal form. This motivates us to parametrize U (s) in the following way, We have decomposed U (s) in terms of R 3 ×GL(2, R) blocks. We could have allowed for off-diagonal 7 Now U when acting on A via eq.(2.28) is not unitary! Rather what happens is that this U can be mapped to a unitary operator which then acts on the wavefunction as we will discuss. This is a notational issue in [14] that we will assume does not create too much confusion. elements. But it is evident that at the end of the evolution those terms have to vanish as the final state is also in the block diagonal form, so allowing those off-diagonal terms will only increase the path length and hence the complexity. Now GL(2, R) can be written as R×SL(2, R). We also note that the first block in (2.21) is diagonal. It is expected that in the normal mode coordinates the quadratic part of the target state (2.11) is always diagonal. As argued previously [14], this induces a flat metric. Keeping this is in mind we can set, x 1 = x 2 = 0. For the rest of the components we choose the following parametrization, x 0 = exp(y 3 (s)) cos(τ 3 (s))cosh(ρ 3 (s)),x 1 = exp(y 3 (s)) sin(τ 3 (s))cosh(ρ 3 (s)), Then the metric in (2.31) becomes,

(2.34)
A suitable functional for this case is (2.35) Using (2.34) we find (2.36) The problem of finding the shortest path will then be mapped to the problem of finding geodesic in GL(5, R) group manifold. So we have to find the geodesic coming from extremizing (2.36) and evaluate (2.36) on this geodesic. The boundary conditions from (2.28) are (2.38) Now we proceed to solve the geodesic equations. These equations can be solved by first finding the conserved charges. As our metric is nothing but the tensor product of R 4 and SL(2, R) matrices, we essentially use the results of [14]. Then we get after using (2.37) and (2.38), Also using similar arguments as [14] we can set θ 0 is just a constant independent of s. We have trivially, Collecting all these together finally our complexity functional evaluates to, We can evaluate U (s) = exp M s on this solution which will give us the optimal circuit. From there we can also identify the unitary operators which are the building blocks of this optimal circuit. by suitably choosing the parameters α, β, γ, δ, ζ, τ, κ and µ. Below we identify the M 's with the corresponding iO's and give the details in the Appendix (A). (2.45) Also we note that, it is evident from this analysis that we can recast the unitaryM by the linear combinations of only those operators which scale the coordinates. This is not surprising given the fact, that we have chosen a non Gaussian reference state which already contains quartic terms apart from the usual quadratic terms. Hence we can reproduce the target state by simply scaling all these terms. Now we end this section by making some comments.

Comments
• First of all we note that unlike for the free theory [14], we cannot just choose the reference state as a product of Gaussians. We have to allow x 4 0 and x 4 1 term in the reference state. If we try to set λ 0 = 0, our boundary conditions will be ill-defined as evident from (2.38), speciallyỹ(1). This is also tied to the fact that we are setting up the problem using the machinery of GL(N, R) group which required us to have the determinant of A(s) to be non zero. More specifically, given the choice of the basis (2.19), if we would chosen our reference state as Gaussian state then the matrix A(s = 0) will have zero determinant because of the absence of quartic terms in the wavefunction. Then the relation (2.28) will not hold as the conjugating by U preserves the determinant of both the reference and target matrix A. So this forces us to make λ 0 = 0. Note that preparation of non-Gaussian states is a hard problem and only partial results exist in the quantum information literature [82].
Since A(s = 1) and A(s = 0) are real symmetric matrices that commute, therefore they both can be expressed as diagonal matrices in a common basis with corresponding eigenvalues as diagonal entries. This diagonalization is brought about by the action of a unitary matrix R via a similarity transformation given by R.A(s = 0/1).R T . More details of this can be found in [14]. It can also be shown that, under this transformation the metric remains invariant and hence the complexity. Then the complexity is simply given in terms of the eigenvalues of the matrix A of the target state. This is expected because once both the reference and target matrices are diagonal in a common basis all that we need are scaling gates corresponding to this diagonal basis to take us from the reference to the target state. This is because any non diagonal entries introduced by the entangling gates corresponding to the new diagonal basis, have to be nullified using more entangling gates before reaching the target state. Thus the use of entangling gates would only add more to the number of gates required and hence would not correspond to minimal complexity. Having said that it can easily be seen that the number of scaling gates say corresponding to the i th diagonal element required to reach from the reference to the target state is 1 2 where λ i is the i th diagonal element of the target matrix and r i is the i th diagonal element of the reference matrix . Plugging this result into the complexity functional F 2 results in equation (2.46). These expressions can be easily generalized for arbitrary lattice size in arbitrary dimensions.
• Now two of the eigenvalues λ 1 , λ 2 coming from the quadratic part of the wavefunction, are of the form a 1,2 + λ b 12 . The other three eigenvalues λ 3,4,5 coming from the quartic part of the wavefunction are of the O(λ). Taking logarithm of these O(λ) eigenvalues gives log(λ) terms in (2.46). This will make the λ → 0 limit ill-defined. To avoid this problem we can choose λ 0 to be proportional to λ such that the λ dependence inside the logarithm cancels out. Further we would want that (2.46) can be expanded perturbatively in λ.
Since we are working in perturbation theory, we would expect to recover the free result by taking λ → 0. The λ → 0 limit of (2.43) is subtle. We would have expect that in λ → 0 limit we would recover the free theory result which is given by But this is not the case. There are actually two problems: -The third term in the expression for d in (2.43) does not have any counterpart in the free theory, but still it seems that in λ → 0 limit it doesn't vanish.
-Also in the fourth term in the expression for d in (2.43) we get log(λ). This make the The second problem can be easily cured by making λ 0 ∝ λ. This will also give that, in λ → 0 limit A(s = 0) will reduce to the product of Gaussians. The first problem is harder to solve. We could envisage having a smooth λ → 0 limit by choosing G IJ differently so that the appropriate components pertaining to the second block are proportional to λ leading to However, this makes the procedure of determining the complexity somewhat ad hoc and introduces a plethora of possible circuits. Also note that apart from the gates corresponding to the generators M 11 , M 22 , the rest of the gates are complicated and hence must be difficult to "manufacture." Hence it makes sense to consider a somewhat different problem where instead of the target state ψ 0,0 we will use the approximate target state ψ 0,0 given in (2.13). This will solve the following problem: Given the gates corresponding to M 11 , M 22 find the circuit complexity to go from a Gaussian reference state to the approximate ground state given in (2.13). This essentially means that in (2.49) we drop the terms proportional to A. This also essentially makes the geodesic problem identical to [14]. In what follows, when we compute complexity expressions from the first block or unambiguous block, this is what we are doing. We will then consider the terms proportional to A which means introducing more complicated gates as in eq.(2.45). This will need us to choose suitable penalty factors that keeps the calculation perturbative.
3 Generalizing to the N -oscillator case Now we generalize previous analysis for N coupled oscillators. The Hamiltonian takes the following form, We then do a discrete Fourier transformation to go to the normal mode coordinates.
The Hamiltonian in terms of these variable becomes, Here we have used the following facts, Also keeping in mind ω = m we define, The ground state wavefunction is (which is generalization of (2.11)), where, , , .
Again we should keep in mind that (3.6) is only valid upto O(λ). We will expand the complexity upto O(λ). (3.6) can be recast into the following form, As seen from equation (3.6), the ground state has the sum of following type of terms within the exponential : e ,x ixjxkxl . Given this, the choice of basis is not unique. But there is a minimal choice such that, the basis for The total number of terms in this basis are N + N (N −1) and for large N it grows as N 2 . For more discussions on the counting of the basis refer to the Appendix (B). Given this choice of basis the matrix A takes a block diagonal form Now again there are several comments are in order.
• The elements of the block A 1 constitutes −2× coefficient of terms of typex 2 a ,x ixj within the exponential in the target state. The form of the matrix A 1 is unique once we fix the target state (3.6). We will at times refer to the block A 1 as the 'unambiguous' block. The reason for this will be clear from the discussion below.
• The elements of the block A 2 constitutes −2× coefficient of terms of typẽ within the exponential in the target state. This block (A 2 ) is not uniquely fixed even after the target state is fixed to be (3.6). This is because elements corresponding tox 2 ax 2 b in the target state, can be either put in the diagonal entry of the matrix A corresponding to the basis elementx axb or as an off diagonal element corresponding to the basis entriesx 2 a andx 2 b . This sort of ambiguity also arises in the entries corresponding tox ixjxkxl . These are always off diagonal, but can be put in the target matrix in more than one way. In the most general target matrix such elements are distributed among all possible entries of the A 2 matrix, such that the sum of all entries add up to the -2× coefficient under consideration. Hence we will sometimes refer to the block A 2 as the 'ambiguous block'.
• Given the freedom in choice of A 2 we now make another choice for this rearrangement inside A 2 such that the determinant of A 2 is always nonzero i.e all the eigenvalues of this matrix is nonzero. This is absolutely necessary for our geodesic analysis for which the determinant of target and reference matrix has to be non zero.
Next we have to choose the reference state. Again as described in the previous section, a desirable property of the reference state is that it should not contain any entanglement in the original coordinates. So we generalize the reference state mentioned in (2.17) for arbitrary N in the following way.
This again can be recast in the following way (after going to the normal mode coordinate), where, Dimension of the identity matrix I k×k is same as that of the dimension of A (2) . It is evident that A s=0 also decomposes in terms of an unambiguous block and an ambiguous block. Again there will be all those ambiguities regarding the rearrangements of the elements inside this ambiguous block of the reference state as discussed previously in the context of the reference state. We also make a choice such that the determinant of A s=0 becomes non zero as in the N = 2 analysis. Armed with this reference state let us now discuss the complexity functional. We will use the following type of functional: (3.13) Then the complexity will be given, Note that this is more general than the functional (F 2 ) that we used in the previous section. As in the two oscillator case, the target (A s=1 ) and reference matrix (A s=0 ) commute with each other and they can be simultaneously diagonalized. So again the complexity will be given by the ratio of the eigenvalues of the reference and target state as stated below 8 , are the eigenvalues coming from the 'unambiguous' block and λ (2) j are the eigenvalues coming from the 'ambiguous' block. h j are some numbers. The number of λ (2) j is the same as the dimension of the unambiguous block and for large N it grows as N 2 . For κ = 1 the complexity functional (3.15) can be rewritten in the following way.
All the individual eigenvalues coming from the unambiguous blocks are positive both for the target and reference states. So the first ratio is automatically positive. But not all eigenvalues coming 8 For arbitrary N, F 2 will be given by, from the ambiguous blocks are positive. Some of them turns out to be negative. As we know there are several ambiguities inside this block. However, we cannot find any choice for the arbitrary parameters such that all the eigenvalues coming from this block will always be positive for all values of N. Similar problem persists for the ambiguous part of the reference block. But we can always check that for a given N, there are always some choices for the rearrangement inside this ambiguous block such that this ratio det(A (2)s=1 ) det(A (2)s=0 ) is always positive. Hence the C κ=1 is well defined. In fact one of the Schatten norms F p with p = 2 is also well defined for our case. For our case that is just √ C κ=1 . It seems that at this moment only these two measures are the only two measures that are well defined for our case. We will use C κ=1 only for all the subsequent discussions for simplicity. It will be an interesting problem to investigate the other measure but we will leave it for future investigations. Next we evaluate this complexity functional. As evident from (3.15) this complexity has two pieces.
where, C κ comes from the unambiguous piece and C κ comes from ambiguous piece including the extra penalty factor A. In the next section we evaluate these two pieces separately.

Evaluation of complexity functional and continuous limit
Before we proceed to compute the complexity functional we like to reinstate the factor of M that we have set to one from (2.5) onwards. With this factor reinstated, the Hamiltonian takes the following form, Now the overall factor of 1 M doesn't change the form of the ground state wavefunction. Only now there is a nontrivial factor of M 2 infront of the x 2 and x 4 part of the Hamiltonian. It will just scales various quantities, In light of this we have the general formula for λ (1) i k as shown below, , N : eveñ , N : odd i k goes from 0 to N − 1 for all k from 1 to d − 1. N denotes the number of lattice points in each of the spatial directions. Then the d − 1 dimensional spatial volume is given by L d−1 = (N δ) d−1 .
Using (3.19) we get the contribution to the complexity coming from the unambiguous block.
From now on we will focus on the κ = 1 case. Also, we will take the continuous limit i.e N → ∞ and δ → 0 such that N δ is finite. We are mostly interested in extracting the leading divergent and finite terms of C κ=1 . From now onwards we will write every expression in terms ofλ = 24 δ d λ and V = (N δ) d−1 , instead of λ, N and Ω = 1 δ . Also we quote here an useful relation which is the generalization of (3.4) for arbitrary d.
with each of the i k goes from 0 to N − 1 for all k.

d=2
Using the general method described in Appendix (C), we arrive at the total expression for the complexity given below: Here E(k) is the elliptic function of the second kind. This captures the exact δ, m dependence upto leading order inλ. We expand in terms of δ and that gives, where a 1 = 0, c 1 = 1/4, c 3 = − 1 96 , f 1 = 1 π , f 1,log = − 1 8 π , f 0 = 0.02. We note that, in comparison with the numerically computed free theory part-expression (E.12) in [14], the exact free theory result has all terms except the constant a 0 and the c 1 V m log m ω ref terms.

d=3
Using the results in Appendix (C), we find that for d = 3 : where, a 2 = 0.29, b 2 = 0.075, c 2 = −0.04 and f 1 = 0.16. Here we have retained only the leading 1 (mδ) term which arises from the integral (C.6) for the interacting sum. Note that there is a log(mδ) factor multiplying V m 2 contrasted to the fact that there is no such logarithmic term multiplying (V m) term in (3.23). Furthermore, unlike [14], the free theory result that our analysis gives is manifestly proportional to volume and only in the interacting part is there a breakdown in V scaling, since it is proportional to V 1/2 . Notice that compared to d = 2, the interaction part has an extra 1/δ dependence suggesting that as d increases, for fixedλ complexity will increase.

General d
For arbitrary d, as we argue in Appendix (C), we have (3.25) Using the results in Appendix (C) we find a 3 = 0.41, a 4 = 0.49, a 5 = 0.55. In principle we can also fix the b i , c i 's for general dimensions, but we will not attempt to do it here. Also we note from

d=4 −
The general method outlined in Appendix (C) enables us to extract information for any dimension, not just integer dimensions. For instance, using the integral form for the free theory sum, we can extract a d−1 as a function of d. This is shown in Figure (1). We can also evaluate the interacting sum as a function ofm ≡ mδ for various dimensions. This is illustrated in Figure (2). The plot is consistent with our findings above-for instance, it shows a divergence form = 0 for d < 4 while it approaches a fixed value for d > 4. Thus we can systematically study the epsilon expansion. At leading order sinceλ * = 16π 2 3 at the fixed point, the result is somewhat trivial since we can replace the f 0 , f 1 by the d = 4 values. However this procedure will prove to be useful when one computes the next order in perturbation.

Complexity in terms of renormalized parameters
Up to now we have given expressions involving the bare parameters. In order to extract physics, we will need to rewrite the expressions above in terms of the renormalized quantities following [83,84]. For the mass we have [84], m R is the renormalized mass and λ R is the renormalized coupling defined at zero momentum. A running renormalized coupling can also be defined at finite momentum µ and since we are interested only in leading order in the coupling, this amounts to simply replacing λ R by λ R (µ). Here, .
For d ≥ 3 we get, (3.29) Note that for d = 4 there is an extra log term. Also we tabulate values of C 0 and C 2 for various dimensions. From this we note that C 0 and C 2 are always positive for all d. At this order we simply havê λ 0 = λ R , where λ R is the renormalized coupling. Now given these two expressions, we get the following.
In both the expressions above, we have indicated the dominant terms in the small δ limit, keeping λ R , m R finite. For d ≥ 4, b 2 = 0. Also note that the genuine extra contributions that came from the first block which were proportional to f 0 , f 1 give a fractional dependence on the volume and are subleading in the large volume limit. Now it is obvious that since the λ R dependence at leading order in V is proportional to δ 4−d , perturbation theory will break down for d > 4. This is expected from the RG picture where d > 4 has a stable Gaussian fixed point but an unstable Wilson-Fisher fixed point. In order to isolate the effect of the interaction we can define (for d > 3) (3.32) In other words we are asking what is the change in complexity when we go from the free theory with mass parameter m R to the interacting theory with the same mass parameter. At the Wilson-Fisher fixed point in the epsilon expansion λ R * = 16π 2 3 and to leading order we use c 2 = −0.02, C 0 ≈ 0.15, the 4-dimensional value. This means that interaction has slightly increased the complexity at the fixed point compared to what happens at the Gaussian fixed point. For both fixed points m R = 0. The sign of c 2 will turn out to be important when discussing the consequences from the flow equation.

Comments about C
(2) κ and structure of penalty factor As discussed around equation (3.15) there is a second contribution coming from the eigenvalues of "ambiguous" block (λ (2) j ). As the name suggests and also from the earlier discussions there are many ambiguities. Nonetheless we will discuss the general structure of the contributions coming from this block to the complexity expression (C Here g(ω i ) has the dimension of m.
Since we are interested in the divergence of the leading term that contributes to complexity, we will assume a simpler form of the eigenvalues (based on above points), Now we have to keep in mind that,ω Now the Dimension of the block A 2 grows as N 2 for large N. But there are only N number ofω i . So given the form of λ (2) j , either in the equation (3.34) or in (3.33), there will be only N number of λ (2) j eigenvalues each with degeneracy N. This argument can be extended straightforwardly for any dimensions d.

Cost function and the penalty factor
We will use the cost function as mentioned in (3.15). We will concentrate on C (2) κ=1 . We will now make the following choice for the penalty factor 10 A as mentioned in (3.15): (

3.35)
A should be dimensionless and for the moment µ and ν are arbitrary but integers. But we will soon make choice based on some physical arguments. With this choice we get, In the original basis the reference matrix is chosen to be the one with eigenvalues λ 0ωref coming from the "ambiguous" block of the target state. This ensures that in the diagonal basis, the eigenvalues are of the form h i k λ 0ωref , where h i k are some constant numbers. We absorb this h i k into the b i k . Also, to make theλ → 0 limit well defined, we assume as discussed previously, λ 0 to be proportional to λ R so that the λ R dependence inside the logarithm will cancel out. Given this, and the form of A, we can easily see that λ R → 0 limit is well defined. Now this leaves us with two possibilities for choosing λ 0 .
Putting all these pieces together (for general dimensions) and reinstating all the necessary factors of δ we get, Now we use the fact these λ (2) j eigenvalues are degenerate. For each i k where k = 1, · · · d − 1, there are N of these eigenvalues with degeneracy N. Using this fact we get, Also going from (3.38) to (3.39) we have ignored the modulus assuming that, the individual terms are positive. Finally we get, where, λ 0 can be set to either one of the expressions given in (3.37). The first sum above yields a factor of V δ d−1 and it is being multiplied effectively by a factor log(δ).(upto some suitable factors to make it dimensionless inside the logarithm.) We had already dealt with the second sum in the previous section, as it has appeared in λ (1) i contribution. We will just use the expression given in (C.15) for that. Now if we focus on the first sum in (3.40), we see that after performing the sum it gives a N d−1 factor. The second sum in (3.40) always grows as N d−1 as evident from (C.15). We demand that the leading volume dependence coming from the first sum in (3.40) will utmost be of the same order as free theory i.e V. If assume this we can now choose the following, (3.41) Alternatively we can demand that the δ dependence can be δ 6−2d like the O(λ R ) contribution of C (1) κ=1 . Then we could have chosen, We can generalize this argument for the penalty factor order by order in higher order in λ R and we put some more details in the Appendix (D).

Flow equations
Here we would like to consider the flow equations for ∆C κ=1 defined via with an infinitesimal change in b, namely b = 1 + db which leads to λ R = λ R + dλ R . Then straightforwardly (for µ = 1, ν = −d for the second block, so that it is subleading for large V ), we find that up to linear order in λ R , These equations are similar to the flow equations for λ R namely dλ R db = (4 − d)λ R + O(λ 2 R ) (which has also been used in deriving the form above). For d > 4 we conclude that the flow for ∆C should be back towards ∆C = 0, which is what we get when we turn off the coupling. The opposite happens for d < 4. Quite pleasingly, this conclusion agrees with the RG picture (see e.g., [85]).

Comparison with Holography
Here we briefly make a qualitative comparison with the holographic results. From (3.25) it is evident that the free theory part is always proportional to spatial volume (V ). This is same as that of the expression for complexity coming from holographic calculation. But interestingly from (3.25) we note that the O(λ) terms are proportional to fractional power of the spatial volume (except for d = 2 where there is no V dependence in the O(λ) correction). Current holographic proposals will never produce such terms. If we consider the complexity equals volume conjecture, then using the results of [86] we can easily see that terms with fractional power of volume will never occur. Complexity equals action gives rise to a logarithmic enhancement of the volume divergence but still it doesn't give rise to fractional volume [40]. Also the occurrence of this fractional power is independent of the choice of the penalty factor (3.35) as this feature shows up in the Gaussian part itself irrespective of whether we are suppressing the contribution from the second block or not. Moreover from the holographic side [45,87], the corrections due to a relevant deformation using existing holographic proposals have been considered leading to V dependence in the complexity expression. Similar type of results for the complexity for d = 2 has been derived using path-integral approach by considering relevant and marginal deformation of conformal field theory in [79]. In light of that, we should emphasize here that, in the path-integral approach (also in the context of holography) one starts from a conformal field theory and then considers perturbation around that. Now we will have a conformal field theory only at the fixed points. In our setup we can only access these fixed point perturbatively, i.e., in the epsilon expansion, hence the comparison has to be made with great care.

Using circuit complexity for inferring fixed points
In this section we will briefly generalize the discussion above for C (1) N κ for a theory with N scalars. We want to study this theory to see if complexity arguments can be used to infer the existence of fixed points. We will consider the following Hamiltonian: Here theλ 1 term has O(N ) symmetry while theλ 2 term breaks this symmetry. There are interesting fixed points that this theory allows for, which have been studied in the epsilon expansion [83]. Next we solve the ground state wavefunction perturbatively inλ 1 andλ 2 as before. Then considering only the contribution from the first block and concentrating first onλ 1 correction, we get the following eigenvalues, where for technical simplification we are considering N odd only. Here a index run from 1 to N counting the number of components of vector field and each of these i k as usual runs from 0 to N − 1 for every a. Basically now we have N number of λ (1) i k eigenvalues compared to the N = 1 case. We also define λ 1 =λ 1 δ −d 24 . The first interaction term above i.e., the b independent term, same as the eigenvalues for the λ 1 φ 4 theory, comes from the λ 1 X a ( n) 4 term in the Hamiltonian, while the second set of terms involving sum over the b index arise from the 2 λ 1 X a ( n) 2 X b ( n) 2 terms in the Hamiltonian. Correspondingly, the complexity C (1) N κ is given by, For the case of N = 1 the b sum does not contribute and the a sum contributes one term, which leads to the same complexity as the λ φ 4 case. Then finally from (4.3) we get, where C f ree is the complexity for the free theory and C int is the interaction term sans the coupling constant for the λ φ 4 theory corresponding to N = 1. The additionalλ 2 φ 4 a 24 interaction piece of the Hamiltonian contributes the following additional terms to the eigenvalues (4.2).
, N : even , N : odd where λ 2 =λ 2 δ −d 24 . We then have the following expression for complexity of the theory given by (4.1),

Connection with RG
Now from (4.6), we see the following. First in the subsequent discussion will use renormalized coupling {λ R 1 , λ R 2 } instead of {λ 1 ,λ 2 }. When we express the bare mass contribution from the free part in terms of the renormalized mass we will get a term in the complexity that is proportional Furthermore, the coupling constants in the interaction part appear in the combination (λ R 1 (N +2)(N +1) 6 + N λ R 2 ). Since the weights of the λ R 1 , λ R 2 are different, it is clear that when N is large, there will be a larger contribution 11 from λ R 1 compared to λ R 2 . Then it is clear that there must be some intermediate N = N c when there is a crossover between which term dominates. In fact in the epsilon expansion the beta functions admit the following fixed points at leading order [83] ( The last corresponds to the Gaussian fixed point. The first corresponds to the Ising fixed point (decoupled theory of N scalars), the third to the Heisenberg/O(N ) fixed point. The second corresponds to the cubic anisotropic fixed point where both couplings are turned on. Note here that in the N large limit, λ R 1 * → 0. This is also to be expected from the circuit complexity expression where it is less favourable to turn λ R 1 on in this limit. From the complexity expression proportional to C f ree , for λ R 1 = λ R 2 3 (N +2) , and for C int , for λ R 1 = λ R 2

6N
(N +1)(N +2) in the plane of couplings, it is "equally favourable" to turn on both couplings. Hence it could be expected that one is allowed to perturbatively turn on both couplings in order to go away from the Gaussian, Ising and O(N ) fixed points. While this stops short of proving the existence of the cubic anistotropic fixed point using the complexity perspective, it hints at its existence (of course for d < 4 since perturbative circuit complexity breaks down for d > 4 even here).

An alternative construction
Here we briefly discuss an alternative construction of the circuit, and we will also estimate the complexity 12 . The advantage of the considerations so far is that we were able to give a meaningful geometrical interpretation to the complexity calculation. We were able to provide an optimal circuit (minimizing C κ=1 functional) which represents the wavefunction upto O(λ). A by product of our calculation was that we were able to minimize the circuit depth and get an answer for the complexity. The downside of this calculation was that there are several ambiguities that enter and we had to make a choice for the penalty factors to make theλ → 0 limit well defined. Now we want to present an alternate calculation which will help us to avoid the problems associated with the non-uniqueness of the O(λ) block. The flip-side of this alternative approach is that we do not know how to give it a geometrical interpretation, leaving this issue for future work. Below we discuss our construction for N = 2 case for d = 2.
We start with Gaussian reference state by setting λ 0 = 0 in (2.17). We consider first the scaling operatorsÕ For this discussion, imagine λ R1 ∼ λ R2 to be of similar magnitude. 12 This is also being considered in [88]. and act on the reference state with them. This gives Then we can set, Here is a small parameter and a 1 , a 2 are defined in (2.12). Then we consider the following operators,Õ 3 =λ These two operators will be responsible for fixing the coefficients ofx 4 , so their coefficient will be of the O(λ 0 ). We will use this fact to approximate our calculation. We use the following unitary, Similarly, forÕ 4 we have the corresponding coefficient α 4 . p 0 and x 0 are arbitrary constants with appropriate dimensions and the ratio p 0 x 0 is positive definite. Here, Now we consider the following two operators, Then we act on the wavefunction with these two operators: 1 + i α 5 p 0 x 0Õ 5 and 1 + i α 6 p 0 x 0Õ 6 and try to fix α 5 , α 6 . However we can determine only one of them at this moment. We determine α 5 without loss of any generality, but we assume here that α 6 is of the order (λ 0 ), although we cannot uniquely determine it at this stage.
• We note that the first line of (5.11) is same as that of C κ=1 in (3.20) for N = 2 and d = 2. This comes from the Gaussian part. This can be readily extended for arbitrary number of oscillator (N ) and the conclusion that we drawn regarding the RG flow in the previous section will remain same.
• We notice that this same set of operators has been used in the proposal of interacting cMERA for φ 4 theory recently [89,90]. Working perturbatively in O(λ) we can work out the algebra satisfied by these operators.
From these one can easily see that these operators can have at least 4 dimensional representations. However given a Gaussian reference state currently we are not able to utilize this algebra to geometrize this problem which would have helped us to achieve "minimal" circuit depth. Moreover, we can use these matrix representations coming from (5.13) and try to take a linear combination of them so that we can reproduceM given in (2.44) coming from the geodesic analysis. We can see that there are no such linear combinations that allow us to reproduceM given in (2.44).
• It would be nice from the tensor network point of view to test whether the cMERA proposal in [89,90] based on these gates achieve minimal complexity or not. Since the Gaussian part of the expression will be the same as our earlier analysis, we expect that the connection we have found with the RG perspective will continue to hold. The interesting question here is, if it will allow us a more refined understanding, in the sense of being able to detect fixed points.
• There is another interesting possibility. Instead of definingÕ 3 ,Õ 4 ,Õ 5 withλ as in (5.4) and (5.7) we could have pushed thisλ inside α 3 , α 4 and α 5 . Then we will have, (5.14) Then D(U ) will have smoothλ → 0 limit without having to choose penalty factors proportional to O(λ). But in this case unlike what is shown in (5.13), the algebra between these operators does not close. We leave these issues for future investigations.

Discussions
We have shown that there is a connection between circuit complexity and renormalization group flows. Namely, we found using the scaling equations for circuit complexity that, for d > 4 the perturbative calculation for circuit complexity breaks down. Put differently, it becomes unfavourable to turn on a φ 4 coupling for d > 4. This conforms with the RG picture, that it is the Gaussian fixed point that is the stable fixed point for d > 4. We also saw that it is possible to argue the existence of the cubic anisotropic fixed point for N scalars. There are several interesting questions to explore: • One of the interesting conclusions from our calculations was that turning on interactions lead to fractional dependence on the volume. While in the large volume limit, these would be subleading, it still warrants the question how holographic calculations would see such fractional dependence. This conclusion is unambiguous and is unchanged by what is happening to the "ambiguous" block.
• Our construction needed us to fix certain ambiguities. In order to make the circuit complexity calculation perturbative, we were forced to make some choices which translated into making the cost functional depend on the coupling. This on its own is not an essential drawback in what we have done, since even in holographic entanglement entropy calculations, the entropy functional depends on the coupling. Nevertheless, we do not have a first principle way of fixing these ambiguities we encountered.
• Our final circuit complexity formulas involved an integral transform. While this is true for the free theory and leading order in perturbation, it is easy to see that it must be true at higher orders in perturbation as well, since the lattice sums can be handled similarly. Especially after picking up the correct s-residue in eq.(C.6) we are left with a Laplace transform over the Schwinger parameter t. A question that arises is if there is a physical interpretation that can be given to this variable.
• In section 5, we showed that there is an alternative set of gates which could be used for circuit complexity calculation which will enable us to use a Gaussian reference state as opposed to a nearly Gaussian reference state we have used in this paper. It will be interesting to develop the geometric picture further for these set of gates. It will also be interesting to develop the ideas in this paper for other interacting theories involving fermions and gauge fields. The idea would be to use circuit complexity to say which theories are reasonable theories from a quantum computer's point of view-any circuit that is prohibitive in allowing a calculation would not be allowed for instance. After all, eventually we would like to understand what is so special about the standard model from this view point.

B
Counting of dimension of the target matrix and the basis elements for arbitrary N : In this appendix, we discuss in detail the counting of the number of components in the basis. From (3.6) and (3.7), it is evident that we have to allow for the quadratic term in the basis v. In general the basis will of the following form, v = x 0 · · ·x N ,x 2 0 · · ·x 2 N ,x ixj , (B.1) where i = j and both i and j run from 0 to N − 1. Now we give an explicit counting for the number of terms in the basis. For the odd N case the minimal basis required to span the target state will always be: Here we have simply denotedx i 's by i andx ixj by i, j. The total dimension (D) of the matrix for the N odd case: For the even N case the basis will include terms: {0, 1, 2, 3, ...., N − 1} → "unambiguous" part For the "ambiguous" part: The set of basis will include the term ab if either one of the conditions given below are met: The number of terms that satisfy these conditions equals N (N +2) 4 , if N is divisible by 4 and otherwise Therefore the dimension (D) of the target matrix A T in the N even case is: To summarize, the number of basis elements for the linear part of the basis is N and the number of quadratic elements in the basis grows as N 2 . As far as we can see, this is the minimal way of extending the basis to include the quadratic terms and this is sufficient to produce all the terms in the wavefunction upto O(λ).

C A general method to do the lattice sums
We are dealing with lattice (multiple) sums involving powers and logarithms of m 2 +4Ω 2 d−1 k=1 sin 2 πi k N . Similar sums occur regularly in lattice field theory. The technique we will focus on is similar to what can be found in [83], except that we will generalize this approach to handle any power and logarithms. This is best handled by making use of an integral transform and writing where for evaluating log(y + x) we choose f (s) = 1/s 2 and for (y + x) a we choose f (s) = 1/(s − a). The contour C encircles the pole. Then we make use of the Schwinger parametrization to write For our sums x = 4Ω 2 d−1 k=1 sin 2 πi k N and y = m 2 . Using 14 where (a) b = Γ(a + b)/Γ(a) is the Pochhammer symbol, we find where I 0 is a modified Bessel function of the first kind. Notice that in the form of the RHS, we can choose d to be fractional as well! Thus we now have a way to compute complexity in the epsilon expansion, say d = 4 − . So in all we now have We now move to a dimensionless variable using t → t/Ω 2 , which gives where in the last line E[k] is the elliptic function of the 2nd kind. For d = 3, explicit expressions can be found as well although these are quite lengthy to quote here. Using these expressions we can expand around the lattice cutoff δ = 1/Ω. For the logarithmic sum we find Expanding this around small δ we get as the leading term, 14) The higher dimensional cases can presumably be done analogously yielding analytic results involving hypergeometric functions of the kind d F d−1 . However, we will not attempt to give a general expression here. Using the d = 2, 3 examples we will resort to general expansions determining the coefficients numerically from the integral expressions. Given the analytical result for d = 2, 3 it is not hard to get the general result for arbitrary d which takes the following form, For the O(λ) term we proceed as follows.
(C. 16) In the second step we have made the change of variables to the dimensionless u = t Ω 2 ; (where, Ω = 1/δ). We also notice that the the integral above is nothing but the Laplace transformation of the function f (u) = √ u (e −2u I 0 (2u)) d−1 .
We will define F{f (u)} := ∞ 0 due −(mδ) 2 u √ u (e −2u I 0 (2u)) d−1 . Now if we set (mδ) = 0 in F{f (u)}, then we observe that: F{f (u)} diverges for d < 4, indicating a 1/(mδ) n>0 behavior. A series expansion of the exact d = 2, 3 results suggest that the dependence on mδ is (mδ) d−4 . Also from the plots shown in (...) F{f (u)} vs (mδ), for d > 4 it can be easily seen that it has a finite value at (mδ) = 0. These plots indicate that the fit for any d should behave as follows, The methods developed in this appendix will prove useful for higher order perturbations, as well as other analogous sums.

D Structure of the Complexity beyond O(λ)
Here, we will discuss higher order in λ corrections to complexity focusing on the structure of the second block. Let us first consider O(λ 2 ) correction to the target matrix. In addition to the correction to the first block, this introduces O(λ 2 ) corrections to the second block and a third block at O(λ 2 ) comprising elements corresponding to coefficients of terms in the O(λ 2 ) block in the target state and O(λ 2 ) cross terms between the second and the third block due to coefficients ofx a 6 ,x 2 dx 4 e , (x ixjxk ) 2 terms in the target state. In each block we only consider the leading order in λ. Henceforth, we focus on the O(λ 2 ) block (A 3 ) and neglect these crossterms.
In d = 2 dimensions, eigenvalues coming from this block (A 3 ) takes the following general form: . (D.1) k runs over the dimensions of A 3 and f k 's are the quadratic polynomials of the frequenciesω a where a, b = 0, · · · N − 1. . Since we are only interested in the divergence structure of the leading contribution to complexity hence we ignore certain terms and consider the more simpler form for the eigenvalues given by, In the general O(λ g ) case we have g + 1 blocks in the target matrix with the (g + 1) th block containing coefficient of terms whose power ofx i sums up to 4g. At this order (order of λ = g) these terms are of purely O(λ g ). The other diagonal blocks-say the j th block, j ∈ {0, 1, ..., g} have leading terms of O(λ (j−1)) with subleading terms upto O(λ g ). The cross terms between the j th block and the (j + 1) th block will contain terms of at least O(λ j ). In each diagonal block we consider only the leading term. With this the eigenvalues of the general j th block is of the general form, Conditions on α 1 ,α 2 ,..,α k ,..,α (j−1) 1 f k (ω α 1 ,ω α 2 , ..,ω α k , ..ω α (j−1) ) (D.3) f k 's in general a j − 1 degree polynomial of the frequenciesω a , where a = 0, · · · N − 1.We consider f (ω α 1 ,ω α 2 , ..,ω α k , ..ω α (j−1) ) =ω α 1 ,ω α 2 , ..,ω α k , ..ω α (j−1) . Since we are only interested in the contribution of the leading divergent term to complexity, therefore we ignore certain terms arising as a result of the conditional sum and consider the much more simpler set of eigenvalues given by, We have counted that dimensions of each of these blocks are N 2 . Now as explained in the main text, each of these N 2 number of eigenvalues has a degeneracy of N. This argument can be easily generalized for arbitrary dimensions d. The contribution to the complexity coming from these blocks are given by (assuming the cost function C κ=1 ) sin 2 π i k N .

(D.5)
We used the following form of the penalty factor, Also, as explained in the main text we can chose for the reference,