Cluster expansion for ground states of local Hamiltonians

A central problem in many-body quantum physics is the determination of the ground state of a thermodynamically large physical system. We construct a cluster expansion for ground states of local Hamiltonians, which naturally incorporates physical requirements inherited by locality as conditions on its cluster amplitudes. Applying a diagrammatic technique we derive the relation of these amplitudes to thermodynamic quantities and local observables. Moreover we derive a set of functional equations that determine the cluster amplitudes for a general Hamiltonian, verify the consistency with perturbation theory and discuss non-perturbative approaches. Lastly we verify the persistence of locality features of the cluster expansion under unitary evolution with a local Hamiltonian and provide applications to out-of-equilibrium problems: a simplified proof of equilibration to the GGE and a cumulant expansion for the statistics of work, for an interacting-to-free quench.


Introduction
One of the main goals of quantum many body physics is to determine the ground state of physically interesting models. Whenever the ground state wavefunction is known, all other ground state properties and observables of interest can be derived. In practice in order to extract physically important information about the ground state, it is valuable to know it in a form that reflects its general physical properties and allows a systematic expansion of observables from more to less relevant information. For example, Renormalization Group methods are supposed to capture the universal ground state properties of physical systems in the thermodynamic limit. One of the most fundamental physical principles that characterize physical systems is locality and this is clearly reflected on the properties of their ground state. Indeed physical Hamiltonians are expected to involve only local couplings or interactions between near different points and this imposes constrains on the form of their ground state. In particular, ground states of local Hamiltonians possess an extensive number of particles, energy and other thermodynamic quantities; relative fluctuations of such extensive quantities about their average values decay rapidly with the system size; local observables, like multi-point correlations, in general have a well-defined limit for large system sizes. These properties are deeply linked to the extensivity of a suitably defined 'free energy' from which all other thermodynamic quantities and local observables can be calculated.
The extensive part of this quantity can be expressed conveniently as a cluster expansion which, roughly speaking, is an expansion in terms of 'joint cumulants' or connected multi-point correlation functions in the state. Cluster expansions have been used extensively in statistical physics, both in classical and quantum systems, typically for the study of the partition function of thermal states [1][2][3]. They have also played an important role in the constructive approach to quantum field theory [4,5] and systematic applications have been developed in various fields, from nuclear physics [6] to quantum optics [7]. Although most studies focused mainly on the partition function, it is also possible to define such expansions directly for the quantum states themselves, which is the approach adopted in the present work.
This approach has the advantage of expressing the state in an exact form, yet and more importantly organizing its information content in order of decreasing physical relevance, while clearly demonstrating its locality properties at any order. It also provides the possibility to study the thermodynamic properties of overlaps between different states. Such overlaps are useful in quantum information theory where they are known as fidelity between two states, which has been shown to manifest criticality by exhibiting non-analytic behavior when one of the two states approaches a critical point [8][9][10][11][12][13]. It is also useful in out-of-equilibrium problems where a quantum system undergoes a rapid change in some parameter, a process known as quantum quench [14]. Quantum quenches have attracted significant interest partially due to experimental implementation in cold atom systems [15][16][17][18][19][20][21][22] and due to the possibility of studying thermalization and more generally equilibration in closed quantum systems [23][24][25][26][27][28][29][30][31][32] (for a review [33]). In this context, overlaps between the initial state and post-quench excited states are relevant in the study of outof-equilibrium quantum thermodynamics as they capture information about quantum fluctuation relations [34] and in particular the statistics of work done by such quenches [35,36]. As we will see, the cluster expansion provides an exact spectral expansion of the state, allowing the systematic derivation of thermodynamic quantities, like the quench work statistics, as well as their quantum fluctuations. We suppose the existence of fields ψ † ( x) and ψ( x) in a d dimensional space that satisfy canonical bosonic commutation rules We suppose the existence of a unique translational invariant vacuum |0 annihilated by ψ( x) and assume that the Hilbert space can be described acting with ψ † on the vacuum. In this perspective, the field ψ † ( x) has the meaning of creating a particle in position x. A general cluster expansion of a quantum state | is of the form The functions W c n ( x 1 , x 2 , ..., x n ) will be called cluster amplitudes and are related to the correlations between different points. The simplest examples of states written in this form are the coherent and squeezed coherent states, in which the expansion terminates after one or two terms respectively. In general any state that is gaussian in terms of the field operators ψ † ( x) can be written in the above form with only the first two cluster terms Indeed only states of this form satisfy Wick's theorem, i.e. their multi-point correlation functions can be decomposed into combinations of two-point and one point correlation functions in all possible ways, or equivalently, connected correlation functions between more than two points vanish. Higher order terms in the above expansion give rise to non-zero multi-point connected correlation functions. In this respect the cluster expansion is clearly analogous to the cumulant expansion of probability theory.
Even though (2) is written in the continuum limit, the general construction can be applied to lattice systems as well, once we substitute the spatial integrals with sums over the lattice sites. The above expansion can be equivalently written in terms of the Fourier transforms of the fields ψ † ( x) and ψ( x), respectively a † k and a k , that create and annihilate particles with definite momentum k. In determining the ground state of an interacting Hamiltonian, the a † k and a k operators can be chosen to diagonalize the free part of the Hamiltonian, so that (2) is an expansion of the interacting ground state in terms of the eigenmodes of the free Hamiltonian, but other choices are in principle allowed and sometimes more suitable. Also note that, although we will focus on bosonic systems, our method could also be extended to systems of fermions or other particle excitations obeying different statistics which will be subject of future studies. The general construction can be applied to condensed matter systems and relativistic quantum field theories; for explicit calculations we mainly focus on the latter context.
As we will see below, the locality of the underlying Hamiltonian imposes general thermodynamic conditions on their ground states that do not rely on the details of the particular system. We call states that satisfy these properties local quantum states.
The cluster amplitudes for the ground state of a specific model can be determined by solving a set of functional equations that are equivalent to the time-independent Schrödinger equation and similarly time evolution can be studied by a set of time-dependent equations. Both sets of equations are consistent with the locality properties at all orders and in particular the evolution under a local Hamiltonian preserves these properties. Apart from the potential of exact description of ground states, the cluster expansion is also valuable in designing approximation schemes, to a large extend independently of the parameter or quality that is considered small in a particular approximation. As we will see, the cluster expansion reproduces efficiently results from standard perturbation theory. More importantly, it is also suitable for non-perturbative methods, like self-consistent approximation or variational methods, but it could be also used to study out of equilibrium problems. Lastly, it can be useful in designing efficient numerical methods [37][38][39][40]. In this context the cluster expansion provides a convenient way to express any ground state that automatically satisfies locality requirements. In general any approximation based on suitable truncation of the cluster expansion captures correctly the qualitative features related to locality and gives increasingly better quantitative results as the truncation order increases. In the next section we present the structure of the paper and a detailed summary of the main results.

Summary of results
The paper is organized in three parts. The first part presents the construction of the cluster expansion in a general abstract way and explains its main properties based on the physical requirements related to locality. The second part shows how the cluster expansion for ground states of general physical models can be derived by means of truncation schemes and recursive approximations, but also perturbative methods with which its consistency is demonstrated. The third part discuss the time-dependent formulation and applies results of the previous parts to out-ofequilibrium problems, especially quantum quenches. Below we summarize the main results of each section:

Construction and properties of the cluster expansion (Sec. 3-4 and their subsections)
(a) We identify the main consequences of locality in physical systems to the structure of their ground state. One of them is the cluster decomposition principle, i.e. the principle that connected multi-point correlation functions decay with the distance between the points. This property is valid for infinite sized systems. Considering instead large finite systems with local or short-range interactions, we expect that their macroscopic properties are described by a (suitably defined) extensive logarithmic partition function log Z.
More generally for such systems we expect that, in the thermodynamic limit, thermodynamic quantities exhibit vanishing fluctuations and correlations functions of local observables are well-defined. We use these general physical requirements to identify a class of local quantum states in which ground states of local Hamiltonians belong. (b) We define the partition function Z[ ] of a quantum state | through its norm, as compared to the vacuum of the bosonic fields over which we expand it and which is supposed to have non-zero overlap with | . (c) We construct the cluster expansion for a pure state, whose information are completely encoded in suitable functions called cluster amplitudes. (d) We construct suitable graphical tools, similar to Feynman diagrams, to compute Z[ ] and the correlation functions in terms of the cluster amplitudes W c n . (e) We obtain a sufficient condition in order for the cluster amplitudes W c n to describe a local quantum state at any order in the graphical expansion, specifically they must decay sufficiently fast with the distance. (f) Through the relation between W c n and the correlators of ψ and ψ † we show that the cluster decomposition property is automatically satisfied at any order of the graphical expansion when the cluster amplitudes W c n satisfy the above condition.

Derivation of cluster expansion for ground states (Sec. 5 and its subsections)
(a) We derive the cluster amplitudes W c n for the ground state of a general class of physical models by means of a perturbative approach. This provides a check of the general statements of the previous section and gives a natural truncation of the cluster expansion at each perturbation order. (b) We derive general functional equations for the cluster amplitudes of ground states.
These constitute the time-independent Schrödinger equation for the cluster expansion of the ground state wavefunction. (c) We demonstrate the application of the above to the case of a free quantum field theory, where we easily derive the known 'squeezed vacuum' form of the ground state, as well as to the interacting φ 4 theory. We verify the consistency with standard perturbation theory. (d) We comment on possible approximation schemes to describe ground states through the new tool of the cluster expansion. In particular the exact information encoded in the cluster expansion of local quantum states opens up the possibility of variational approaches, which could give insight beyond standard perturbation theory.

Applications to out of equilibrium thermodynamics (Sec. 6 and its subsections)
(a) We verify that, when a local quantum state evolves according to a local Hamiltonian, it remains a local quantum state during the evolution, that is to say the general properties of the cluster amplitudes of this class of states are not spoiled along the time evolution. (b) The evolution of quantum systems are studied in the Schrödinger picture through the time-dependent Schrödinger equation expressed in terms of the cluster amplitudes. This could give insight in new approximation schemes to describe systems out of equilibrium. (c) In the context of quantum quenches, it has been earlier shown that initial states that are ground states of free or interacting models and evolve under massive free Hamiltonians, tend for large times to some generalized sort of equilibrium, described by the so-called Generalized Gibbs Ensemble (GGE), which is expected to be valid more generally for the evolution under an integrable Hamiltonian. This can be shown much easier by using the general information about local quantum states. (d) A direct consequence of the cluster expansion properties is also that in the case of quenches where the initial state consists solely of pairs of particles of opposite momenta (which is a very common case) the initial state is in fact a squeezed state. This simplifies dramatically the study of its time evolution. (e) Starting from the knowledge of the cluster expansion, we compute the statistics of the work done during a quantum quench from an interacting to a free quantum field theory.

General construction of cluster expansion
This section is mainly divided in two parts: after having introduced the cluster expansion for general states, in Section 3.1 we describe the thermodynamic properties that a ground state of a local or short range Hamiltonian is supposed to possess. Through these properties, in Section 3.2 we will introduce the class of local quantum states and anticipate the main restrictions imposed by locality requirements on their cluster expansion, leaving the technical details to Section 4. In order to construct the cluster expansion of a general state, we consider a system of local bosonic fields ψ( x) defined in a finite box of length L with periodic boundary conditions and we will later take the limit L → ∞. We define creation and annihilation operators of the bosonic fields in momentum space through the Fourier transform of the fields in coordinate space ψ( x): The momenta k are quantized and the bosonic operators a k satisfy canonical commutation rules where δ k, q is the Kronecker delta over all momentum components and n i are integers. Any state | can be constructed by acting with the ψ † ( x) operators on the vacuum of the theory, therefore we can expand it as The integrations from now on are meant to be in the finite box of length L. In the following, we will suppose 0| = 0, so we can normalize the state with respect the reference state |0 , that is to say 0| = W 0 = 1. It can be simply shown that any state such that 0| = 0 can be written in the exponential form: where the W c n is the connected part of the W n function, defined in a way analogous to the cumulant expansion. The first few terms are: The consistency check of (6) and (7) is a straightforward calculation. From now on, we will focus on states invariant under coordinate translations, thus each W c n must be translationally invariant. It is useful to express the above expansion also in momentum space, which is simply done by a Fourier transform of the W c n functions and of the ψ † fields: where K n is linked to the Fourier transform of W c n in the infinite volume limit: In the limit of infinite system size L → ∞ we can write: where now the commutator of a( k) is a Dirac delta function instead of a Kronecker delta, i.e.
[a( k), a † ( q)] = (2π) d δ d ( k − q). Notice that in the above expression we have explicitly separated the Dirac delta function which is due to the translational invariance. Also notice that any explicit dependence on L has disappeared. The cluster expansion can be written for all those states with non-zero overlap with the bare vacuum, but ground states of local Hamiltonians are not arbitrary states and satisfy precise thermodynamic requirements: the next section is devoted to describe these properties and to the introduction of the local quantum states.

Thermodynamic properties of ground states and local quantum states
As we already pointed out in the introduction, ground states of local Hamiltonians are not random states of the Hilbert space, but satisfy certain conditions following from locality. Below we enlist some common manifestations of locality.

Extensive thermodynamic quantities
Ground states of extended systems are characterized by well-defined extensive thermodynamic quantities. Extensivity of the energy, for example, means that the expectation value of the system's Hamiltonian increases linearly with the system's volume V = L d in the limit L → ∞. The same behavior is also expected for other thermodynamic quantities that correspond to spatial integrals of local density operators over the whole system where Q(x) is considered as 'local' if it only involves operators acting on a neighborhood of the point x. In a lattice model this would mean operators acting on the site x and those within a finite radius around it. In the continuum limit, this is equivalent to considering operators that may include spatial derivatives of any finite order.

Cluster decomposition and vanishing relative fluctuations
Another crucial property of such states is that relative fluctuations of thermodynamic quantities vanish in the thermodynamic limit L → ∞. This property is the basis of an important ingredient of standard statistical physics, the equivalence of statistical ensembles. It is also related to locality, more precisely to the cluster decomposition principle, which is another fundamental requirement that ground states of local Hamiltonians are supposed to satisfy. It states that connected correlation functions between observables localized at large spatial separations must vanish e.g. for two point correlations This principle originates from quantum field theory in infinite space and expresses the requirement that quantum measurements at distant spatial points must be independent from each other. In the context of statistical physics, where we are interested in the behavior of systems of finite size L in the thermodynamic limit, it refers to distances R that tend to infinity not faster than the system size L i.e.
where the subscript 'L' means that the expectation value is calculated in a system of finite size L. The vanishing of relative fluctuations of extensive quantities in the thermodynamic limit is a direct consequence of the above relation. Indeed due to (14) the variance of Q as defined by (12) scales slower than L 2 and therefore the relative variance decays in the thermodynamic limit

Extensivity of 'free energy'
In thermal ensembles, thermodynamic quantities are in general derived as logarithmic derivatives of the partition function Z th with respect to intensive parameters [41]. For example, for a thermal state of temperature β −1 the partition function is Z th = Tr e −βH where H is the system's Hamiltonian, and the extensive energy is H = −∂(log Z th )/∂β. It is useful to introduce the concept of partition function also for a generic pure state | . We will suppose that the state | has non-zero overlap with the vacuum |0 defined as the reference state in Section 1. Then we define its partition function as the norm when the normalization is fixed in reference to the vacuum of the theory It should be mentioned that the above normalization presupposes a suitable UV regularization. Indeed in most cases in relativistic QFT the above overlap, when both states are normalized to one, is UV divergent in perturbation theory. We will see this fact also in our framework and see that, assuming the normalization 0| = 1 with an explicit UV cut-off, the norm | presents UV divergences when the cut-off is taken to infinity. The aim of the present work is to study general thermodynamical properties of pure states that are linked to the infrared behavior of the theory and are not related to their UV behavior. UV properties are, in some sense, less general than the infrared ones: UV divergences and the proper renormalization scheme depend both on the dimensionality and on the specific model, but different systems share the same general thermodynamic properties. Moreover, if the field theory is the continuous description of a lattice system, its lattice space provides a natural cut off: actually, in the presence of a lattice, our construction can be still applied provided we replace the integrals in the coordinate space with the proper summations. Of course, a proper treatment of the UV singularities is necessary to compute physical expectation values in actual models and this problem deserves further investigation, but it is beyond the scope of the present work that aims to a description of general thermodynamic properties. Because of the above, in the following we will assume that all the expressions are UV finite by means of a suitable cut off. As a general remark, note that a proper renormalization scheme must guarantee finite physical quantities (i.e. correlation functions), but the overlaps in the renormalized theory could be still UV ill-defined when we take the cutoff to infinity, since they are not physical observables.
In the following, we define log Z[ ] as the 'free energy' of the state. This definition for the free energy has some similarities with the usual thermal one and in particular it is expected to be an extensive quantity. We can justify this statement when | is a non-degenerate ground state of some local Hamiltonian H and the vacuum |0 is the ground state of another local Hamiltonian. For example, if H is made by a free part plus an interaction, |0 could be the ground state of the free part of the Hamiltonian, but this is not the only possibility. In this case, the free energy log Z[ ] becomes a well known object in the study of the Casimir effect [42] known as surface free energy. Consider a new quantity Z 0 so defined: In the same way that a thermal partition function of a quantum system can be seen as the partition function of a classical system in the d + 1 dimensional cylinder of circumference β, in the present case Z 0 can be seen as that of a slab geometry. The width of this slab is given by β and the boundary conditions on the slab are fixed by |0 , that can be interpreted as a boundary state [36]. The free energy log Z 0 is the sum of several terms [42]. The leading term for large β is extensive in the whole volume of the classical system, therefore proportional to the volume of the d dimensional quantum system and to the length β of the extra dimension: this is the bulk contribution to the free energy. A second term due to the boundaries is also present and it is usually called surface free energy: this quantity scales with the size of the boundaries, but it is not affected by their distance. Passing from the quantum to the classical system, the size of the boundary is simply the volume of the d dimensional quantum system: the surface free energy is independent from β, but it is extensive in the volume of the quantum system. Subleading terms in β → ∞ are generically present, but we are not interested in them. We are going to see that the free energy log Z[ ] is nothing else than the surface free energy in the classical system, therefore it must be extensive in the volume of the quantum system. In order to do so, we let β → ∞ in (18): in the zero temperature limit, the thermal ensemble becomes a projector on the ground state.
Above, E G is the energy of the ground state | . If the theory presents a gap between the ground state and the first excited state, further corrections in (19) are exponentially suppressed. We explicitly insert the normalization factor | because we are using the non-standard normalization 0| = 1. Thanks to (19) we can immediately evaluate the zero temperature limit of Z 0 : Above, further corrections are exponentially damped when β → ∞. Looking at (20) we can exactly recognize the contributions described before. The term βE G is proportional both to the length of the classical extra dimension and to the volume of the d dimensional quantum system, since the ground state energy is an extensive quantity. Therefore, −βE G is the bulk free energy of the classical system on the slab. The next term is independent from β and therefore it must be recognized as the surface free energy, extensive in the volume of the quantum system: apart from a minus sign, it is exactly the free energy of the state defined through its partition function (16).

Well-defined correlation functions in the thermodynamic limit
Another property commonly required in quantum field theory is that correlation functions have a well-defined limit when the volume becomes infinite. In fact cluster expansion methods have been developed in the past primarily in order to check the validity of this requirement [4]. In this context the system size L plays the role of an infrared regularization and in the end all physically meaningful quantities must be independent from it as L → ∞. In many physical contexts, this property may be considered as a consequence of cluster decomposition, in the following sense. Cluster decomposition means that regions of the system very far apart are completely uncorrelated to each other. A small increase of the system size means physically that we modify the system only at its boundaries and therefore a measurement performed deep in the bulk should not be affected by this small change. As a result, bulk correlations should be independent from the system size when this is already very large. However we will recognize this property as a separate one in order to avoid restrictions coming from such more special considerations.

Local quantum states
Motivated by the general arguments of the previous section, we define as local quantum states the states that share these thermodynamic properties. As we saw the various properties are not all independent but instead they are related to each other. We therefore identify three general properties which are more fundamental, in the sense that other properties, like the extensivity of thermodynamic quantities, vanishing of relative fluctuations and ensemble equivalence in the thermodynamic limit, can be deduced from them. Also we can consider these three properties as rather independent from each other, in the sense that one does not necessarily imply the other, except in more special settings.

Local quantum states and their properties:
We denote a state | as a local quantum state if the following properties are satisfied: • Property 1: The correlators of the fundamental local fields ψ, ψ † have a well-defined L → ∞ limit, that is, in this limit they are independent from the system size. • Property 2: The cluster decomposition property for the ψ , ψ † fields holds. • Property 3: The 'free energy' log Z[ ] is an extensive quantity. Property 2 is the cluster decomposition property expressed in terms of the fields ψ: if ψ is a physical field itself, then this requirement is obvious. We will see at the beginning of Section 5 how to properly define ψ in the context of common relativistic field theories.
For the sake of clarity, we anticipate the connection of Properties 1, 2, 3 with the form of cluster amplitudes of local quantum states. In Section 4 we introduce a diagrammatic technique for the expansion of the free energy and correlation functions in terms of cluster amplitudes and then we find a sufficient condition to ensure that any graph of this expansion is consistent with Properties 1, 2 and 3. We require the cluster amplitudes to be independent of the system size and to decay fast enough such that: The above condition must be interpreted as follows: any UV singularities in W c n must be regularized by means of an explicit UV cut off, in such a way that W c n does not have any such singularity. In this case the bound M n depends explicitly on the cut off, as it can be seen in the special case j = n of (21): Once W c n becomes bounded, condition (21) becomes a constraint on the decay of the cluster amplitude at large distances: it must decay fast enough to have that the function we obtain integrating out some coordinates in |W c n | is still bounded. Of course, if W c n has some singularities, after we remove the UV cut off (22) is no longer valid. In this case the diagrams of Section 4 are not guaranteed to be finite any longer, but eventual divergences are not due to the thermodynamic limit L → ∞, rather to UV singularities of the cluster amplitudes.
Condition (21) has some non-trivial consequences in Fourier space: The Dirac delta in the definition (10) of K n has been eliminated by fixing one of the coordinates of W c n and integrating on the others. Using (21) with j = 1 we have: Note that (21) is a sufficient condition that ensures Property 1, 2, 3 at any order of the expansion, but the latter could not converge: in principle there could be local quantum states whose cluster amplitudes do sot satisfy (21). Nevertheless, in Section 5.1 we provide an explicit computation of cluster amplitudes of ground states in perturbation theory. At any order in perturbation theory, K n are non-singular functions of their momenta. Moreover, in Section 5.2 non-singular functional equations for the cluster amplitudes of ground states are derived: these evidences strongly suggest that (21) could be even a necessary condition to have Property 1, 2 and 3. We note that in the case of states with a fixed number of particles N the properties 1, 2, 3 are violated at first sight. For example, the cluster property is spoiled: if we measure a many-body observable and find that all particles are in a subset of space then obviously further measurements would find no particle anywhere else, independently of the distance, i.e. the measurements are strongly dependent. In such systems with fixed number of particles, the correct definition of the thermodynamic limit is that both the number of particles N and the system size L tend to infinity so that the density N/L remains constant. With this definition, any observable that measures a finite number of particles would satisfy the cluster property and our analysis is again applicable. In Appendix A we outline how to apply our formalism in systems with a finite number of particles through a familiar example, the prototypical model of the interacting Bose Gas in d dimensions.

Diagrammatic technique
In this section we present a diagrammatic technique, analogous to the Feynman diagrams, which is our main tool for the study of the cluster expansion. Only a basic exposition is presented here, leaving further considerations to Appendix C. In Section 4.1 we use this formalism to study the effects of locality requirements on the cluster expansion of local quantum states and derive the results we outlined in the previous section.
The fundamental quantities we would like to compute are the free energy and the correlators of the fields ψ , ψ † . These quantities, apart from the particular case of gaussian states that have only W c 1,2 = 0 (see Appendix B for the explicit calculations), cannot be computed in a closed form and we will have a sum of complicated terms that involve products and integrations which can be expressed conveniently as graphs. We first construct the graphical rules to evaluate the norm | ; next we generalize them to include also correlators. As can be guessed, the 'interaction' vertices in our diagrams will represent the W c n functions. When we compute | we have W c n coming from the ket and [W c n ] * coming from the bra, therefore we need two different kinds of vertices. To distinguish them, instead of using different vertex notation, we attach an arrow to the lines and the convention is that a vertex with n outward pointing arrows is associated to W c n , while a vertex with n inward pointing arrows is associated to [W c n ] * (Figs. 1-2). To the edges of each leg is associated a coordinate.
The contraction of the fields ψ and ψ † in the bra and ket means that we are inserting delta functions on the coordinates; we can indicate this by joining two legs of the vertices. Notice that we can contract only fields that come from the ket with fields that come from the bra, which is the reason why we use arrows: each vertex has either only ingoing arrows or only outgoing arrows; mixed configurations are not allowed. This recipe ensures that we are contracting the right fields and each vertex representing a term of the ket will be joined with a vertex representing a term of the bra, as it should. When we want to compute Z[ ], all the W c n and [W c n ] * are contracted with each other. Thus the graphs contributing to Z[ ] have no external legs: all the arrows start from some vertex and finish in another. To obtain the contribution we are looking for, we should integrate over all the variables of the internal legs. Moreover there are some additional symmetry factors. Since we must perform all the possible contractions, in principle we can get several times the same graph from different contractions and we must take care of this degeneracy. The factorial terms of (7) are there to automatically take care of most of the possibilities, but often they produce an over-counting and each contribution must be divided by the number of permutations of vertices that do not change the graph. Then we must also divide by the number of permutations of internal lines that leave the graph unchanged. We give some examples in Fig. 3. Exactly as happens in the usual Feynman diagrams theory [44], the symmetry factors are such that Z[ ] can be exponentiated in terms of the connected graphs: The connected graphs are simply the graphs that cannot be split into two or more disconnected pieces. Taking the logarithm of (25) we discover that the free energy can be computed as the sum over the connected blocks: Once we understood the graphical rules to compute the free energy, it is easy to generalize them to the correlation functions. Thus suppose we want to evaluate a generic correlator: When we perform the contractions in the numerator, we must take into account also the possible contractions with the fields of the correlators. To do so we introduce the possibility of having external legs, that is to say, arrows that do not link two vertices. Also in this case the arrow notation is pretty useful: since a ψ field can be contracted only with an element of the ket, an arrow departing from a vertex will be associated to a contraction with the ψ field. Similarly an arrow that ends on a vertex will be associated to ψ † . The coordinates associated to external legs are fixed by the field in the correlator and we should not integrate over them. The remaining rules are the same as for Z[ ]: a simple example is given in Fig. 4.
With these rules, it is quite easy to formally compute the numerator of (27), moreover, as happens in the usual Feynman diagrams [44], all the blocks in the graphs of the numerator of (27) that do not contain external legs can be exponentiated, giving exactly | . Thus, at the end of the story, correlators such as (27) can be evaluated simply summing the contributions of all the graphs in which each block has at least one external leg.
Diagrammatic rules to compute the free energy (coordinate space): 1. Draw a connected graph, that is to say, it does not contain any disconnected pieces. The vertices are connected with arrows and each vertex has either all arrows outgoing or all ingoing; mixed configurations are not allowed. External legs are not allowed. 2. The contribution of each graph can be computed as follows: assign to a vertex with n outgoing arrows the function W c n , to a vertex with n incoming arrows assign [W c n ] * . Take the product of all the vertices functions. Vertexes joined with an arrow have a coordinate in common. Integrate over all the free coordinates. 3. Divide by the symmetry factor of the graph, that is to say, the number of permutations of vertices and legs that leave the graph unchanged. 4. Sum over all possible connected graphs.
Diagrammatic rules to compute the correlators (coordinate space): 1. Draw a graph: disconnected graphs are allowed, but each connected block must have at least one external leg. Assign a coordinate to each external leg. The remaining rules are the same as for the free energy. Notice that, if by chance we want to compute the connected part of a correlator, we must consider only connected graphs. 2. The correlator to be computed is determined by the external legs. Suppose we have n outgoing external legs with coordinates y i and m incoming external legs with coordinates x i , then we are computing Integrate over all the free coordinates, that is to say, the coordinates associated to internal legs. 4. Divide by the symmetry factor obtained exactly as in the free energy case. Notice that the external legs are fixed and they cannot be interchanged. 5. Sum over all the graphs so constructed. This means summing over the distinct topological configurations, but also on the possible permutations of the external momenta that change the graph.
Equivalent rules hold in the momentum space. They are valid in the thermodynamic limit where we can substitute summations in the momentum space with integrals. Translational invariance of W c n implies a conservation law of the momenta at each vertex.
Diagrammatic rules to compute the free energy (momentum space): 1. Draw the same graphs as in the coordinate space case. 2. The contribution of each graph can be computed as follows: assign to a vertex with n outgoing arrows the function K n , to a vertex with n incoming arrows assign [K n ] * . Take the product of all the vertices functions. Vertexes joined with an arrow have one momentum in common. 3. Impose the conservation of momenta at each vertex. Integrate over all the remaining free momenta, the integration measure for a momentum k is d d k/(2π) d . 4. Divide by the symmetry factor of the graph, that is to say, the number of permutations of vertices and legs that leave the graph itself. Multiply the obtained value by L d : this factor comes from the invariance under translation of the entire graph. 5. Sum over all possible connected graphs.

Diagrammatic rules to compute the correlators (momentum space):
1. Draw a graph exactly as in the coordinate space case. The interpretation of vertices and arrows is the same as in the computation of the free energy in momentum space. Assign a momentum to each external leg. For simplicity consider connected graphs: the value of disconnected graphs will be simply the sum of its connected parts. 2. Suppose we have n outgoing external legs with momenta k i and m incoming external legs with momenta q i , then we are computing The L factor is necessary because of the normalization of (4). Impose the global conservation law of momenta: Integrate over all the free momenta, that is to say, the momenta associated to internal legs. A momentum k has as integration measure d d k/(2π) d . 4. Divide by the symmetry factor obtained exactly as in the free energy case. Notice that the external legs are fixed and they cannot be exchanged between them. 5. Sum over all the graphs so constructed. This means summing over the distinct topological configurations, but also on the possible ways of permutations of the external momenta that change the graph.
A couple of examples of the computations in momentum space are given in Fig. 5. The rules in coordinate space hold without any particular hypotheses on W c n , in particular we do not need to suppose (21). We will use them in Section 4.1 to explore the consequences of Properties 1, 2 and 3 on the analytic properties of W c n . Instead the rules in the momentum space can be more suitable to some explicit calculations, as we will see in Section 5 and Section 6. Notice that we can also generalize these rules to compute overlaps between cluster expansions.

Diagrammatic rules to compute overlaps:
Consider | and | two cluster expansions and we want to compute | , or even insert operators between these states. The rules we explained in the case of a single state keep always distinct the vertices of the bra and the vertices of the ket by means of the direction of the arrows. If we want to differentiate the two we must simply associate to the vertices with outgoing arrows the cluster amplitudes of | and to the vertices with ingoing arrows the functions describing |, once we take the complex conjugate, of course. All the remaining rules are exactly the same. Each graph with N vertices can be thought as a subgraph with N − 1 vertices eventually attached to a single vertex that has always at least one external leg. In this figure is drawn a typical example, the number of legs of the figure has not a precise meaning and want to represent a typical situation, also the nature of the selected vertex can change and have only ingoing arrows instead of outgoing. It could be also that the vertex and the subgraphs with N − 1 are not really attached and form two disconnected pieces, what really matters is that we can always select a vertex with at least one external leg.

Thermodynamic properties of the graphical expansion
In this section we will analyze the relations between the Properties 1, 2, 3 and the features of the cluster amplitudes. In principle we would like to find necessary and sufficient conditions on the cluster amplitudes in order to have the desired thermodynamic properties, but we are able to properly address only one implication and provide a sufficient condition for the validity of Properties 1, 2, 3 at each order of the expansion. Assuming (21) i.e.
and that W c n has no explicit dependence on the system size is sufficient in order to satisfy Properties 1, 2, 3 at any order of the expansion. It could be questioned whether or not the graphical expansion converges and the thermodynamic properties of each graph imply Properties 1, 2 and 3 on the true correlators and free energy. If the expansion does not converge any longer, it could be that Properties 1, 2 and 3 are still valid, but not respected individually by each graph. Indeed, as we argue in Appendix D, in critical systems the expansion becomes divergent, even though Properties 1, 2, 3 are expected to be still valid on the correlators and free energy; nevertheless we argue that (28) should be still valid.

Property 1
We now show that (28) guarantees that any graph that appears in the expansion of correlation functions is well defined in the infinite size limit. Note that a generic graph for the correlators is nothing but a complicated convolution of the W c n functions, therefore we need to study such an operation. The proof is by induction on the number of vertices that contribute to a graph: in each graph with N vertices we can identify a subgraph with N − 1 vertices linked to a single W c n vertex (or its complex conjugate) with some external legs attached to it (Fig. 6).
Then, some of the external legs should depart from W c n . Suppose { x 1 , .., x j } are associated with external legs of W c n , then the remaining x i coordinates will be associated to external legs departing from the subgraph with N − 1 vertices. The internal legs connecting the vertex with the subgraph will be denoted with coordinates y i : in Fig. 6 we drew the case of only one leg connecting them, but there could be more than one or even none. We denote with G N −1 ( y 1 , .., y n−j , x j +1 , .., x l ) the value of the subgraph with N − 1 vertices. The function G N is nothing but a convolution between W c n and G N−1 : Above, we inserted an unimportant symmetry factor S, that depends on the details of the connection between W c n and the subgraph. Now, we show that if G N−1 is uniformly bounded, then provided (28), G N exists and it is also uniformly bounded. Note that the graphs with only one vertex are simply the W c n functions that are uniformly bounded because of (28). Then, by induction we can say that each graph of the correlators is well-defined in the infinite size limit. Therefore, suppose |G N−1 | < M with M an unknown, but finite number, then: Because of (28) the integral in the above expression is finite and uniformly bounded for any of the remaining x coordinates. Then G N exists, because the integral in (29) is absolutely convergent and moreover G N uniformly bounded. This guarantees that each graph in the expansion of the correlators is finite in the L → ∞ limit.

Property 2
We now show that (28) is enough to guarantee Property 2 (i.e. cluster property) at any order of the expansion, therefore we will show that each connected graph decays to zero whenever the distance of two of its coordinates is sent to infinity. Actually, we will prove a stronger statement: any connected graph satisfies a similar condition to (28). Consider G N ( x 1 , ..., x l ) a generic connected graph formed by N vertices, then where the integral is meant to be over the coordinates in the subset A. The convergence of the integral holds for any bipartition of the set of coordinates such that B has at least one element, therefore there is always at least one coordinate that is not integrated. Of course the convergence of (31) implies that the value of the connected graph decays to zero whenever the distance of two coordinates is sent to infinity, apart from a possible zero-measure subset. The proof of (31) is by induction: note that any connected graph can always be seen as two smaller connected graphs joined together (apart the graphs with only one vertex formed by only W c n or the complex conjugate vertex). Therefore, assume G N can be divided in two connected graphs GÑ and G N−Ñ with Ñ < N; by re-labeling the coordinates we can assume that the first j coordinates x 1 , .., x j belong to GÑ . The value of G N can be written similarly to (29) where S is an unimportant symmetry factor, the coordinates y are associated to the internal lines that connect the two subgraphs: since the graph G N is connected, the two subgraphs must be linked through at least one arrow. We can now write the inequality: ., x j , y 1 , .., y n−j )| |G N−Ñ ( y 1 , .., y n−j , x j +1 , .., x l )| (33) Since the integration domain A is never the whole set of coordinates x, at least one among the two graphs GÑ and G N−Ñ possesses a x coordinate that is not in A. Without loss of generality we can suppose that at least one of the x coordinates in GÑ is not an integration variable, then consider G N−Ñ . Because of the inductive step, G N−Ñ satisfies (31): for some constant M. Note that in the above we are considering only the integration over the x coordinates that are both in A and variables of G N−Ñ , instead the y coordinates are not integrated. Note that we can use (31) because, since G N is connected, there is at least one internal leg that connects GÑ and G N−Ñ , therefore there is at least one coordinate y. Using this inequality in (33) we have: Now we can use that GÑ satisfies (31), since in the above integration GÑ possesses at least one x coordinate that is not an integration variable. Using (31) on GÑ it follows that also G N fulfills (31). Therefore we have shown the inductive step and this concludes the proof. Note that we can also say something about the Fourier transform of G N . At the end of Section 3.2 we showed that the Fourier transform of the cluster amplitudes that satisfy (28) are non-singular, analogously the condition (31) implies that the Fourier transform of G N has no singularity (after we factorize the overall δ function caused by the translation invariance). This simply means that if (28) is valid then none of the graphs of the correlators computed in the momentum space has any singularity.

Property 3
The proof of this implication is simple and relies on our graphical expansion. We start with an interesting mapping between the graphs that appear in the free energy and the graphs appearing in the two point correlator ψ † ( x)ψ( y) . Pictorially it is very simple: consider a graph of the two point correlator ψ † ( x)ψ( y) . Its external legs must be an outgoing arrow and an ingoing one. If we join these two legs, we get a graph for the free energy. We can also proceed in the reverse direction taking a graph of log Z[ ] and cutting an internal line to get a graph of ψ † ( x)ψ( y) . This is shown in Fig. 7. It should not be forgotten that this mapping does not properly take in account the symmetry factors that can change during the gluing procedure, but since we are merely interested in the system size dependence of each graph it is sufficient for our purpose. This correspondence is enough to show that any graph of log Z[ ] is extensive. Because of the invariance under translation, each graph of the two point correlator can depend only on the difference between the two points, thus the integration in Fig. 7 becomes simply L d G(0, 0). Now, using the fact that each graph of the correlators does not depend explicitly on the system size, we have that G(0, 0) is independent of L, thus each contribution to the free energy is extensive. It must be said that G(0, 0) may be a divergent quantity when the UV cut off is removed, but this does not change the extensivity since the latter refers only to the scaling with the system size, i.e. the infrared cut off.

Explicit calculation of the cluster amplitudes of ground states
In the previous sections we have studied the general properties that a cluster expansion of a local state should respect. Now we shall link our abstract discussion to some explicit models. Even if the general construction can be applied also to condensed matter systems with a finite number of particles (Appendix A), from now on we will focus mainly on relativistic theories, in which we avoid the subtlety of dealing with a finite number of particles. In particular we will consider a real bosonic field φ whose dynamics is governed by a relativistic Lagrangian with local self-interactions.
Above, β n are the coupling constants of the theory, λ is a parameter to tune the strength of the interaction, that is inserted to make eventual perturbative calculations more straightforward. As we stated in the introduction and further investigated in Section 3.1, the ground state (or ground states, if many of them are present) of such a theory is supposed to satisfy some locality requirements encoded in Properties 1, 2 and 3. Even if most of the content of this section can be extended to the case of many degenerated ground states, in the following we will imagine β n to be tuned so that there is only one ground state. Our physical properties have been expressed in terms of some abstract fields ψ ψ † : the first step is to correctly identify these fields. In the canonical quantization scheme, the field φ is divided in terms of the modes that diagonalize the free Hamiltonian. In finite volume this decomposition is expressed by: Where a k and its conjugate follow standard commutation rules. The natural thing to do is to define the ψ field from the a k field. Therefore we invert relation (4) and use it as the definition of ψ: Proceeding in this way, the cluster expansion in the momentum space (9) (11), will be simply an expansion in the Fock space of the theory constructed acting with the creation operators of the free theory on the bare vacuum. The Hamiltonian derived from (36) is: Above, we ignored the zero point energy by normal ordering the free part. A cautionary remark must be made about Properties 1 and 2: the locality requirements are written in terms of the fields ψ ψ † , but the physical fields that should respect them are φ and its conjugated momentum. It turns out that, at least in the massive case m = 0, it is pretty simple to show that if ψ satisfies Properties 1 and 2, then also φ and its conjugated momentum satisfy them. Also the inverse implication is true. As a matter of fact, combining (37) and (38) we get: The kernel is exponentially decaying D( x − y) ∼ e −m| x− y| from which it is straightforward to verify our claim. From our reasoning, ground states of these theories should have a cluster expansion as (9), with regular K n functions. In Section 5.1 we verify that the ground state possesses indeed a proper cluster expansion through standard perturbative calculations and describe the suitable method to calculate K n in a perturbative way. At least in the massive case, such a computation confirms the absence of singularities in K n at each order in perturbation theory. In the massless case (m = 0) the naive perturbation theory does not work and we attempt a nonperturbative check of the absence of singularities in Appendix E, confirming that at least K 2 has no singularities. In Section 5.2 a set of exact functional equations for the cluster amplitudes K n is presented: they provide a consistency check about the analytic properties of K n , reproduce the results of standard perturbation theory of Section 5.1 and give access to other approximation schemes or even to some exact predictions, as in Appendix E.

Perturbative computation with Feynman graphs
This section is devoted to the construction of a perturbative method to compute the ground state of an Hamiltonian such as (39). In the following we will suppose that the ground state is unique and 'perturbatively close' to |0 , that is to say the ground state in the limit λ → 0 reduces to |0 . These assumptions give constraints on the allowed form of the interaction, but it is still a rather general theory. The Hamiltonian of the system can be split in a free and an interacting part: Given the ground state |G , we can project it on a state with n particles at definite position x i and find the wavefunction W n : Now we should extract from both sides their connected part, which is the cluster amplitude that enters in (7). Going to Fourier transform, we find: The L factors come from (9) and are such that when L → ∞ the right hand side becomes L independent. The Kronecker delta is due to the translation invariance of the state. Instead, the physical requirement of Property 1, 2 and 3 are completely encoded in the analytic properties of K n . In particular there should be no extra δ factors, and no other singularities are expected: we will check this claim in the next computations. For simpler notation, it is better to compute the complex conjugate of (43): The next steps will bring the right side of the above expression in a form more suitable to perturbation theory. To compute the right term side we can use the same trick as in (19) to obtain the ground state. As a matter of fact, applying a thermal ensemble on the bare vacuum and then taking the zero temperature limit we are able to extract the ground state.
Above, we explicitly used our normalization 0|G = 1. We can use this trick to compute the right side of (44).
Above, N G = G|G . Instead of using the Hamiltonian H , it is more suitable to switch to the interaction picture defined as: where T exp is the τ ordered exponential and V I is the potential in the interacting picture: Using the fact that 0|H 0 = 0 we arrive at: To proceed further, it is useful to define the field φ in the interaction picture and in Fourier space: so that we can write: We can express (49) in terms of φ I : where : : stands for normal order. Notice that from the relation above we can read: Therefore we arrive at: Notice the strong resemblance of the expression above with the starting point of standard adiabatic perturbation theory [45], but there are some important differences. Apart from being in euclidean time, notice that the integral runs from 0 to β, not from −β to β as usual. This fact has strong consequences in the definition of the diagrams, as we will see. Now we should expand the exponential in (54) and perform the Wick contractions. To do this, it is useful to express the time ordered propagator via an auxiliary momentum along the euclidean time direction: Now we can evaluate (54) via Feynman diagrams. The rules of the game are quite the same as those we encounter in standard textbooks [41,[44][45][46] to compute correlators, with some crucial differences. Since the derivation of these Feynman rules is a standard computation, we merely stress what changes in comparison with the standard references and then state the rules. We can anticipate that, as usual, the denominator of (54) has the effect of eliminating the disconnected diagrams with no external legs in the numerator. Moreover, relation (44) tells us that, in order to obtain K n , we should compute only the connected part of (54) which, at the end of the story, will be represented by connected Feynman graphs. Notice that, since we are working in an 'euclidean time' formulation, the propagators have no poles (at least in the massive case) and we do not need to introduce any pole prescription. So far nothing is different from usual Feynman graphs, but here is the crucial difference: we do not have conservation of energy at each vertex. Conservation of energy in [45] comes out from the integration over the whole time axis in the analogue of (54), while in our case the integration is over the half line (0, β → ∞) and this has crucial consequences. Imagine that we expand (54) and contract the fields; then using (55) we can separately integrate along the time direction of each vertex and we encounter terms of the form: The last equality is true in the distribution sense, in the limit β → ∞, σ i depends on the contraction we are doing, but for now their value is not important. Notice that, if in (56) instead of integrating between (0, β) we integrate between (−β, β), then we would get a δ-function instead of (56): this is what happens in [45]. Having in mind these differences from [45], we can write down the Feynman rules: for simplicity we implement the L → ∞ limit substituting summation in the momentum space with the proper integrations. To avoid confusion with the graphs of Section 4 the propagators will be indicated by dotted arrows (Fig. 8) and the vertices as simple points.
Diagrammatic rules to compute K * n : 1. To compute (54) draw a graph with n external outgoing arrows and assign (k 0 i , k i ) i=1,...,n to them. It is important to draw only outgoing arrows, otherwise we would spoil the vertices convention. If we are interested only in K n , we draw only connected graphs in order to extract the connected part of (54). If instead we want the full expression (54), disconnected graphs are allowed, but each disconnected block must have at least one external leg. Each term in (54) has a global L −d(n−2)/2 factor that exactly cancels with the same factor in (44), as it should. Therefore, we can ignore it. 2. Each vertex with s external lines labeled by s momenta {(q 0 j , q j )} j ∈{1,s} counts as: with → 0 and σ i is chosen −1 if the arrow is ingoing and +1 otherwise. 3. For each vertex of s legs we must impose a momentum conservation law s i=1 k i = 0 where σ i is chosen as at the previous point. Notice that the conservation refers to the spatial momentum while the k 0 j momenta remain completely free. 4. Integrate over all the free momenta with integration measure for a component q j equal to dq j /(2π). Notice that we must also integrate over the zero component of the momenta of the external legs: only the spatial external momenta are fixed. 5. Divide each graph by the corresponding symmetry factor, that is the number of permutations of vertices and of internal legs that leave the graph unchanged. 6. Imposing the conservation laws at each vertices we will end up with a global conservation law over the external momenta, as expected.
7. Sum over all the connected graphs to obtain K * n ( k 1 , .., k n )/ n i=1 2E( k i ). Notice that a connected graph has only a global conservation law n i=1 k i = 0, therefore the momenta are not conserved in 'separate groups'. Since in (43) we already took into account the global conservation law, K n does not contain any additional δ factor, as expected from our locality constraints.
These rules are best suited for findings the functions K n as power series in the coupling constant. As an example, we focus on the φ 4 theory: Fig. 9. Graphs for K 4 and K 2 at first order in λ in φ 4 theory. Their values are given in (59) (60).
We can readily see that the graphs that contribute to K n at first order in λ are those of Fig. 9.
We can easily perform the integration over the momenta along the time direction to get: All other K n with n > 4 are of order higher than O(λ). All K n with odd n vanish due to the symmetry of the interaction under φ → −φ. We can easily count the minimum number of vertices needed in order to construct the graphs that contribute to the computation of K 2n : they correspond to a connected graph with no closed loops, i.e. a tree-like graph with 2n external legs. Thus the order of K n in perturbation theory is: More generally, for an interaction of the form φ j , there is a topology relation between the number of external legs n, the number of the interaction vertices V and the number of internal loops L: This can be easily derived from the Euler's polyhedral formula and can be found in [45]. Using this formula it is not difficult to pin down the first order in λ for each K n , since it is simply the number of vertices of the graph. We get K n ∼ λ (n+2L min −2)/(j −2) where L min is the minimal number of loops needed to draw the graph with n external legs. Notice that, if m is even, then we cannot construct any graph with an odd number of external legs. Truncation in some order of λ means truncation in the functions K n : at order O(λ α ), only K n with n ≤ n max = (j − 2)α + 2 (65) are non-vanishing. As a self-consistency check of this diagrammatic technique, it can be verified that the φ field correlators computed using the above results for K n and the rules of Section 4, agree with the zero temperature correlators derived from standard perturbation theory. This is the content of Appendix F. As expected, some of the diagrams may involve integrals that diverge for large momenta, as for example in (62): this is due to the ultraviolet divergences of the theory and in order to make sense of them we should implement a renormalization scheme, which must be compatible with the computation of the field correlators (see for example [45]). This procedure is implemented with the insertion of the proper renormalization counterterms in the Lagrangian with the net effect of modifying the definition of the coupling constants. As we already stated, we will not proceed along this path and we rather assume that integrals are finite because of a UV cutoff. Notice that, as long as m = 0, we do not have any singularity in each diagram we can draw, therefore K n have no singular behavior: this is a non-trivial check of the condition (21).

Exact functional equations from the Schrödinger equation
The method of Section 5.1 is a perturbative method mostly suited for an expansion of K n in terms of λ, even thought a careful treatment of analogous expansion could lead to nonperturbative results, for example the mass renormalization in quantum field theory [45]. To obtain such information at least a partial summation of the graphs is needed, but this is not an easy task, mostly because the of lack of energy conservation at the vertices. Actually, there exists another way to implement a non-perturbative study of the coefficients K n : we can derive a set of exact equations that K n must satisfy in order for |G to be the ground state of the Hamiltonian. It is worth to anticipate that the general structure of these equations is consistent with the expected analytic properties of the K n functions. Moreover these equations allow us to find the λ expansion of K n in a systematic way, avoiding the graphical technique of Section 5.1, but consistently with it. Moreover, if in Section 5.1 we must ask that the ground state is unique in order to extract it from the vacuum through (46), the equations we are going to use do not use this fact, therefore can be used in principle to study systems with degenerated ground states: of course, the hypothesis in which (2) can be written must still be valid, therefore we need a non-zero overlap of the state with the free vacuum of the theory.
We start with the observation that the ground state is simply an eigenstate of H : To solve the eigenstate equation we can directly impose a solution |G in the form of the cluster expansion (7), so that we would obtain a set of equations for the unknown functions K n . Since each translational invariant state can be written as (7), among the solutions of (66) there are all the translational invariant eigenstates, but since we are interested in the ground state we can look for a solution only in the set of local quantum states with the precise analytic properties of the cluster amplitudes we identified in Section 3. At this point we could have more than one solution: it would be interesting to study which eigenstates different from the ground state are local quantum states. In principle, to select the ground state we should compute the energy of all of them and choose the lowest, but if we make in the same hypothesis as in Section 5.1 we can isolate the correct ground state with the requirement that, for λ → 0, it becomes the ground state of the free theory, i.e. the vacuum.
As before we assume an Hamiltonian of the form (39), but without loss of generality it is easier to suppose that it is normal-ordered, defining new coupling constants β n : In the above : : means that we consider the normal-ordered expression, with all the creation operators on the left of the annihilation operators. The action of H on |G can be computed using the observation that a k , because of the commutation rules, acts as a differentiation with respect to a † k . Since |G is an exponential state, acting on it with as many operators a k i we want, we will always obtain a state of the form |G , with an operator containing only a † operators. For example, if we consider a − k |G we have: With this, we can see immediately that: Dealing with the interaction term of the Hamiltonian is more involved, because we do not have only a derivative: The combinatorics of such an expression is difficult, because |G contains all the cluster amplitudes K n and having more than one derivative links them in a non-trivial way. In Section 5.2.1 we show a graphical representation to deal with this complicated combinatorics. After this cumbersome computation, we can write H |G as: In the second line we have explicitly written the simplest part of (71) that does not involve derivatives and all the remaining action of the interaction is contained in F n . This is a complicated object involving sums and convolutions of K with the interaction terms as we are going to explain in Section 5.2.1. Now we should impose equation (66) and we obtain: Note that we used a factor L d to put the ground state energy inside the summation: what is contained in square brackets is now well-defined in the L → ∞ limit, since the ground state energy is extensive. From these equations we can extract directly equations for the cluster amplitudes K n : Notice that, if equations (74) are satisfied, then also (73) is true. It is also possible to show the inverse using suitable test states. The argument is this: applying 0| to (73), all the terms with n > 0 will automatically cancel; only the n = 0 term survives.
Since 0|G = 1, we have that such an equation is satisfied only if the equation n = 0 of (74) is satisfied. Now, move on with another test state containing one particle, that is to say 0|a k : applying it on the left of (73), knowing that the term n = 0 is zero and using the commutation rules, we have that only the term with n = 1 is not annihilated by this test state. This time, we would find that the equation (74) with n = 1 must be satisfied. At this point using a test state with two particles, we will find n = 2 out of (74), then proceed step by step increasing the number of particles: with this idea, it is immediate to see by induction the equivalence of (73) and (74). Notice that in (74) there are no extra δ-like contributions, which would be inconsistent (there are no extra δ hidden in F, as it will be clear in Section 5.2.1). Moreover, since there are no explicit singularities in the equations (at least for m = 0), we should expect a smooth behavior for K, without any singularity. Some care must be taken in the massless case in which E( k) = | k| and therefore poles arise; we discuss about this fact in Appendix E. These features of the equations can be regarded as a consistency check for the analytic properties of the cluster expansion. It should also be explicitly said that these equations could also been used to compute the ground state energy, simply taking the n = 0 term in (74): In the above, F 0 is completely fixed once we determine K n through the other equations: this relation is particularly important if we have more than one solution of (74), since in order to find the ground state we must choose the solution that minimizes E G . Notice that starting from the solution for λ = 0, these equations can be solved recursively to obtain an expansion of K in powers of λ: the first order solution can be obtained simply neglecting F.
In this way we have that K is of first order in λ and λF is of order at least λ 2 , since F contains itself at least one K. Inserting (77) in F we can update the right hand side of (74) and find new K corrected up to λ 2 . We can then iterate the procedure: this is a systematic way to obtain exactly the same perturbative expansion as Section 5.1. As a matter of fact in the φ 4 model, once we normal order the Hamiltonian obtained from (58) we get the coupling constants: The other coupling constants β n are zero. Inserting (78) in (77) we obtain again (61) and (62), as expected.

Diagrammatic representation of the functional equations
We now want to present a simple graphical way to deal with the combinatorics of (71) and compute F n of (72). We invite the reader to check its consistency with the explicit computation of some terms through (71). We still use arrows to connect vertices and in this graphical approach we need two types of vertices, one to represent K n and one to represent β n : the vertices related to K n are simple dots, instead β n is represented by a circle (Fig. 10). The following rules hold in the thermodynamic limit, where we can replace sums with integrals.
Diagrammatic rules to compute F n : 1. A graph associated to F n is made by one and only one circle representing the interaction vertex and a number of dots representing the K n amplitudes. 2. Each dot must be connected to the circle by an arrow, which automatically implies that the graph is connected. Arrows can only point outwards from dots, but at the circle they can be either ingoing or outgoing. Following this recipe, the dots related to K can never be directly connected. 3. To get F n ( k 1 , .., k n ) we must have n external legs with assigned momenta { k i } i=1,..,n .
All the external arrows must point outwards: this is necessary to respect the sign convention we adopt. 4. Vertices joined with an arrow have the same momenta; an internal arrow with momentum q enters as q in the respective K. 5. At each vertex, dot or circle, we must impose the momentum conservation. Notice that regarding the β circle we must take in account also the direction of the arrows. (For example in Fig. 10 the requirement is n i=1 p i = m i=1 q i .) 6. After we impose all the conservation laws we have a global conservation law on the momenta and no other restriction over them, as it should. We must integrate over all remaining free momenta. A free momentum q has integration measure d d q/(2π) d . 7. Divide by the symmetry factor, that is to say, the number of permutations of vertices and internal lines that leave the graph (topologically) unchanged. External lines are fixed ant they cannot be exchanged. 8. Sum over all possible graphs. This means sum over all the possible topologies and over all the possible permutations of the external momenta that change the graph. To avoid confusion, once we fix the topology of the graph we will always mean that we are also summing over all non-equivalent permutations of external legs. Fig. 10. With a blank circle we denote the interaction vertex between the Hamiltonian and the cluster amplitudes, its value is n+m . Fig. 11. Example: this graph enters in the computation of F 3 , when the coupling constant β 4 is present, therefore it appears in the φ 4 theory.
As an example, we are going to see in detail the calculation of Fig. 11, that contributes to F 3 . We assign the external momenta and arrived at point 5 we have: The free momentum is q and it must be integrated out; the symmetry factor is 1/2 because we can exchange the two internal lines. At point 7 we get: Now we must sum over the distinct permutations of the external momenta, therefore after point 8 we get: Now the quantity f Sym 3 ( k 1 , k 2 , k 3 ) is completely symmetric in its argument and it is the contribution for F 3 ( k 1 , k 2 , k 3 ) generated by such a graph. The expressions for F n are quite involved and usually analytically intractable; there is only one case in which we can analytically solve them and it is the case of a gaussian state. This is something expected, since we know that the gaussian state can be explicitly solved by other methods (Appendix B) and we can use this special case as a check for our method. The quadratic case is such that in (67) we have only β 2 and β 0 different from zero, without loss of generality we can redefine the λ parameter in order to have β 2 = 1. With this choice, the Lagrangian we are considering is simply Fig. 12. Graphs representing the two contributions to F 2 keeping only K 2 = 0, to obtain F 2 we must sum them.
while the normal ordered Hamiltonian is: β 0 does not matter for what follows, since an additive constant does not change the ground state, but only its energy. In Appendix B we study such an Hamiltonian by standard methods and find its ground state is a squeezed state with K 2 given by (B.9); here we want to obtain the same result through our equations. From our equations, it is not difficult to understand that only K 2 will be non-zero. Consider the perturbative solution described in Section 5.2: at first order the only non-zero term is K 2 , as it is clear from (77). Using now this solution in the next iterative step, we should compute F n with the first order solution (77). Using the diagrammatic rules to compute F n we see immediately that, thanks to the fact that only K 2 is not zero at first order, only F 2 and F 0 are not zero at the same order. At the next order, F 2 corrects K 2 as it is clear from (74), but all K n =2 remain zero: repeating the same reasoning we see all K n =2 remain zero at any perturbative order. Thanks to this reasoning, we know that the solution we are looking for possesses only K 2 different from zero. Using this fact we can now write the exact equations that K 2 must satisfy, then exactly solve them. From the topology of the graphs it is clear that only F 2 will be non-zero and has two different contributions (Fig. 12). Using the rules we explained, we can write F 2 as: Notice that in the first graph of Fig. 12 we get two different graphs if we exchange the momenta of the two external legs, but the value is unchanged: this gives a multiplicative factor of 2 to the linear term in K 2 , correctly counted in the expression above. Now we should insert this quantity in (74) and we finally get the equations for K 2 : These have two solutions: We can identify the correct solution requiring that, for λ → 0, the ground state becomes the vacuum, therefore we must have K 2 → 0 and so we must choose the plus sign. With a little bit of algebra, the expression can be rearranged as: This is exactly expression (B.9) of Appendix B, as it should. As we anticipated at the beginning of this section, we have found also another solution to the eigenstate equations, but note that the state we would have obtained choosing the minus sign in (86) is such that |K 2 | ≥ 1 and this is an unphysical state, as we explain in Appendix B.
After this check in the gaussian case we would like to use these equations to study a real interacting field theory, like the φ 4 model. The equations in this case become involved and we do not attempt an analytical solution of them, nevertheless they are suited for different approximation schemes, as we are going to discuss, and we can even acquire precious, non-perturbative information through them, as we will do in Appendix E.

Approximation schemes
In this brief section we want to sketch some possible approximation schemes that can be used to get the cluster amplitudes of a ground state, having in mind the example of a relativistic field theory (such as (36)): we do not want to be quantitative and simply comment on the possible validity and problems of these approximated methods, leaving their development for further studies.

Naive perturbation theory:
This is the simplest method that we can think of: we can use the perturbation theory of Section 5.1 or, equivalently, solve recursively in λ the functional equations of Section 5.2. Solving up to order O(λ α ) we get an approximation to K n that provides a natural truncation in the cluster expansion, as commented at the end of Section 5.1. The appealing fact is that such an approximate method is straightforward to perform, but we know that often naive perturbation theory lacks of predictive power and we need to go beyond it. This truncation scheme, from the viewpoint of the functional equations, can be generalized also to the dynamical case: in Section 6 we will see that the evolution of a local quantum state with a local Hamiltonian follows a simple generalization of the functional equations.

Perturbation theory starting from Hartree-Fock approximation:
This method represents the first improvement of the naive perturbation theory: we focus on the φ 4 theory as an example, but can be done in general. The idea is simply to separated in the action a quadratic part that captures as much information as possible, in a consistent way.
To do this, we add and subtract terms in the φ 4 Lagrangian (58) and write: This is exactly the same Lagrangian as (58), in which we added and subtracted φ 2 in a proper way: now, in square brackets, it appears a quadratic Lagrangian with a shifted mass. The unknown expectation value φ 2 must be computed selfconsistently. With this trick we have a better control on the interaction term: now the small parameter is not the bare λ parameter anymore, but rather the variance of the square of the field (φ 2 − φ 2 ) 2 . The crudest approximation is obtained neglecting the latter term in the action and is known as Hartree-Fock approximation [47]. Since the latter term would be zero if the ground state were gaussian, such a method is valid as long as the gaussian state is a good approximation. To proceed further we should write the eigenmodes of φ using the shifted mass; we call ã k these new operators which are related to the old a k and a † k simply through a linear transformation, known as Bogoliubov rotation.
In the above, the shifted mass m is defined as: The Bogoliubov transformation in terms of these parameters is: The Hilbert space will be constructed out of a new vacuum |0 that is annihilated by all the ã k . Then all the construction of Section 5 can be repeated simply replacing a k and m with the new operators and mass. We can use Section 5.1 or Section 5.2 to perform naive perturbation theory as explained before, but notice that such an expansion is not perturbative in λ due to the shifted mass: this approximation goes beyond naive perturbation theory.
After we obtain the proper truncated approximation of K n , the expectation value φ 2 must be computed in a self-consistent way. We can formally compute φ 2 with the methods of Section 4 and Appendix C: the graphical expansion of φ 2 must be truncated to the proper order, consistently with the approximation of K n . In the latter expression φ 2 would enter as a parameter and therefore we would have obtained a set of self-consistent equations for φ 2 that have to be solved. Such an approximation scheme has the nice feature of going beyond the perturbative approach; moreover through the generalization of the functional equations presented in Section 6 it can be also applied to study the time evolution of local states. Unfortunately, this method has two main drawbacks: the self-consistency equations for φ 2 will become more and more involved as soon as we proceed further in perturbation theory in the ã k basis. Moreover, the cluster expansion will be expressed in terms of ã k : if, by chance, we need the cluster expansion in terms of a k (for example, in Section 6.2 we focus on quenches from an interacting to a free Hamiltonian and the time evolution is trivial in terms of the a k operators) it is not an easy task to switch from one language to another.

Variational approach:
The above methods use the already found equations. Here instead we want to describe a way to derive new equations from a variational principle. The ground state is the state of minimum energy, therefore it is the state that minimizes the expectation value H . As explained in standard books [48], one way to approximate the ground state is to minimize the energy over a selected set of trial states. Within these trial states, the state with the lowest energy is the best approximation for the ground state. At this point, we can use our knowledge about the generic properties of a ground state: it is a local quantum state and we know the form of its cluster expansion. Let | be a generic local quantum state with cluster expansion (2); the ground state |G is identified by the set of functional equations: In the above, H is the expectation value of the Hamiltonian evaluated on the local quantum state | , the derivative is a functional derivative where W c n and [W c n ] * must be thought of as independent quantities. Equations (92) (that can be written also in momentum space with K n as well) express only the stationarity of the energy with respect to variations of the state: if more solutions are present, in order to find the ground state we should keep the solution with the lowest energy. So far if we were able to perform the above procedure, we would obtain the exact ground state, since it belongs to the set of local quantum states: the approximation comes out because we cannot compute H in an exact way. For example, we can focus on the general Hamiltonian (39): computing its average reduces to computing multipoint correlators of a k and a † k , that we can compute only through a graphical expansion. Therefore the best we can do is to approximate H through a truncation of the set of its graphs. Actually, there could be different meaningful truncation schemes; the easiest one we can imagine is to implement corrections to the gaussian state: if the gaussian state is a good approximation, then K n>2 must be small. Using the rules of Appendix C we can organize the graphs of H in terms of how many K n>2 vertices they contain and truncate the expansion. Notice that we can exactly take into account the two point correlators through the resummation of one particle irreducible graphs: of course, the sum of irreducible graphs appearing in (18) must be truncated as well. Using this approximation for H we would find an approximation to the equations (92): once we find the solution, we should check self-consistently that K n>2 are actually small, as well as the neglected terms in H . Such a method has the appealing feature that it can give non-perturbative information; moreover if by chance we know some further properties of K n we can use them from the beginning in our trial local state. As a drawback we must understand the proper way to truncate the expansion H and to solve the set of complicated equations coming from its variation. Another appealing fact is that such a variational procedure can be applied also to systems whose ground state does not have a truly cluster expansion, but the latter represents only a faithful thermodynamic description of it as we explained at the end of Section 3.2: such a variational method could improve our knowledge of the ground states of many interesting condensed matter systems, for example the interacting Bose gas [43].

Dynamics of the cluster expansion and application to quenches
So far we have studied only static properties of local quantum states; we shall now turn to the study of dynamics. In this section we ask ourselves what happens to a system that evolves with a local Hamiltonian, provided it is initialized in a local quantum state, which is the case in many out of equilibrium situations.
More specifically, we consider a thermodynamically large closed system lying on the ground state | of some local Hamiltonian H 0 . At time t = 0 the Hamiltonian changes abruptly from H 0 to H and the system is let to evolve for a long time under the new Hamiltonian. This process has been termed "quantum quench" [14]. Even though the system remains always in a pure state, for large times expectation values of local observables may become stationary. This is a consequence of the interference of incoherent quasiparticles produced at the time of the quench and coming from distant points. Any finite subsystem of the infinite system exhibits such stationary behavior, meaning that the reduced density matrix of any subsystem becomes stationary. The complement of the subsystem can therefore be seen as a 'bath' in contact to the subsystem, which drives it to equilibrium.
Such equilibration is not however necessarily thermal. If the evolution is constrained by extra conserved quantities, apart from the total energy, then the system is conjectured to equilibrate to a 'Generalised Gibbs Ensemble' (GGE) [23] which maximizes the entropy subject to all conserved quantities [49,50]. This is the case when the evolution is done under an integrable Hamiltonian, which is characterized by an infinite set of local conserved charges. Free (i.e. gaussian) Hamiltonians are the simplest examples of integrable models. Non-trivial integrable models are the one-dimensional Bethe Ansatz solvable models. The lack of thermalization in a one-dimensional integrable system was demonstrated in a seminar experiment in ultra-cold atoms [15], while the experimental observation of a GGE was recently achieved [22]. The GGE hypothesis has been successfully tested in a large number of quenches with integrable evolution, while in some cases 'quasi-local' conserved charges [51][52][53] must be included apart from the local ones, in order for the GGE to describe correctly the large time limit [32,54,55].
The role of clustering properties of the initial state in equilibration has been first recognized in [27,56] (see also the more recent [57,58]) in the context of lattice systems, where it is shown that the system equilibrates locally to a gaussian ensemble. Cluster expansions have been applied to quantum quenches in continuous models in [59,60] where it is shown that any quench in a bosonic field theory (both in one and in higher dimensions) from any initial Hamiltonian to a massive free Hamiltonian leads to the GGE, as a consequence of the cluster decomposition property of the initial state. This is not so in the case of massless evolution and interacting initial Hamiltonian in one dimension [61] where more memory of the initial state survives at large times. Recently, the cluster property has been used to rigorously study the relaxation of a general lattice system after a quantum quench and the locality properties of the charges that should be included in the GGE [62]. Numerical studies based on cluster expansions have also been fruitful, even in the case of evolution under interacting Hamiltonians [63].
In most of the above studies the calculation is done at the level of local observables. Here instead we present an analysis of the evolution of the quantum state itself. We derive timedependent equations for the cluster amplitudes for the time evolution under a general Hamiltonian and solve them in the case of a gaussian Hamiltonian before and after the quench. The starting point of our general analysis is to understand what happens with the locality-related Properties 1, 2 and 3 of Section 3.1 when we evolve the system with an Hamiltonian H that is supposed to have local or short range interactions. As we will show, evolution under such Hamiltonians preserves the local character of local quantum states. Before we present the mathematical argument, we first give a physically intuitive picture. We can start analyzing Property 2, i.e. the cluster decomposition property. At the initial time, points at large distance to each other are completely uncorrelated. Once we start to evolve the system, excitations will start to travel from one point to another, eventually building up correlations between them. It is plausible that far distant points require a large time to be correlated, since we must wait until an excitation starting from one point reaches the other one. This should be true if the interaction is short ranged: long range Hamiltonians correlate distant points instantaneously. This statement can be made more precise if the evolution has a light-cone form, which is certainly true in relativistic field theories but also in short-range lattice models (Lieb-Robinson bound [64]) (Fig. 13). When this reasoning is correct, it means that two points at infinite distance cannot be correlated at any finite time, therefore the cluster property is preserved during the time evolution. Note that this has already been established in [62] (Theorem 6.3) where it is shown that the cluster property of a quantum state is preserved by the evolution under a local Hamiltonian. We will see below how we can naturally reach this conclusion also in our formalism. For the same reasons, also Property 1 i.e. the finiteness of correlation functions in the thermodynamic limit, should not be spoiled: a point deep in the bulk cannot be affected by what happens at the boundaries placed at infinite distance, at least not at finite time. These considerations suggest that a local quantum state, when evolved with a short range Hamiltonian, remains a local quantum state, therefore we expect that: In other words, the effect of the evolution is simply to change K n in time, without spoiling their analytic properties. Note that the zero-order cluster amplitude K 0 (t) must be included, since it is not guaranteed that the time evolution preserves our conventional normalization 0| = 1. The notation is chosen in such a way that we can absorb K 0 in the cluster expansion including n = 0 in the summation: we explicitly put a factor L d because we will see in a while that, with this choice, K 0 is L independent.
The above statement, which is based on physical arguments, can be made much more solid and, in order to be quantitative, we will focus on systems that evolve with an Hamiltonian of the same form (67) we used in Section 5.2. Remember that the cluster expansion can be written for any state such that 0| t = 0, therefore (93) can be written without supposing that | t is a local quantum state. What we should do is to write the equations that K n must satisfy during the time evolution and, assuming that the analytic properties of the initial condition are those of a local quantum state, to check whether they are respected also at any later time. To do this, we can write the time dependent Schrödinger equation: We already studied the time independent version of the Schrödinger equation in Section 5.2; including time derivatives to obtain the time dependent version is straightforward and with the same steps we arrive at: Now F is time dependent because the cluster amplitudes K n that enter in its definition evolve in time. Notice that the time independent equation (74) is readily obtained from the above if we look for eigenvectors: since eigenvectors evolve with phases e −iEt , therefore we should require from (95) ∂ t K n>0 = 0 and K 0 = −iE/L d and (74) is obtained.
Even though these equations are rather involved, we see that at least in the massive case m = 0, they respect the expected analytical properties of K n for a local state: if the initial conditions for K n are non-singular functions, then in the equations of motion do not appear any singularities, which implies that K n itself has no singularity at any finite time t . Additional care should be taken when we consider the massless case m = 0 due to the presence of singularities in (95): also in the static case we had this problem and in Appendix E we see how the singularity disappears, but we do not attempt a similar analysis in the time dependent case.
Notice that equations (95) are valid even when the Hamiltonian is taken time dependent, simply allowing for time dependent β n : this has important consequences regarding generic out of equilibrium problems, where the interaction can be varied in time. It does not matter how we change the coupling constants of the theory: the state preserves its form with time dependent K n .
Equations (95) are non-trivial, but they present two simple special cases that can be explicitly solved: • Time evolution governed by the free Hamiltonian (λ = 0). In absence of interaction the equations become trivial and we simply have: This is obvious, since if we take a generic cluster expansion and let it evolve freely, we simply have to multiply each creation operator by an oscillating phase factor a † k → e −iE( k)t a † k : in Schrödinger picture, these phases are attached to the K n rather than a † . • Quadratic Hamiltonian and gaussian initial state (K n =2 = 0). In this case the interaction preserves the gaussianity and the equations are exactly the same as (85) with the addition of the time derivative part: where for brevity we suppressed the explicit dependence on time and momenta. If β 2 is chosen time independent, we can rescale λ to have β 2 = 1. The solution can be readily written as: where: Notice that if we search for time independent solutions, we must impose C( k) = 0: this ensures K 2 is constant in time and it is exactly equal to (87), as it should.
We can briefly comment on the result above, keeping in mind the discussion of Appendix B. As we show there, |K 2 | ≤ 1 is a necessary requirement in order to have a meaningful gaussian initial state. Note that E( k) is the energy of the modes of the full Hamiltonian with λ = 0. Therefore, in order to have a meaningful theory we should impose E( k) to be a real quantity, which implies α 2 > 1. Then, with this hypothesis, it is a matter of simple algebraic manipulations to show that |K 2 ( k, − k, 0)| ≤ 1 implies |C| < 1. Knowing that |C| < 1 we see that K 2 ( k, − k, t) never acquires singularities during the time evolution, and moreover |K 2 ( k, − k, t)| ≤ 1, as it should.
It must be said that the evolution of K 2 can be obtained also through the method of Appendix B: (B.11) remains valid also if we replace the initial K 2 with the evolved counterpart and the initial correlators aa † S with the evolved ones aa † S (t). At this point from (B.11) we see that we can find K 2 through the correlators: Then the time dependent correlators can be evaluated using the suitable Bogoliubov transformation (B.6) that diagonalizes the full Hamiltonian and have a trivial evolution. We do not report these calculations since, even though straightforward, they are lengthy, but the result would match (98) as it should. We will not proceed further in the general discussion of these equations and for the time being we are going to see two possible applications of this method to quench problems.

Relaxation to GGE in quenches from interacting to free theories
As it is clear from the above, the cluster expansion is especially convenient for the study of quench dynamics in which the evolution is done according to a Hamiltonian diagonal in terms of the creation and annihilation operators a † k , a k . This can be achieved if the evolution Hamiltonian is gaussian (free) in terms of local fields and we choose to expand the initial state as a cluster expansion in terms of the corresponding creation and annihilation operators. In such an interacting to free quench, our approach gives a particularly simple proof of validity of the GGE conjecture.
The original GGE conjecture [23] explicitly claims that in thermodynamically large systems and for a large class of initial states | , including (but not necessarily limited to) ground states of local Hamiltonians, the evolution under an integrable Hamiltonian is such that the large time value of any local observable O equilibrates to a stationary point, whose value can be computed through a Generalised Gibbs Ensemble density matrix ρ GGE that maximizes the entropy subject to all local conserved charges Q n that characterize the integrable system. Regarding these charges, it is common a little abuse of notation: we talk about local charges as conserved quantities that can be written as an integral over all the space of some point wise operator, constructed out of the field and its derivatives. In a more quantitative way, the GGE conjecture says: where ρ GGE is defined as and the Lagrange multipliers λ n are defined by the condition for all the local charges Q n . By local observable O we mean any observable defined within a finite size subsystem, i.e. including also multi-point correlation functions. Thus the GGE conjecture equivalently claims that the reduced density matrix ρ A of any finite size subsystem A at t → ∞ is equal to the reduced density matrix corresponding to the GGE Notice that integrable field theories possess an infinite set of local conserved charges, apart from the Hamiltonian itself, therefore the GGE prediction preserves lot of information about the initial state, but it is still economic. As a matter of fact, any quantum system possesses the maximal set of conserved quantities constructed as the projectors on the eigenstates of the Hamiltonian. The size of this set increases exponentially in the system size, but the GGE conjecture states that only local conserved charges matter in the long time limit: this latter set increases only linearly in the system size. Locality of the charges is a desirable feature, since statistical ensembles are expected to satisfy locality requirements which are linked, as we explained, to the extensivity of thermodynamic quantities, but using only these conserved quantities in the construction of the GGE has been proven to be a too tight constraint: quasi-local charges must also be used [53]. Quasi-local charges are a relaxation of the locality constraint, that is to say such charges, instead of being integrals of a point wise operator, are constructed from short range operators integrated on the whole space. The interacting-to-free quench is the simplest situation in which the post-quench Hamiltonian is integrable and constitutes an important benchmark both to test the GGE prediction itself, but also to study the locality properties of the conserved charges that enter in its construction. It has been proven in [59,60] that in relativistic quenches in free theories, the expectation values of local observables indeed equilibrates and the right GGE density matrix can be written in terms of the occupation modes a † k a k : The GGE density matrix is gaussian and closely resembles a thermal ensemble, but each mode equilibrates at a different, momentum dependent, temperature. In this case, (104) is rephrased as: This result can be rephrased as follows: computing large time correlators of the relativistic field φ and its conjugate field π , only the two point connected correlators survive and, once we have expressed the latter in terms of the mode operators, only a † k a k gives a non-vanishing contribution to the observable. The key information that was used in [59,60] to obtain this result, is that the initial state respects the cluster property: this assumption was used to keep under control the analytic properties of the correlator in momentum space, similarly to what we will see in a while. However long calculations were necessary in order to use directly the φ and π operators, while we know that the cluster property, in a massive theory, can be expressed in terms of ψ and ψ † as well. Testing the GGE prediction on these operators is very simple, consider for example the connected two point correlator function for the field ψ: since the evolution is free, we can obtain the evolved fields in Heisenberg representation simply replacing a k → e −iE( k)t a k .
At this point in the limit t → ∞ we can apply the stationary phase method: if a k a − k c is smooth enough, then the integral becomes zero at t → ∞. In this context, smooth enough means that we require that a k a − k c and the other higher multipoint connected correlators have not extra δ singularities: if we would had a delta placed somewhere, then even at large times we would retain an oscillating time dependent term. As long as the cluster property for the ψ and ψ † fields works, multipoint correlators in the momentum space have only one δ singularity on the overall momentum, due to translational invariance in the coordinate space. In principle, from the cluster property (Property 2), other singularities in the momentum space are allowed, but they must be mild enough to have that, transforming back in the coordinate space, the correlator is well-defined in the infinite size limit (Property 1). Having under control the singularities, we can straightforwardly apply the stationary phase argument and (108) vanishes. With the same steps, it is easy to see that all the higher multipoint connected correlators vanish when t → ∞, instead one out of the two point correlators and the one point correlators survive: Notice that from our rules of Section 4 a 0 has embedded in it a L d/2 factor, therefore ψ is actually L independent, as it should. On these correlators we cannot apply the stationary phase argument: the delta structure in momentum space of (109) makes it time independent, instead in (110) we have no momentum to integrate over. If we have ψ = 0 because of symmetries of the initial state (for example if | is the ground state of the φ 4 model this is the case), then the long time behavior of local observables constructed with ψ and ψ † can be completely described in terms of a time independent gaussian ensemble ruled by the two point correlator in the momentum space, i.e. the GGE construction applied to the free case. Actually, there is a subtlety when we study relativistic field theories: in these theories, the ψ field is somewhat artificial and the real physical fields are the bosonic field φ and its conjugated momentum. As we stated at the beginning of Section 5, we can switch from the ψ field to φ with (40): as long as m = 0 it is easy to show that cluster property for ψ and ψ † holds if and only if it holds for φ and the conjugated momentum, this is thanks to the exponential damping of the kernel in (40). This implies that, as long as we require the cluster property on the physical fields, we have it on ψ and ψ † , therefore we have the same control on the singularities of the correlators in the momentum space. At this point, we can study the long time limit of connected correlators of the field φ exactly as we have just done for the ψ fields. Their expression, thanks to (37), can be written as the Fourier transform of multipoint correlators of a k , with some additional E( k) factors, but these in the massive case do not introduce any additional singularity and the stationary phase argument can be applied. This is no longer true in the massless case, but also the fact that cluster property of the field φ guarantees cluster property for ψ does not hold anymore, since the kernel in (40) decays too slowly to guarantee it. It can happen that, even if we have the cluster property on the physical field φ, this is spoiled on the ψ fields: therefore the argument we used before to study the δ singularities is spoiled and we do not expect the GGE to describe the large time behavior. This is consistent with [59], where relativistic field theories are studied and the fact of having massive excitations was crucial to have the GGE prediction realized.

Spectral expansion and work statistics
Among many other possible applications of the cluster expansion to quantum quench problems, a simple and straightforward application is the computation of the work statistics. The work statistics i.e. the probability distribution of the work done in the system by a quantum quench, is a quantity of central interest for several reasons. Firstly, the study of work statistics is important for the verification of quantum fluctuation relations which constitute the thermodynamic laws that govern non-equilibrium physics [65]. Second, it turns out to be related to other independently studied quantities like the Loschmidt echo and the Casimir force [35,36]. Through its relation to the latter quantity it can be shown that it exhibits universal characteristics [36,66]. Third, it can be used to probe the spectral properties of the pre-quench and post-quench system. Indeed the work statistics is nothing but the probabilities of the transitions from the initial state to any possible post-quench excitation. As such its general properties can be easily deduced. The minimum work level corresponds to the transition from the pre-quench ground state to the postquench ground state, which appears as an isolated δ-peak. Above it and starting from a threshold related to the energy gap of the lowest possible excitations there exists a series of other peaks or edge singularities that correspond to any isolated bound states or multi-particle excitations respectively that are present in the post-quench system [36]. From the type of the edge singularity we can extract information about the nature of the quasiparticles that describe the corresponding excitations [66,67]. By classifying the various excitations in increasing order of particle number, the cluster expansion provides directly a spectral decomposition of the transition channels and a characterization of the corresponding peaks and edge singularities of the work statistics.
In general, supposing our system is initialized in a ground state |G of some physical Hamiltonian, the work statistics P (W ) is simply defined as: Where |ψ n are the eigenstates of the post-quench Hamiltonian H with eigenvalues E n , E G is the energy of the pre-quench state. We introduce explicitly the norm N G = G|G to match the notation used so far for cluster expansions. Actually, instead of computing P (W ), it is often easier to compute its generating function g(t): Once we have it, we can compute the moments W n or the cumulants κ n : Notice that once we know all the cumulants we can reconstruct the moments, simply comparing the analytic expansion of the two expressions above. In particular, we will focus on the computation of cumulants.
The appealing feature of g(t) is that it acquires the physical meaning of a particular overlap, called Loschmidt echo [35,36]. Inserting (111) in (112) we simply get: The Loschmidt Echo for a state |G is defined to be the last quantity in (114), where H 0 is the pre-quench Hamiltonian. From the above, we can trivially notice that the momenta of the work distribution are associated with the expectation value of the powers of the post quench Hamiltonian on the initial state: Of course, the cumulants κ n are nothing else than the connected part of G|(H − E G ) n |G G . Momenta and cumulants can be computed with standard perturbative methods, since they are simply expectation values of some observables (powers of the post quench Hamiltonian). Nevertheless through the cluster expansion we can directly compute the generating function g(t), since it can be expressed as an overlap of local quantum states (114). The logarithm of g(t) is often more convenient for computations: The method of Section 4 is pretty suitable for evaluating the above expression. If |G has the structure of a local quantum state and we define |G t = e −iH t |G , from the discussion of Section 6 we expect also |G t to be a local quantum state, whose cluster amplitudes are changing in time. Therefore, we can also apply our graphical tool to evaluate the logarithm of the overlap log G|G t , as we explained at the end of Section 4, once we know the proper cluster expansion. There are different possible ways to compute the ground state energy, for example we can use (76). Notice that log g(t) is certainly an extensive quantity: this implies that all cumulants of the energy are extensive, as expected from the discussion of Section 3.1.
To simplify our presentation we perform a computation in the simpler case in which H is free: in this case the K n factors evolve with simple phases (96) and the problem reduces to determining the right initial conditions and dealing with our graphical expansion. In this section we will suppose the state initialized as the ground state of an Hamiltonian like (67), in particular we focus on the exactly solvable case of a quadratic Hamiltonian and then on the prototypical example of φ 4 theory. Actually, it is quite easy to have exact information about P (W ) when W is tuned very close to −E G , simply inserting the lowest eigenstates of the free Hamiltonian, constructed acting with a † on the bare vacuum, whose energy is chosen as reference energy (therefore zero). Increasing W we have no contribution to (111) until we reach the energy −E G : this is simply the gap between the pre-quench ground state and the post-quench ground state. In both cases, when 0 < W + E G < 2m we cannot excite any particle, but as soon as W + E G reaches the threshold 2m a pair of opposite momentum excitations can be produced and P (W ) is no longer zero. Then another threshold is reached at W + E G = 4m: in this case we can produce two such pairs or, in the φ 4 theory for example, also a cluster of 4 particles represented by K 4 . Then other thresholds are reached for any even multiple of the particle mass. As long as the odd amplitudes K vanish, these are the only allowed thresholds; otherwise other thresholds for odd multiples of m arise too.
Explicitly expressing (111) in terms of the eigenstates of the free theory we have: Notice that, since the energy of each particle is at least m, we can keep only the terms such that lm ≤ W + E G . The overlaps that appear are simply the wavefunctions with l particles, whose connected part is K l (43): therefore, a truncation in the sum of (117) up to l, implies that only K n with n ≤ l will enter in the computation of the work statistics. Therefore they do not enter as an infinite series as it happens for correlators, but as a finite sum, even though rather involved. In principle, the lowest terms of P (W ) can be easily calculated once we know K n . For example, if the odd cluster amplitudes are zero as it happens for the ground state of the φ 4 model, the contribution up to W + E G < 4m is: From the above expression it is easy to study the behavior of the work statistics at the first threshold W + E G = 2m: in this case, only the value of K 2 (supposed to be not zero) at zero momentum matters. Once we replaced the summation with the proper integration, we easily find: above, is the Heaviside Theta function, d is the d dimensional solid angle. It must be noticed that this and the other thresholds are immediately washed away by the thermodynamic limit, since N G is exponentially divergent in the system size. We can start with the simple case of a quadratic pre-quench Hamiltonian: The post-quench free Hamiltonian is simply obtained setting λ = 0. As reported in Section 5.2.1 and Appendix B, the initial state is a gaussian state with K 2 given by (87): it remains only to apply our graphical rules. This case has already been studied in literature with other methods [36], but nevertheless it is interesting to derive it again in our framework and use it as a check.
Since only K 2 is non-zero, the graphs that contribute to log G|G t are simply closed chains of 2n vertices, n vertices of the ket and n of the bra, as shown in Fig. 14. Using the rules of Section 4 we have that such a graph contributes as: Fig. 14. Graphs for log G|G t when a squeezed state is considered: they are made by closed chains of 2n vertices, the permutations that leave the graph unchanged are given by the dihedral group D n that counts 2n elements.
With S 2n the proper symmetry factor. Looking at Fig. 14 we can easily see that a graph with 2n vertices has 2n equivalent permutations, therefore we can write: Since we have |K 2 | ≤ 1 we can sum the series: Then, inserting this expression in (116) we get: The last step is to compute E G . There are many standard ways of doing this, but it is rather simple to obtain it from (76). Once we evaluate correctly F 0 with the rules of Section 5.2 we get: We have now the exact expression of log g(t), once we insert K 2 from (87) of Section 5.2.1. The expression for the Loschmidt echo agrees with the result of [36], as it should. Now that we exactly solved the gaussian case, we can focus on the φ 4 model: here the situation is much more complicated since, unlike before, the cluster amplitudes K n are non-vanishing for any even n and this makes impossible to explicitly sum all the graphs as we did before. To obtain an approximate result we can perform a naive perturbative expansion in λ, using the results (61) and (62) that tell us K 2,4 ∼ λ. With this information, the lowest order contributions in λ to log G|G t are given by the two graphs of Fig. 15 and they are O(λ 2 ): Consistently, also E G must be evaluated at the same order. Using (76) and keeping only O(λ 2 ) in F 0 we get: The expression for β 0,2,4 in the φ 4 case can be found in (78). At this point, putting all the pieces together we can find log g(t) up to order O(λ 2 ): the expression is rather long, but notice that through (113) we can easily evaluate the cumulants, simply deriving log g(t) around zero.
Once we insert (61) and (62) we get: The cumulants are all extensive quantities, as they should. Notice that the above integrals are UV divergent and a UV cut-off must be used. This means in particular that they are nonuniversal quantities. Instead universality is expected in the new edge singularity associated with the 4-particle cluster at the threshold W = 4m. In fact as we explained these edge singularities are low-energy effects and their exponents are essentially universal. This expression for the cumulants can be compared with those obtained from usual perturbation theory, since the momenta are nothing else than expectation values of powers of the free Hamiltonian on the φ 4 ground state (113) and the expression match, as it should. From the O(λ 2 ) expression of G|G t we can also study the O(λ 2 ) the thresholds W + E G = 2m and W + E G = 4m. The general expression for the threshold at 2m has already been reported in (119), therefore we need simply to plug in it the expression for K 2 at first order in λ (62). To study the threshold at 4m we write P (W ) as Fourier transform of the generating function g(t). Close to the threshold W + E G = 4m we have: Consider the first term in the above: its threshold is at 2m, therefore it is smooth in correspondence with the 4 particle threshold. The second term is not smooth at the four particle threshold and can be analyzed similarly to (118) giving: The O(λ 2 ) expression for K 4 is reported in (61) |K 4 (0, 0, 0, 0)| 2 = λ 2 4(2m) 3 (132) and Q is a numerical factor. We have therefore found that the introduction of φ 4 interactions in the initial state produced an extra edge singularity at W = 4m that was absent in the quadratic case. This is a signature of the presence of clusters of four particles and as we see it is characterized by a different exponent 3 2 d − 1 which was derived by means of an essentially phase-space calculation, while the presence of only pairs is characterized by the exponents d − 1 and d 2 − 1. This can be seen considering in full generality the threshold, without restricting to the O(λ 2 ) contributions. The fact that only K 4 contributes to the 4 particles threshold is true only at O(λ 2 ): the threshold at 4m in general depends on both K 2 and K 4 as it can be easily understood from (117), but these terms are of higher order in λ.

Conclusions and outlook
In this work we studied the effects of the general physical assumptions of locality on the analytic structure of a cluster expansion for local quantum states, in particular we find a wide class of states that satisfy the cluster properties for local observables (Section 3). We denoted such states as local quantum states and their cluster expansion is studied. Since ground states of physical Hamiltonians are known to satisfy the cluster property, we proposed to look for them among the general class of states we found: this ansatz was confirmed by perturbation theory at all orders (Section 5.1). The knowledge of the analytic structure of local quantum states allows us to study ground states directly from the Schrödinger equation (Section 5.2), changing the usual point of view of field theory. With the knowledge of the ground state itself, rather than the usual information about the expectation values computed on it, we have the possibility of studying physical quantities that are not directly linked to correlators of fields and we sketch the study of the spectral decomposition and the work statistics.
The knowledge of the exact analytic form of the ground state opens up the possibility of implementing efficient approximation schemes over it, as we commented in Section 5.3. Using the cluster expansion we can naturally go beyond the well known gaussian approximation, increasing our understanding of ground states. Moreover, using approximation schemes on the cluster expansion guarantees that any approximation we are implementing will produce a state whose qualitative features do not differ from those of a ground state: for example, cluster property will be always satisfied. This is quite important both on a conceptual level and on a practical one, since any approximation that gives a state qualitatively different from a ground state is doomed to be a very poor approximation of it. In particular, having exact information about the state gives us the possibility of implementing variational approximation schemes: these could permit to go beyond perturbation theory and such a program will be subject of future studies.
As we already commented, we expect that such an approach could find interesting applications also in out-of-equilibrium problems: in Section 6 we studied the behavior under time evolution of local quantum states, showing that evolution under local Hamiltonians does not spoil their physical properties and the evolved state remains a local quantum state. The cluster expansion permits us to explicitly write the time dependent Schrödinger equation for the state: this is an important change of perspective, since usually the time evolution is explored in the Heisenberg picture, where are the operators that evolve.
This work is far from being exhaustive and it has been mainly thought to present the basics of the cluster expansion for ground states, giving in the meanwhile as many tools as possible to handle such a construction. During the exposition we tried to explicitly comment over all the possibilities of such a construction and they will be subject of future investigations.
The operators follow standard commutation rules [ ( x), † ( y)] = δ( x − y) and the Hilbert space is a Fock space constructed acting with † on a vacuum |0 annihilated by . The number of particles is fixed: Note that for d ≥ 2 the ground state of (A.1) exhibits condensation of the zero mode, that is macroscopically occupied. The idea is to define |0 as the condensate state and define ψ † as operators that move particles from this reference state to other states, in such a way that the number of particles (A.2) is surely conserved. Let α k be the modes of the field , then following [43,68] we define the reference state |0 as: The particles in |0 are all in the condensate at zero momentum. Now, as mode operators a k define: where n 0 is the density of particles in the condensate and it must be determined self-consistently. The a † operator does not create a particle, rather excites it outside of the condensate |0 . It is possible to show [43] that in the thermodynamic limit a k can be regarded as a canonical bosonic operator in a weak sense, when we compute correlators on the ground state. As a matter of fact: when the above is inserted in a correlator the first term is O(L 0 ), since α † 0 α 0 = L d n 0 , while the extra correction is suppressed in L → ∞ and we obtain the usual commutation rules [a q , a † k ] = δ k, q . Further details of the limit can be found in [43]. The ψ operators are defined as the Fourier transform of a k .
From the expectation value of (A.2) we get the self-consistency equations for the condensate density n 0 .
The final step is to take the thermodynamic limit of (A.1) in terms of ψ and obtain an effective Hamiltonian directly defined in the thermodynamic limit L → ∞ [68]: Note that the above is a thermodynamic limit of the Hamiltonian valid only in a weak sense, nevertheless it is enough to study the thermodynamic properties of the ground state. The ψ † operators create localized excitations, therefore it is plausible that they satisfy the cluster property. Moreover, the number of excitations is not fixed by the above Hamiltonian. The methods of Section 5 can be easily adapted to the Hamiltonian (A.8) and used to compute the cluster expansion of the ground state of (A.1) in terms of ψ.

Appendix B. Gaussian states and quadratic Hamiltonians
A well known type of state in quantum field theories is a translational invariant squeezed state that is actually a particular case of (9) obtained keeping only K 1 and K 2 , but we will focus on the case in which K 1 is absent. Therefore: Because of its simplicity many quantities (correlators, free energy..) can be computed exactly and we often use this state as a check and a guide for our more general approach, therefore we want to enlist some of its properties. A squeezed state naturally arises considering ground states of Hamiltonians with a perturbation quadratic in the fields. As an example, consider the simple Lagrangian: If we consider λφ 2 as a perturbation, then we will proceed defining creation and annihilation operators based on the mass m: Actually, we could have proceeded also in a smarter way, interpreting (B.2) as a mass shift: Now we have a simple free theory and we can diagonalize its Hamiltonian through new bosonic operators b and b † such that: Therefore the ground state |S (with this notation we anticipate it will be a squeezed state) is simply the vacuum of the b operators, thus it is completely identified by the requirement b k |S = 0. This last relation permits us to express |S in terms of the a operators and their vacuum state |0 : as a matter of fact comparing the expressions (B.3) and (B.5) and the analogous expressions of the conjugated momenta of the field φ we can express the b operators in terms of a. Where: With these relations we can readily obtain an equation for the ground state in terms of the a operators: The solution of this last equation is a squeezed state (B.1) with: We can use (B.6) also to compute correlators of the a fields on the ground state, simply inverting (B.6). For example: With similar computation we can arrive at the following result for the correlators on a squeezed state as (B.1). In particular the expression below is valid for each K 2 , not necessary in the form (B.9) and also complex K 2 are allowed: Since the squeezed state (B.1) is gaussian with zero mean, all the correlators can be computed once we know the two point correlators above. Notice that, since our squeezed state is translational invariant, all the correlators with a non-zero total momentum must be zero. Notice that something odd happens in (B.11) if we let |K( k, − k)| 2 > 1: in that case the correlator a k a † k S becomes negative. This situation is impossible to reach when K 2 describes the ground state of (B.2) since (B.9) implies |K 2 | ≤ 1, nevertheless we can learn something important studying such a state. Having a negative correlator is non-sense, since: The inner product is positively defined and we must have S |S ≥ 0, therefore a k a † k S ≥ 0. We can solve this absurd fact if we notice: The two expressions are identical whenever |K 2 | < 1, but they are rather different as soon as this bound is violated. We can see that the right expression for a k a † k S is indeed the series rather than the fraction. In this way whatever is the value of K 2 we always have a k a † k S ≥ 0 as it should and when |K 2 | > 1 the correlator does not become negative, but it rather becomes a positive divergent quantity. From this expression it is clear that the correct interpretation of the bound |K 2 ( k, − k)| < 1 is as the radius of convergence of the geometric series: this is a very useful observation we are going to use in Appendix D. Notice that if the inequality is saturated, we have a divergence in the correlation function and this is a signal of criticality, as we will see in Appendix D. For the moment, notice that if we are studying the ground state of a Lagrangian such as (B.2), then for λ > −m 2 the condition |K 2 | < 1 is always met, instead when λ = −m 2 we have K 2 (0, 0) = 1 and the Lagrangian (B.4) becomes massless M = 0, therefore it is critical. In the remaining case λ < −m 2 we have that the theory described by (B.2) is ill defined because M 2 < 0, therefore the Hamiltonian is no more bounded from below and has no more a ground state.

Appendix C. Partial resummation in the two point correlator
In Section 4 we have seen the basic graphical tool to compute the correlators and the free energy, but apart from formal calculations, we can really evaluate only a finite number of graphs and truncate the expansion. In this section we are going to perform a partial resummation of graphs that can give us some insight on important, non-perturbative phenomena, for example we are going to use these results to correctly insert critical systems in our framework (Appendix D). In Feynman graphs a central role is played by the one particle irreducible graphs, for example the renormalization of the masses is to due their resummation in the perturbative series [45]. Also in our case we can define irreducible graphs and in this section we explain how they can be summed, exactly as it happens in the usual Feynman graphs [45]. There will be some complications due to the presence of two kinds of vertices, but the ideas remain the same. The best way to introduce the irreducible graphs is looking at the two point connected correlator in momentum space: It is immediate to notice that the graphs that contribute to the connected correlators are all and only the connected graphs, among them we can define the irreducible graphs. A graph is called reducible if after we cut one internal line, the external legs become disconnected, irreducible otherwise. Some examples are given in Fig. 16. The irreducible graphs can be used as building blocks to construct general graphs of the two point connected correlators, simply joining them. This is quite obvious, since if a graph is reducible then it can be divided in two blocks, if one of them is reducible it can be divided again: we can proceed until we have divided the graph in its irreducible parts. It must be stressed that the presence of two different kinds of vertices means some difficulties, as a matter of fact in the correlator a † a c do not enter contributions only from its irreducible graphs, but also from the irreducible graphs of aa c and so on. With a little thought, it is possible to organize this combinatorial problem in a simple manner, if we define a 2 × 2 matrix G ir ( k), whose entries are defined as: 17. Graphical representation of the square of G ir ( k), graphically the products of the matrix entries are simply the graphs obtained joining the arrows. In the square of G ir appear all the graphs of the two point correlator made of two irreducible pieces.
Taking the square of such a matrix is equivalent to consider the contributions of all the graphs that are made of two irreducible pieces. The key to understand this fact is to notice that for graphs with two external legs the value of the graphs obtained by joining two of them is simply the product of the two contributions. With this in mind consider Fig. 17, in each entry we can put the value of an irreducible graph: when we perform the matrix product, the product of the values of graphs is simply the value of the graph obtained joining them. Therefore the square of such a matrix contains the contributions of all the graphs constructed joining two of the irreducible graphs. If instead of considering only an irreducible graph for each entry we consider the sum over all the possible graphs, then the square produces the sum of the values of all the possible graphs made by two irreducible pieces. In a similar way the nth power produces all the possible graphs made of n irreducible graphs. Notice that each entry of the nth power is still associated to the same correlator of the same entry of G ir , for example [G ir ( k)] n 1,1 is a contribution to a † k a k c . It is not difficult to convince ourselves of the last statements after some specific examples are considered. At this point we can do the last step to find the two point correlators: we have simply to sum over all the possible graphs. This means that we have to sum over all the possible powers of G ir , therefore we have a geometric series.
The last equality holds only if we are within the radius of convergence of the series and if G ir is too big, the correlators diverge: in Appendix B we discuss the specific case of a squeezed state, but the same discussion can be easily generalized. We can write this expression in an even more compact form if, instead of the normal order, we reverse it and use a k a † k c = 1 + a † k a k c : 18. Irreducible graphs for squeezed state. Fig. 19. Dressed propagators can be represented with two arrows and their value changes with their direction: with this notation we can recycle the old graphical rules and is still true that each vertex has all ingoing or all outgoing arrows.
In the very simple case of a squeezed state (K n =2 = 0) we can exactly compute G ir , since we have only the two irreducible graphs of Fig. 18, therefore: If we plug this expression in (C.7) we match exactly the result (B.11) given in Appendix B, as it should. In the general case there are many contributions to the irreducible graphs, but all the other contributions have one or more K n≥2 vertices, therefore the series of irreducible graphs must be approximated with a proper truncation.
Actually we can proceed even further: consider any graph for a generic correlator, then we can construct other graphs from it replacing the 'bare' arrows with some irreducible graphs. This leads us to the concept of dressed propagator [45] and to a definition of new graphical rules to compute propagators in which we have explicitly taken care of all the possible insertions of irreducible graphs. In the usual Feynman graphs [45] there is only one kind of vertex, thus only one type of dressed propagator, here instead we have more. To recycle as far as we can the old rules, we represent the different propagators as in Fig. 19: differently from the arrows of Section 4 these dressed propagators have also a weight, determined by the connected two point correlators. Notice that the single arrow has the same momentum at the two extrema, instead the double arrows exchange their sign: as before the arrows departing from a vertex are all ingoing or outgoing. The presence of double arrows allows a vertex also to self-interact; moreover a graph can contain single arrows: for example, it is clear from the definition that the two point correlators are graphs made of a single arrow. Through these rules we do not compute correlators in the normal order, but rather in the opposite order: all the creation operator a † are on the right, the destruction operators a on the left. This comes out from the fact that the first arrow of Fig. 19 brings a contribution a k a † k c . Since we have already summed on the one particle irreducible graphs, in the graphs with the dressed propagators we must avoid all the graphs in which we can find subgraphs with two external legs: in Fig. 20 there are some examples of graphs to be avoided. All the other rules (conservation of momenta, integration on internal lines, symmetry factors) are the same as before, in Fig. 21 there are some allowed graphs with their values.

Appendix D. The cluster expansion at criticality
There is an important point that should be clarified and regards the presence of criticality. As we stated in Section 4.1, asking certain decaying properties for the cluster amplitudes (21) leads to the conclusion that any graph in the momentum space does not have any singularity. On Fig. 20. When we draw graphs with the dressed propagator we must take care of having already summed on the irreducible graphs, this means that we cannot draw subgraphs with only two external legs. Therefore, neither of these graphs can be drawn: in the second, in the dotted lines, there is the subgraph with two external legs. the other hand, we know in advance that the two point correlators in critical systems typically exhibit divergences at zero momentum. Of course, even in critical systems Properties 1, 2 and 3 are expected to be valid. There are two possible ways to acquire the singularity in the two point correlator: the amplitudes W c n do not decay fast enough (21) or the expansions we wrote in Section 4.1 do not converge anymore. Note that it is very unlike that we can drop the decay condition about the cluster amplitudes: as a matter of fact (21) guarantees that all the graphs we compute do not have infrared divergences. If infrared divergences arise in some graphs in the expansion of correlators, in order to have the latter independent on the volume we would need cancellations of the infrared divergences among the graphs. On the other hand, it is simpler to imagine that the expansions do not converge anymore: this is exactly what happens in the gaussian case of Appendix B. In particular consider the expression (B.13) a k a † k S = (1 − |K 2 ( k, − k)| 2 ) −1 . In that case it is clear that the singularities of the two point correlator are not associated with singularities of |K 2 |, but rather with the points where it equals 1: if we look at the geometric series, those are the points where it becomes divergent. Then, rephrasing the last sentence, we can say that a singularity in the two point correlator signals a problem in the convergence of a series. Actually, this interpretation is valid also in the most general case: consider the formal expression for the two point correlators (C. 7) found in Appendix C, in terms of the one particle irreducible graphs.
The expression (C.7) is the natural generalization of (B.13) and gives to the singularities in the two point correlators the same interpretation of the gaussian case. To have a singular correlator (C.7) we should ask det(1 − G ir ( k)) = 0: this condition is associated to the lack of convergence of the geometric series (C.6) and det(1 − G ir ( k)) = 0 can be fulfilled even if we ask for (21), therefore any graph that appears in G ir ( k) is not singular. We can conclude that at criticality the state still has a cluster expansion whose amplitudes are expected to satisfy (21), but the graphical expansions of the correlators are not expected to be convergent any longer: in particular, we are at the boundary of the radius of convergence of (C.6).

Appendix E. Zero bare mass, failure of perturbation theory and non-perturbative study
As we discussed in Section 3.1 and Section 4.1, the K n cluster amplitudes of quantum local states should have no singularities in order to guarantee our locality properties and the perturbative calculation of Section 5.1 confirms this behavior, at least in the massive case (m = 0). The troubles arise when we let m → 0: in this case the expressions (61) and (62) become singular because of poles at zero momenta. This is something not expected and we have a contradiction between the perturbative result and the supposed analytic structure of K n , on the other hand near the poles the perturbation theory is not reliable anymore. Perturbation theory makes sense if each contribution is small, but near a pole, regardless how small λ is, the contributions (61), (62) are always diverging quantities. This argument makes us suspect that the poles in (61), (62) are not really poles of K 2 and K 4 , but rather signal a problem with the perturbation theory. To understand better this feature we should go beyond the perturbative regime and this is not a simple task. We will study directly the functional equations (74): we have seen in Section 5.2 that (61) and (62) are simply the O(λ) solutions of (74). Before doing so, we should see what happens in the exactly solvable case of the gaussian state, therefore consider (B.9) of Appendix B. In particular we can extract the first order in λ: We can see that when m = 0, the first term of the expansion in λ becomes singular at zero momentum, on the other hand the exact solution is finite in k = 0, in particular K 2 (0, 0) = −1 regardless the value of λ. Therefore, at least in the gaussian case, we see that the picture discussed above is correct and K 2 has no singularities. In an interacting theory like the φ 4 model we cannot solve exactly the functional equations (74), nevertheless we can attempt a consistency check: we will suppose that the cluster amplitudes K n are not singular, then check this is consistent with (74). In particular, we will check consistently that K 2 does not have singularities at zero momentum, where the perturbation theory fails. Such a consistency check is not exhaustive, since we will see only the absence of singularities in K 2 : it would be important to check also the other K n and, even better, to study directly the cluster amplitudes in a non-perturbative way. Nevertheless, because of the involved form of (74), even a consistency check for K 2 is remarkable. In order to proceed further, we need to use the equation (74) with n = 2 and study the features of F 2 ( k, − k): this quantity can be computed summing the contributions of the proper graphs, as explained in Section 5.2.1. As we are going to see, these graphs can be organized in terms of their singular behavior in the momentum k. The rule is this: each outgoing arrow that departs directly from the 'circles' vertices carries a singularity ∼ | k| 1/2 , the same singularity is carried by the outgoing Fig. 22. Graphs for F 2 constructed using only K 2 and K 4 . arrows associated with K 2 vertices. Other outgoing arrows have no singularities associated with them. This counting of singularities holds in dimensions higher than one (d > 1) and we will suppose that this is our case. To understand this fact it is useful to see some explicit examples, then the general case will be obvious: in Fig. 22 are drawn all the graphs that contribute to F 2 constructed using only K 2 and K 4 . Consider for example the first graph in the upper left corner of Fig. 22, its contribution can be readily written as: As promised, it has a divergence as | k| −1 , that is to say ∼ | k| −1/2 for each of its external legs. Instead, the graph immediately below of it contribute as: Again, we have a singularity ∼ | k| −1 , that is to say ∼ | k| −1/2 for each of its external legs, as we said. Now, look instead to the contribution of the first graph from above, in the middle column of Fig. 22: our counting argument asserts the whole graphs should have a singularity ∼ | k| −1/2 . Its contribution is: This value diverges as | k| −1/2 . As a matter of fact, if d > 1 and the term E( k) is omitted, the remain integration does not diverge when k → 0. It is not difficult to repeat the check also for the remaining graphs of Fig. 22 and other graphs constructed with general cluster amplitudes K n . The reason why the counting of divergences is so simple is that integrations over closed loops remove singularities. Now we should write down the equations for K 2 : their exact form does not matter, what is important is taking care explicitly of the pole structure and the presence of K 2 attached to external legs. We will see in a while that the equations (74) with n = 2 can be organized as: Where ξ 0 , ξ 1 and ξ 2 are non-singular functions of k. This statement can be easily seen if we organize the graphs in a proper way. For example, consider the contribution of the first and second graphs from the above, in the second column of Fig. 22: the value of the second graph is simply K 2 ( k, − k) times the value of the first. Therefore, once we have summed these contributions we get 1 + K 2 ( k, − k) times the value of the first graph: this value has exactly the | k| −1/2 factor we need. In a similar way and being careful with the symmetry factors, we can organize also the other graphs in the same way and arrive to (E.5): we invite the reader to check some additional explicit examples. Without a full knowledge of K 2 and K 4 it cannot be excluded that both ξ 1 and ξ 2 become zero at k = 0, but if at least one among ξ 1 and ξ 2 is not zero at zero momentum we conclude: Therefore, K 2 (0, 0) = −1. Notice that K 2 is not singular at zero momentum, moreover we have obtained also a simple non-perturbative result. If the bare mass is zero (m = 0) and the dimension is larger than one (d > 1), then we must have K 2 (0, 0) = −1. Notice that the value of K 2 (0, 0) is not only non-perturbative, but even independent from λ, exactly as happened in the gaussian case (E.1).

Appendix F. Consistency check for correlators in φ 4 theory
In Section 4 we developed a method to obtain correlation functions once the cluster expansion is known, while in Section 5 we saw how it is possible to compute the cluster expansion for ground states, at least perturbatively, and we gave the first order computation for the φ 4 theory. A simple self-consistency test is to check if the correlators of the fields computed from this method match with the standard perturbative computation, i.e. the usual Feynman graphs [44]. This is what we are going to do here, therefore in this section we deal with a φ 4 Lagrangian (58). The field correlators in the interacting ground state can be computed extracting the latter from the free ground state, i.e. the vacuum, using (19). We choose to evaluate the normal ordered correlations to make the combinatorics easier, but the calculations can be done also without this requirement. Using the same notations of Section 5.1 and evaluating the correlators in momentum space, we have: This expression resembles very much (46), but notice that in (46) the exponential of the Hamiltonian was only on the left and not also on the right of the product of fields. Therefore, proceeding exactly as we did in Section 5.1 we arrive at the expression: where we used that, by definition, φ I ( k, 0) = φ( k). Notice that the only difference with (54) is that the domain of integration of the exponential is (−β, β) instead of (0, β). From the expression above, it is not difficult to develop the formalism of Feynman graphs [44] (that is different from Section 5.1 because of the different domain of integration of the euclidean time), but since we are interested only in the first order in λ, it is easier to expand directly the exponentials. The correlators can be easily computed using (55). For example, for the four-point correlator we get: The integrations can be easily performed giving: In a similar way we can also compute the two-point correlator. All other correlators vanish at this perturbation order.
We will now compute the same object, but using the cluster expansion (9) of |G with the expressions (61) and (62) for the functions K n . The graphical rules of Section 4 allow us to find the expectation values of the fields a k , a † k that are related to φ through the expression (50), therefore the four point correlator on the ground state is: At first order, only the graphs in Fig. 23 contribute, therefore among the four point correlators of the fields a and a † only aaaa G and a † a † a † a † G are of order λ, while all others are of higher order. From Section 4 we can read the value of the graphs in Fig. 23: Thus, at first order in λ we get: Plugging the result for K 4 of equation (61) into the expression above we get exactly the same result as in (F.4), as it should be. With similar computation also the two point correlator can be computed and the two first order computations match, all the other correlators vanish at this order, as they should.