# Quantum Cluster Variational Method and Message Passing Algorithms Revisited

###### Abstract

We present a general framework to study quantum disordered systems in the context of the Kikuchi’s Cluster Variational Method (CVM). The method relies in the solution of message passing-like equations for single instances or in the iterative solution of complex population dynamic algorithms for an average case scenario. We first show how a standard application of the Kikuchi’s Cluster Variational Method can be easily translated to message passing equations for specific instances of the disordered system. We then present an “ad-hoc” extension of these equations to a population dynamic algorithm representing an average case scenario. At the Bethe level, these equations are equivalent to the dynamic population equations that can be derived from a proper Cavity Ansatz. However, at the plaquette approximation, the interpretation is more subtle and we discuss it taking also into account previous results in classical disordered models. Moreover, we develop a formalism to properly deal with the average case scenario using a Replica-Symmetric ansatz within this CVM for quantum disordered systems. Finally, we present and discuss numerical solutions of the different approximations for the Quantum Transverse Ising model and the Quantum Random Field Ising model in two dimensional lattices.

## I Introduction

Exact solutions of problems involving many interacting particles in finite dimensional systems are very difficult to find. This is particularly true in specific complex situations where disorder is present, like in the classical Edward-Anderson model, or the quantum Anderson transition Anderson (1958). When dealing with classical systems a very successful approach in the last few years has been the well known Cavity Method Mézard and Parisi (2001, 2003); Kabashima (2003), which is exact for models designed on a tree or a random graph. Moreover, it is possible to show that the method corresponds to the Bethe approximation of the free energy for a model defined in a finite dimensional latticeYedidia et al. (2005a); Pelizzola (2005), and that it is intimately connected with message passing algorithms in single instances of a specific problemKschischang et al. (2001); Braunstein et al. (2003); Mézard and Zecchina (2002); Montanari et al. (2008).

Although the success of the method in many models is undisputed Mézard and Montanari (2009), it took some time to understand how to improve over this Bethe approximation for a finite dimensional disordered system Rizzo et al. (2010); Yedidia et al. (2005a, b); Chertkov and Chernyak (2006a, b); Pelizzola (2005); Kappen (2002); Zhou et al. (2012); Xiao and Zhou (2011). The idea behind most of these improvements is based on a Cluster Variational MethodMorán-López and Sánchez (1995), that applied to specific instances of the problem usually leads to what is called Generalized Belief Propagation (GBP) algorithm Domínguez et al. (2011); Lage-Castellanos et al. (2011a, b, 2013). However, it is also possible to use the same approach to the replicated free energy and then to choose the RS ansatz or the more general Parisi’s hierarchical ansatz to send the number of replicas to zero Rizzo et al. (2010); Lage-Castellanos et al. (2013); Ferraro et al. (2016). Remarkably, also in this more general approximation one encounters a connection between average case predictions and the behavior of message passing algorithms in single instances of the same problem Lage-Castellanos et al. (2013); Ferraro et al. (2016); Lage-Castellanos et al. (2014); Dominguez et al. (2015).

For quantum disordered models such a comprehension is still lacking, and even the simpler quantum models on random graphs require effort and are actively studied. A clear breakthrough in that direction was the use in Krzakala et al. (2008) of the path-integral representation of quantum spin models that was suitable to derive a closed equation where the histories of the spins constitute the proper variables to iterate. As in classical systems the method is exact in trees, however, it is much more demanding from the computational point of view. Because of its clear connection with the Cavity Method in classical systems it is usually coined Quantum Cavity Method. Somewhat similar in spirit is the approach followed in Ioffe and Mézard (2010a, b); Dimitrova and Mézard (2011) where the quantum model is also defined in tree-like structures but is properly parametrized to simplify the solution of the corresponding closed equation. The approach can also be understood as a further approximation within the Quantum Cavity Method presented by Krzakala et al. (2008).

Moreover, again within a Bethe approximation in Ramezanpour (2012); Biazzo and Ramezanpour (2013, 2014) the authors proposed a message-passing algorithm to compute the Hamiltonian expectation of a quantum disordered system. The technique rests on the use of an appropriate trial wave function and the connection of quantum expectations to average quantities in a classical system with both local and global interactions. On the other hand, starting from a more algorithmic point of view the authors in Evans and Stephens (2008); Leifer and Poulin (2008); Poulin and Bilgin (2008); Bilgin and Poulin (2010); Poulin and Hastings (2011); Jouneghani et al. (2014) proposed different versions of what they called Quantum Belief Propagation (QBP) algorithm suitable also to be used in single instances of disordered systems with tree-like topology. Here again approximations are unavoidable to guarantee a polynomial performance of the algorithm.

These approaches, however, rest on Bethe like models or Bethe like approximations to finite dimensional systems and are difficult to extend to more general scenarios. Fortunately, already Morita and Tanaka Morita (1957a, b); Tanaka et al. (1994); Katori and Suzuki (1994); Morita and Tanaka (1994) developed a general approach to derive a set of self-consistent equations within a Cluster Variational Method formalism. However, as far as we know these equations were mainly studied in homogeneous systems and only formally presented for disordered models.

In this work we first take profit of this derivation and the parametrization done in Dimitrova and Mézard (2011) to study single instance implementations of message passing algorithms for quantum disordered systems at both the Bethe level and the plaquette approximation. Moreover, inspired by previous results on classical models Rizzo et al. (2010); Lage-Castellanos et al. (2013); Dominguez et al. (2015), we propose a set of fixed point equations to describe these quantum models without the specification of the instance and compare both approaches. At the Bethe level these equations are a generalization of the dynamic population equations already derived in Morita (1957a, b). At the plaquette level, to our knowledge, this is the first time that these equations are presented for quantum models. Finally we present how to properly generalize the cluster variational method for quantum disordered models within the RS-ansatz in an average case scenario. These approaches were tested in the Quantum Ising model in a transverse field and in the Quantum Ising model with a random field.

The rest of the work is organized in the following way. In the next section we present the models used to test the different approximations. This will make clear the kind of models of interest and help to simplify the notation below. Then we will present the Cluster Variational Method in the form derived by Morita and Tanaka Morita (1957a, b); Tanaka et al. (1994); Katori and Suzuki (1994); Morita and Tanaka (1994) in their seminal works but using a modern notation and connecting it to message passing equations. We do this for the Bethe and the plaquette approximation. Then we show how to derive similar equations but taking into account the average over the disorder. Finally we present and discuss the numerical results of the different approximations for the two models under study.

## Ii Models

The Quantum Ising Model in a transverse field is one of the most basic models displaying quantum phase transitions. Under this name we identify a general class of models described by the Hamiltonian:

(1) |

where the choice of the interacting pairs defines the underlying lattice. The properties of this system obviously depend on the topological structure as well as on the distribution of interaction constants and local fields . Popular choices for the fields and interactions include the ferromagnetic case, where is uniform and positive everywhere and the local field is also uniform. More interesting is the so called Random Field Ising Model (RFIM) where interactions are homogeneous too but local fields fluctuate from site to site according to a given distribution .

In contrast to the classical version, the local fields applied perpendicularly to the easy direction of the material included in Eq. (1) result in non-commuting terms that make the Hamiltonian difficult to diagonalize. The basis with a definite direction of each spin is no longer an eigenstate of the system, meaning that transverse fields introduce quantum fluctuations that can destroy any long range order for sufficiently strong field intensities.

An exact solution for this kind of systems has been found only for special cases, mainly mean field and fully connected models or very particular tree-like topologies such as a random latticeKrzakala et al. (2008). For the latter, the Quantum Cavity Method is the most successful technique, but solutions for finite dimensional lattices remain elusive, and this Quantum Cavity Method represents just an approximation at the Bethe level of the actual problemIoffe and Mézard (2010a, b); Dimitrova and Mézard (2011). In finite dimensional lattices a quantum extension of the CVM is suitable because of the geometric and periodic character of the structure. We follow this line of thought in the next sections to obtain a description of both, the ferromagnetic and RFIM at the level of the Bethe (pairwise) and Kikuchi (plaquette) free energy approximations for a bi-dimensional square lattice. Extensions to other lattice models should follow straightforwardly.

## Iii Quantum Cluster Variational Method: Single Instances and Message Passing

The Cluster Variational Method Kikuchi (1951); Pelizzola (2005); Yedidia et al. (2005a); Morita (1957a, b); Tanaka et al. (1994); Katori and Suzuki (1994); Morita and Tanaka (1994); Morán-López and Sánchez (1995) relies on a constrained minimization procedure of a region-based free energy functional. As a result, it is possible to get estimates for the local probability distributions of the system under study; distributions that are otherwise hard to obtain because of the well known difficulty of tracing the complete distribution over the exponential number of allowable states in the large limit.

For quantum problems the equivalent of local probability distributions are the projections of the full system density matrix on the subset of variables of each region. With these local density operators it is possible to write a region free energy and repeat essentially the same classical minimization.

### iii.1 Bethe approximation

In the specific case of pairwise interactions at the Bethe level, the variational parameters of the problem are the pair and site density operators, denoted and respectively. The Bethe free energy is thus written in terms of these operators as a sum of weighted contributions of all pairs and sites of the interaction network or lattice:

(2) |

where

(3) |

can be regarded as the free energy of the pair and the other term represents the contribution of the individual sites:

(4) |

The prefactors and take integer values such that the contribution of each variable to the total free energy is counted only once. These are the two first terms in a cumulant expansion of the total free energy Tanaka et al. (1994). For the 2D square lattice and . The minimization of is performed under a set of constraints enforcing normalization and consistency of the set of operators and :

(5) | |||||

We can use now Lagrange multipliers to put together the objective function and the restrictions (5). Up to this point there is no particular relevance in the choice of the basis for the Hilbert state space. In fact, the Lagrange function for this problem can be nicely written in a basis-independent way:

In this last equation (III.1) the different ’s are real numbers and are Hermitian operators acting on the single site Hilbert spaces. The stationarity condition for is obtained by setting to zero the linear part in of the increment

The resulting equations are the operator version of the belief propagation (BP) Yedidia et al. (2005a) equations for the local distributions in terms of the Lagrange multipliers:

(6) | |||||

(7) |

It is customary to make the linear transformation:

(8) |

where the sum includes all links containing , denoted as , except the link itself. This substitution gives the more familiar BP-like structure:

(9) | |||||

(10) |

The operator represents the effective interaction of the spin with its neighbor , both forming the link . Since the site-site interaction in (1) is directed along the direction we chose to parametrize this operator as . This form resembles the projected cavity solution in Dimitrova and Mézard (2011). In this work the authors make a recursive construction in a tree to obtain a self-consistent equation for a cavity field. In our parametrization, (without a hat) can be interpreted again as a kind of cavity magnetic field, in the sense that it stands for the interaction of spin with the portion of the network growing in the direction of link .

The set of fields can be determined by a fixed point iteration after plugging (9) and (10) into the consistency conditions in (5). The procedure is analogous to the use of the BP algorithm in the classical case. This time though, we have a set of coupled operator equations:

(11) |

where the (local) partition functions are fixed from the normalization condition . Spin , for example, will have the following normalization:

(12) |

Then, instead of working in (11) directly with operators that are hard to translate into actual numerical values it is convenient to write the consistency between them by matching their moments. This is, we make use of the relation a relation that is sufficient for our purposes of finding .

The algorithmic procedure to solve for all the values in a given lattice is the following. First, select at random a link region . Then, using all the and acting on each one of the spins from outside , form the local density operator and find the magnetization of spin , . Now focus on the expression for . All cavity fields in the exponent were used in the link equation except precisely , the cavity field of link on spin . Its value is now found numerically as the root of the equation . This step is repeated many times on different links and spins pairs until the value of changes less than a certain tolerance everywhere in the lattice. The process is schematically explained in Fig.1. In thinner red lines appear the messages taken at step to find a new value (green, thick line) at step . The reader familiar with the BP algorithm may notice that in this case we need to keep the cavity fields on during the calculation, whereas in the classical case they can be canceled out in both sides of (11) due to the commutation properties of the effective Hamiltonian.

### iii.2 Kikuchi approximation

For finite dimensional lattices it is important to take into account explicitly the existence of short loops, which are completely disregarded by the Bethe approximation Yedidia et al. (2005a). A sound improvement of the Bethe choice of regions could be obtained by including larger regions into the free energy. The simplest generalization for a square lattice is precisely the inclusion of plaquette regions, each formed by the four spins of the elementary cell. The Kikuchi free energy will have an extra term with respect to (2), comprising the contribution of all plaquettes:

(13) |

The free energy of a plaquette is defined similarly to (3) and (4) by means of a plaquette density operator . The prefactors take the values , and in this case. Constrained minimization of is technically similar to the Bethe case except for the use of some extra Lagrange multipliers that enforce marginalization of plaquettes over link distributions. These new multipliers are Hermitian operators acting on the two-spins Hilbert space corresponding to the spins in each link. They represent an effective -directed interaction of the spins in the link with the rest of the lattice. Hence, we write them in the following way:

(14) |

This is not the most general expression for an operator in the product space of two spins. This choice is based on the fact that the spin to spin interaction lies always in the OX direction. The value of the Lagrange parameters should be obtained from the marginalization conditions (5) and the equivalent for plaquettes. For the local density operators of each region we get:

(15) | |||||

(16) | |||||

(17) |

where are indexes corresponding to plaquette, link and site regions respectively. For a region the set contains all its sub-regions and the set is populated with all the regions to which belongs. For example, for the link , contains the two spin regions and , whereas includes all the plaquettes intersecting on .

The fixed point iterations are again performed via moment matching. This is, the magnetization predicted by must be consistent with and this last operator must also produce the same magnetization and correlation as the parent plaquettes :

(18) | |||||

(19) | |||||

(20) | |||||

(21) |

The set (18)-(21) gives the solution for the cavity fields implicitly. For future use it is convenient to formally define the functions that would give this fields in a explicit way:

(22) | |||||

The algorithmic details for the implementation of the fixed point iteration are very similar to the Bethe case. We take a random plaquette with the corresponding external cavity fields at step (see Fig. (2)) and evaluate the RHS of (19) and (21) for every connected pair and single spin in the plaquette. Then, by a fixed point iteration, we make the LHS of the same equations be consistent with the plaquette prediction. The consistency equation (18) is used too at this step. This way we find the internal fields of the plaquette at step in terms of the external ones at .

Once the approximated marginal distributions are known, the local contributions to the free energy can easily be obtained through the normalization factors of (15)-(17):

(23) |

Finally let us discuss the connection of the formalism to the classical GBP. The classical results are obtained effortlessly from (15)-(21) when the transverse field is zero. In these cases, all terms in the exponentials conmute and the cavity fields that appear in both sides of the consistency equations can be cancelled out. When conmutation is important though, the terms do not cancel and most be tacking into account.

## Iv Quantum Cluster Variational Method: Average case scenario

In the previous sections we discussed how to deal with single instances of disordered systems withing a Quantum Cluster Variational method. More frequently in physics one expects to be able to average over the disorder right from the beginning. Inspired by the replica-CVM methodology introduced in Rizzo et al. (2010) for the classical version we now try to perform the average case calculations for disordered Ising quantum models.

The replica trick is a general framework that in principle allows to account for different degrees of complexity in the structure of the state space. Here we present the solution of the problem within the Replica Symmetric (RS) approximation that assumes that a single state dominates the thermodynamics of the problem. Within this approximation we studied the RFIM at the level of the Bethe and Kikuchi approximations starting from (1) with being i.i.d random variables in the interval . The extension to more general ansatz follows directly from the work Rizzo et al. (2010) and the approach presented here.

As a starting point let us consider an alternative approach to find the CVM free energy of a single instance. With a given realization of the disorder one way to find the region free energy of the system is to minimize a variational expression that is equivalent to the Lagrange function discussed in previous sections:

(24) | |||||

(25) |

Here is the subset of cavity fields appearing in the expressions for , the normalization constant of each region distribution. The minimization process makes this set of fields depend in principle on all the external parameters :

(26) |

The average free energy density is now defined in the thermodynamic limit as the free energy per spin after averaging over the distribution of :

(27) |

Here stands for the Bethe or Kikuchi approximations and defined in (2) and (13) or any other valid region based free energy approximation. Putting (26) into (27) we see that the problem now reduces to averaging the logarithm of the local partition functions:

(28) |

This is not an easy task, mainly because we do not know the analytic dependence of the cavity fields on the realization of the disorder; remember that the cavity fields are found via a fixed point iteration. To overcome this problem let us take a step back, return to the variational character of the expression (24) and define:

(29) |

The difference between (28) and (29) is that in the latter the fields are free parameters to be optimized. Notice also that each region depends only on the external fields acting locally. Since the cavity fields are independent parameters (to be fixed later by a minimization process) the average over all the set reduces to only the fields in the region . Another useful manipulation is to split the sum over according to the kind of region . Here is an index that goes over the types of regions used in the approximation and labels different regions of the same kind. Also, let us define as the number of regions of type in the system. This way we arrive to an expression that has the form of an average over the disorder and the cavity fields:

(30) | |||||

(31) |

Notice that the sum in (31) is not extensive anymore; for the Bethe and Kikuchi approximation it has 2 and 3 terms respectively. This is now a functional on the cavity field distribution and at the end of the calculations it should be minimized. For each region type the set includes a different number of link-to-spin and plaquette-to-link fields; in order to make this explicit it is convenient to use the notation .

The RS anzats implies that the cavity fields in the expression for are described by a certain distribution :

(32) |

This is the same expression that is obtained by the explicit replica symmetric calculations in Rizzo et al. (2010). Symmetry breaking considerations may define an average (32) that include distributions of distributions in the case of 1RSB or more levels in general.

A further assumption in our RS calculation is that the joint is factored out in terms of simpler distributions. In general we consider that each in the terms is independently distributed with probability and the same for the fields in , distributed according to . This simplifications are justified for Bethe lattices but for finite dimensional problems remain an approximation. The disordered fields are distributed independently so their distribution factorizes too.

Our task now is to find the and such that the functional reaches a minimum. We have to solve the stationarity conditions:

(33) | |||||

(34) |

### iv.1 Bethe approximation

Let us start by writing down with detail the Bethe average case calculation. The (variational) average free energy for a network with fixed connectivity is, using (31):

(35) |

The averages in (35) are over the external field and the cavity fields. The explicit expressions for each term are:

(36) |

and

(37) |

The partition function for the link and spin regions are as usual defined as the trace of the Boltzmann factor for the corresponding effective hamiltonian:

(38) | |||||

(39) | |||||

(40) | |||||

(41) |

The important quantity in (36) and (37) are the distributions and . In the classical version of this calculation one assumes that cavity fields are uncorrelated and these distributions factorize. This is literally true for random networks and only an approximation for lattices with short loops. The factorized forms we consider are:

(42) |

and

(43) |

The next step is to plug everything into the free energy and use the stationarity condition . A simple functional derivative shows that the minimization implies that:

(44) |

where the explicit dependence on means that this variable is not averaged out. Differentiating both sides we get a nice relation between the average magnetization predicted by the link and the spin terms:

(45) |

Multiplying both sides by and integrating to average out also :

(46) |

where averages on left and right hand sides are done using expressions similar to (37) and (36) respectively:

(47) |

and

(48) |

We have used the shorthand and to increase readability of the expressions.

Let us focus for a moment on the magnetization functions and . In an actual lattice, if both are referred to the same spin they must have the same value. This is a consequence of the consistency relations (18). The function depends on fields, of which are also arguments of . The extra field can be obtained from the condition . This is an implicit equation that we can solve. Formally this solution is written as:

(49) |

Using the previous definition we can now relate the link magnetization to the spin one:

(50) |

Putting (46), (47), (48) and (50) together we get an expression that allows the determination of by means of a population dynamics scheme:

(51) | |||||

The above equation has the following interpretation: the LHS represents the average magnetization of a spin obtained by sampling the cavity fields from their distribution. The RHS, on the other hand, represents also the average magnetization of a spin but calculated by taking fields from their distribution and the other fixed by the delta function to be consistent with the link to spin marginalization.

From (51) we can obtain a numerical approximation for using a population dynamics method. The idea is to represent by a list of field values, that is, one typical sample of N values taken independently from . In order to obtain the right we perform a sampling process that simulates (51). First, two groups of fields are selected randomly from the list. Consider that each set acts on one of the spins of an hypothetic link region and find the magnetization of the two spins. In this step, the disordered external fields need to be sampled too. Once we have the magnetization predicted by the link for the spin, say , we demand the spin region predict the same magnetization as the link. This fixes the total cavity field on the spin. Finally, substracting the values initially sampled from the total cavity field we get the effective cavity field as represented in (49). This value is returned to the list of fields in a random position. The convergence of this procedure is monitored by following the evolution of the first two moments of the list. Once the process has converged, the distribution satisfying (51) is obtained as a histogram of the sample list. Algorithm (1) shown below summarizes the method to follow:

### iv.2 Kikuchi approximation

The formalism for the plaquette approximation is essentially the same. We will assume for definiteness a 2D configuration but the results are easy to extend to more dimensions or other kind of lattices, for example, triangular ones. The intensive variational free energy according to (31) is:

(52) |

Below we include for clarity the resulting equations for all the for a spin region , a link and a plaquette :

(53) |

(54) | |||||

(55) | |||||

To lighten the formulas above we have used again the convention and . Notice that the field probability distribution of each region is factored in terms of single and . The real interactions are put together with the cavity ones into an effective hamiltonian that includes all the terms of the Bethe case plus the plaquette-to-link fields:

Now we use the stationarity conditions (34) to obtain the relation between the first and second average moments (i.e. the magnetization and correlation) predicted by each region:

(56) | |||||

(57) |

From the first equality in (56) and repeating the steps for the Bethe case we get an expression for the distribution:

(58) | |||||

In (58) the function is the one which gives the effective field that makes the magnetization predicted by the spin consistent with the magnetization predicted by the link. Compared to the Bethe case it now includes the dependence on the plaquette-to-link messages .

Let us work now with the rightmost equation of (56) or, equivalently, with (57). These involve the consistency between plaquettes and links. It is to be noted that the single-instance update equations based on a plaquette to link marginalization alone do not determine the value of the fields . Instead, only the sum and is completely specified. As a a consequence, equations (56) or (57) will not give us an expression for . Alternatively, we get an equation for the joint distribution of the correlation field and the sum of the magnetization fields and . This distribution is defined as the convolution of the original and ,

(59) |

and obeys an equation that is structurally very similar to (58):

(60) | |||||

In the above equation, the symbol in the RHS stands for all the cavity fields acting on the plaquette from neighboring regions. It includes
also the local magnetic field on each spin. The symbol includes a subset of ; just those fields acting on the link .
The functions , and give the effective correlation and magnetization fields on link due to the interactions in plaquette .
^{1}^{1}1The expressions , and
are again only formal representations of the result of the self-consistent determination of the fields.

## V Numerical Results

### v.1 Quantum Transverse Ising Model

For an homogeneous system , in an homogeneous transverse field , the numeric solution of the update equations in single instances simplifies significantly. The model at equilibrium should be described by a fixed point iteration of a single combination of the parameters , and . Therefore we can take only one pair of plaquette-to-link and link-to-spin marginalization equations and iterate them recursively until the combination reaches a fixed point. To simplify the notation we also drop specific spatial indexes and write only .

Our description below makes emphasis on the plaquette approximation considering that the Bethe case is extensively presented in a vast literature. All the same, the reader interested only on the Bethe approximation can formaly put and to zero and iterate only the link-to-spin marginalization condition.

The sequence of steps is the following. First, the triad is initialized to some arbitrary real values. Then these fields are used to evaluate the moments of of an imaginary plaquette:

(61) | |||||

(62) |

Calculations are made only for one spin and one link for symmetry reasons. As a consequence of (18) and (19) a new value for the field is generated such that spin magnetization equals . Explicitly, the equation that must be solved is:

where is the connectivity of the spin in 2D. After, new values and are obtained from the LHS of (19) and (21). Their value must make the magnetization and correlation from consistent with (61) and (62). We have:

(63) | |||||

(64) |

and then we would have to solve:

(65) | |||||

(66) |

In the equations above the value of from the previous step is used when solving for . This algorithm is repeated until stability is reached i.e. until the variation of the field values drop below certain prefixed threshold. For each temperature and/or field it is convenient to use as initial values the results obtained for a nearby point in the phase diagram. This improves the convergence speed significantly.

To compare this result with actual message passing equations in single instances we studied a 16x16 square lattice with homogeneous field and periodic boundary conditions iterating (18)-(21) starting from random initial conditions and following a random update scheme until convergence (SI). Since the system is homogeneous, there is no need of running a large number of instances nor using a large system; field values tend to the same value everywhere. The only difference between samples would be the initial conditions and the actual random update order. In this case we averaged results for 10 different initial field configurations.

Although the population dynamics (PD) solution to the problem is introduced properly in the next section, since it is mainly relevant for disordered systems, its application to this model is shown here for completeness. Broadly speaking, it is similar to the FP but focuses on the stability of a population of cavity fields instead of a single set of values. The population of field values is supposed to represent the distribution of fields for the average case scenario. In the homogeneous system populations will be represented by a single value, this is, distributions are delta-shaped around the fixed point fields.

The numerical results obtained for the three methods are shown in Fig.(3). In this figure we present the H-T phase diagram of the Quantum Ising model in a transvere magnetic field for the Bethe and Kikuchi approximations.

From Fig.(3) we observe that in all the approximations we get a line dividing a paramagnetic solution where the spontaneous magnetization in the direction is zero from another region where long range order dominates. For low transverse field the transition temperature coincides with the classical case prediction. We also find that above a critical value of the external field quantum fluctuations destroy the possibility of ferromagnetic order at any temperature. The results for the Bethe case are a lot less noisy that for the plaquette approximation. It is interesting that in the latter all methods find a region for intermediate values of the external field where there is a gap of non-convergence between the paramagnetic region an the ferromagnetic one. It is not clear to us whether it is a numerical problem or if it is an intrinsic property of the approximation.

In Fig.(4) we show vertical cuts of the phase diagram taken at and . The behavior of the longitudinal magnetization is qualitatively equivalent to the classical ferromagnetic case. In the transversal direction the system presents always a magnetization in the same direction of the applied field. In these plot we see again that for small fields the critical temperature depends weakly on , taking values very close to the classical one. As the external fields increases, the system needs to lower the temperature to establish the long range order. This can be done up to a certain critical field, above which quantum fluctuations destroy the possibility of an ordered phase.

### v.2 Quantum Transverse Ising model in a Random Field

To numerically approach the average case scenario at the Kikuchi level is not as simple as in the Bethe case a result that is already known from classical modelsRizzo et al. (2010). The problem relies in the proper equations (58) and (60) that do not define a simple closed equation for and to obtain it obtaining by deconvolving and is numerically challenging. Further complications may also arise from the fact that need not to be positive definite. Nonetheless, one can still study the average properties of the update equations in the Kikuchi approximation using population dynamics. The result, though, will be a solution of (58) and (60) only in the paramagnetic regime. In the case of distributions with permanent magnetization, it constitutes only an heuristic tool.

In a 2D square lattice, a plaquette has four other neighbor plaquettes with which it shares a link. From each of these regions it interacts via one plaquette to link triplet and two link-to-spin fields , , five fields in total. Given all those external messages one can, using the update equations, find the fields inside the plaquette. In order to keep as much information as possible we define a population representing the joint distribution of the five values mentioned before. Following the scheme for the Bethe case, we sample the surroundings of the plaquette, calculate a new set of messages and return it to the population. Once the population stabilizes, all the observables can be found by sampling repeatedly the resulting distribution.

In Fig.(5) we compare the results of using the population dynamics algorithm and single instance simulations for the RFIM. Similar to the ordered case, the phase diagram of the magnetization is divided in two regions, para and ferromagnetic. The classical limits of low fields are in agreement with the previously known results and of course with the corresponding values in Fig.(3). For the SI calculations, 100 samples of a 32x32 square lattice with periodic boundary conditions are averaged. For high values convergence is an issue for both PD and SI simulations. Also, in this region the longitudinal magnetization is rather small in the ferromagnetic region. We did not managed to observe a critical value as in the ordered model.

The shape of the field distributions on the lattice changes for the para or ferromagnetic phase. In the paramagnetic region the magnetization fields, and distribute as delta functions around zero, see for example Fig.((a)a). The correlation fields are also well centered around a given value for high temperatures, see Fig.((b)b). On the other hand, inside the ferromagnetic phase, due to the heterogeneous local fields in the direction, we observe that all distributions spread suggesting the possible existence of a glassy phase.

## Vi Conclusions

In this work, we first re-derived the equations for the Cluster Variational Method for models involving quantum phase transitions. Starting from a variational expression for a region based free energy we managed to find approximations to local probability distributions. The minimization of the region free energy is somewhat hindered by the quantum nature of the hamiltonian and the non-conmutativity of the operators appearing on it. As a consequence the cavity fields of the classical models transform in our approach into hermitian operators, parametrized by Pauli matrices. We then approximate the problem transforming these equations for operators into an approximate set of equations for the parameters describing the density operators defining the variational method.

This quantum-CVM is a good framework for studying finite dimensional models. For ordered systems the standard approach exploits the translational symmetry and reduce the problem to the determination of a handful of parameters a technique very well known in the literature. On the other hand, for disordered models we were able to transform the consistency relations imposed between overlapping regions into proper message passing equations that can be treated in polinomial time. We showed by studying the Quantum Ising model in a transverse uniform external field, that both approaches are equivalent when disorder is absent. When disorder is present, like in the Quantum Ising model in a transverse random external field, the message passing equations derived here become nevertheless a very efficient computational approach to study the properties of the model.

In a more general setting, in this work we also presented a version of the CVM for quantum models within an average case scenario, i.e. where the average over the disorder is done without specifically treating single instances. Although this generalization translates into a very complex set of population dynamic equations between operators, we can approximate them through complex populations of physically sound parameters, here magnetization and correlations that can be solved using a variation of standard techniques. The results of all the approaches were compared studying the Quantum Ising model in a transverse random external field.

## References

- Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958), URL https://link.aps.org/doi/10.1103/PhysRev.109.1492.
- Mézard and Parisi (2001) M. Mézard and G. Parisi, Eur. Phys. J. B 20, 217 (2001).
- Mézard and Parisi (2003) M. Mézard and G. Parisi, J. Stat. Phys. 111, 1 (2003).
- Kabashima (2003) Y. Kabashima, Journal of the Physical Society of Japan 72, 1645 (2003), URL http://jpsj.ipap.jp/link?JPSJ/72/1645/.
- Yedidia et al. (2005a) J. Yedidia, W. T. Freeman, and Y. Weiss, IT-IEEE 51, 2282 (2005a).
- Pelizzola (2005) A. Pelizzola, J. Phys. A 38, R309 (2005).
- Kschischang et al. (2001) F. R. Kschischang, B. J. Fret, and H.-A. Loeliger, IEEE Transf. Inf, Theory 47, 498 (2001).
- Braunstein et al. (2003) A. Braunstein, R. Mulet, A. Pagnani, M. Weigt, and R. Zecchina, Phys. Rev. E 68, 036702 (2003).
- Mézard and Zecchina (2002) M. Mézard and R. Zecchina, Phys. Rev. E 66, 056126 (2002).
- Montanari et al. (2008) A. Montanari, F. Ricci-Tersenghi, and G. Semerjian, J. Stat. Mech. p. P04004 (2008).
- Mézard and Montanari (2009) M. Mézard and A. Montanari, Information, Physics and Computation (Oxford University Press, Cambridge, 2009).
- Rizzo et al. (2010) T. Rizzo, A. Lage-Castellanos, R. Mulet, and F. Ricci-Tersenghi, J. Stat. Phys. 139, 375 (2010).
- Yedidia et al. (2005b) J. Yedidia, W. T. Freeman, and Y. Weiss, IT-IEEE 51, 2282 (2005b).
- Chertkov and Chernyak (2006a) M. Chertkov and V. Y. Chernyak, Phys. Rev. E 73, 065102 (2006a).
- Chertkov and Chernyak (2006b) M. Chertkov and V. Y. Chernyak, J. Stat. Mech. p. P06009 (2006b).
- Kappen (2002) H. J. Kappen, in In Modeling Bio-medical signals (World Scientific, 2002), pp. 3–16.
- Zhou et al. (2012) H. Zhou, C. Wang, J.-Q. Xiao, and Z. Bi, J. Stat. Mech. L12001 (2012).
- Xiao and Zhou (2011) J.-Q. Xiao and H. Zhou, J. Phys. A 32, 425001 (2011).
- Morán-López and Sánchez (1995) J. L. Morán-López and J. Sánchez, Theory and Applications of the Cluster Variation and Path Probability Methods (Plenum Press, 1995).
- Domínguez et al. (2011) E. Domínguez, A. Lage, R. Mulet, F. Ricci-Tersenghi, and T. Rizzo, J. Stat. Mech. p. P12007 (2011).
- Lage-Castellanos et al. (2011a) A. Lage-Castellanos, R. Mulet, F. Ricci-Tersenghi, and T. Rizzo, Phys. Rev. E 84, 046706 (2011a).
- Lage-Castellanos et al. (2011b) A. Lage-Castellanos, R. Mulet, F. Ricci-Tersenghi, and T. Rizzo, Journal of Statistical Mechanics 2011, P12007 (2011b).
- Lage-Castellanos et al. (2013) A. Lage-Castellanos, R. Mulet, and F. Ricci-Tersenghi, Journal of Physics A: Mathematical and Theoretical Physics 46, 135001 (2013).
- Ferraro et al. (2016) G. D. Ferraro, C. Wang, H.-J. Zhou, and E. Aurell, Journal of Statistical Mechanics: Theory and Experiment 2016, 073305 (2016), URL http://stacks.iop.org/1742-5468/2016/i=7/a=073305.
- Lage-Castellanos et al. (2014) A. Lage-Castellanos, R. Mulet, and F. Ricci-Tersenghi, Europhysics Letters 107, 57011 (2014).
- Dominguez et al. (2015) E. Dominguez, A. Lage-Castellanos, and R. Mulet, Journal of Statistical Mechanics 2015, P07003 (2015).
- Krzakala et al. (2008) F. Krzakala, A. Rosso, G. Semerjian, and F. Zamponi, Physical Review B 78, 134428 (2008).
- Ioffe and Mézard (2010a) L. B. Ioffe and M. Mézard, Physical Review Letters 105, 037001 (2010a).
- Ioffe and Mézard (2010b) L. B. Ioffe and M. Mézard, Physical Review E 82, 184534 (2010b).
- Dimitrova and Mézard (2011) O. Dimitrova and M. Mézard, Journal of Statistical Mechanics: Theory and Experiment 2011, P01020 (2011), URL http://stacks.iop.org/1742-5468/2011/i=01/a=P01020.
- Ramezanpour (2012) A. Ramezanpour, Physical Review B 85, 125131 (2012).
- Biazzo and Ramezanpour (2013) I. Biazzo and A. Ramezanpour, Journal of Statistical Mechanics: Theory and Experiment 2013, P04011 (2013).
- Biazzo and Ramezanpour (2014) I. Biazzo and A. Ramezanpour, Physical Review E 89, 062137 (2014).
- Evans and Stephens (2008) Z. W. E. Evans and A. M. Stephens, Physical Review A 78, 062317 (2008).
- Leifer and Poulin (2008) M. Leifer and D. Poulin, Annals of Physics 323, 1899 (2008).
- Poulin and Bilgin (2008) D. Poulin and E. Bilgin, Physical Review A 77, 052318 (2008).
- Bilgin and Poulin (2010) E. Bilgin and D. Poulin, Phys. Rev. B 81, 054106 (2010), URL https://link.aps.org/doi/10.1103/PhysRevB.81.054106.
- Poulin and Hastings (2011) D. Poulin and M. B. Hastings, Physical Review Letters 106, 080403 (2011).
- Jouneghani et al. (2014) F. G. Jouneghani, M. Babazadeh, D. Salami, and H. Movla, arxiv:1409.2048.v1 (2014).
- Morita (1957a) T. Morita, Journal of the Physical Society of Japan 12, 1081 (1957a).
- Morita (1957b) T. Morita, Journal of the Physical Society of Japan 12, 754 (1957b).
- Tanaka et al. (1994) T. Tanaka, K. Hirose, and K. Kurati, Progress of Theoretical Physics Supplement 115, 41 (1994).
- Katori and Suzuki (1994) M. Katori and M. Suzuki, Progress of Theoretical Physics Supplement 115, 83 (1994).
- Morita and Tanaka (1994) T. Morita and T. Tanaka, Progress of Theoretical Physics 92, 1081 (1994).
- Kikuchi (1951) R. Kikuchi, Phys. Rev. 81, 988 (1951).