{'fraction_non_alphanumeric': 0.035529178397747635, 'fraction_numerical': 0.011837727156385975, 'mean_word_length': 4.2404224997904265, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the path integral representation for quantum spin models and its application to the quantum cavity method and to Monte Carlo simulations

The cavity method is a well established technique for solving classical spin models on sparse random graphs (mean-field models with finite connectivity). Laumann et al. [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 1 . Another example is dynamical mean-field theory 2 , 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 N 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 3 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 N ) 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 4 . Cayley tree models have however pathological surface effects, and the theory of finitely-connected mean-field frustrated systems is better defined on random graphs 5,6 , 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 7,8,9 . Finite-connectivity models are technically much more involved than their fully-connected counterparts. The replica method 10 that has been first developed to solve the Sherrington-Kirkpatrick model becomes less practical in this setting 11 , and the alternative cavity method turned out to be more useful 6 .
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 12 . More generally the issue of the nature of the quantum phase transitions at zero temperature 13 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 14,15,16,17,18,19,20 , with a few exceptions that appeared in the last year 21,22,23,24 . In a very interesting contribution 21 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 23 , who developed a similar approach in the framework of the so-called static approximation 14,16,18 . 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) 25,26,27 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 28 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 B, 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 B 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.
From a classical energy function of N Ising spins, E(σ 1 , . . . , σ N ), one can construct an operator E = E(σ z 1 , . . . , σ z N ), 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, Our goal is then to compute the quantum statistical mechanics properties at inverse temperature β, i.e. the partition function Z and the average of observables (operators) O, defined by A well-known way of tackling such problems is to transform them into an extended Ising model by using the Suzuki-Trotter formula 29 , as summarized in the following lines : In the two last lines σ Ns+1 = σ 1 . For a finite value of the number of Suzuki-Trotter "slices" N s , the problem has thus become one of N × N s Ising spins, each of the σ i being promoted to a ring (σ 1 i , . . . , σ Ns i ) with nearest neighbor ferromagnetic interactions along the "discrete imaginary time" α axis (with periodic boundary conditions). The original interactions E 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 σ i = (σ 1 i , . . . , σ Ns i ) is the configuration of the ring of Ising spins at site i and σ = (σ 1 , . . . , σ Ns ) is the full configuration of the N × N s spins. We can thus introduce a probability measure on the N × N s Ising spins, such that the normalization constant Z Ns reduces to the partition function Z in the N s → ∞ limit (in the following we shall sometimes keep implicit the dependence on N s ). To write in a compact way this last equation we have defined the average of independent copies of the classical energy on the various slices, and the ferromagnetic interaction along the imaginary time axis induced by the transverse field (we use σ Ns+1 = σ 1 ). One can easily show that the average value of observables can be obtained in this formalism as 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 O diagonal in the {|σ } basis, i.e. that can be written as O(σ z 1 , . . . , σ z N ): A non-diagonal observable we shall study in the following is the transverse magnetization σ x i (written here for an arbitrary site i), which can be computed as

B. The continuous imaginary time limit
To recover the truly quantum properties of the model one has to perform the limit N s → ∞. The basic degrees of freedom σ i which were the configurations of a ring of Ising spins (σ 1 i , . . . , σ Ns i ) then becomes piecewise constant functions σ i (t) ∈ {−1, 1} of an imaginary time parameter t, the discrete coordinate α ∈ [1, N s ] being mapped to t ∈ [0, β] with the correspondence t = βα/N s . 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 30 easier to formulate than Feynman path-integrals for continuous coordinates 31 , and can be given a rigorous mathematical content 32,33,34,35 . Note that these continuous time trajectories can be easily represented in the memory of a computer, as the trajectory of site i is fully specified by σ i (t = 0) and the times at which the spin flips. Actually numerous continuous time quantum Monte Carlo algorithms do exist, see for instance 36,37,38,39 .
The rest of the paper will crucially rely on the procedure developed in Sec. II B 2 and II B 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 σ i , 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 σ i , or of the trajectory σ i (t) 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 σ i . This is precisely the issue we shall tackle in Sec. II B 2 and II B 3, after having recalled in Sec. II B 1 the well-established path-integral representation of a spin-1/2.

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 (B and h respectively): The diagonalization of the order 2 matrix hσ z + Bσ x easily leads to The path-integral representation of this propagator reads 30,38 Each term of these expressions corresponds to a spin trajectory that changes value at times t 1 < t 2 < . . . ; it is weighted by a factor B raised to the number of such discontinuities, and by exp[h λ 0 σ(t)]; 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 N s → ∞ limit 68 . Alternatively one can notice that Eqs. (12),(13) coincide with (11) at λ = 0 and that they obey the same set of first order linear differential equations, which implies that they coincide for all values of λ.

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 t ∈ [0, λ] to {−1, +1}, conditioned on their initial (σ(t = 0)) and final (σ(t = λ)) values. More explicitly, for instance for σ(t = 0) = σ(t = λ) = σ, the probability of a trajectory with 2n flips at times in the infinitesimal intervals [t j , t j + dt j ], with t 1 < · · · < t 2n , is defined to be 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: 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 σ(0) = σ(λ) = σ is either constant on the whole time interval, or made of a constant part upto time u, followed by a jump to −σ and a second part of the trajectory representative of W (−σ → σ, h, λ − u). Similarly the second one expresses the necessity for a trajectory from σ to −σ to have at least one discontinuity at a given time u, followed by a trajectory accounting for W (−σ → −σ, h, λ − u). These equalities can be proven either from the path-integral representation of (12), (13), or from the explicit expressions of W 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 σ(0) = σ, σ(λ) = σ ′ : draw a random variable u ∈ [0, λ] with density proportional to e σhu W (−σ → −σ, h, λ − u) (see below in eq.(18) for some details on how to perform this step) set σ(t) = σ upto time u call the same procedure to generate a trajectory from −σ to −σ on the remaining interval of length λ − u • if σ = σ ′ -with probability e σhλ /W (σ → σ, h, λ), set σ(t) = σ on the whole time interval, and exit the procedure otherwise, * draw a random variable u ∈ [0, λ] with density proportional to e σhu W (−σ → σ, h, λ − u) * set σ(t) = σ upto time u * call the same procedure to generate a trajectory from −σ to σ on the remaining interval of length λ − u In order to draw u ∈ [0, λ] with a density proportional to e σhu W (−σ → −σ, h, λ − u), we compute its cumulative distribution, FIG. 2: Definition of the effective field trajectory.
A simple way to draw u amounts to draw G uniformly at random on [0, 1], and to invert the above expression to obtain u(G). One can proceed similarly for the generation with the density proportional to e σhu W (−σ → σ, h, λ − u), which involves another cumulative distribution G.

Generating trajectories for a piecewise constant longitudinal field
In the previous subsection we considered the particular case of a constant longitudinal field h. Let us now address the general case of a piecewise constant h = h(t) on the interval of imaginary time [0, β], with the following definitions illustrated in Fig. 2: we shall call p the number of times it changes value between t = 0 and t = β, 0 = t (0) ≤ t (1) ≤ · · · ≤ t (p) ≤ t (p+1) = β the times of these changes, λ (i) = t (i+1) − t (i) the length of these intervals for i ∈ [0, p], and finally h (i) 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, and to generate spin trajectories according to the corresponding weights. The computation of the partition function can be performed by inserting p + 1 representations of the identity in (19), corresponding to the spin values at the imaginary times t (i) where h(t) is discontinuous: with σ p+1 = σ 0 . Given the trajectory h this computation is easily performed, necessiting only the multiplication of the p matrices of order 2 defined in (10). The sampling of the spin trajectory σ(t) on the interval [0, β] is done as follows. The values (σ 0 , . . . , σ p ) of the spin at times (t (0) , . . . , t (p) ) are generated 69 according to the probability Z(σ 0 , . . . , σ p |h)/Z(h). Then for each interval i ∈ [0, p] a spin trajectory from σ i to σ i+1 is generated according to the procedure of Sec. II B 2, the longitudinal field being constantly equal to h (i) on this interval of time. Finally the p + 1 trajectories are concatenated to obtain the full trajectory from t = 0 to t = β. Let us emphasize that the path integral representation of the imaginary time evolution is well known in the literature 30,31,32,33,34,35 ; 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 36,37,38,39 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 T ) or quantum (increasing B) fluctuations destroy this order outside a region delimitated by a critical line in the (B, T ) 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.

A. The quantum cavity method treatment
The quantum Bethe lattice ferromagnet is defined by the Hamiltonian where the first sum runs over the edges of a random l + 1-regular graph of N vertices 40 . This means that the graph of interactions is uniformly chosen among all graphs on N vertices for which each vertex has the same number l + 1 of neighbors. These graphs have the good properties to realize finite size Bethe lattices: the set of vertices at distance 70 smaller than a given cutoff d from an arbitrarily chosen vertex is, with a probability which goes to one when N diverges with d held fixed, a regular tree of connectivity l + 1. 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 N s 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 µ i→j (σ i ) as the probability law of the configuration of the ring σ i when the interaction with its neighbor j has been removed from the graph. If we denote by ∂i the set of vertices neighbors of i, the recurrence equations for these distributions are: where z i→j is a normalization constant, and we introduced for two arbitrary imaginary time dependent quantities a · b = Ns α=1 a α b α /N s . A graphical representation of this equation is given in Fig. 3 in the case l = 2. The site i has then two neighbors, k 1 and k 2 ; the distributions µ k1,2→i encode the effect on k 1,2 of the part of the tree that does not involve i; this is represented by the dashed line in figure 3. In absence of site i the sites k 1 and k 2 are decoupled (as we assume the graph to be a tree there are no paths between k 1 and k 2 that do not go across i) and hence independent. Therefore the distribution of σ i in absence of site j 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 41,42 . 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" µ i→j , 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 σ i in the complete tree is given in terms of the messages sent by its neighbors, with z i a normalization constant, while the joint law for two neighboring sites i and j reads Moreover the free-energy per site f can be expressed in terms of the normalization constants of these laws, 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 µ i→j are all equal to a single law, that we shall denote in the following η(σ), which is seen from (22) to satisfy 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 43 . 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 B = 0. In such a case the weight w forces all spins along the ring to take the same value, hence where h is solution of The classical model thus exhibits a continuous paramagnetic to ferromagnetic transition when the inverse temperature β crosses its critical value, β c = arctanh(1/l). We now come back to the general B case and rewrite explictly the cavity method prediction for the free-energy per site, that is obtained from Eq. (25) by substituting µ i→j = η: and for the probability law on one site and two adjacent sites, Note that the expression (29) is variational, in the sense that its derivative with respect to η vanishes on the solution of equation (26).

B. 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 N s Ising spins, its complete characterization should involve 2 Ns − 1 real numbers. Such a direct representation becomes very soon impossible to handle when N s grows, and in particular in the continuous time limit N s → ∞. We shall however see that an alternative representation allows to bypass this difficulty. First we define a probability law p(σ|h) on the configurations of a ring of spins σ by Z(h) ensuring the normalization of p(σ|h). These definitions allow to rewrite Eq. (26) as Suppose now that one has an estimation of η(σ) given by a representative weighted sample of N traj elements {σ i }, that is such that the weights a i add up to one. A new estimation of η (i.e. a new set of configurations σ ′ i and weights a ′ i ) 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 N traj representants of η, one repeats in an independent way the steps: • generate a configuration σ ′ i according to the law p(σ|h), and set a ′ i = Z(h) Once the N traj new elements have been generated one just has to multiply the weights a ′ i 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 N s → ∞ 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 p(σ|h) corresponds exactly to the procedure we defined in Sec. II B (we show explicitely this correspondence for the normalization Z(h) 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 1,6 . 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 N s → ∞ limit of the formulas (9),(10), which yields 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 j(σ) as the number of discontinuities of σ(t) on the interval [0, β].

C. Continuous-time results
We present now numerical results for l = 2; the behavior is qualitatively the same for any l. We followed the method of resolution presented above, using N traj = 10 4 trajectories. For a fixed temperature, we initialized the population at B = 0; 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 B with a step dB = 10 −2 and for each step we let the population equilibrate for 10 2 iterations and averaged the observables over 10 3 subsequent steps. We checked that the weights a i defined in (33) remain quite uniformly distributed over the population at all investigated temperatures and transverse fields.
In figure 4 we plot the magnetizations m z and m x as functions of the transverse field B for three different temperatures T < T c (B = 0) = 1/arctanh(1/2) = 1.820 . . .. At these low temperatures, the system undergoes a continuous phase transition at a critical value B c (T ) of the transverse field. The transition is characterized by the vanishing of m z and by a jump in the derivative of m x . Unfortunately obtaining reliable values of m z close to B c is a difficult numerical task due to critical fluctuations that cause strong finite size effects in the size of the population N traj . Therefore we located the transition by the jump in the derivative of m x ; we fitted m x by linear laws close to B c from above and below, and determined their intersection.
The resulting phase diagram is reported in fig. 5. We found that for T 0.3 the temperature dependence of all observables is very weak and B c ∼ 2.232, that is a reliable estimate for the T → 0 limit and is in excellent agreement with the value determined in 24 . The scaling of B c for T → 0 is compatible with the essential singularity that is found in the Curie-Weiss model, see Appendix A and in particular Eq. (A8). The very weak (practically unobservable) temperature dependence of the free energy below T = 0.3 make us confident that the entropic term is negligible at these low temperatures and the free energy for T = 0.1, plotted in figure 6 is representative of the ground state energy e GS . The latter has a weak singularity (discontinuity of the second derivative) at B c which is not observable with our numerical precision, but is easily observed by a direct computation of m x = − deGS dB , see figure 4. Overall, our results are in very good quantitative agreement with the ones obtained at T = 0 in 24 by a matrix product state description, except for the value of the exponent β characterizing the vanishing of the magnetization close to B c , m z ∼ (B c − B) β . Even if we do not have very precise data close to B c due to finite population-size effects, our data are compatible with the mean-field exponent β = 0.5 at all investigated temperatures, while in 24 a slightly smaller value of β ∼ 0.41 is reported at T = 0. We will further comment on this discrepancy in section IV C.

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

Resolution at finite Ns
At finite N s , the distribution η(σ) can be encoded with 2 Ns − 1 real numbers (e.g. the probabilities of each of the 2 Ns 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 N s ). This method is extremely precise but computationally very heavy, as the time to solve the equation scales exponentially in N s . We are then limited to N s ≤ 13; the computation for the largest N s = 13 took two weeks on a standard workstation, while we estimated the one for N s = 14 to take at least two months.
Therefore, in order to study the approach to the limit N s → ∞, we resorted to the population method described in Sec. III B. We represented η(σ) by a population of discrete time trajectories σ = (σ 1 , · · · , σ Ns ) and solved Eq. (26) by iteration following the procedure detailed in section III B without taking the limit N s → ∞. Similarly to the continuous time case, we used N traj = 10 5 and performed 10 2 iterations to achieve convergence, after which data have been collected along 10 3 iterations. The computation time is now polynomial in N s so we can go to much larger values (at the price of numerical precision due to finite N traj ). Note that the continuous time computation is much more efficient. In the discrete time case, at low B and large T , the spins σ α are constant along large time intervals and the information encoded in a discrete trajectory is redundant. On the contrary, at low T and large B the method at finite N s is obviously incorrect because of the discretization. One could think to adapt N s in order to find the optimal value for a given T ,B; 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 N s 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 N s ≤ 13 and by the population method for N s > 13. We do not report detailed results for the magnetization as a function of B; 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 T = 1.0 and B = 1.5; the second at very low temperature T = 0.1 and B = 1.8. In figure 7 we plot the longitudinal magnetization m z for different values of N s . As expected, finite N s effects are much stronger at low temperature. Note that for large N s the leading correction is ∝ N s −2 , as in the Curie-Weiss model (see Appendix A). For T = 1 deviations from the leading term are very small and N s ∼ 10 is enough to get a fair estimate of the true m z . On the contrary, for T = 0.1 deviations are very large and N s 256 is needed.
In figure 8 we report the critical temperature as a function of B for different values of N s together with the continuous time result. Note that in the limit of large transverse field, each of the N s time slices becomes independent of the others. The model reduces to N s copies of the classical system (B = 0) with a ferromagnetic coupling rescaled by 1/N s . Then in this limit the critical temperature is given by the classical one divided by N s , T c (N s , B → ∞) = 1/[N s arctanh(1/l)]. 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-1/2 models, the finite N s approximation might be useful in order to understand the behavior at intermediate temperatures, and might also give quantitatively accurate results far from T = 0. 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.

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 43 ). 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 η(σ): 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 14,15,16,18 and has been recently applied to the random k-satisfiability in a transverse field 23 . The present derivation based on the cavity method gives back the same results originally derived in 23 , 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 23 . Equivalently we can rewrite the equation above as where with a slight abuse of notation we used p for the distribution of both the effective field h and its associated magnetization m = tanh h. 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 p(h). From Eq. (36) it follows that, in the limit N s → ∞: where in the last line we defined the function The following explicit expression for w s (m) can be found in 23 (where it was called e −βu0(m) ): with I 1 the modified Bessel function of the first kind. For completeness we provide a proof of this result in App. C. From the equations above one obtains the free energy (29) as a functional of p(m), that has then to be minimized. Remarkably, the resulting expression for the free energy is equivalent to the cavity free energy for a classical system whose variables are the continuous m i ∈ [−1, 1] and whose Gibbs measure is defined by i.e. it is obtained from the quantum measure (5) by replacing the quantum operators by their average and the transverse field term by its average as defined in (38). The analogy with a classical system, or a direct differentiation of the free energy with respect to p(m), allows to write a cavity equation for p(m) which is the analog of (26): Given the structure of w s (m) it turns out that p(m) can be decomposed as a + δ(m − 1) + a − δ(m + 1) + p(m), where p has its support strictly between −1 and 1. We solved Eq. (41) by iteration, using a discretized representation of the regular part p and keeping explicitly the weights a ± . Note that at B = 0 the first term in (39) vanishes. Then a ± = (1 ± tanh(βh))/2 where h is the solution of Eq. (28): the static approximation is obviously exact in the classical case. For a given temperature T < T c (B = 0), we initialised p(m) on the classical solution and increased gradually the transverse field B, at each step iterating Eq. (41) until convergence (typically after ∼ 10 2 iterations with a discretization step dm ∼ 10 −3 ). Eq. (41) can be further simplified in the limit T → 0. As I 1 (x) ∼ e x for large x, the Dirac distributions can be neglected in (39) for any B > 0, and at the leading exponential order w s (m) ∼ e βB √ 1−m 2 . We can thus use the simplified expressions p(m) ∼ e βζ(m) , z l ∼ e −βe l at this order, and hence obtain from (41): where the normalization e l must be determined in such a way that max m ζ(m) = 0. Like Eq.  In summary, the static approximation is a reliable variational tool to study qualitatively the phase diagram of the system. Remarkably, it can be solved down to T = 0 (the solution at T = 0 being easier than for finite T ).

A. Methods and algorithms
The quantum Monte Carlo method is an important tool that is widely used to study quantum statistical physics and quantum phase transitions, especially in the context of lattice models 48 . In this section, we discuss the application of the ideas discussed in this paper to quantum Monte Carlo simulations.

A generic quantum heat bath Monte Carlo scheme
The procedure we discussed in sec. II B to generate the "continuous time" spin configurations can be actually used directly as a continuous time quantum heat-bath algorithm. Once the procedure that generates the new continuous time configuration given the local field trajectory is available, the implementation of the Monte Carlo simulation is rather straightforward: one just randomly picks a site, computes the local field trajectory due to its neighbouring spins, and generates a new imaginary-time trajectory of this spin according to the rules discussed in sec. II B.
The advantage of this method is that one can apply it to any discrete spin models on any kind of lattice. In the case of ferromagnetic non frustrated interactions, there exist of course a number of algorithms that allow to considerably speed up the simulation and this is the subject of the next section. However, for highly disordered and frustrated systems such as spin glasses, for instance, no such generic cluster algorithm exists even at the classical level (although interesting progresses are being made, see 49,50,51,52 ) and most classical simulations 53 rely on the Metropolis or the heat bath algorithm and the parallel tempering technique 54 . In the case of quantum spin glasses, to the best of our knowledge, only the usual discrete time Suzuki-Trotter decomposition has been tried 55,56,57 . Another important case where no loop algorithm is known is the one of multi-spin interactions, where the problems are defined on a factor graph (see section V), which is very common in the context of constraint satisfaction problems such as satisfiability or in coding theory. There is thus a clear need of an efficient generic heat-bath strategy when no cluster or loop algorithm exists, that would result in a substantial gain of simulation time.
In Fig. 10, we show the first application of our continuous time quantum heat-bath algorithm to the case of the ferromagnetic model on a regular random graph of connectivity 3, as studied in the rest of the present paper. The results of our algorithm converge fastly, with respect to the number of Monte Carlo sweeps, to the asymptotic ones for very large samples (i.e. N = 10 5 spins). In the paramagnetic phase, the convergence is even faster than for the loop algorithm discussed in the next section. The application of these ideas and methods to more complex models is an interesting direction of study.
Since we are dealing with a ferromagnetic, non frustrated model, we expect the usual cluster and loop algorithms 36,39,58 to allow for a significant speed up of the simulation with respect to the heat-bath procedure, even on random graphs. We have thus implemented the algorithms devised by Rieger and Kawashima in 38 , which is an adaptation of the classical Swendsen-Wang 59 algorithm to quantum systems in continuous time. The results are compared with the heat-bath method in Fig.10. Indeed, the loop algorithm performs very well and better than the heat-bath method close and below the critical point, and we thus have used this method to compare the results of the cavity approach with finite size instances.

B. Cavity method versus simulation
How does the cavity predictions compare with numerical simulations? To answer this question, we have performed quantum Monte Carlo simulations on random regular graphs of N spins and connectivity 3. For large N we expect the results to be self-averaging and thus we always consider a single instance, and do not perform averages of many realizations. The results for two different temperature T = 1 and T = 0.1 are shown in Fig.11 where we plot the longitudinal and transverse magnetization as function of B for different sizes from N = 64 to N = 2048. The agreement with the asymptotic cavity result is perfect (apart from finite size effects in m z for B ≈ B c , see next paragraph): this demonstrates the correctness of the approach we have developed in the present paper.

C. Finite size scaling and critical exponent
The transition between the ferromagnetic and the paramagnetic phase is of second order; below the threshold B c (T ) the longitudinal magnetization grows with an exponent β, i.e. m z ∝ (B c − B) β , in the thermodynamic limit. We have used finite size scaling techniques 60,61 to analyze our data in the neighborhood of the transition and to check the mean-field value of the exponent β = 1/2. Let us first briefly recall the basic idea of the finite size scaling method and the way in which it has to be amended for mean-field models.
In a generic infinite size d-dimensional model, in the vicinity of a second order phase transition driven by a parameter denoted B, the correlation length diverges as ξ ∝ |B − B c | −ν . For a system of finite extent L the observables depends on the size through a scaling function of the ratio L/ξ. This has to be corrected for dimensions d larger than the upper critical dimension d u of the considered universality class, or when the model lacks any underlying finite dimensional structure. In that case the scaling function is found 62,63,64,65,66 to depend on the size N (total number of degrees of freedom, equivalent to L d of a d-dimensional model) through N 1/(duν) (B − B c ), where ν takes its mean-field value in the universality class under investigation. More explicitly the scaling forms of the longitudinal magnetization and of the Binder cumulant g = 1 We present in Fig. 12 the results of such an analysis for our Monte-Carlo data obtained at temperature T = 1. At this positive temperature the critical behavior of a quantum d-dimensional system is equivalent to a classical d-dimensional Ising model 13,29 , which implies d u = 4 and ν = 1/2. Using the cavity method predictions B c (T = 1) = 1.98 and β = 1/2 we found a very good data collapse (see Fig. 12), which confirms the validity of these values of β and B c . To close this section, let us discuss the value of the critical exponent β at T = 0. We believe that because of the mean-field nature of the Bethe lattice model, β keeps the same value 1/2 at positive and zero temperatures. This is indeed what happens for the Curie-Weiss model, see Appendix A. A simple argument is the following. For ferromagnetic models the Suzuki-Trotter formalism and numerical simulations suggest that the critical behaviour of the quantum d-dimensional model at zero temperature corresponds to the one of the classical model in d + 1 dimensions 13,29,38 . Hence if a model behaves in a mean-field way at positive temperature it should also do so at zero temperature, having formally gained one more spatial dimension. To give further credit to our theses we performed quantum Monte-Carlo simulations at very small temperatures, T = 0.01. The plot of Fig. 13 shows again a very good collapse with the value β = 1/2, using the zero temperature critical value of the transverse field extracted from the cavity computation, B c = 2.232. This collapse has been obtained with d u = 3; indeed at such a low temperature we are well inside the low temperature regime where the transition is truly quantum, hence the shift of the upper critical dimension to account for the imaginary time supplementary dimension.
A different value of the zero temperature exponent, β = 0.41, has been obtained in 24 using the matrix product state ansatz for the description of the groundstate. Though we cannot rule out the possibility that a crossover occurs at temperatures even lower than T = 0.01, we believe that β = 1/2 down to T = 0, in accordance with the mean-field nature of the Bethe lattice, and that the different value reported in 24 might be due to the truncation of the matrix product state ansatz.

V. THE GENERIC REPLICA SYMMETRIC QUANTUM CAVITY METHOD
In this final Section we give a more generic description of the replica symmetric quantum cavity method. As the main ideas should have already been conveyed by the example of the ferromagnet on the random regular graphs we shall mainly emphasize the differences and complications that arise in more general models.
Let us consider the class of Ising spin models whose classical energy E(σ) can be decomposed into the sum of M B c = 2.232 , T = 0.01 interaction terms a = 1, . . . , M , each of them depending on a finite number of spins, In this expression ∂a is the subset of spin indices the a-th interaction effectively depends on, σ ∂a = (σ i : i ∈ ∂a) is a shorthand for the configuration of those spins, and J a denotes coupling constants that might appear in the definition of the a-th interaction term. This decomposition is conveniently represented as a factor graph 41 , i.e. a graph with two kind of vertices (see Fig. 14 while the transverse weights w(σ i ) are local in the spin indices i. As we did on the ferromagnet example we first state the solution of this extended model on a tree. The difference is that we have now two kind of "messages", from variables to interactions and vice versa. Let us denote µ a→i (σ i ) the probability law of σ i in the model corresponding to the factor graph where all interactions but a have been removed in the neighborhood of i, and similarly µ i→a (σ i ) for the factor graph with only a removed from ∂i. These quantities obey the following recursion equations, on all edges of the factor graph, the various constants z being normalization factors. The reader will easily verify that these equations reduce to Eq. (22) in the ferromagnetic case with ε(σ i , σ j ) = −σ i σ j ; in this case we could eliminate one type of message (from interactions to variables), as the interactions were only pairwise.
We shall now consider ensembles of random factor graphs, denoting E[·] the expectation over the distribution of the graphs and coupling constants. We assume all interaction nodes to involve a fixed number k of variables (the ferromagnet corresponded to k = 2), while the degree distribution of the variables is specified by a probability law q d over the positive integers (all random graphs verifying these constraints are equiprobable in the ensemble; the lack FIG. 14: An example of a factor graph of a finite dimensional a priori structure is the origin of the mean-field character of this family of models). We shall denote γk = d dq d the average variable degree. Moreover the M = γN coupling constants J a are drawn in an identical, independent way for each of the interactions. Suppose the recursion equations (46) are solved on a factor graph sampled at random from the ensemble under consideration, and that an edge from a variable to an interaction, call it i → a, is chosen uniformly at random. The probability law µ i→a it bears is itself a random variable, denoted η in the following, due to the random choices of the graph and of its coupling constants. This random variable can be related to its neighboring equivalents by the local relations (46). Suppose that i has d adjacent interactions apart from a (see Fig. 15 for an illustration). Then one has, from (46), where the η a,i are the d(k − 1) probability laws that determine η. The hypothesis of the replica symmetric cavity method is to assume that the above form is correct in spite of the graph not being globally a tree, and that the η a,i are independent identically distributed copies of the random variable η, which is written in formula as the symbol d = denoting identity in distribution of random variables, and the function f being an abbreviation of the right-hand-side of (47). Note that the connectivity random variable d is not distributed according to q d but rather q d = (d + 1)q d+1 /γk: as one uniformly choose a random edge of the graph, and not a random site, this sampling favours larger connectivity variables. For the ferromagnet on the random regular graph one had q d = δ d,l+1 while q d = δ d,l .
The physical observables of the system can be computed from the solution of this distributional equation. The thermodynamic limit of the free-energy per site is found to be where the expectations of the left-hand-side are over the choice of the random factor graphs, and for the right-handside over independent copies of the random variable η and the coupling constants J, while in the first term d is drawn from the law q d . This is a generalization of the expression (29) found for the ferromagnet. Similarly the average marginal of the law (5) for an arbitrary site reads with d drawn from q d , while the average marginal law of the k variables in an arbitrary interaction is the factors 1/z in these last two equations ensuring their normalization (cf. Eq. (30) for the regular ferromagnet).
Let us now discuss a method of resolution of the distributional equation (48), consisting in encoding the distribution of η by a sample (or population) of a large number N of representatives η i . From an arbitrarily initialized population, one applies iteration steps according to Eq. (48): one draws an integer d from q d , d coupling constants J 1 , . . . , J d and d(k − 1) indices in [1, N ]. The right-hand-side of Eq. (48) is then computed from the corresponding d(k − 1) representants of η randomly chosen in the population, and the resulting η is used to replace one discarded element of the population. The iteration of these steps brings the population close to a fixed point of the distributional equation; then physical observables like the average free-energy (see Eq. (49)) can be computed, evaluating the expectations over the random variable η as a sampling from the approximate representation provided by the finite population. Compared with the classical replica symmetric cavity method the difficulty is that each of the η i in this population is itself a probability distribution over configurations of spin rings, or in the continuous limit over spin trajectories. It is however possible to apply the trick explained in the simpler case of the ferromagnet, and encode each of the η i as a sample of N traj spin trajectories. Let us introduce some further notations in order to give an explicit form of the updating procedure of the population of η i 's.
The dependence of ε(·, J) on one of the Ising spins can always be parameterized through two functions u and v as ε(σ, σ 1 , . . . , σ k−1 , J) = − σ u(σ 1 , . . . , σ k−1 , J) + v(σ 1 , . . . , σ k−1 , J) , u playing the role of an effective magnetic field acting on σ. For N s -fold replicated variables we define u(σ 1 , . . . , σ k−1 , J) = (u(σ 1 1 , . . . , σ 1 k−1 , J), . . . , u(σ Ns 1 , . . . , σ Ns k−1 , J) , We can thus apply the sampling procedure explained in Sec. III B, the only change being the different definition of the effective field trajectory h, and the inclusion of a contribution arising from the function v in the weights of the generated spin trajectories. Let us emphasize that at variance with the classical cavity method, here we have to deal with a population of population of spin trajectories (of total size N × N traj ) already at the replica symmetric level.

VI. CONCLUSIONS
The explicit procedure for the construction of continuous time spin trajectories presented in Sec. II B allowed us to make progresses both on the analytical side, with an improvement over the discretized time cavity method of 21 and the static approximation of 23 , and on the numerical simulations side, with a generic quantum Monte Carlo procedure for spin-1/2 models in transverse field. The ferromagnetic model we studied in this article is only the simplest of a large family that can be tackled, at the price of an heavier computational cost, with the same methods. We believe that its study was in any case worthwhile from a methodological point of view: its simplicity permitted a complete resolution of the quantum cavity equations, which (i) showed a perfect agreement with Monte Carlo simulations, hence giving credit to the conjecture that the replica symmetric cavity method leads to exact results in the thermodynamic limit for sparse mean-field ferromagnetic systems, as was proved for classical models in 43 ; (ii) allowed to test quantitatively some approximate treatments (finite N s and static approximation).
The more interesting cases we plan to address in the future by mean of the cavity method will involve fluctuating connectivities, in order to study the role of these local fluctuations on the critical behaviour of the models, glassy phases at low temperatures (a particularly motivating case will be the regular multispin ferromagnet studied at the classical level in 47 ), and models related to quantum computing issues, as the random k-satisfiability model in a transverse field 23 . Note that already the replica symmetric treatment of these models will involve a population of populations of spin trajectories, as explained in Sec. V, which will make an exact treatment of the one step replica symmetry breaking version of the cavity method extremely challenging. We hope however that the better control of the approximative treatments we gained on the simple ferromagnet will help to devise appropriate approximate strategies to handle these cases.
At the same time, we proposed a heat bath Monte Carlo strategy that should be applicable to general spin-1/2 models. The performance of this algorithm should be tested on frustrated models for which cluster algorithms are not easy to implement, and we expect that in such cases the heat bath method could give interesting results.
The strategy developed above become inefficient at very low temperatures, because the spins jump many times along a trajectory and the amount of information needed to encode η(σ) by a population of trajectories becomes very large. Therefore an important directions of research would be to look for possible simplifications of the formalism in the limit β → ∞. Moreover, one could use the recursion equations (46) on a given sample, thus constructing a quantum belief propagation algorithm 67 . Finally, it is worth to note that the procedure described in section II B, on which our results are based, is in principle not restricted to spin-1/2 models. It should be possible to generalize it to any model built on discrete degrees of freedom 30 , for instance hard-core bosons or higher spin models.
where the scaling of the coupling constant is chosen appropriately to make the thermodynamic limit well-defined. At variance with the Bethe lattice model here each spin interacts with all others. Applying the Suzuki-Trotter decomposition described in Sec. II A leads for a finite number N s of Suzuki-Trotter slices to the following expression of the partition function: We then perform N s Hubbard-Stratanovitch transformations to disentangle the quadratic terms and obtain Z = βJN 2πN s This can be further simplified if the N s → ∞ limit is performed afterwards and yields for the free-energy per site Using the variational character of this expression the longitudinal and transverse magnetizations (i.e. m z and m x ) can be obtained by taking the explicit derivatives with respect to h and B respectively, where m is taken as the solution of the saddle-point equation. The latter is easily found to imply that m = m z at the saddle-point. In the following we take J = 1 to simplify the notations. In absence of the transverse field (B = 0) one recovers the classical Curie-Weiss model, with m = m cl (β, h) solution of the traditional equation m = tanh(β(m + h)). The ferromagnetic transition is signaled by the appearance of a non-trivial solution at h = 0, which is possible for small enough temperatures, i.e. for β > β c (B = 0) = 1. Consider or equivalently, The bottom panel of Fig. 16 shows the evolution, as a function of the transverse field, of the longitudinal and transverse magnetizations at a temperature smaller than the classical critical temperature T c (B = 0) = 1. The transverse magnetization is continuous but its derivative has a finite jump at the transition. One can indeed show from Eq. (A7) that m x = B in the ferromagnetic phase, while m x = tanh(βB) for the paramagnetic one. Another thermodynamic quantity easily computed for the Curie-Weiss model is the ground-state energy per spin, e gs (B), obtained from (A5) in the limit β → ∞. One finds which shows that e gs (B) and its first derivative are continuous at the transition, while its second derivative has a finite jump.
Let us finally argue about the scaling of the finite N s corrections, reconsidering the step we took between Eqs. (A4) and (A5). At the next-to-leading order one would have obtained One can then check explicitely that the eigenvalues of the matrix in the exponential are not modified at the order N s −1 , hence the finite N s correction on the observables should be of order N s −2 , as we observed in the study of the Bethe lattice ferromagnet.
In the last step we have taken the N s → ∞ limit to obtain the expression of Eq. (19). In this appendix we derive the expression (39) for the average transverse weight in the static approximation. Rewriting the sum over σ of (38) as a sum over paths in the continuous time limit, one obtains dx 2n δ (2x 1 − 2x 2 + · · · − 2x n + 1 − mσ) , where in the second line we have isolated the contribution of the constant trajectories, and changed variables from t i to x i = t i /β. Let us concentrate on the term corresponding to a given value of n > 0. We perform a further change of variables, setting y i = x i − x i−1 (with x 0 = 0), the intervals between the reduced time of flips in the trajectory. The Jacobian of this change of variables being 1, the integral over x 1 , . . . , x 2n can be rewritten as 1 0 dy 1 . . . dy 2n 1 1(y 1 + y 2 + · · · + y 2n ≤ 1) δ(2y 2 + 2y 4 + · · · + 2y 2n − (1 − mσ)) (C2) Here we have used 1 1(·) as the indicator function of an event, S o = y 1 + y 3 + · · · + y 2n−1 and S e = y 2 + · · · + y 2n the sum of the odd/even y i 's, and ρ n (S) as the density of the distribution of the sum of n independent random variables uniformly distributed on [0, 1]. It is easy to prove by recurrence that for S ∈ [0, 1] one has ρ n (S) = S n−1 /(n − 1)!. Collecting these facts together leads to which is indeed equal to (39) thanks to the series expansion of the Bessel function.