PDE/statistical mechanics duality: relation between Guerra's interpolated $p$-spin ferromagnets and the Burgers hierarchy

We examine the duality relating the equilibrium dynamics of the mean-field $p$-spin ferromagnets at finite size in the Guerra's interpolation scheme and the Burgers hierarchy. In particular, we prove that - for fixed $p$ - the expectation value of the order parameter on the first side w.r.t. the generalized partition function satisfies the $p-1$-th element in the aforementioned class of nonlinear equations. In the light of this duality, we interpret the phase transitions in the thermodynamic limit of the statistical mechanics model with the development of shock waves in the PDE side. We also obtain the solutions for the $p$-spin ferromagnets at fixed $N$, allowing us to easily generate specific solutions of the corresponding equation in the Burgers hierarchy. Finally, we obtain an effective description of the finite $N$ equilibrium dynamics of the $p=2$ model with some standard tools in PDE side.


Introduction
During the last decades, Statistical Mechanics of disordered systems has acquired a prominent role in describing complex phenomena and emerging properties of system with highly non-trivial dynamics. In particular, one of the major success is the possibility to analyze the relaxation and the equilibrium dynamics of spin-glass models [62,75], i.e. systems of simple degrees of freedom (the spins, adopting the physics jargon) whose interaction are frustrated (or, in other words, spins in the system can compete with each others) and known only through their probability distribution. Random and frustrated interactions are responsible for a rich phenomenology (e.g. the existence of multiple timescales for the relaxation toward the equilibrium), which is due to the complex structure of the free-energy landscape (see for example [62]). The analysis of spin-glass equilibrium dynamics is commonly accomplished by computing the so-called quenched free-energy 1 in the thermodynamic limit (whose analysis is motivated by its self-averaging property [53,54]) in terms of few order parameters (in particular, the so-called Edward-Anderson overlap [40]) and tracing the phase diagram in the parameter space. Remarkably, Statistical Mechanics of disordered systems has gained a renewed interest with the advent of Deep Learning [60], since it offers the ideal tool for investigating information processing in neural networks (see for example [2,3,17,19,23,35,38,43,56,58,71]). In particular, Statistical Mechanics commonly deals with the Machine Retrieval regime, i.e. the long-term relaxation of neural networks which are fed with some inputs to be classified according to some previously learnt features. Focusing on the methods, the solution of such complex systems can be found by means of the straightforward (but unfortunately non-rigorous) replica trick approach. At the same time, rigorous and mathematically transparent approaches can be developed, see for 1 With quenched free-energy, we mean the expectation value w.r.t. to the intrinsic disorder (the couplings) is performed after taking the logarithm, i.e. the intensive free energy is f Q N (β) = −(βN ) −1 E log Z N (β), to be compared with its annealed counterpart, i.e. f A N (β) = −(βN ) −1 log EZ N (β). Thus, the disorder in the two versions of the free energy is treated in a substantial different way. The relevant quantity for spin-glass models is the former, since it deals with the equilibrium dynamics of the system at fixed couplings realization. example [13,47,49,51,67,68,70,[80][81][82]. For our concerns, the most important one is the Guerra's interpolation scheme, which makes use of sum rules for the computation of relevant quantities [3-5, 14, 24, 50] or taking benefit of mechanical analogy or relation with PDEs [1,9,[19][20][21][22]. In the latter approach, the quenched free-energy can be interpolated in an effective mean field scheme, in which the interpolation parameters can be interpreted as coordinates of a fictious spacetime (t, x). In this case, it is possible to use the entire set of PDE technology in order to solve the thermodynamics of the model, also for finite N . Even if such interpolation methods were initially developed in the context of spin-glass, they are quite general, and by this they can be easily brought to similar (and possibily simpler) models, such as p-spin ferromagnets. The latters are spin systems in which the spins interact with coupling of order p and cooperatively (i.e. they are not spin-glasses). In this paper, we study the relation of Guerra interpolated partition functions of p-spin ferromagnets and prove that the expectation values of the order parameter (the global magnetization) w.r.t. to the associated Boltzmann-Gibbs measure is in 1-1 correspondence with the elements of the Burgers hierarchy [59,66,74,83]. This duality can be explored in both sense: for instance, we can generate specific solutions of the Burgers hierarchy by exploiting the finite size solution of the p-spin systems; on the other hand, we can achieve informations about the thermodynamics of the ferromagnets from the PDE side.
Despite being interesting by itself, this duality could be a promising tool for investigating the principles of information processing in biological phenomena. Indeed, as firstly suggested in [84], Statistical Mechanics turned out to be an effective tool in the description of universal phenomena in biological systems, see also [7,12,16,29,57,63,72,78]. For instance, kinetics of biochemical reactions can be framed in a purely statistical mechanics picture in terms of the Curie-Weiss ferromagnetic model, see [8,11]. However, since in standard statistical mechanics one considers the thermodynamic limit N → ∞, this scenario is a good approximation only for systems with large size. As opposite to this picture, recently the research in biochemical information processing has focused on phenomena involving system with small size [73], for which the statistical mechanics scenario loses its predictive power. Thus, for such practical application one has to consider the thermodynamics at finite size N . On the other side, it has been established [31] that the binding of spins in systems such as long chain molecules (such as proteins) can takes place with multisite interactions (i.e. p > 2) beyond the pairwise scenario, see also [39]. In this picture, our duality can thus provide rigorous mathematical tools to investigate these interesting biological situations.
The paper is organized as follows. In Section 2, we provide the basic tools of statistical mechanics for p-spin ferromagnetic models, also discussing the thermodynamic solution of the free-energy (both with physical and mathematically rigorous methods, i.e. Guerra's interpolation scheme). In Section 3, we provide the fundamental notions about the Burgers hierarchy and the reduction to linear PDEs through Cole-Hopf transforms. In Section 4, we prove the duality between the two sides and some results (such as the interpretation of phase transition in the ferromagnets as the development of shock waves in the PDE side and the search of solutions of the Burgers hierarchy by means of finite size solutions). We also give a description of p = 2 ferromagnetic model with the tools offered by the Burgers equation.

A cursory look at mean-field p-spin ferromagnetic models
In this Section, we will introduce the fundamental notions about the p-spin ferromagnetic models. In particular, we will discuss the solution of the model, both with purely physics arguments and Guerra's interpolation schemes. We will also discuss the existence of the thermodynamic limit for the free-energy of the models. Let us start by introducing the following Definition 1. Let σ ∈ Ω N = {−1, +1} N be a general configuration of the system at finite size N , and let J > 0 be the interaction strength. The Hamilton function of the p-spin ferromagnetic model is defined as (2.1) Remark 1. The normalization factor 1/N p−1 is inserted in order to ensure the linear extensivity of the Hamilton function, meaning that the energy of the system scales linearly in its volume: where N (σ) (the energy per site of the model associated to the configuration σ) is finite in the infinite-size limit N → ∞.
Remark 2. Traditionally, the Hamilton function of a p-spin ferromagnetic model is given by However, the difference between the two formulations is negligible in the thermodynamic limit. Indeed, by noticing that , holding in N → ∞ limit (see for example [33]), it is easy to see that thus the thermodynamics of the systems within the two frameworks are equivalent (and they differ only up to a trivial rescaling of the temperature by a factor p!).
Definition 2. Let β ∈ R + the level of thermal noise (i.e. the inverse temperature β = T −1 ). Then, the partition function of the ferromagnetic p-spin model is defined as

2)
where σ ≡ σ∈Σ N is the sum over all possible configurations of the system. The Boltzmann factor corresponding to the partition function (2.2) is defined Definition 3. The global magnetization of the ferromagnetic p-spin model is defined as Remark 3. Since the model is ferromagnetic, this is the only order parameters we need to fully describe the equilibrium of the model.

Remark 4. The expression of the Hamilton function in terms of the global magnetization is
Remark 5. Notice that, as usual, the ferromagnetic strength J has the only effect of rescaling the thermal noise in the system. Therefore, without loss of generality, we can set J = 1. In this way, we will denote H N (σ) ≡ H N (σ|J = 1) and Z N (β) ≡ Z N (β, J = 1).
Definition 4. Given a function F (σ) of the spins in the system, its expectation value is defined as Remark 6. We stress that, for odd p, the Hamilton function is not invariant under gauge transformations σ → −σ. This means that, as opposite to the even p cases, there are no gauge-equivalent solutions.
Definition 5. The intensive statistical pressure A N (β) of the system is defined as Remark 7. The intensive statistical pressure is related to the usual (intensive) Helmholtz free energy f N (β) as A N (β) = −βf N (β). Let us denote with P (σ) the Boltzmann-Gibbs distribution of the system, which is of course defined as Given the energy per site at fixed system configuration we can write down the equality which, a part for a factor −β, is the relation between the (intensive) Helmholtz free energy and the expectation values of the relevant thermodynamic observables for the system. Indeed, the first contribution is nothing but the entropy per site, since, according to the definition (2.6), we have Clearly, S N [P ] = − σ P (σ) log P (σ) is the Shannon entropy associated to the probability distribution P (σ). Thus, Eq. (2.8) exactly paints the relation between the intensive pressure A N (β) and the intensive Helmoltz free energy f N (β).
For our concerns, it is relevant the statistical pressure in the thermodynamic limit The solution of the p-spin ferromagnetic models can be straightforwardly obtained by purely statistical tools. The expression of the expectation value of the energy per site follows directly from (2.5), and reads Regarding the entropy per site, due to the mean-field nature of the model we make the assumption that the equilibrium probability distribution to observe the system in a given configuration can be factorized as product of probabilities of independent sites P (σ) = N i=1 P (σ i ), and the spin orientation is driven by the global magnetization. Then, in the thermodynamic limit we can write wherem is the thermodynamic value of the magnetization: Remark 8. When working at finite size N , the computation of the expectation value of the energy per site requires evaluating correlation functions of the form ω[σ i1 σ i2 . . . σ ip ], which is in general a non-trivial task. However, due to the mean-field nature of the ferromagnetic model, simplifications occur in the thermodynamic limit. Indeed, as remarked above, in this limit the probability distribution in the configuration space factorizes, i.e. P (σ) = i P (σ i ), so that the evaluation of the aforementioned correlation functions becomes trivial. Likewise, it is possible to achieve the same result by assuming the self-averaging property for the global magnetization. In simple words, we require the fluctuations of the order parameter w.r.t. its equilibrium valuem to vanish as N → ∞. This is translated in mathematical terms by requiring that the probability distribution of the global magnetization converges to a Dirac-delta distribution which is centered around the equilbrium value, i.e. lim β almost everywhere, with m being the global magnetization in the thermodynamic limit. The selfaveraging assumption is a reasonable hypothesis for ferromagnetic systems, see for example [18], as the number of pure states does not depend on the system size. As a consequence, the expectation value of a general function F of the global magnetization can be simply evaluated: The previous remark implies that the expectation value of the energy per site in the thermodynamic limit can be evaluated as Regarding the entropy contribution, the mean-field assumption (2.10) allows us to discard correlations between the spins and then reduces the problem to one-body computations. Then, with this expression of the probability distribution we can write down the entropy using the Shannon prescription: which is valid in the thermodynamic limit. Then, putting everything together we get (2.11) By imposing the extremality condition for the intensive pressure ∂mA(β) = 0, we have βpm p−1 − arctanh(m) = 0 ⇒m = tanh(βpm p−1 ), (2.12) which is exactly the self-consistency equations for the p-spin ferromagnetic model [18].

Existence of the thermodynamic limit
In the previous Section, it was tacitly supposed the existence of the thermodynamic limit for the p-spin ferromagnetic models as described by the partition function (2.2). For the sake of simplicity, we rigorously prove the existence of the thermodynamic limit for even p (however, this results can be carried out also for models with Hamilton functions which are polynomial in the global magnetization m N (σ), therefore including also the odd p case, see for example [26]). The key idea is to separate the ferromagnetic models with N interacting spins in two distinct non-interacting subsystems respectively with N 1 and N 2 spins, such that N = N 1 + N 2 . Then, we introduce an interpolating partition function with the following Definition 6. Let σ i = ±1 for i = 1, . . . , N be the binary spins of the model and t ∈ [0, 1] an interpolating parameter. Let us build two separate, non-interacting subsystems by collecting respectively the spins σ i with i = 1, . . . , N 1 and σ i for i = N 1 + 1, . . . N 1 + N 2 , with corresponding order parameters Then, the interpolating partition function is defined as 14) The Boltzmann factor associated to the partition function (2.14) is The statistical pressure of the model is defined as Remark 9. We stress that, within the Guerra's interpolating framework, we do not divide the system of size N in two subsystems with resp. N 1 and N 2 spins. Indeed, in order for this decomposition to hold in the previous interpretation, one has to discard mutual interaction (the surface term) between the two subsystems (at least in the thermodynamic limit), which is only possible in finite-dimensional (i.e. non-fully connected) models. In fully connected spin systems, the surface and volume terms are of the same order, thus the former cannot be neglected. In order to avoid this, the Guerra's generalized partition function (2.21) works by interpolating between the original model and two independent systems with sizes N 1 and N 2 , such that N 1 + N 2 = N . In this way, surface terms are not present.

Remark 10.
The global magnetization of the composite system is a convex linear combination of the order parameters of each component, i.e.
where ρ = N 1 /N . Further, the interpolating free energy satisfies the boundary conditions where of course A N (β) is the same as (2.7).
Definition 7. Given a function F (σ) of the spins in the system, its expectation value for the interpolating system (2.14) is defined as Having introduced the central quantities, we are now in position to state the following Theorem 1. The thermodynamic limit of the model (2.2) for even p exists and it is given by where inf stands for the infimum.
Proof. First of all, we show that, for fixed β ∈ R + , the intensive pressure A N (β) is bounded from below for each N ∈ N. To see this, we write down its explicit expression: Thus, we can bound the intensive pressure as is finite. Next, we consider the intensive pressure associated to the interpolating partition function (2.14). Its t-derivative is Provided that p is even, the mapping x → x p is convex, meaning that ∀ρ ∈ [0, 1] and for all values of m 1 (σ) and m 2 (σ). This directly implies that thus the interpolated intensive pressure is a decreasing function w.r.t. the interpolating parameter t. As a straightforward consequence, we have In other words, the sequence {N A N (β)} N is sub-additive, and by virtue of the Fekete's lemma, we easily get lim This proves our statement.

Solution via Guerra's interpolating scheme
The solution of the model can be equivalently carried out via the interpolating techniques [] originally developed in the context of spin-glass as alternative to the replica trick route. The key idea of the method is to introduce a generalized partition function interpolating between the original model and an effective field scenario (i.e. this limit is a 1-body model). Then, we introduce Definition 8. Let t ∈ [0, 1] an interpolating parameter. Then, the Guerra's generalized partition function is defined as where ψ is a constant to be set a posteriori. The Boltzmann factor associated to this partition function is The associated generalized intensive statistical pressure is therefore The key idea of the interpolating program is resumed in the Proposition 1. The intensive free energy of the model in the thermodynamic limit is obtained by means of the sum rule Proof. The proof is a straightforward application of the fundamental theorem of integral calculus.
Definition 9. Given a function F (σ) of the spins in the system, its expectation value w.r.t. to the partition function (2.21) is defined as Let us call againm the thermodynamic limit of the value of the global magnetization at the equilibrium, i.e.m = lim Remark 11. The equilibrium valuem of the global magnetization in the thermodynamic limit can can depend on the interpolating parameter t. However, it is possible to show that in the thermodynamic limit, the free parameter ψ can be conveniently chosen in order form to be independent on t almost everywhere (see next Proposition). To do this, we shall again assume the self-averaging property of the order parameter β almost everywhere and for all t ∈ [0, 1]. This again implies that the variance of the global magnetization vanishes in thermodynamic limit, i.e.
Further, we make the assumption that the variance of the order parameter scales as N −1 , i.e. 27) or, in other words, fluctuations around the thermodynamic value of the order parameter scales as 1/ √ N . In the physics jargon, this is equivalent to require that the magnetic susceptibility diverges only at the critical point, which is reasonable for ferromagnetic systems. The assumption (2.27) will be very useful in Proposition 2, since we can introduce a function ∆(σ) accounting for fluctuations around the thermodynamic value of the global magnetization: for sufficiently large N and β almost everywhere. The function ∆(σ) has zero mean (ω t [∆(σ)] = 0, which trivially follows from taking the expectation value of both the terms in the expansion) and finite variance ω t [∆(σ) 2 ], even in the thermodynamic limit. In other words, by restricting ourselves to the pure state with positive magnetization (without loss of generality), we can express the global magnetization of a general configuration σ in terms of fluctuations (the second term on the right side in (2.28)) around the expectation value ω t [m N ]. The two sides of (2.28) are compatible, since with the expansion (2.28) it is easy to show that for sufficiently large N , thus satisfying the assumption (2.27).
Now, we have all of the ingredients needed to prove the following Proposition 2. It is possible to suitably choose the parameter ψ ∈ R such that the expectation value in the thermodynamic limit of the global magnetization is independent on the interpolating parameter almost everywhere, i.e. dm dt = 0 a.e. (2.29) Proof. First of all, we assume the decomposition (2.28), which we report here for the sake of completeness: The t-derivative of the expectation value of the global magnetization is clearly (2.31) Adopting the expression (2.30), we can now compute this quantity term by term at the non-trivial order in the 1/N expansion. In particular We are interested in considering only the contributions up to the 1/N order (the subleading terms will vanish in the thermodynamic limit). Then where R 1 (m N (σ)) accounts for the rest of the expansion, i.e.
The leading contribution in R 1 is of order N −3/2 , thus the whole quantity will not contribute to (2.31) in the limit N → ∞. Since the fluctuations have zero mean, the second term in (2.32) identically vanishes, leaving us only with In a similar fashion, it is easy to prove that where also in this case R 2 (m N (σ)) accounts for the subleading corrections scaling at least as N −3/2 for large N . Finally, it is clear that where we defined Q(m N (σ)) = R 1 (m N (σ)) − R 2 (m N (σ)) whose leading contribution scales itself as N −3/2 in the large N limit (thus, N Q(m N (σ) scales as N −1/2 ), and therefore it is negligible in the N → ∞ limit). The r.h.s. of the last line in (2.37) is well-defined in the thermodynamic limit, thus -callingm = lim N →∞ ω t [m N [σ]], we have (recall that the variable ∆ has finite variance in the N → ∞ limit) Here, we used the fact that N Q → 0 in the limit N → ∞, since the leading contribution of the quantity Q is of order N −3/2 . This means that the sequence { dωt[m N (σ)] dt } N converges almost everywhere to the r.h.s. of the previous equation, so by virtue of Egorov's theorem [41], it is almost uniformly convergent. As a consequence, the relation holds almost everywhere, which means Then, it is simple to note that the requirement dm dt = 0 can be consistently fulfilled almost everywhere by choosing ψ = βpm p−1 , which proves our assertion.
Proposition 3. The expectation value of the p-th power of the global magnetization can be written in terms of the centered momenta with degree lower or equal to p.
Proof. The proof is a straightforward application of binomial theorem. Indeed Remark 12. We stress that we can also extract the lower two terms from the sum, in order to get or in other words With all of these ingredients in our hand, we can prove the following Theorem 2. The thermodynamic limit of the intensive statistical pressure for the model (2.2) is given by Proof. First of all, we compute the derivative of the statistical pressure. In this case, we have We notice that, recalling our choice ψ = βpm p−1 , we have so that we can apply the Rem. 12. Indeed, using the relation (2.40), we easily get We stress that, by assuming the self-averaging property (2.26) of the order parameter, the sum in r.h.s. clearly vanishes in the thermodynamic limit (since centered momenta will disappear in the N → ∞ limit), therefore leaving us only with which is independent on t (due to Prop. 2). Thus, its t-integration is trivial. On the other side, the t = 0 is easy to handle with, since it is a 1-body computation. Indeed which leads to A N (β, t = 0) = 1 N log 2 N cosh ψ = log 2 + log cosh(βpm p−1 ), (2.43) where we recalled our choice ψ = βpm p−1 . Now, using the sum rule (2.23), we obtain the thermodynamic limit of the intensive pressure A(β) = β(1 − p)m p + log 2 + log cosh(βpm p−1 ), (2.44) as claimed. This is in agreement with the results coming from purely statistical mechanics arguments (2.12).

Few words on Burgers hierarchy
The Burgers equation [25,30,87] is one of the most studied nonlinear evolutive equations, and it can be written in the form where α is the viscosity parameter and (. . and so on. This equation is known to emerge as a 1 + 1-dimensional reduction of Navier-Stokes equations for an incompressible fluid in absence of pressure gradient [86,87]. The most important peculiarity of this equation is that it can be linearized in the heat equation via Cole-Hopf transform [34,55]. Furthermore, it is one of the simplest models describing the development and propagation of shock waves. A related and well-studied equation is the Sharma-Tasso-Olver (STO) equation [66,74,83], which can be written as This equation is known to be integrable (in particular, it presents infinitely many symmetries, a bi-Hamiltonian structure, solitary wave solution and an infinite number of conservation laws [66,74,83,85]). These two equations are the lowest elements of the so-called Burgers hierarchy, which can be presented in the form The next two higher equations in the hierarchy are respectively It is clear that the complexity of the elements in the Burgers hierarchy dramatically increases with the index n. However, all of these equations share the same property of Burgers equation, see also [66].

Remark 13. Notice that the Burgers hierarchy is also commonly written in the form
This is related to the form given in (3.3) is given by applying on the former the transformation x → δx with δ = 1 α .
Proof. The proof of this Lemma works in the same way of Lemma 1 in [59]. However, in order to make this Section self-contained, we reported here for the sake of completeness. The proof works by induction, so let us prove it for n = 1 first: Assuming now that the identity holds for n, we will prove it for n + 1. Indeed: Using the thesis for n, the last member of the equation is which proves our assertion.
Lemma 2. The following identity holds: Proof. The proof works by straightforward computation. Indeed

∂ ∂t
Assuming that the function Ψ is analytic in x and t, then Ψ x,t = Ψ t,x , leading to ∂ ∂t (3.14) We are now in position to state the following Proof. First, we perform the Cole-Hopf transform By plugging it into the Burgers hierarchy (3.3) we have Now, using Lemmas 1 and 2, we have By setting the argument of the derivative to zero and assuming that Ψ = 0 for all x and t, we finally get the linear equations Ψ t + α n Ψ n+1,x = 0.

Guerra's mechanical scheme and relation with the Burgers hierarchy
This final Section is devoted to prove the connection between the Guerra's interpolated partition function of the p-spin ferromagnets and the Burgers hierarchy. Before proceeding, it is worth to mention that relations between PDEs and statistical models have been extensively analyzed in the literature, for instance by Ellis and Newman [42] and Bogolyubov and co-workers [27,28]. More recently, PDE methods have gained a consistent role in the analysis of statistical spin models, especially for disordered systems, see for example [20,22,32,46] and references therein.
Having introduced both the players in the duality, we are now in position to prove it. We will start by defining the generalized quantities in the Guerra's interpolation scheme.

Definition 10. The Guerra's generalized partition function is defined as
with associated Boltzmann factor The intensive statistical pressure of the model is Remark 14. We interpret the variable t and x respectively as temporal and spatial coordinates in a 1+1-dimensional space. This interpretation will be clear in a moment.
Remark 15. We stress that the original p-spin model (2.2) (without external fields) is recovered with the choice t = −β and x = 0. The inclusion of a (uniform) magnetic field is reproduced by setting x = h = 0, so that the present framework will still work.
Definition 11. Given a function F (σ) of the spins, its expectation value for the interpolating system (4.1) is defined as We assume that the function A N (t, x) is a differentiable function w.r.t. the space-time coordinates. We can therefore compute its derivatives. In particular, we have In order to find differential equations, we also need higher spatial derivative of the generalized statistical pressure (or equivalently, derivatives of the magnetization expectation value by virtue of (4.5)). It is trivial to note that differentiating the expectation value of the global magnetization would generate expectation values of polynomial in the magnetization itself. This is due to the fact that the x-derivative should increase of a unity the power. This is clear, for example, by considering the first derivative of the expectation value of the magnetization. Indeed, it is clear that A similar relation holds for the expectation value of a generic power of the magnetization: , then we have the following Proposition 5. The functions u q (t, x) satisfy the recurrence relation Proof. The proof of this Proposition follows from a trivial rearrangement of Eq. (4.6).
Practically, we generate the expectation value by repeatedly applying the operator Now, taking q + 1 = p and recalling that u Further, taking the x-derivative of the entire equation, commuting the derivatives ∂ t and ∂ x acting on A N (t, x) and recalling that ∂ x A N (t, x) = u 1 (t, x) (which we simply call u(t, x) for the sake of simplicity), we arrive to state the following Theorem 4. For each p ≥ 2, the expectation value of the global magnetization u(t, x) of the p-spin model described by the Guerra's generalized partition function (4.1) satisfies the equation of the Burgers hierarchy ∂u(t, x) ∂t where α = N −1 is the viscosity parameter.
Remark 16. An alternative route for the proof of the duality is to start with the partition function Z N (t, x) and compute the space-time derivatives. It is easy to see that In this setup, the relation between the expectation value of the magnetization ω t,x (m N ) and the partition function Z N (t, x) is precisely the Cole-Hopf transform, then the duality is easily understood.

The inviscid limit and gradient catastrophe
From the point of view the duality with the Burgers hierarchy, the thermodynamic limit corresponds to the inviscid limit of the Burgers hierarchy α → 0. Then, the whole class of non-linear equations dramatically simplifies, so that we get In order to solve this equation, we also need the initial profile of the solution u 0 (x) = u(t = 0, x). Again, this is a trivial computation, since the initial profile is a 1-body problem. Indeed, we have Therefore, we can state the following Proposition 6. The solution of the self-consistency equation in the thermodynamic limit of the model (4.1) is given by the solution of the initial value problem We stress that such a (first order) system describes the motion of traveling waves in 1 + 1 dimensions with effective velocity v(t, x) = p u(t, x) p−1 . Then, as standard in this case, we can look for solution in implicit form as . Now, recalling that the original p-spin model is recovered with the choice x = 0 and t = −β and that u(−β, 0) = lim N →∞ ω −β,0 (m N (σ)) =m in the inviscid limit (here, we dropped the dependency on t and x, but we again assume the self-averaging property of the order parameter), we easily get the selfconsistency equationm = tanh(βpm p−1 ), (4.11) in perfect agreement with previous results (2.12) and (2.45).
Remark 17. For each p, the model equilbrium dynamics (as resumed in (4.11)) undergoes an ergodicity breaking transition. However, this phase transition is of second order (in the standard Erhenfest classification, i.e. the second derivative of the free energy is discontinuous) only for p = 2, while for all p > 2 the phase transition is of first order (i.e. the first derivative of the free energy is discontinuous).
On the Burgers' side, the inviscid limit has the peculiarity of the appearance of the gradient catastrophe and the related development of shock waves. Indeed, it is easy to understand that these two phenomenons (i.e. gradient catastrophe and ergodicity breaking) are equivalent in this mapping. To understand this, let us analyze the characteristic curves of the system (4.10), given by the system while u is constant along the characteristic curves, i.e. du/dt = 0. The characteristic curves are straight lines which can be parametrized in terms of a quantity ξ as x ξ (t) = ξ + F (ξ)t, where ξ clearly is the x value at t = 0 and F (ξ) = v(u 0 (ξ)) = p tanh(ξ) p−1 . It is well known that the gradient of the solution of (4.10) diverges as t c (ξ) = −1/F (ξ), where the characteristic lines start to cross each other.
Remark 18. By inspecting at its plot, we see that, for even p, the function F (ξ) is only positive, meaning that the gradient catastrophe only takes place for t < 0. This is consistent with our approach, since t = −β and β ∈ R + . On the other hand, rigorously the solution of (4.10) is defined for t > 0. This is not a problem, since the solution can be analytical continuated for negative values of t just before the shock (i.e. for t > t c ). A sketch of the family of characteristic curves is depicted in Fig. 1, left panel.
Since we ultimate want to set x = 0 and t = −β, we should consider only the characteristics for which the gradient catastrophe holds at x = 0. In order to ensure this, we have to choose xξ(t c ) = 0 for someξ, meaning thatξ = −F (ξ)t c . By using the formula t c (ξ) = −1/F (ξ), we have thatξ = F (ξ)/F (ξ). In other words, with the last formula we find for the values of ξ for which the gradient explodes in the position x = 0, then with the standard gradient formula we can compute the time at which the solution develops a shock wave. Recalling that t = −β = − 1 T , we arrive at the following Proposition 7. The critical temperature for the ergodicity breaking phase transition can be identified resolving the following system: where F (ξ) = p tanh(ξ) p−1 .
The results are reported in Fig. 1, right panel.

Burgers hierarchy solution from p-spin thermodynamics
In this Section, we will take benefit of the duality in order to find explicit solutions of the Burgers hierarchy with non-vanishing viscosity α and initial profile u 0 (x) = tanh x. To this aim, we should find explicit expressions for the expectation value of the global magnetization at finite N . In this case, the global magnetization m N (σ) can only take discrete values. To understand this, we can start with a system configuration in which σ i = 1 for all i = 1, . . . , N , whose corresponding magnetization is trivially m N (σ) = 1. All of the other values can be obtained by progressively  N ). Every time we flip a spin, there will be a net difference in the value of the magnetization whose magnitude is 2/N . Therefore, the possible values of the magnetization are With these ingredients, we can write the generalized partition function (4.1) at finite N as (4.14) Now, using the basic relation ω t, .

(4.15)
In the last equation, we used the subscripts p to distinguish between different elements of the Burgers hierarchy and N to stress that the corresponding viscosity parameter is fixed as α = 1/N . As is clear from (4.15), the solutions we can build by explicit use of the duality are rational functions in which both the numerator and denominator are linear combinations (with coefficients not depending on x and t) of exponential waves of the form e Ax+Bt . Other solutions sharing this structure can be found in [59]. To conclude this analysis, we provide some specific examples of solutions of the initial value problem for fixed N and p.
For p = N = 2, the function is solution of the initial value problem For p = 2, N = 3, the function 1 + e 6x + 3e is solution of the initial value problem For p = 2, N = 4, the function is solution of the initial value problem .
We can also vary the value of the order p of the interactions. Indeed, for p = 3 and N = 2, the function u 3,2 (t, x) = −e 2t + e 2x e 2t + e 2x , is solution of the initial value problem Finally, for p = 4 and N = 2, the function is solution of the initial value problem

A representation for finite-size p = 2 solution through Burgers duality
As mentioned above, the interesting feature in the duality between p-spin ferromagnets and the Burgers' hierarchy is that the analysis of the statistical model can be reduced to the study of solutions of linear PDEs for finite N by means of the Cole-Hopf transform. By reversing the duality, we can give a representation of the finite size solution of such spin systems by means of purely PDE methods. In particular, for p = 2 we can take advantage of the heat-kernel technology to find a description of the Curie-Weiss model by deriving an effective self-consistency equation for the order parameter. However, we stress that, for these simple systems, the Burgers duality route is not needed, as the solution is easily found even at finite size N . Despite this, it is worth to analyze how the whole connection works in both senses.
Recall that the Guerra's generalized partition function is with associated free energy The spatial derivative u(t, x) = ∂ x A N (t, x) satisfies the Burgers equation The solution of (interpolated) Curie-Weiss model is equivalent to search the solution of the Burgers equation with initial profile u 0 (x) = u(t = 0, x) = tanh x. Using the Cole-Hopf transformation u(t, x) = 1 N (log Ψ) x , the problem is reduced to the heat equation By using the previous definitions, we can identify A N (t, x) = 1 N log Ψ(t, x), so that the Ψ function is nothing but the Guerra's generalized partition function (4.16).
The heat equation can be solved via the heat kernel technology, so that the general solution is given by where Ψ 0 is the initial profile of the Cauchy problem, and Remark 19. We stress that solutions of the form (4.20) are well-defined for t < 0, due to the "wrong" sign of the temporal derivative in the heat equation (4.19). However, this is coherent with our setup, since the link with the thermodynamic model is achieved with t = −β with β ∈ R + . Now, since u(t, x) = 1 N (log Ψ) x and the initial profile of the solution of Burgers equation is u 0 (x) = tanh x, we immediately have that Ψ 0 (x) = cosh N (x), so that, according to the duality, we have and then the argument of the logarithm is an integral representation of the partition function. In order to make contact with the thermodynamic picture of the Curie-Weiss ferromagnet, we should match the order parameter in terms of the relevant variable on the Burgers side. To do this, we use the fact that ω t,x (m N ) = ∂ x A N (t, x), so we take the spatial derivative of the free energy (4.22). Then, we have . (4.23) By performing the transformation y = −2tȳ + x, we thus have Then, it is proved the following Proposition 8. The (finite-size) expectation value of the global magnetization for the (interpolated) Curie-Weiss model is equivalent to the first moment of the random variableȳ distributed according to the probability distribution where C is a normalization constant.
Remark 20. Since the probability distribution (4.25) is the same associated to the partition function (which in this setup is represented by the Ψ function), we can exactly identify the global magnetization with the random variableȳ.
In terms of theȳ variable, the intensive pressure is (4.22) is We can now take the spatial derivative of the intensive pressure. To do this, we consider the identity Thus, we have By denoting with · t,x the average with respect to the probability distribution P t,x (ȳ) and comparing Eqs. where · = · −β,0 . The substantial difference lies in the fact that in this case the average is performed at finite N , so that the r.h.s. is not expressed as a function of the expectation value of the magnetization, since the probability distribution is smooth (or equivalently, it is not peaked on the equilibrium value of the magnetization). Indeed, we checked that the probability distribution (4.25) at x = 0 and t = −β exhibits the behaviour we expect from the Curie-Weiss picture, see Fig.  3. In particular, below the critical temperature, the probability distribution displays two different peaks (which are related by the symmetry transformation y → −y) becoming sharper and sharper as N increase, ultimately tending to two Dirac deltas. Above the critical temperature, there is a single peak centered at y = 0, mimicking the fact that the system has lost its ferromagnetic behavior.
The comparison between the explicit solution (4.15) and the prediction of (4.24) (for which we checked the equality (4.28)), is reported in Fig. 2.

Conclusions and further developments
Dualities are always powerful tools in mathematical and theoretical physics, since they allow a two-faced investigation of apparently unrelated fields with new mathematical tools, often making easier to derive non-trivial results. In this paper, we examined the relation between the equilibrium dynamics of the p-spin ferromagnetic models and the Burgers hierarchy; in particular, the expectation value of the order parameter in the first side is identified with the solution of the initial value profile of the Burgers hierarchy. We also present some examples of application of the duality on both sides. The methods here developed can in principle be applied to other spin models. Indeed, particularly interesting extensions of the present work would be the application of the PDE-statistical mechanics duality to p-spin systems with random external field [15,61], as well as to diluted and finite-connected ferromagnets [10,44,64,65]. Even more interesting would be the application of PDE-statistcal mechanics duality to disordered system (i.e. spin glass) models with interactions of order p. In this case, as is easily understood, the situation is much harder than the simple p-spin ferromagnets, due to the intrinsic complexity of such systems. Since the structure of pure state is strongly dependent on the realization of the internal disorder and the system size, the structure of Guerra's generalized partition function has to be consistently adapted. For the Sherrington-Kirkpatrick (SK) model, it reads [52] Z and the associated intensive pressure is defined as A N (t, x) = 1 N E log Z N (t, x). Here, both the couplings J ij and the effective external fieldsJ i are normally distributed, i.e. J ij ,J i ∼ N (0, 1) for all i, j = 1, . . . , N , and E stands for the average w.r.t. these random variables. The presence of the square root of the "space-time" coordinates t and x is motivated by an extensive use of the Wick theorem in the computations, so that the derivatives of the intensive pressure does not explicitly depend on t and x. However, for spin-glass models the present methods can only be useful in the thermodynamic limit. This is due to the fact that generally a comprehensive description of spinglass models is possible in the limit N → ∞, where it is possible to use rigorous methods (such as TLCs theorems or concentration inequalities [54,79]) to achieve the thermodynamic solution. On the technical side, even when working at a replica symmetric level, the mapping of the spin-glass model to mechanical systems explicitly depends on the fluctuation of the order parameter (the replica overlap) w.r.t. its thermodynamic value (see for example [1] for transport-like PDEs for the SK and Hopfield models), thus in general we can take benefit of the PDE-statistical mechanics duality only in the thermodynamic limit (where fluctuations vanish). Further, replica symmetry breaking within the Guerra's interpolating framework can be conveniently addressed enriching the interpolation structure. In particular, for the SK model, the K-step RSB equilibrium dynamics is achieved in terms of the generalized partition function where now x = (x (1) , x (2) , . . . , x (K) ) and again J ij ,J (a) i ∼ N (0, 1) for all i, j = 1, . . . , N and a = 1, . . . , K, see [21]. The intensive pressure of the model is now obtained as A N (t, x) = 1 N E 0 log Z 0 (t, x), where Z k−1 (t, x) θ k = E k (Z k (t, x) θ k ) for each k = 1, . . . , K. Here, E k is the expectation value w.r.t. the random variablesJ (k) i , while E 0 is the average w.r.t. the coupling realizations J ij . Finally, the θ k parameters quantify the intensity of each peak in the N → ∞ K-RSB overlap distribution [1,21]. This implies that the K-RSB approximation of spin-glass equilibrium dynamics can in principle be mapped to K + 1-dimensional mechanical systems. Also, the application of the whole technology to equilibrium dynamics of neural networks would be highly desirable. To conclude, we strongly believe that PDE-statistical mechanics duality can be, at least in the thermodynamic limit, a very useful tool to investigate properties at the equilibrium of spin-glass models for each interaction order p [33,36,37,45,69,76,77]. We leave this discussion open for future works.