On the path integral representation for quantum spin models and its application to the quantum cavity method and to Monte Carlo simulations

# On the path integral representation for quantum spin models and its application to the quantum cavity method and to Monte Carlo simulations

Florent Krzakala Laboratoire PCT, Unité Mixte de Recherche (UMR 7083 Gulliver) du CNRS et de l’ESPCI ParisTech, 10 rue Vauquelin, 75231 Paris, France    Alberto Rosso LPTMS, Unité Mixte de Recherche (UMR8626) du CNRS et de l’Université Paris-Sud, Bât. 100, Université Paris-Sud, 91405 Orsay Cedex, France    Guilhem Semerjian    Francesco Zamponi LPTENS, Unité Mixte de Recherche (UMR 8549) du CNRS et de l’ENS, associée à l’UPMC Univ Paris 06, 24 Rue Lhomond, 75231 Paris Cedex 05, France.
June 30, 2019
###### Abstract

The cavity method is a well established technique for solving classical spin models on sparse random graphs (mean-field models with finite connectivity). Laumann, Scardicchio and Sondhi [arXiv:0706.4391] proposed recently an extension of this method to quantum spin-1/2 models in a transverse field, using a discretized Suzuki-Trotter imaginary time formalism. Here we show how to take analytically the continuous imaginary time limit. Our main technical contribution is an explicit procedure to generate the spin trajectories in a path integral representation of the imaginary time dynamics. As a side result we also show how this procedure can be used in simple heat-bath like Monte Carlo simulations of generic quantum spin models. The replica symmetric continuous time quantum cavity method is formulated for a wide class of models, and applied as a simple example on the Bethe lattice ferromagnet in a transverse field. The results of the methods are confronted with various approximation schemes in this particular case. On this system we performed quantum Monte Carlo simulations that confirm the exactness of the cavity method in the thermodynamic limit.

## I Introduction

Mean-field approximations are often useful first steps to unveil the physical content of realistic models. This is all the more true when exact solutions are probably impossible to obtain in the finite-dimensional setting, in particular when quenched disorder and/or quantum effects have to be taken into account, as for instance in the case of Anderson localization popu_first (). Another example is dynamical mean-field theory DMFT (), that has been a very fertile approach to the problem of strongly correlated fermions. It can be sometimes preferable to study mean-field theories not by making an approximation to a finite-dimensional model, but rather by formulating a model which is mean-field by nature, this allowing in particular to state the results in a mathematically clearer way. The simplest of such examples is the Curie-Weiss model of ferromagnetism, in which classical Ising spins all interact attractively which each other, with a coupling constant scaling inversely with the size of the system to ensure a well-defined thermodynamic limit. The equivalent for quenched disordered systems is the Sherrington-Kirkpatrick model SK () of a spin glass, where again all spins interact weakly with other, yet with coupling constants of random signs.

The mean-field character of the above mentioned models arises from their infinite connectivity (in the thermodynamic limit). There exists however another class of models, which are still mean-field yet keep a finite connectivity, each of the degrees of freedom they possess interacting only with a finite (with respect to ) number of neighbors. For ferromagnetic models they can be obtained by the Cayley tree construction, where one draws an infinite regular tree and studies the magnetization of the root site Baxter (). Cayley tree models have however pathological surface effects, and the theory of finitely-connected mean-field frustrated systems is better defined on random graphs VB (); MePa (), of fixed or fluctuating connectivity. Classical models of spins on such random structures have been the subject of extensive study in the last decade. These works were motivated on the one hand by their somehow more physically realistic features, namely the finite connectivity, and also because of their strong relationship with issues originating from computer science, namely the understanding of phase transitions in random constraint satisfaction problems MoZe (); MePaZe (); pnas (). Finite-connectivity models are technically much more involved than their fully-connected counterparts. The replica method Beyond () that has been first developed to solve the Sherrington-Kirkpatrick model becomes less practical in this setting replica_diluted (), and the alternative cavity method turned out to be more useful MePa ().

The interplay between quenched disorder and quantum fluctuations can lead to a very rich phenomenology, and in particular the properties of the glass phase found at low temperatures in classical models can be qualitatively modified when a transverse field acts on the system Giulio (). More generally the issue of the nature of the quantum phase transitions at zero temperature Sachdev () in presence of disorder is a very rich one. In the context of mean-field theory this point has been mainly studied in fully-connected models BrayMooreQSK (); Goldschmidt (); Nish_stat (); Georges (); Leticia (); Chamon01 (); MI08 (), with a few exceptions that appeared in the last year BT07 (); qc_first (); qksat (); qbl (). In a very interesting contribution qc_first () Laumann, Scardicchio and Sondhi made a first step in extending the cavity method to quantum spin models in a transverse field, and in this paper we shall develop further this idea by solving a discretization problem which plagued their proposal. Let us also mention here the work of Knysh and Smelyanskiy qksat (), who developed a similar approach in the framework of the so-called static approximation BrayMooreQSK (); Nish_stat (); Leticia (). The motivations for this line of work is two-fold. On the one hand one can expect an even richer physical behavior of finitely-connected quantum models with respect to the fully-connected ones. The possible fluctuations in the local geometry, and some notion of distance which was absent in fully-connected models opens the way to a more complex phenomenology. On the other hand one can aim at a better understanding of some issues of quantum computing, and in particular on the use of quantum annealing (or adiabatic algorithm) diego (); annealing (); adiabatic () to solve random constraint satisfaction problems. These quantum algorithms do indeed rely on the application of a transverse field on spin models that have been extensively studied at the classical level with the cavity method. Computing the location and nature of the phase transitions qrem () encountered along the annealing path (as the transverse field is progressively turned off) might give some informations on the behaviour of these quantum algorithms themselves.

The remaining of the article is organized as follows. In Sec. II we recall the Suzuki-Trotter approach to spin-1/2 models in a transverse field, and develop our main technical contribution in Sec. II.2, where we show how to actually build the spin trajectories of the path integral representation of the imaginary time evolution operator. Section III is then devoted to the study of a very simple example of finitely-connected quantum model, namely the Bethe lattice ferromagnet. We first explain the continuous time quantum cavity treatment of this model, before presenting the results of the method and confronting them with some approximate approaches. In Sec. IV we present numerical results of Monte Carlo simulations we performed for this model, and show how the computations of Sec. II.2 can be turned in a simple and versatile quantum Monte Carlo method. The generic formalism of the quantum cavity method is developed in Section V; we hope this order of presentation, and the inclusion of a fully worked-out example before the general case, will ease the reading of this work. We finally draw our conclusions and put forward perspectives for future work in Sec. VI. Some technical details are deferred to a series of Appendices.

## Ii Path integral representation for spin-1/2 models

### ii.1 Spin models in a transverse field: Suzuki-Trotter formalism

Let us consider the Hilbert space spanned by the orthonormal basis of kets , where denotes a configuration of Ising spins, . This space can be viewed as the tensorial product of spins , with operators and , whose action on the base vectors is defined by

 σzi|σ––⟩ = σi|σ––⟩ , (1) σxi|σ––⟩ = |σ1,…,σi−1,−σi,σi+1,…,σN⟩ .

From a classical energy function of Ising spins, , one can construct an operator , diagonal in the basis. The Hamiltonian operators investigated in this paper are obtained from such a classical energy by the addition of a transverse field,

 ˆH=ˆE−BN∑i=1σxi ,with B≥0 . (2)

Our goal is then to compute the quantum statistical mechanics properties at inverse temperature , i.e. the partition function and the average of observables (operators) , defined by

 Z=Tr(e−βˆH) ,⟨ˆO⟩=Tr(ˆO e−βˆH)Tr(e−βˆH) . (3)

A well-known way of tackling such problems is to transform them into an extended Ising model by using the Suzuki-Trotter formula suzuki (), as summarized in the following lines :

 Z = Tr((e−βNsˆE+βNsB∑Ni=1σxi)Ns) = limNs→∞Tr((e−βNsˆEeβNsB∑Ni=1σxi)Ns) = limNs→∞∑σ––1,…,σ––NsNs∏α=1 ⟨–σα|e−βNsˆE eβNsB∑Ni=1σxi|σ––α+1⟩ = limNs→∞∑σ––1,…,σ––NsNs∏α=1e−βNsE(σ––α)∏i,α⟨σαi|eβNsBσx|σα+1i⟩ .

In the two last lines . For a finite value of the number of Suzuki-Trotter “slices” , the problem has thus become one of Ising spins, each of the being promoted to a ring with nearest neighbor ferromagnetic interactions along the “discrete imaginary time” axis (with periodic boundary conditions). The original interactions acts indentically and independently on each of the configurations . For notational convenience we shall use bold symbols for quantities that depend on the slice , for instance is the configuration of the ring of Ising spins at site and is the full configuration of the spins. We can thus introduce a probability measure on the Ising spins,

 μ(σ––) = 1ZNse−β˜E(σ––)N∏i=1w(σi) , (5) ZNs = ∑σ––e−β˜E(σ––)N∏i=1w(σi) ,

such that the normalization constant reduces to the partition function in the limit (in the following we shall sometimes keep implicit the dependence on ). To write in a compact way this last equation we have defined

 ˜E(σ––)=1NsNs∑α=1E(σ––α) , (6)

the average of independent copies of the classical energy on the various slices, and

 w(σ)=Ns∏α=1⟨σα|eβNsBσx|σα+1⟩=Ns∏α=1(cosh(βBNs)δσα,σα+1+sinh(βBNs)δσα,−σα+1) , (7)

the ferromagnetic interaction along the imaginary time axis induced by the transverse field (we use ). One can easily show that the average value of observables can be obtained in this formalism as

 ⟨ˆO⟩=∑σ––μ(–σ)⟨σ––α|ˆOe−βNsˆH|σ––α+1⟩⟨σ––α|e−βNsˆH|σ––α+1⟩ , (8)

where the slice number is here arbitrary, thanks to the cyclic invariance around the discrete imaginary time axis. This can be simplified further for observables diagonal in the basis, i.e. that can be written as :

 ⟨ˆO⟩=∑σ––μ(–σ)O(σ––α)=∑σ––μ(σ––)1NsNs∑α=1O(σ––α) . (9)

A non-diagonal observable we shall study in the following is the transverse magnetization (written here for an arbitrary site ), which can be computed as

 ⟨σxi⟩ = ∑σ––μ(σ––)1NsNs∑α=1⟨σαi|σxeβNsBσx|σα+1i⟩⟨σαi|eβNsBσx|σα+1i⟩=∑σ––μ(σ––)1NsNs∑α=1(tanh(βBNs))σαiσα+1i .

### ii.2 The continuous imaginary time limit

To recover the truly quantum properties of the model one has to perform the limit . The basic degrees of freedom which were the configurations of a ring of Ising spins then becomes piecewise constant functions of an imaginary time parameter , the discrete coordinate being mapped to with the correspondence . In this limit the sum over in the expression (5) of the partition function is naturally interpreted as a path-integral. The discreteness of the spin degrees of freedom actually make such a path-integral representation Farhi_Gutmann () easier to formulate than Feynman path-integrals for continuous coordinates Feynman53 (), and can be given a rigorous mathematical content Ginibre (); Gallavotti (); Aiz (); Ioffe (). Note that these continuous time trajectories can be easily represented in the memory of a computer, as the trajectory of site is fully specified by and the times at which the spin flips. Actually numerous continuous time quantum Monte Carlo algorithms do exist, see for instance Beard_Wiese (); worm (); qcluster (); evertz ().

The rest of the paper will crucially rely on the procedure developed in Sec. II.2.2 and II.2.3. Though it will also be useful for analytical purposes, it is more intuitively motivated by the following simulational consideration. Maybe the simplest way to ensure the detailed balance condition in a Monte Carlo simulation which aims at sampling an arbitrary measure is to perform transitions from the current configuration to a configuration obtained by replacing the value of a randomly chosen degree of freedom , by a random value drawn from the measure conditioned on all other degrees of freedom. This procedure is known in classical simulations as the heat-bath, or Glauber algorithm. Its equivalent in quantum simulations consists in drawing a new configuration of the ring , or of the trajectory in the continuous imaginary time, according to the equilibrium measure induced by the spin trajectories of all other sites. A moment of thought reveals that this boils down to study the evolution of a single spin-1/2 in the presence of a constant transverse field and a piecewise constant longitudinal field, the latter being the effective field induced by the rest of the system on . This is precisely the issue we shall tackle in Sec. II.2.2 and II.2.3, after having recalled in Sec. II.2.1 the well-established path-integral representation of a spin-1/2.

#### ii.2.1 The path integral representation of a single spin in constant fields

Let us define the propagator for the evolution during an interval of imaginary time of a spin in constant transverse and longitudinal fields ( and respectively):

 W(σ→σ′,h,λ)=⟨σ|eλ(hσz+Bσx)|σ′⟩ . (10)

The diagonalization of the order matrix easily leads to

 W(σ→σ′,h,λ)=⎧⎪⎨⎪⎩cosh(λ√B2+h2)+σh√B2+h2sinh(λ√B2+h2)%if σ=σ′B√B2+h2sinh(λ√B2+h2)if σ=−σ′ . (11)

The path-integral representation of this propagator reads Farhi_Gutmann (); qcluster ()

 W(σ→σ,h,λ) = ∞∑n=0B2n∫λ0dt1∫λt1dt2…∫λt2n−1dt2nexp[σh(2t1−2t2+⋯−2t2n+λ)] , (12) W(σ→−σ,h,λ) = ∞∑n=0B2n+1∫λ0dt1∫λt1dt2…∫λt2ndt2n+1exp[σh(2t1−2t2+⋯+2t2n+1−λ)] . (13)

Each term of these expressions corresponds to a spin trajectory that changes value at times ; it is weighted by a factor raised to the number of such discontinuities, and by ; a spin trajectory with identical (resp. opposite) initial and final value has to jump an even (resp. odd) number of times. There are two ways to convince oneself of the correctness of this result. Applying the Suzuki-Trotter formalism to this single spin problem leads to such a weight in the limit111The weight becomes raised to the number of discontinuities in the spin trajectory, and the factors are absorbed in the change of variables from discrete to continuous time.. Alternatively one can notice that Eqs. (12),(13) coincide with (11) at and that they obey the same set of first order linear differential equations,

 ∂λ∂λW(σ→σ′,h,λ) = σ′h W(σ→σ′,h,λ)+B W(σ→−σ′,h,λ) , (14)

which implies that they coincide for all values of .

#### ii.2.2 Generating trajectories for a constant longitudinal field

The above expressions (12),(13) can be interpreted as the normalizing constants of probability measures on the set of piecewise constant functions from to , conditioned on their initial () and final () values. More explicitly, for instance for , the probability of a trajectory with flips at times in the infinitesimal intervals , with , is defined to be

 1W(σ→σ,h,λ)B2neσh(2t1−2t2+⋯−2t2n+λ)dt1…dt2n . (15)

Our goal is now to construct a procedure for actually sampling from these probability measures, that is constructing spin trajectories according to these weights. We shall do this by exploiting the following two identities:

 W(σ→σ,h,λ) = eσhλ+B∫λ0du eσhu W(−σ→σ,h,λ−u) , (16) W(σ→−σ,h,λ) = B∫λ0du eσhu W(−σ→−σ,h,λ−u) . (17)

The path-integral interpretation of these relations, more easily conveyed by the drawing of Fig. 1, is as follows. For the first one, it means that a spin trajectory starting and ending at the same value is either constant on the whole time interval, or made of a constant part upto time , followed by a jump to and a second part of the trajectory representative of . Similarly the second one expresses the necessity for a trajectory from to to have at least one discontinuity at a given time , followed by a trajectory accounting for . These equalities can be proven either from the path-integral representation of (12),(13), or from the explicit expressions of given in Eq. (11). As a consequence of (16),(17) one obtains the following recursive procedure to draw a spin trajectory for a constant longitudinal field on a time interval of length , constrained to , :

• if

• draw a random variable with density proportional to (see below in eq.(18) for some details on how to perform this step)

• set upto time

• call the same procedure to generate a trajectory from to on the remaining interval of length

• if

• with probability , set on the whole time interval, and exit the procedure

• otherwise,

• draw a random variable with density proportional to

• set upto time

• call the same procedure to generate a trajectory from to on the remaining interval of length

In order to draw with a density proportional to , we compute its cumulative distribution,

 G(u) = ∫u0dt eσht W(−σ→−σ,h,λ−t)∫λ0dt eσht W(−σ→−σ,h,λ−t)=1−eσhusinh((λ−u)√B2+h2)sinh(λ√B2+h2) . (18)

A simple way to draw amounts to draw uniformly at random on , and to invert the above expression to obtain . One can proceed similarly for the generation with the density proportional to , which involves another cumulative distribution .

#### ii.2.3 Generating trajectories for a piecewise constant longitudinal field

In the previous subsection we considered the particular case of a constant longitudinal field . Let us now address the general case of a piecewise constant on the interval of imaginary time , with the following definitions illustrated in Fig. 2: we shall call the number of times it changes value between and , the times of these changes, the length of these intervals for , and finally the values the field takes in each of these intervals.

We shall have to compute the partition function of a spin acted on by such a field,

 Z(h)=Tr(p∏i=0eλ(i)(h(i)σz+Bσx)) , (19)

and to generate spin trajectories according to the corresponding weights. The computation of the partition function can be performed by inserting representations of the identity in (19), corresponding to the spin values at the imaginary times where is discontinuous:

 Z(h) = ∑σ0,…,σpZ(σ0,…,σp|h) , Z(σ0,…,σp|h) = p∏i=0W(σi→σi+1,h(i),λ(i)) , (20)

with . Given the trajectory this computation is easily performed, necessiting only the multiplication of the matrices of order defined in (10). The sampling of the spin trajectory on the interval is done as follows. The values of the spin at times are generated222This can be done with a number of operations of order . Let us define . First draw with probability , then according to , and so on and so forth. according to the probability . Then for each interval a spin trajectory from to is generated according to the procedure of Sec. II.2.2, the longitudinal field being constantly equal to on this interval of time. Finally the trajectories are concatenated to obtain the full trajectory from to .

Let us emphasize that the path integral representation of the imaginary time evolution is well known in the literature Feynman53 (); Farhi_Gutmann (); Ginibre (); Gallavotti (); Aiz (); Ioffe (); however we could not find in previous works such an explicit sampling procedure for generating the spin trajectories. Actually, as far as we know, all continuous time quantum Monte Carlo algorithms Beard_Wiese (); worm (); qcluster (); evertz () do not proceed in an heat-bath way, generating “from scratch” a new spin trajectory conditioned on the local effective field, but rather constructing the spin update using the current configuration itself.

## Iii The quantum Bethe lattice ferromagnet

The remainder of this article will be devoted to the study of the simplest of the finitely connected models that can be handled by the quantum cavity method, namely the transverse field quantum spin-1/2 ferromagnet on the Bethe lattice (more precisely on a random regular graph). The physical properties of such a model are very intuitive: at low temperature and transverse field the model is ferromagnetically ordered, with a positive spontaneous longitudinal magnetization. Thermal (increasing ) or quantum (increasing ) fluctuations destroy this order outside a region delimitated by a critical line in the plane, that ends up in a quantum critical point at zero temperature.

An even simpler model displaying these features is the (fully connected) quantum Curie-Weiss model, whose solution we recall in Appendix A. Both of them are of a mean-field nature, and should share most of their qualitative properties, yet the Bethe lattice model is quantitatively different, and technically more involved because of its finite connectivity.

### iii.1 The quantum cavity method treatment

The quantum Bethe lattice ferromagnet is defined by the Hamiltonian

 ˆH=−∑i−jσziσzj−BN∑i=1σxi , (21)

where the first sum runs over the edges of a random -regular graph of vertices rgraphs (). This means that the graph of interactions is uniformly chosen among all graphs on vertices for which each vertex has the same number of neighbors. These graphs have the good properties to realize finite size Bethe lattices: the set of vertices at distance333This is the graph theoretic distance between two vertices and , defined as the length of a shortest path of adjacent vertices joining and along the graph. smaller than a given cutoff from an arbitrarily chosen vertex is, with a probability which goes to one when diverges with held fixed, a regular tree of connectivity . On the other hand such a graph has no surface, contrarily to the usual Cayley tree, and the problem of the boundary conditions for frustrated models (that will be encompassed by the generic treatment of Sec. V) is absent with such a definition.

The hypotheses of the cavity method are more simply explained assuming first that the graph of interactions is actually a finite tree, and that the quantum aspects of the problem have been handled by a finite number of Suzuki-Trotter slices. In such a case it is easy to solve the model exactly by taking benefit of the natural recursive structure of a tree: one breaks the graph of interaction into subtrees that are then glued together. Let us explain this with more precise formulae, defining as the probability law of the configuration of the ring when the interaction with its neighbor has been removed from the graph. If we denote by the set of vertices neighbors of , the recurrence equations for these distributions are:

 μi→j(σi)=1zi→jw(σi)∏k∈∂i∖j∑σkμk→i(σk)eβσi⋅σk , (22)

where is a normalization constant, and we introduced for two arbitrary imaginary time dependent quantities . A graphical representation of this equation is given in Fig. 3 in the case . The site has then two neighbors, and ; the distributions encode the effect on of the part of the tree that does not involve ; this is represented by the dashed line in figure 3. In absence of site the sites and are decoupled (as we assume the graph to be a tree there are no paths between and that do not go across ) and hence independent. Therefore the distribution of in absence of site is given by (22). A proof of such recursion equations and a discussion of the connections with the Bethe-Peierls approximation and the Belief Propagation algorithm can be found in fgraph (); Yedidia (). On a given tree this set of equations (one for each directed edge of the graph) has a unique solution, that is a set of “messages” , that can be efficiently determined by sweeping the edges from the leaves towards the center of the tree. From them all the relevant thermodynamic quantities can be computed. For instance the probability law of the configuration of the ring in the complete tree is given in terms of the messages sent by its neighbors,

 μ(σi)=1ziw(σi)∏k∈∂i∑σkμk→i(σk)eβσi⋅σk , (23)

with a normalization constant, while the joint law for two neighboring sites and reads

 μ(σi,σj)=1zi−jμi→j(σi)μj→i(σj)eβσi⋅σj . (24)

Moreover the free-energy per site can be expressed in terms of the normalization constants of these laws,

 −βf=1NlnZ=1NN∑i=1lnzi−1N∑i−jlnzi−j , (25)

where the second sum runs over the (undirected) edges of the tree.

The above derivation was exact because we assumed the graph of interactions to be a tree. The scope of the cavity method is to extend these results to random graphs which are only locally tree-like, in the precise sense explained above. In its simplest version, called replica symmetric for historical reasons, one assumes the existence of a single pure state in the configuration space of the random graph model. This implies a decay of correlations at large distance in the graph, hence the effect of the long loops neglected in the tree derivation amounts to create a self-consistent boundary condition which traduces the absence of a surface in the random graph. As in the present model all sites have the same neighborhood (there is no fluctuation neither in the connectivity nor in the intensity of the interaction couplings), one has to look for a solution of (22) where the messages are all equal to a single law, that we shall denote in the following , which is seen from (22) to satisfy

 η(σ)=1zlw(σ)(∑σ′η(σ′)eβσ⋅σ′)l . (26)

The assumption on the unicity of the pure state can fail for two kind of reasons. In ferromagnetic models, as the one considered now, there is an ordered phase at low temperature and transverse field in which two pure states coexist because of the up/down symmetry of the model. This is not a serious limitation of the method, in the following it will be kept understood that an infinitesimal longitudinal field is applied to the system in order to select one of the two pure states. In fact the exactness of the cavity method for classical ferromagnetic models on random graphs has been recently proven rigorously ferro_rigorous (). A much more serious problem, that shall not be discussed in this paper, arises in frustrated models, when an exponential number of pure states proliferate at low temperatures; the replica-symmetry breaking version of the cavity method is then required to solve the problem.

Let us first write the solution of the equation (26) in the classical situation, for . In such a case the weight forces all spins along the ring to take the same value, hence

 η(σ) = 1+tanh(βh)2(Ns∏α=1δσα,1)+1−tanh(βh)2(Ns∏α=1δσα,−1) , (27)

where is solution of

 h=lβarctanh(tanh(β)tanh(βh)) . (28)

The classical model thus exhibits a continuous paramagnetic to ferromagnetic transition when the inverse temperature crosses its critical value, .

We now come back to the general case and rewrite explictly the cavity method prediction for the free-energy per site, that is obtained from Eq. (25) by substituting :

 −βf = ln⎛⎝∑σw(σ)(∑σ′η(σ′)eβσ⋅σ′)l+1⎞⎠−l+12ln⎛⎝∑σ,σ′η(σ)η(σ′)eβσ⋅σ′⎞⎠ , (29)

and for the probability law on one site and two adjacent sites,

 μ(σ) = 1zl+1w(σ)(∑σ′η(σ′)eβσ⋅σ′)l+1 , μ(σ,σ′) = 1zedgeη(σ)η(σ′)eβσ⋅σ′ . (30)

Note that the expression (29) is variational, in the sense that its derivative with respect to vanishes on the solution of equation (26).

### iii.2 A convenient representation of η(σ)

We turn now to the problem of the determination of the probability law solution of Eq. (26). As this is a law on the probability of Ising spins, its complete characterization should involve real numbers. Such a direct representation becomes very soon impossible to handle when grows, and in particular in the continuous time limit . We shall however see that an alternative representation allows to bypass this difficulty. First we define a probability law on the configurations of a ring of spins by

 p(σ|h)=1Z(h) w(σ) eβσ⋅h ,Z(h)=∑σw(σ) eβσ⋅h , (31)

ensuring the normalization of . These definitions allow to rewrite Eq. (26) as

 η(σ)=∑σ1,…,σlη(σ1)…η(σl) p(σ|σ1+⋯+σl) Z(σ1+⋯+σl)zl . (32)

Suppose now that one has an estimation of given by a representative weighted sample of elements , that is

 η(σ)=Ntraj∑i=1ai δ(σ−σi) , (33)

such that the weights add up to one. A new estimation of (i.e. a new set of configurations and weights ) can then be obtained by plugging this estimation in the left-hand-side of Eq. (32), which leads to the following procedure. To generate each of the new representants of , one repeats in an independent way the steps:

• draw independently integers in with probability

• set

• generate a configuration according to the law , and set

Once the new elements have been generated one just has to multiply the weights by a global normalization factor, and the new estimates can be again plugged in the right-hand-side of Eq. (32) to approach its fixed point solution.

At this point the continuous time limit can be taken without any difficulty of principle: as long as is finite the spin trajectories have only a finite number of discontinuities, and can then be easily encoded with a finite set of reals corresponding to the time they change values. Moreover one can easily show that in this limit drawing a spin trajectory from the law corresponds exactly to the procedure we defined in Sec. II.2 (we show explicitely this correspondence for the normalization in App. B).

This representation of the probability law by a sample of representative elements is a widespread technique in the field of disordered systems popu_first (); MePa (). We warn however the reader accustomed to the classical cavity method that we did not use it exactly in the usual way, as the classical equivalent of is a single real number, not a probability law.

Before closing this section let us discuss the computation of the physical observables in the continuous limit, taking as representative examples the longitudinal and transverse magnetizations. These can be obtained by taking the limit of the formulas (9),(II.1), which yields

 mz = limN→∞1NN∑i=1⟨σzi⟩=∑σμ(σ)1β∫β0dt σ(t) , mx = limN→∞1NN∑i=1⟨σxi⟩=∑σμ(σ)1βBj(σ) , (34)

where in the right-hand-side is to be computed from Eq. (30), the sums over are understood as sums over continuous time trajectories, and we have defined as the number of discontinuities of on the interval .

### iii.3 Continuous-time results

We present now numerical results for ; the behavior is qualitatively the same for any . We followed the method of resolution presented above, using trajectories. For a fixed temperature, we initialized the population at ; then each trajectory is constant and its value is chosen in such a way that the average magnetization is equal to the classical magnetization that can be obtained from the solution of Eq. (28). In this way we select one of the two possible ferromagnetic states, that we can follow by increasing gradually the transverse field. We increased with a step and for each step we let the population equilibrate for iterations and averaged the observables over subsequent steps. We checked that the weights defined in (33) remain quite uniformly distributed over the population at all investigated temperatures and transverse fields.

In figure 4 we plot the magnetizations and as functions of the transverse field for three different temperatures . At these low temperatures, the system undergoes a continuous phase transition at a critical value of the transverse field. The transition is characterized by the vanishing of and by a jump in the derivative of . Unfortunately obtaining reliable values of close to is a difficult numerical task due to critical fluctuations that cause strong finite size effects in the size of the population . Therefore we located the transition by the jump in the derivative of ; we fitted by linear laws close to from above and below, and determined their intersection.

The resulting phase diagram is reported in fig. 5. We found that for the temperature dependence of all observables is very weak and , that is a reliable estimate for the limit and is in excellent agreement with the value determined in qbl (). The scaling of for is compatible with the essential singularity that is found in the Curie-Weiss model, see Appendix A and in particular Eq. (63). The very weak (practically unobservable) temperature dependence of the free energy below make us confident that the entropic term is negligible at these low temperatures and the free energy for , plotted in figure 6 is representative of the ground state energy . The latter has a weak singularity (discontinuity of the second derivative) at which is not observable with our numerical precision, but is easily observed by a direct computation of , see figure 4.

Overall, our results are in very good quantitative agreement with the ones obtained at in qbl () by a matrix product state description, except for the value of the exponent characterizing the vanishing of the magnetization close to , . Even if we do not have very precise data close to due to finite population-size effects, our data are compatible with the mean-field exponent at all investigated temperatures, while in qbl () a slightly smaller value of is reported at . We will further comment on this discrepancy in section IV.3.

### iii.4 Comparison with approximate treatments

In this section we compare the previous results with approximated solutions of Eq. (26). First we consider the finite case, then we consider variational approximations to the solution.

#### iii.4.1 Resolution at finite Ns

At finite , the distribution can be encoded with real numbers (e.g. the probabilities of each of the configuration, with their sum constrained to be 1). Then Eq. (26) can be rewritten as a fixed point equation for these numbers, and the solution can be found by simple iteration (below we refer to this procedure as “exact solution” for finite ). This method is extremely precise but computationally very heavy, as the time to solve the equation scales exponentially in . We are then limited to ; the computation for the largest took two weeks on a standard workstation, while we estimated the one for to take at least two months.

Therefore, in order to study the approach to the limit , we resorted to the population method described in Sec. III.2. We represented by a population of discrete time trajectories and solved Eq. (26) by iteration following the procedure detailed in section III.2 without taking the limit . Similarly to the continuous time case, we used and performed iterations to achieve convergence, after which data have been collected along iterations. The computation time is now polynomial in so we can go to much larger values (at the price of numerical precision due to finite ). Note that the continuous time computation is much more efficient. In the discrete time case, at low and large , the spins are constant along large time intervals and the information encoded in a discrete trajectory is redundant. On the contrary, at low and large the method at finite is obviously incorrect because of the discretization. One could think to adapt in order to find the optimal value for a given ,; however this is most naturally done in the continuous-time framework where the number of spin flips is the natural variable describing a trajectory. Therefore we think that the finite population method is useful only for the illustrative purposes of this section.

Data reported in this section and in figures 7, 8 have been obtained by exact solution for and by the population method for . We do not report detailed results for the magnetization as a function of ; as expected, the agreement with the continuous time computation is very good at high temperature, and becomes poorer and poorer on lowering the temperature. We focused on two points in the phase diagram; the first at intermediate temperature and ; the second at very low temperature and . In figure 7 we plot the longitudinal magnetization for different values of . As expected, finite effects are much stronger at low temperature. Note that for large the leading correction is , as in the Curie-Weiss model (see Appendix A). For deviations from the leading term are very small and is enough to get a fair estimate of the true . On the contrary, for deviations are very large and is needed.

In figure 8 we report the critical temperature as a function of for different values of together with the continuous time result. Note that in the limit of large transverse field, each of the time slices becomes independent of the others. The model reduces to copies of the classical system () with a ferromagnetic coupling rescaled by . Then in this limit the critical temperature is given by the classical one divided by , . Below this temperature the system is always in the ferromagnetic phase, so that the quantum phase transition cannot be studied within this approximation.

For generic spin- models, the finite approximation might be useful in order to understand the behavior at intermediate temperatures, and might also give quantitatively accurate results far from . This might be useful in cases where more complicated cavity treatments (such as the so-called one-step replica symmetry breaking) are necessary. Still the continuous time method appears to us preferable as it remains reliable down to very low temperatures.

#### iii.4.2 The static approximation

A different strategy to obtain approximate solutions of these models is to use the variational property of the Bethe free energy (29): the free energy of the model is the minimum of the free energy defined by Eq. (29) over the function (in the classical case the correctness of this procedure has been proven in ferro_rigorous ()). Then one can propose a variational form for and minimize the free energy with respect to the free parameters to obtain an upper bound to the true free energy of the problem.

A popular variational form is the so-called static approximation, that in our context amounts to the following ansatz for :

 η(σ)=η[1β∫β0dtσ(t)]=η[1NsNs∑α=1σα] , (35)

i.e. the probability of a trajectory depends only on its average spin value along the imaginary time. This approximation has been widely used in computations based on the replica method for fully connected models BrayMooreQSK (); Nish_stat (); Goldschmidt (); Leticia () and has been recently applied to the random -satisfiability in a transverse field qksat (). The present derivation based on the cavity method gives back the same results originally derived in qksat (), but is simpler because the use of replicas is avoided. Here we discuss only the case of the ferromagnet on a regular graph but the equations can be easily generalized to more complicated cases where fluctuations of the couplings and/or the connectivity are present qksat ().

Equivalently we can rewrite the equation above as

 η(σ) = ∫∞−∞dhp(h)Ns∏α=1[ehσα2coshh]=∫1−1dmp(m)Ns∏α=1[1+mσα2] , (36)

where with a slight abuse of notation we used for the distribution of both the effective field and its associated magnetization . This expression makes evident that the assumption of the static approximation is that the spins in a trajectory are uncorrelated along the imaginary time, and subject to a common field extracted from the distribution .

From Eq. (36) it follows that, in the limit :

 (37)

where in the last line we defined the function

 ws(m)=∑σw(σ)δ(m−1NsNs∑α=1σα) . (38)

The following explicit expression for can be found in qksat () (where it was called ):

 ws(m)=βB√1−m2