Evolution of complexity following a quantum quench in free field theory

Using a recent proposal of circuit complexity in quantum field theories introduced by Jefferson and Myers, we compute the time evolution of the complexity following a smooth mass quench characterized by a time scale $\delta t$ in a free scalar field theory. We show that the dynamics has two distinct phases, namely an early regime of approximately linear evolution followed by a saturation phase characterized by oscillations around a mean value. The behavior is similar to previous conjectures for the complexity growth in chaotic and holographic systems, although here we have found that the complexity may grow or decrease depending on whether the quench increases or decreases the mass, and also that the time scale for saturation of the complexity is of order $\delta t$ (not parametrically larger).


Introduction
Entanglement, the most distinct trait of quantum mechanics, has proven to be a key concept in many different areas of physics. In particular, it has recently become a standard term in the high-energy physics vocabulary due mainly to the AdS/CFT correspondence, the Ryu-Takayanagi formula [1] and its generalization [2]. Many previous works have explored these connections (see, e.g., [3] for a review) with the aim of elucidating the basic mechanism of the AdS/CFT correspondence and space-time emergence.
Some questions remain unanswered though, in particular how the behavior of the interior of black holes (which is a region often inaccessible to Ryu-Takayanagi surfaces) in asymptotically AdS space-times is encoded in the dual CFT. In particular, by a classical gravity computation, it is possible to show that the volume of the Einstein-Rosen bridge connecting the two sides of an eternal AdS black hole grows linearly in time, which raises the question of what quantity this corresponds to in the dual CFT [4].
As we will explain below, it has been proposed that the computational complexity in the field theory, or complexity for short, could be the quantity whose growth is the dual to that of the Einstein-Rosen bridge in the gravity side. Complexity is a concept originally introduced in the field of classical computation [5], though it is also extensively studied in quantum information. The basic idea is to consider some reference state |Ω and some target state |ψ , as well as a set of quantum operators or gates, O I . The complexity of |ψ is then the length of the minimal circuit connecting it to |Ω by successive applications of the operators O I . To properly define "length"and "minimal", we need a notion of distance, or cost, associated to a given circuit, that in principle we are free to define. The definitions of cost and the choice of O I should be made as to capture the intuition that the gates are easy to apply. For instance, the operators should act on only a few degrees of freedom at a time. It is also possible to study the complexity of quantum operators themselves [6] instead of quantum states, but in this paper we focus on the latter.
Two conjectures which aim to connect the complexity of the boundary field theory to geometric constructs in the bulk space-time have been proposed, the so called "complexity = action" [7] and "complexity = volume" [8] proposals. In these works, it was argued that these objects in the bulk indeed grow as expected. The bulk computation of the evolution of holographic complexity was put on firmer grounds by [9], and it has also been subject of many recent works [10,11,12,13,14,15,16,17,18,19]. For instance, in [20], the holographic complexity following a quantum quench was calculated by modeling the quench with an AdS-Vaidya spacetime and it was found that the complexity evolves approximately linearly.
To truly check the proposed correspondence, it is necessary to have a precise definition of complexity in quantum field theory in the first place. That was only put forward very recently [21,22,23,24,25], where the authors defined and calculated the ground state complexity of free quantum field theories. Their approach was inspired by previous work by Nielsen [26] where the task of finding the optimal circuit connecting two quantum states was turned into a geometric problem of finding geodesics in a curved space. Some similarities with the general relativistic computation of holographic complexity were reported in [21], and [22], even though the theories with classical gravity duals are known to be very different, i.e., strongly coupled and with a large number of degrees of freedom.
The study of complexity in quantum field theory might also prove useful for subjects unrelated to holography. For instance, in [27], a lower bound on the complexity of an ensemble of operators in lattice systems was obtained in terms of out-of-time order correlators. So, in this sense quantum chaos and complexity seem to be quantitatively related and it would be interesting to check if something similar happens in continuum QFT's. It has also been conjectured [28] that a universal second law of complexity should hold for operators in quantum chaotic systems (analogous to the second law of thermodynamics).
In the present work, we use the definition of [21] and, following a construction similar to the one done in [29] and [30] in the context of cMERA, we calculate the complexity of circuits preparing the instantaneous state after a quantum quench and study its time evolution. Our basic aim here is to explore the definition [21] in a dynamical, non-equilibrium, setting, to compare our results with holographic calculations, and to use them as a benchmark for future studies on interacting theories. The calculation becomes particularly simple if the so called F 2 complexity is used, since in this case the problem of finding the optimal circuit geometrizes to that of computing geodesics in a Riemannian curved space. We shall deal with mass quenches in free bosonic theory. By choosing to work with free field quenches, we keep the problem computationally tractable while still allowing for non-trivial dynamics. In this case, one should not expect thermalization in the usual sense of Gibbs ensemble, since the free scalar field is an integrable system. Nevertheless, it was shown in [31] that many such integrable systems do thermalize in a broader sense, approaching the so called Generalized Gibbs Ensemble (GGE). In fact, it has been shown in [32] that after the quench the entanglement entropy in momentum-space agrees with the prediction from a GGE, with the density matrix after the quench being explicitly computed in terms of the quench data.
The paper is organized as follows. We begin with a brief review of smooth mass quenches in free scalar theories in Section 2. We then discuss the definition of circuit complexity in QFT as proposed in [21] in Section 3. In Section 4, we introduce our circuit by defining the reference and target states, as well as the elementary gates. The complexity is then calculated in Section 5, the final result being given by (74). The final remarks appear in Section 6. In Appendix A, we prove the equivalence of Schrödinger and Heisenberg picture calculations of the complexity, which is an important step of the calculation.
Note added: After the paper appeared online we have been informed of an upcoming paper [33] containing a very similar analysis.

Mass quenches in free scalar fields
We begin with a quick review of the approach introduced in [34] (the reader is referred to [32] for details) to study mass quenches in a free bosonic quantum field theory. The action describing the system is where the mass profile m(t) asymptotes to constants m in ≡ m(−∞) and m out ≡ m(+∞) at early and late times, respectively, so that asymptotic states exist before and after the quench. We shall focus on d = 2 for simplicity. This problem is equivalent to a standard quantum scalar field of constant mass m 0 placed in a cosmological background and therefore can be understood using intuition from quantum field theory in curved spacetimes [32]. We focus on the particular choice a smoothened version of a step function connecting the initial and final masses where the nontrivial time dependence happens roughly within the scale δt, since in this case exact solutions are known for the classical mode functions.
There are two types of mode functions, the so-called in-and out-modes φ in/out k (t, x) = (2π) −1/2 e i kx χ in/out k (t), depending on how we choose to fix the boundary conditions. They are given by where 2 F 1 [...] is the hypergeometric function and The in-mode satisfies lim t→−∞ χ in k (t) ∼ e −i ωint and therefore is adapted to observers at early times (at t = −∞, before the quench starts), while the out-mode satisfies lim t→∞ χ out k (t) ∼ e −i ωoutt and is adapted to post-quench observers at t = +∞.
Canonical quantization of the field φ then proceeds in the usual way. Both in and out-modes are equally good and the field can be expanded using either of them, i.e., To each set of modes (in/out) one can associate a corresponding vacuum state |0 in/out , namely, as well as particle excitations by acting with a † in/out k . They are of course different, being adapted either to early or late-time measurements only, though both correct. In fact, the two definitions are connected by a Bogoliubov transformation that is block-diagonal with respect to k (i.e., it mixes only opposite momentum modes) of the type where and |A k | 2 − |B k | 2 = 1. From this it is straightforward to check that the in and out vacua are indeed physically different: |0 in , which contains no particles for an initial observer, is populated by |B k | 2 particles of the out type (and vice-versa), which is the phenomenon of particle production by the quench.

Complexity in quantum field theory
We start by reviewing the definition of circuit complexity in quantum field theory as proposed in [21]. As mentioned in the Introduction, to properly define complexity one first needs to specify a reference state |Ω , a target state |ψ , and a set of elementary gates (quantum operators) {O I } to be used to build a quantum circuit connecting the them. The circuit then consists of a continuum set of operations parametrized by some parameter s (we take s ∈ [0, 1] with no loss of generality) as where the set of complex functions Y I (s) characterize a particular circuit roughly by deciding whether or not the gate O I is applied at "level s" of the circuit. As we shall see, in general there will be many distinct such functions preparing the same final state. The path-ordering symbol P here ensures that the circuit is built in the direction of growing s, i.e., that operators are applied sequentially from s = 0 to s = 1. In the following it will be convenient to introduce such that the circuit becomes Next, we need a notion of circuit depth that will allow us to quantify the costs of applying different circuits to accomplish the same task of connecting a pair of reference and target states. The computational complexity of this task (or, equivalently, the complexity of the state |ψ relative to |Ω ) will then be defined as the cost of the optimal circuit amongst all, i.e., the one with minimal depth. Given a circuit, or set of functions {Y I (s)}, its depth can be quantified by any functional of the form with the cost function F satisfying a list of reasonable physical requirements such as smoothness, positivity, and the triangle inequality (see [26] for details). It remains unclear whether there is a preferred choice of F and, for the sake of simplicity, we shall follow [21] and focus here on the so called F 2 distance defined by F = δ IJ Y I (Y J ) * , even though the strategy is exactly the same for any other distance measure. This means that, for our purposes, the circuit depth is defined by and the corresponding complexity is simply The main advantage of this approach, to be discussed in detail below, relies on the fact that the problem of minimizing D F2 [U ] to calculate the complexity can in principle be rephrased as a geometric problem of finding geodesics in a Riemannian manifold. Namely, for specific choices of reference and target states the set of gates O I needed to implement the task closes a Lie algebra and hence the circuits U are elements of the corresponding Lie group. One can then cover the group manifold with coordinates x in such a way that quantum circuits become curves in this manifold; the curves are defined by Y I (s) = Y I (x(s), ∂ s x(s)) and finding the optimal circuit amounts simply to finding the minimal length curve connecting the reference and target states. In particular, using the F 2 cost function above reduces the problem to that of calculating a geodesic in a Riemannian space 1 whose metric is given by In the following we shall design quantum circuits preparing the time-dependent state following a quantum quench from a product state having no entanglement over spatial subregions. We will see that an infinite number of elementary gates is needed, each acting on a particular momentum mode k of the quantum field. 2 These gates will be denoted by O I,k and the circuit by a set of functions Y I k (s). Thus, we will write the circuit as The sum runs over discrete values of k for the field in the compact space and must be replaced by V dk in the continuum limit, where V is the spatial volume factor needed to make the exponent dimensionless. In our case, when it is time to take the continuum limit we shall ignore this factor and set V = 1 for simplicity.

Introducing the circuit 4.1 Reference and target states
In this section, we specify the target and reference states that define our circuit. We shall choose our reference |Ω to be a product state, characterized by the absence of entanglement in any spatial subregion. A nice way to impose this condition (see [30,29,21] for similar work in the context of cMERA tensor networks) is by demanding the state to satisfy Reference: for some arbitrarily chosen mass scale M . It is straightforward to check that this definition implies [29] φ which indeed indicates that the state have no spatial entanglement. Our target state, on the other hand, will be time-dependent and shall be denoted |ψ(t) . We choose it to be the instantaneous state of the system at time t following a mass quench that starts from the vacuum of the pre-quench Hamiltonian at t = −∞, i.e., Target: where E is the time evolution operator. Notice that the states above are defined in the Schrödinger picture, which is why no time dependence shows up in the fields in (17) and the whole time dependence is contained in the target state |ψ(t) . Indeed, in the field representation |{φ(x)} (where π(x) = −i δ δφ(x) ) equation (17) can be seen as a functional differential equation defining the wave functional Ω[φ(x)] ≡ {φ(x)}|Ω , whose solution is a simple Gaussian function with width controlled by M . By discretizing the field into a chain of coupled harmonic oscillators this becomes precisely the Gaussian product state used by the authors of [21] as the reference state with respect to which the free scalar field complexity is evaluated. Similarly, in the absence of a quench, i.e. the limit m in = m out , the time evolution in (19) is trivial and the target state reduces simply to the vacuum of a free massive scalar field, which also corresponds to the choice of target state used in [21]. 3 Nevertheless, since the analysis of the exactly solvable quantum quench model of previous Section is done in the Heisenberg picture, it will prove more convenient to redefine our whole circuit and work in the Heisenberg picture instead. The physical intuition becomes slightly less obvious in this case though, since the roles are interchanged: the time dependence due to the quench now appears in the constraint defining the reference state (since the fields themselves are now time dependent), whereas the target state becomes constant. Anyway, the Heisenberg picture reference and target states are defined by Reference: or equivalently, using (5), (6) and (7), by the constraints where A and B are the Bogoliubov coefficients (8) and we have introduced Presenting the circuit in terms of the out modes is interesting since we are interested in studying time evolution of the complexity and, in particular, its late-time behavior after the quench has finished (a regime where the out modes are the most natural ones to use). For simplicity of notation, hereinafter we drop the superscript out and convention that mode functions χ k and creation/annihilation operators a k , a † k carrying no superscript always refer to out.
In principle it is not quite obvious that the circuit complexity will be the same in both Schrödinger and Heisenberg pictures, so in Appendix A we show explicitly that this is the case. Also, as discussed above, in the limit where m in = m out (no quench) our setup is equivalent to the one studied in [21], so as consistency check we shall be able to recover their results using Heisenberg picture calculations.
Having defined the target and reference states, the next task is choosing a set of gates connecting them. From now on it is always implicit that we are working in the Heisenberg picture.

Choosing the gates
The choice of gates used to build a circuit connecting the reference and target states is ambiguous, since there are not many requirements to constrain the options. In addition to being able to implement the desired task in a unitary way, we only require that the gates should act on the smallest possible number of degrees of freedom at each step, so that they constitute the most basic possible set of operations needed to construct the final state. This captures the intuition that the gates are easy to apply, and that the complexity then measures how difficult it is to complete the task in the sense of counting the minimum number of elementary operations needed. In a lattice system, for instance, this can be achieved by using gates that only act on neighboring lattice sites, i.e., gates comprising a notion of real space locality. But this is not the only option, and indeed the gates we shall introduce are different: they are local in momentum space, in the sense that they only act on specific fixed-momentum subspaces of the full Fock space (recall that the Fock space of free theories factorizes into subspaces of fixed-momentum modes [32]). Clearly these are non-local operations from the point of view of real space, but the important thing is that they also allow us to break the task of connecting two states in the full Fock space into a sequence of operations that act only fixed-momentum subspaces (the minimal number of which defines the complexity). This is, for instance, the way in which cMERA quantum circuits are constructed for free fields [29], where the procedure of acting with the gates has the nice interpretation of introducing entanglement at different energy scales. Ultimately, the definition of the gates is subjective and we often have to resort to choices of practical convenience, but the underlying idea here is to search for universal features of the complexity that do not depend on this specific choice. We shall see in the following that this is precisely what happens in the present work, namely, the results of [21] for the free scalar field complexity (obtained using a different set of gates based on real-space intuition) are exactly recovered using our set of gates in the limit of no-quench.
The target state |ψ(t) that we want to construct is built from the initial vacuum using essentially the quench Hamiltonian (through time evolution). So, to motivate our choice of gates, it is instructive to write the Hamiltonian associated with our model (1) in terms of creation and annihilation operators [32] We see that the Hamiltonian is written as a linear combination of operators O I,k which mix only k and −k modes (here I = 1, 2, 3 or, equivalently, It is straightforward to check that (for each k) these operators close a SU(1, 1) Lie algebra, 6 As we shall see below, using them as fundamental gates is enough to build a quantum circuit connecting the target and reference states introduced in the last section. Hence, they constitute a natural set of operators to construct our circuit, and will be our choice. A similar choice of gates has been used recently in [22] to study a different notion of complexity for quantum field theory states. For future convenience we shall adopt the set of gates {O 1,k , O 2,k , O 3,k }, so from now on whenever a summation over I appears it refers to I = 1, 2, 3.
The quantum circuit is written as where the functions Y I k (s) define the circuit in the same spirit as in [21]. We emphasize here the abuse of language when writing dk, which is to be understood as k for the field in the compact space or as V dk in the continuum limit (V is the spatial volume, which we ignore by taking V ≡ 1). Unitarity of U is guaranteed by taking the functions to be real, Y I k * = Y I k , since the O I,k are all Hermitian. Our task is then to find the functions Y I k (s) that minimize the circuit depth (13) while satisfying the boundary conditions Note that differentiation of (25) with respect to s yields The main task at this point is to isolate the Y I k to find the line element (15), but this is not possible in the present form. A similar problem was faced by the authors of [21] when dealing with circuits that build the ground state of free scalar fields. Their trick in that case was to reformulate the problem in matrix language and focus on how the circuit acted on 2 × 2 matrices defining the ground state, which allowed them to solve for the Y I k . We shall follow the same steps in the next sub-section.

Matrix representation
Let us begin by taking a look at how the transformation U (s) acts on the creation and annihilation operators. We first define their transformed versions by and differentiate with respect to s, i.e., d ds Now we use (25) and notice that, due to the path-ordering symbol P (which puts higher values of s always to the left), the derivative of U yields a factor on the left, For U (s) −1 (defined with anti path-orderingP) we have the opposite and the operator appears on the right, such that d ds The commutator is straightforwardly calculated using the canonical commutation relations for a k , a † k and yields a linear combination of a k and a † k themselves. These become a k (s) and a † k (s) after taking into account the U −1 and U in equation (32), so that one gets d ds where G k (s) is the matrix We immediately recognize G k (s) as an element of the SU(1, 1) algebra in the 2 × 2 matrix representation, i.e., where satisfy the SU(1, 1) algebra To simplify the notation, in the last step we have introduced Equation (33) is formally solved by The crucial point is that M k (s), being the exponential of the Σ I , is itself nothing but a SU(1, 1) transformation.
In other words, U (s) acts on the creation and annihilation operators as where M k (s) is a generic SU(1, 1) transformation that we can parametrize as in terms of some p k , q k satisfying |p k | 2 − |q k | 2 = 1. Hence our whole circuit can be understood as a sequence of This fact is what allows us to express the functions Y I k (s) defining a particular circuit in terms of a particular sequence of SU(1, 1) transformations, solving the problem raised at the end of the previous subsection. Using the fact that the generators Σ I satisfy the orthogonality relation equation (39) can be easily solved for the X I k (s) through differentiation with respect to s, namely, This determines the Y I k (s) associated with a given SU(1, 1) transformation M k (s) via (37).

Boundary conditions for the geodesics
Having isolated the Y I k (s) we now proceed to the second issue raised on the last section, namely, to express the boundary conditions (26) on U (s) as boundary conditions for the corresponding trajectories in the SU(1, 1) manifold.
Notice that the condition (7) can be written using (26b) as Applying U (1) −1 on both sides and using (40) we have and by comparing with the defining constraint (21a) for the reference state |Ω(t) we immediately identify where N k is a proportionality constant. This is a pair of equations that can be solved for p k (1) and q k (1) (the matrix components of M k (1)) in terms of A k , B k , α k , β k (the coefficients characterizing the reference and target states). The boundary conditions (26) are then rephrased as where we included also the last two which correspond simply to the statement that M k (0) = 1. The proportionality constant is fixed by requiring compatibility with the SU(1, 1) normalization |p k (1)| 2 − |q k (1)| 2 = 1, which leads to We see that there is still a phase factor freedom in N k . This is related to an ambiguity in the definition of our reference and target states which we shall discuss in detail in the next section, so for now one can simply take N k to be real without loss of generality and postpone the discussion of the phases to the next section.

Phase factor ambiguity
It is important to notice (see [29] for a similar discussion in the context of cMERA) that there is an ambiguity in the definitions (21) of our reference and target states. If one multiplies any of the constraints by a timedependent phase factor we still get the same physical constraint, i.e., Equivalently, this means that there is the ambiguity of replacing (α k , β k ) → e i η k (α k , β k ) and (A k , B k ) → e iη k (A k , B k ). 7 Since they appear in a mixed way in (47), there is a resulting relative phase ambiguity between the two terms that is nonvanishing. Even though the initial and final states are obviously not affected by this choice of phase, the circuit connecting them in general is and, as a result, the circuit depth associated to a given path can depend on it. In practice, this means that the boundary conditions (47) for p k , q k have to be replaced by (here ξ k ≡η k − η k ) and, when seeking for the optimal circuit to find the complexity, in addition to finding the geodesics associated with the metric (57) we also need to find the choice of relative phase ξ k (t) that minimizes the circuit depth. We will show explicitly that such a choice of ξ k is not the simplest one (ξ k = 0).

Geometrizing the problem
We are now ready to geometrize and solve the complexity problem. The first step is to parametrize the manifold of SU(1, 1) matrices using the set of coordinates {θ k , φ k , χ k } defined by which clearly satisfy |p k (s)| 2 − |q k (s)| 2 = 1. The boundary conditions (50) then can be translated as with χ 0 an arbitrary constant. Here, p * k (1) and q k (1) are to be understood as shorthand for the expressions (50a),(50b).
The functions Y I k (s) characterizing our circuit can be expressed in terms of the coordinates using equations (37), (41), (43) and (51). One finds where a prime is used here to denote d ds . The circuit depth (13) associated with a given curve Y I k (s) ≡ {θ k (s), φ k (s), χ k (s)} in the SU(1, 1) manifold is then given by where we defined the new variables The problem now becomes a geometrical one by noticing that this can be seen as the length of a particle's trajectory in a curved Riemannian manifold with metric (again, we stress that there is an abuse of notation here: the way to make sense of this is to formulate the problem in a compact space so that dk is replaced by a discrete sum over k). The optimal path is simply a geodesic in this geometry, the length of which defines the circuit complexity.
One can see that the metric is block diagonal with respect to the momentum-mode label k, i.e., ds 2 = Λ −Λ dk ds 2 k . It is also important to notice that each block k depends only on the coordinates {θ k , ψ k , φ k } associated with that same momentum k (there is no mode mixing). This means that we can find the curves minimizing the length in each of these momentum subspaces and add them together to find the geodesic in the whole manifold, which can also be checked by simply looking at the geodesic equation. Therefore the main goal now is to calculate the geodesic associated with the metric subject to the boundary conditions (52). Before delving into the calculations to find the geodesics, a comment is in order here. The metric above is diffeomorphic to the one found in [21] restricted to a hypersurface y = const. This should not come as a surprise since our choice of gates {O 1,k , O 2,k , O 3,k } leads to geodesics in the space of SU(1, 1) transformations, while the gates used by the authors of [21] lead to geodesics in the space of GL(2, R) transformations. It is well-known that SU(1, 1) is isomorphic to SL(2, R), a subgroup of GL(2, R), so at the end of the day we should expect to deal with a submanifold of the one appearing in [21].

Symmetries and Killing vectors
The task of finding the geodesics associated with the metric (57) can be simplified by taking advantage of symmetries and their corresponding conserved quantities. By inspection one can check that the metric (57) possesses a trivial Killing vector, k 0 ≡ ∂ v k , but the question remains whether there are others. As shown in [21], to each generator of the symmetry group there must be a corresponding Killing vector. This is because of the way the metric is defined, via the functions Y I k (s) in (43): if we redefine M k (s) → M k (s)M , whereM is a constant SU(1, 1) matrix, the Y I k (and consequently the metric) remain unchanged. So in our case we should expect at least 3 Killing vectors (one for each SU(1, 1) generator). It turns out that none of them happens to be the k 0 above, so there are actually 4 in total. 8 Following [21], let us takeM = e I Σ I with I (I = 1, 2, 3) constants. To first order in I , the statement is that the metric is invariant under To find the corresponding Killing vectors, we assume this change in M k is coming from a diffeomorphism δx j , i.e, and use (42) to find This can be inverted in the (j, I) indices to yield with the Killing vectors (k I ) j given by Indeed, we see that there are as many of them as there are SU(1, 1) generators Σ I . In terms of the coordinates {θ k (s), φ k (s), χ k (s)}, the Killing vectors above read which together with k 0 = ∂ v k discussed before consist of all the isometries of the metric (57).
To each of these Killing vectors one can associate conserved momenta c I = (k I ) i g ij x j which remain constant along the geodesics and can be used to simplify the task of integrating the geodesic equations. Their explicit expressions (apart from overall numerical factors) are the following We now proceed to use these equations to find the geodesics.

Geodesics and the complexity
From (64a) and (64d) it follows that Now, the boundary conditions at s = 0 (see (52)) demand that c 0 = c 3 = 0, otherwise u k , v k (and hence the length of the curve) would be divergent at s = 0. This implies u k = 0, v k = 0, i.e., where in the last step we have expressed the result in terms of the boundary condition for χ k using the definition (55) of u k , v k in terms of the original variables φ k , χ k . By plugging these results into (64b) and (64c) and imposing the boundary condition θ k (0) = 0 we get simply What remains is to determine the constant c 1 in addition to fixing the global phase factor ξ k using the boundary conditions at s = 1. These can be rewritten (using (66) and (67)) as with p k (1), q k (1) given in (50). Equation (68a) allows us to fix the phase ambiguity, namely we have to solve for ξ k . The solution is easily found to be which then implies This completely fixes θ k (1) in (68b) in terms of the target and reference state data, i.e., The remaining constants χ 0 , c 1 are also uniquely determined in terms of these data by equations (68b) and (68c), but the expressions will not be needed here.
Having obtained the geodesic {θ k (s) = θ k (1) s, u k (s) = χ 0 , v k (s) = −χ 0 }, the complexity is then readily found by substituting it into expression (54) to find its length (the minimal circuit depth). Since u k and v k have vanishing derivatives, we see that the complexity is fully determined by the boundary value of θ k (s), namely, Using (72) and (71) one can write this explicitly as This analytic expression for the complexity in terms of the reference and target state data (A k , B k , α k (t), β k (t)) given in (8) and (22), with the out-modes given in (3) is the main result of the present paper. The time evolution due to the quench is encoded in the coefficients α k (t), β k (t) and is highly nontrivial, via hypergeometric functions. For this reason the integral over k (or a discrete summation over k, to be more precise) cannot be carried out analytically, but it is straightforward to evaluate it numerically for specific choices of the parameters characterizing the quench if we impose a finite cutoff Λ. We shall do that in a moment, but first let us analyze some interesting limits.

Static limit: the vacuum complexity
As a consistency check, we can particularize to the static or no quench case that corresponds to m out = m in . Equivalently, this can be seen as the early time limit t → −∞ of the time evolution following the mass quench. In this limit, the construction above gives simply the ground state complexity of a massive scalar field (with mass m in ) previously studied in [21] (remember that our choices of reference and target states in this limit become exactly the same as those of [21]).
The coefficients A k , B k , α k , β k take the very simple form and the complexity (74) reduces to Using the identity we get which is exactly (the continuum version of) the expression for the complexity in a discretized massive scalar field computed in [21]. It is interesting to stress that the choice of elementary gates used in [21] differs significantly from the one used here, yet both results for the complexities agree. This indicates some degree of universality of the complexity with respect to the choice of elementary gates.

Sudden quench limit
Another interesting limit to look at is the infinitely fast or sudden quench limit δt → 0, where the mass profile (2) reduces to a step function. For t < 0 there is essentially no dynamics since m = m in is held fixed and we get simply the vacuum result above C vacuum F2 for the complexity. For t > 0, on the other hand, we have and the complexity (74) reduces to (80) In other words, in this limit the complexity jumps from the constant value (78) at t < 0 to an oscillatory behavior at t > 0 that persists forever. The frequency of oscillation is not easy to guess though, since ω out depends on the momentum variable k which is being integrated over. As we shall see below, this behavior is compatible with the smooth quench results, where the same late time oscillations appear for t δt.

Small and large reference scale M
One can also play with the free scale M that characterizes the reference state and obtain reasonably simple exact expressions for the complexity (74) in the extreme limits where M → 0 and M → ∞. Since M appears only in the functions α k (t), β k (t), taking these limits is straightforward using (22), namely while the difference |α k | 2 − |β k | 2 = Im (χ kχ * k ) is independent of M . This leads to where we have used arccosh(x) ≈ log(2x) (x → ∞) to convert the integrand into a logarithm. The main difference between the two limits appears in the second ratio inside the argument of the log, being the time dependence dictated byχ k (t) as M → 0 and by χ k (t) as M → ∞ as already anticipated by (81). Intermediate scales of M are therefore characterized by an interplay between these two behaviors.

UV behavior
It is important to notice that the ground state complexity (74) is UV divergent. This is easily seen by analyzing the large k behavior of the integrand, namely However, the leading UV divergence happens to be exactly the same as the one of the ground state complexity (76), since ω in ∼ |k| for large k. It is tempting then to subtract this constant divergent piece and introduce the "renormalized" complexity (i.e., the relative complexity of the the instantaneous state |ψ(t) with respect to the ground state one) which is guaranteed to be UV finite during the whole time evolution. By definition, this quantity vanishes at early times since the pre-quench state is nothing but the vacuum of the initial Hamiltonian. In the following we shall study it numerically for different choices of quench parameters as well as reference state data.

Numerical analysis
Having looked at special solvable limits of interest, we now analyze the renormalized complexity ∆C F2 (t) introduced in (84) in its full generality using numerical methods. As discussed above, this object is UV finite unlike the complexity (74) itself. In practice, the analysis shall be carried for a fixed and finite (but large) value for the momentum cutoff Λ. We emphasize once again that the integrand appearing in (74) is known exactly in terms of the reference and target state data, but its k-dependence is highly nontrivial such that doing the integration analytically becomes hopeless. That is why we need to resort to numerical integration, at this very last step of the analysis. The goal is to understand how the dynamics of ∆C F2 (t) is affected by the quench parameters δt, m in , m out (or δm ≡ m out −m in ), which characterize the target state |ψ(t) , as well as by the scale M that characterizes the reference state |Ω . Notice that, with the mass taken to depend on time as defined in (2), the quench effectively starts at t ≈ −δt and ends at t ≈ δt, so the plots below will always focus on this interval (−δt, δt) where a nontrivial time evolution is expected. The relevant results are shown in Figure 1. For numerical convenience we have fixed Λ = 100 and all the plots were made using specific numerical values for the parameters chosen for presentation purposes only, although the qualitative behavior illustrated in each one has been checked to hold generically.
Regardless of the choice of parameters one can see two distinct regimes for the time dependence of the complexity that are illustrated in Figures (a) through (d): an approximate linear evolution for −δt t δt, followed by a saturation phase at t δt, where the complexity approaches a constant with eventual fluctuations around this value. As illustrated in (a) and (b), the slope of the linear phase can be positive or negative (i.e., the complexity can grow or shrink) depending on the sign of δm. Namely, for quenches that increase the mass (δm > 0) the slope is always positive, while for mass-decreasing quenches it is negative. The amplitude of the final state oscillations is influenced not only by the mass difference δm, but also by the initial value m in as one can see by comparing (a) and (b): the former has m in = 1 and the oscillations become more pronounced as δm is increased, while the latter with m in = 6 and decreasing δm shows almost imperceptible oscillations.
The quench rate δt also has significant influence in the time evolution as indicated in (c). First, the slope of the linear phase increases the faster the quench is done. Second, the late time oscillations of the complexity are also sensitive to δt and become more pronounced as δt decreases. This is compatible with the infinitely fast quench limit δt → 0 discussed analytically in Section 5.3.2, which is also shown in the plot, where the complexity jumps (hence having an infinite slope) at t = 0 from a constant value to an oscillatory behavior for t > 0.
In (d) we analyze also the dependence of the complexity of the energy scale M that defines the reference state. To do this, we keep the quench data m in , m out , δt (i.e., the target state) fixed. Interestingly, we see that the renormalized complexity growth rate decreases and becomes even negative as the scale M is increased. In particular, there is a particular scale (M ≈ 3 for the set of parameters used in the plot) for which there is a transition from growth to decrease of the renormalized complexity. At this scale, ∆C F2 (t) remains approximately zero during the whole time evolution, i.e., the quench barely increases the complexity of the initial state. It is important to keep in mind that varying M affects not only the instantaneous complexity C F2 (t), but also the vacuum contribution C vacuum F2 that is being subtracted to yield ∆C F2 (t).

Higher dimensions
The analysis of the complexity above is easily extended to an arbitrary number of spacetime dimensions. Namely, the gate operators now carry a vector index k, i.e. O I,k , but still close the same SU(1, 1) algebra. The reference and target states remain unchanged, since the out modes χ out k (t) are the same in higher dimensions (see e.g. [34]) and therefore the defining constraints (21) do not care about the number of dimensions. The circuit is thus essentially the same and the calculation of geodesics proceeds in a similar way by focusing on a single k to later join them all together. As a result, the complexity in d = D + 1 dimensions takes the same form as in (74) but with where S D−1 is the area of the (D − 1)-sphere. 9 6 Discussion and Conclusion In this work, we have studied the time evolution of the computational complexity as the result of a smooth mass quench in a free bosonic field theory. This model at late times is known to thermalize in the sense of generalized Gibbs ensemble [31,32], so this can be seen as a toy model for studying the evolution of complexity during a thermalization process. The reference state was chosen to be one with no spatial entanglement, while the target state was a time-dependent one obtained from the pre-quench vacuum through time evolution. As elementary gates, we chose the ones appearing in the momentum-space Hamiltonian (23) of the model. Since they close a SU(1, 1) algebra, one can follow the procedure of Jefferson and Myers in [21] to geometrize the problem of finding the complexity as the one of finding a geodesic in the manifold of SU(1, 1) transformations.
The main result of the paper is formula (74) for the time evolution of the complexity C F2 (t) in terms of the data (A k , B k , α k (t), β k (t)) characterizing the reference and target states. In the limit where there is no quench and we have simply a massive scalar field in the vacuum state, the time dependence cancels out and we recover the same expression (78) for the complexity obtained in [21]. This is interesting since the set of elementary gates chosen here is quite different from the ones used in that case.
The complexity itself is well-known to be UV divergent, since there are infinitely many momentum modes (gates) k contributing to the circuit. However, we have shown that the UV divergent piece happens to be time-independent and exactly the same as the one from the initial state complexity C vacuum (the "renormalized" complexity) is by construction UV finite. We then studied the time dependence of this object as a function of the quench data m in , m out , δt as well as the free scale M characterizing the reference state.
We found that the time evolution of the renormalized complexity is marked by two distinct phases, separated by a time scale set by the quench rate δt (or, equivalently, the thermalization scale for the model). For −δt t δt, the complexity grows (or decays, depending on the sign of δm ≡ m out − m in ) approximately in a linear way, with a slope dC dt inversely proportional to δt and directly proportional to δm. Then, at a time of order δt the complexity saturates and starts oscillating around a mean value. We checked that the amplitude of the late time oscillations is sensitive not only to δt and the mass difference δm, but also to the starting point of the quench (m in ) and the reference scale M . Also, for fixed quench data the renormalized complexity shifts from growing to decreasing in time as the scale M is increased, and in particular there is an intermediate scale for which the change in complexity due to the quench is nearly vanishing.
As mentioned before, it has been conjectured that the complexity growth might be dual to the growth of Einstein-Rosen bridges in eternal AdS black hole geometries. Classical gravity computations show that observables characterizing this growth [9] evolve linearly at late times. The holographic complexity was observed to grow linearly as well in AdS-Vaidya model of a quantum quench [20]. Also, a linear growth followed by saturation with oscillations around a mean value had been previously conjectured to be the universal behavior of the complexity for chaotic systems [28], [35]. The time scale for saturation in these systems is believed to be parametrically larger than the thermalization scale. Our results show some similarities with these conjectures and calculations even though we dealt with a free (hence non-chaotic) field theory, with the important difference that the time scale for saturation here is of the order of the thermalization scale δt. For this reason, it seems to be an interesting future prospect to extend the calculations presented here to interacting theories, which are the ones most likely to have a simple classical gravitational dual, and check if the early linear behavior shown here can be pushed to last up to larger time scales. Another interesting prospect in this direction would be to use the κ = 1 measure discussed in [21] instead of the F 2 -distance used here, which has been argued to be more appropriate to get closer to holographic results. scholarship of CAPES -Brazilian Federal Agency for Support and Evaluation of Graduate Education within the Ministry of Education of Brazil, and by the CNPq grant 146086/2015-5. G.C. acknowledges financial support from the Brazilian ministries MCTI and MEC.

A Relation between Schrödinger and Heisenberg picture circuits
As introduced in Section 4, in the Schrödinger picture the target state is time-dependent while the reference state is not. Namely, the target state is |ψ (t) obtained by evolving on time the initial state |0 in while the reference state |Ω is constrained to satisfy where we have chosen, with no loss of generality, the time slice t = 0 for the Schrödinger picture operators. The optimal circuit is |ψ (t) = U Sch optimal (s = 1, t) |Ω ≡ P e i 1 0 ds Y I (t,s )OI |Ω , which can also be written as |ψ (0) = E (t, 0) −1 P e i 1 0 ds Y I (t,s )OI |Ω , where E (t, 0) is the time evolution operator. In the Schrödinger picture, the task is therefore to find a circuit as above satisfying the boundary condition One can write the optimal circuit as |ψ (0) = U Hei optimal (s = 1, t) |Ω (t) = P e i 1 0 ds X I (t,s )OI |Ω (t) = P e i 1 0 ds X I (t,s)O I E (t, 0) −1 |Ω . (A.8) and the task is to find the set of functions X I that satisfy the boundary condition U Hei optimal (s = 0, t) = 1 (A.9) and minimize ds δ IJ X I X J * . The question then is how to relate the two apparently distinct ways to calculate the complexity. To answer that, let us multiply the last condition by E −1 E = 1 on the left to obtain Recall that the Hamiltonian for our problem is a combination of creation and annihilation operators that close a SU(1, 1) algebra, and we have chosen the {O I } to be exactly that set of operators, that is, for some functions f J (t) that can be seen from (23). This means that