1 Introduction

Random Matrix Ensembles appear in relation to a variety of problems in mathematics and physics, often showing intriguing and unexpected connections. Symmetric, Hermitian and Symplectic Ensembles have been introduced to describe the statistics of energy levels of heavy nuclei and complex systems [13, 29, 36]; the Circular Unitary Ensemble and the zeros of the Riemann zeta function, respectively, in the thermodynamic limit and in the far limit on the critical line appear to follow the same statistics [18]; the Hermitian Matrix Ensemble (HME) arises from discrete approximations of Topological Field Theory [10, 11, 20, 37]; the partition function of the orthogonal ensemble appears in the calculation of the generating function for specific subsequences of permutations [4]. In addition to the above-mentioned connections, a remarkable result is the identification of the partition function of Random Matrix Ensembles with particular solutions of nonlinear integrable systems [2, 20, 34, 37]. For example, the partition function for the HME can be identified simultaneously with the \(\tau \)-function of a particular solution of the Toda Lattice hierarchy and of the Kadomtsev–Petviashvili hierarchy [1].

In the present paper, we propose an approach to the study of the Symmetric Matrix Ensemble (SME) based on the method of differential identities developed and effectively applied to a variety of statistical mechanical models (see, e.g. [7, 8, 21, 23]). In the case of SMEs, the method relies on its underlying integrable structure realised by the Pfaff Lattice [5, 6, 17, 31, 32]. In the thermodynamic limit, such integrable structure allows the derivation of a system of partial differential equations (PDEs) of dispersionless type for the order parameters. We note that a first direct connection between the dispersionless limit of the Pfaff hierarchy and the thermodynamic limit of matrix models has been studied in [19]. In this paper, we prove that for the SME with even interactions, the order parameters satisfy an integrable hydrodynamic chain of PDEs. Integrable hydrodynamic chains, of which the moments Benney chain is the prototypical example [9], represent an important class of integrable systems of dispersionless type that has attracted a great deal of interest over the last two decades in relation to their classification, construction of new integrable systems via reductions and associated Hamiltonian structures [15, 24, 26,27,28]. We find, at the best of our knowledge, a new example of integrable hydrodynamic chain and prove its integrability via the method of hydrodynamic reductions [14, 15].

The SME is characterised by the partition function of the form

$$\begin{aligned} Z_n({\mathbf {t}})=\int _{{\mathcal {S}}_n} \text {e}^{-H(M)}\, dM, \end{aligned}$$
(1.1)

where dM is the Haar measure, i.e. \(dM :=\prod _{1\le i\le j \le n} dM_{ij}\), and \({\mathcal {S}}_n\) is the set of real symmetric matrices. The function H(M), chosen as

$$\begin{aligned} H(M) = - {{\,\mathrm{Tr}\,}}\left( -\,\frac{M^2}{2}\,+\,\sum _{k\ge 1} t_k M^k \right) , \end{aligned}$$

is referred to as the Hamiltonian and parameters \({\mathbf {t}} = \{t_{k}\}_{k\ge 1}\) are the coupling constants. This terminology refers to Matrix Models of interest in Quantum Field Theory where H(M) is interpreted as the Hamiltonian of the system and \(t_{k}\) are the coupling constants associated with different degrees of interaction [37]. Hence, the free particle theory corresponds to the case where all coupling constants \(t_{k}\) vanish. Notice that for \({\mathbf {t}} = {\mathbf {0}}\), the expression (1.1) reduces to the partition function of the Gaussian Orthogonal Ensemble (GOE).

Based on a classical result by Weyl [35], observing that H(M) depends on the eigenvalues \(\{ z_{k}\}_{k=1}^{n}\) of M only, the integral (1.1) can be reduced to an integral over the eigenvalues of the form

$$\begin{aligned} Z_n({\mathbf {t}})= C_n \int _{{\mathbb {R}}^n} |\Delta _n(z)| \, \prod _{i=1}^n \text {e}^{- \frac{z^2_i}{2} + \sum _{k\ge 1}t_k\,z_i^k}\,dz_i \, \end{aligned}$$
(1.2)

where \(C_n\) is a constant obtained from the integration over the remaining degrees of freedom, and \(\Delta _n(z)\) denotes the Vandermonde determinant \( \Delta _n(z)=\prod _{1\le i < j \le n}(z_i-z_j)\). A fundamental result by Adler and van Moerbeke [5] establishes that for \(2n \times 2n\) symmetric matrices, the function

$$\begin{aligned} \tau _{2 n}({\mathbf {t}}) := \frac{1}{(2n)! C_{2n}}Z_{2n}({\mathbf {t}}) = \frac{1}{(2n)!}\int _{{\mathbb {R}}^{2n}} |\Delta _{2n}(z)| \, \prod _{i=1}^{2n} \text {e}^{- \frac{z^2_i}{2} + \sum _{k\ge 1}t_k\,z_i^k}\,dz_i \, \end{aligned}$$
(1.3)

is the Pfaffian of the moments matrix

$$\begin{aligned} m_{2n}({\mathbf {t}})=\left( \langle x^{i}, y^{j} \rangle _{\mathbf {t}}\right) _{0 \le i,j < 2n-1} \end{aligned}$$

where \(\langle \,\cdot \,, \cdot \, \rangle _{{\mathbf {t}}}\) is a skew-symmetric scalar product induced by the measure on the SME. More specifically, \(\tau _{2n}({\mathbf {t}}) = pf \left( m_{2n} \right) \) is a particular solution of the Pfaff Lattice, an integrable system arising in relation to the algebra splitting of \(gl (\infty )\) into \(sp (\infty )\) and the algebra of \(2\times 2\) blocks lower triangular matrices [5]. The Pfaff Lattice equations are constructed based on the following unique factorisation of the semi-infinite moments matrix

$$\begin{aligned} m_{\infty }({\mathbf {t}})= \left( Q({\mathbf {t}})^{-1} \right) \,J\, \left( Q({\mathbf {t}})^{-1}\right) ^{\top } \end{aligned}$$
(1.4)

where J is the semi-infinite skew-symmetric matrix such that \(J^2=-I\) and Q is a semi-infinite lower triangular matrix. The Lax matrix of the form

$$\begin{aligned} L({\mathbf {t}})= Q({\mathbf {t}}) \Lambda Q({\mathbf {t}})^{-1}, \, \end{aligned}$$
(1.5)

where \(\Lambda = \{\delta _{i,j-1} \}_{i,j=1}^{\infty }\) is the shift matrix with \(\delta _{i,j}\) the Kronecker delta, satisfies the Lax equation [5]

$$\begin{aligned} \frac{\partial L}{\partial t_{k}} = \left[ \,- (L^{k})_{{\mathfrak {t}}}\,, L \, \right] \,. \end{aligned}$$
(1.6)

The projection \(\left( A \right) _{{\mathfrak {t}}}\) for a given matrix A is defined as follows

$$\begin{aligned} A_{{\mathfrak {t}}}:= A_{-} - J\,(A_{+})^{\top }\,J + \frac{1}{2} \left( A_{0} - J (A_0)^{\top } J, \right) \end{aligned}$$
(1.7)

with \(A_{\pm }\) denoting, respectively, the upper and lower triangular part of A, with all \(2\times 2\) diagonal blocks equal to zero, and \(A_{0}\) the block diagonal part of A with \(2\times 2\) diagonal blocks. The entries of the Lax matrix L depend on the sequence of \(\tau \)-functions \(\{\tau _{2n}({\mathbf {t}})\}_{n\ge 1}\) and their derivatives with respect to the coupling constants \(t_{k}\). The Lax matrix associated with the SME partition function is a solution of the Lax equation with initial condition specified by \(\tau _{2n}({\mathbf {t}})\) and its derivatives evaluated at \({\mathbf {t}} = {\mathbf {0}}\). It is important to note that in this case the integrals are specified by the Gaussian measure and can be evaluated using Selberg’s theorem (see, e.g. [22]). The study of the form of the Lax equation (1.6) and its asymptotic properties in the large n limit provides important information on the generic properties of solutions, as for example their singularities and breaking mechanisms, independently of the particular initial condition. More specifically, we observe that the components of the Lax equation (1.6) can be organised as two coupled systems of ODEs, a double chain in infinite components, of the form

$$\begin{aligned} \begin{aligned} \partial _{t_{k}}{\mathbf {v}}_{n}&= F_{k}[{\mathbf {v}},{\mathbf {w}}] \\ \partial _{t_{k}}{\mathbf {w}}_{n}&= G_{k}[{\mathbf {v}},{\mathbf {w}}] \end{aligned} \end{aligned}$$
(1.8)

where \({\mathbf {v}}_{n}\) and \({\mathbf {w}}_{n}\) are the entries of the Lax matrix L suitably recast in the form of infinite vectors, e.g. \({\mathbf {v}}_{n} = (\dots ,v^{-k}_{n},\dots ,v^{-1}_{n},v^{0}_{n},v^{1}_{n},\dots ,v^{k}_{n},\dots )^{\top }\), each associated with a position n on the lattice. \(F_{k}\) and \(G_{k}\) are nonlocal functions on the lattice, evaluated on specific subsets of sites that depend on the chosen \(t_{k}\)-flow.

We then proceed with the study of the Lax equations for SMEs with even power interactions such that the partition function is of the form

$$\begin{aligned} Z_{2n}({\mathbf {t}})=\int _{{\mathcal {S}}_{2n}} \text {e}^{{{\,\mathrm{Tr}\,}}\left( -\,\frac{M^2}{2}\,+\,\sum _{k \ge 1} t_{2k} M^{2k} \right) }\, dM. \end{aligned}$$
(1.9)

The above choice automatically selects a reduction of the even Pfaff Lattice given by the hierarchy (1.6) restricted to the even times \(t_{2k}\). Hence, the system (1.8) is replaced by a single chain of the form

$$\begin{aligned} \begin{aligned} \partial _{t_{2 k}}{\mathbf {w}}_{n}&= H_{k}[{\mathbf {w}}] \end{aligned} \end{aligned}$$
(1.10)

where, similarly to the more general case (1.8), \(H_{k}[{\mathbf {w}}]\) is a nonlocal function on the lattice. Introducing the variable \(x = \varepsilon n\) and the interpolation function

$$\begin{aligned} {\mathbf {u}}(x) := {\mathbf {w}}\left( \frac{x}{\varepsilon } \right) \end{aligned}$$

such that \({\mathbf {u}}(x\pm j \varepsilon ) = {\mathbf {w}}_{n \pm j}\) for some integer j, the thermodynamic limit of the matrix ensemble, i.e. the limit for \(n \rightarrow \infty \), corresponds to the continuum limit of the reduced even Pfaff hierarchy (1.10) obtained by taking \(\varepsilon \rightarrow 0\) such that \(x = \varepsilon n\) remains finite and \({\mathbf {u}}(x)\) is an infinite component vector field of the continuous variable x. Substituting the interpolating function in (1.10) and expanding in Taylor series for \(\varepsilon \rightarrow 0\), at the leading order, one obtains a hierarchy of compatible partial differential equations. A direct calculation performed for the first flows associated with \(t_{2}\), \(t_{4}\) and \(t_{6}\) shows that the equations so obtained constitute an infinite system of first-order PDEs of hydrodynamic type, referred to as hydrodynamic chains, of the form

$$\begin{aligned} \mathbf{u}_{T_{2 k}} = A^{(2 k)}(\mathbf{u}) \, \mathbf{u}_{x} \end{aligned}$$
(1.11)

where \(T_{2k}\) corresponds to the time variable \(t_{2k}\) suitably rescaled, e.g. \(T_2 = \varepsilon t_2\). \(A^{(2k)}({\mathbf {u}})\) is a sparse matrix where each row contains a finite number of elements depending on a finite number of components of the vector field \({{\mathbf {u}}(x)}\). Infinite matrices of this type are referred to as class C (chain-class) matrices [15]. We conjecture that the form (1.11) holds for all \(k=1,2,\dots \) . We do not have at the moment a proof of this conjecture.

We prove that the hydrodynamic chain so obtained passes the diagonalisability test (see Proposition 5) introduced in [15] and is integrable in the sense of hydrodynamic reductions (Theorem 7) [14, 15]. As the compatibility of the hierarchy (1.11) implies that the matrices \(A^{(2k)}\) commute, it is sufficient to perform the diagonalisability test for the first flow of the hierarchy only. The test establishes that a hydrodynamic chain of class C defined by the infinite matrix \(A({\mathbf {u}}) = \{a^{i}_{j} \}_{i,j=-\infty }^{\infty }\) is diagonalisable if and only if all components of the Haantjes tensor

$$\begin{aligned} H^{i}_{jk} = N^{i}_{pr} \, a^{p}_{j} \, a^{r}_{k} - N^{p}_{jr} \, a^{i}_{p} \, a^{r}_{k} - N^{p}_{rk} \, a^{i}_{p} \, a^{r}_{j} + N^{p}_{jk} \, a^{i}_{r} \, a^{r}_{p}, \end{aligned}$$
(1.12)

where \(N^i_{jk}\) is the Nijenhuis tensor

$$\begin{aligned} N^{i}_{jk} = a^{p}_{j} \, \partial _{u^{p}} a^{i}_{k} -a^{p}_{k} \, \partial _{u^{p}} a^{i}_{j} - a^{i}_{p} \left( \partial _{u^{j}} a^{p}_{k} -\partial _{u^{k}} a^{p}_{j} \right) \,, \end{aligned}$$
(1.13)

vanish identically. The notion of integrability in the sense of hydrodynamic reductions for a hydrodynamic chain extends the similar concept introduced in the context of finite component systems [14] and characterises the chain under consideration as integrable if it admits N-phase solutions of the form \({\mathbf {u}}(R^{1},\dots ,R^{N})\) for any integer N, where \(R^{i} = R^{i}(x,t)\) (Riemann invariants) satisfy the semi-Hamiltonian diagonal system of hydrodynamic type

$$\begin{aligned} R^{i}_{t} = \lambda ^{i}\left( R^{1},\dots ,R^{N}\right) R^{i}_{x} \qquad i=1,\dots ,N. \end{aligned}$$
(1.14)

The time t can be identified, subject to a suitable re-scaling, with any of the “times” \(t_{2k}\) and \(\lambda ^{i}\left( R^{1},\dots ,R^{N}\right) \) are the characteristic speeds of the corresponding flow. The system for Riemann invariants (1.14) is required to fulfill the semi-Hamiltonian property which can be expressed in terms of the following differential constraint on the characteristic speeds

$$\begin{aligned} \partial _{k} \left( \frac{\partial _{j} \lambda ^{i}}{\lambda ^{j} - \lambda ^{i}} \right) = \partial _{j} \left( \frac{\partial _{k} \lambda ^{i}}{\lambda ^{k} - \lambda ^{i}} \right) , \qquad i \ne j \ne k, \end{aligned}$$
(1.15)

with the notation \(\partial _{i} = \partial / \partial R^{i}\). The condition (1.15) guarantees that equation (1.14) constitutes a system of conservation laws [30]. A classical result by Tsarev establishes that the system (1.14) is completely integrable by the generalised hodograph method ([33], see also [12]) and the solution is given by the following equation

$$\begin{aligned} x + \lambda ^{i}\left( R^{1},\dots ,R^{N}\right) t = \mu ^{i}\left( R^{1},\dots ,R^{N}\right) \qquad i =1,\dots ,N \end{aligned}$$

where the functions \(\mu ^{i}\left( R^{1},\dots ,R^{N}\right) \) satisfy the system of linear PDEs of the form

$$\begin{aligned} \frac{\partial _{j} \lambda ^{i}}{\lambda ^{i} - \lambda ^{j}} = \frac{\partial _{j} \mu ^{i}}{\mu ^{i} - \mu ^{j}}. \end{aligned}$$
(1.16)

The solution to the system (1.16) is parametrised via N functions of one variable that can be fixed by the initial conditions on the functions \(R^{i}\).

The paper is organised as follows. In Sect. 2, we review some results regarding SMEs and their relationship with the Pfaff Lattice (see, e.g. [5, 34]). In Sect. 3, we study the structure of the Lax matrix (1.5) and write the explicit evolution equations for the first flow of the Pfaff Lattice. Section 4 is devoted to the SME with even degree interactions and its relation with the even Pfaff hierarchy. In Sect. 5, we study the thermodynamic limit and show that the resulting hydrodynamic chain is diagonalisable and integrable. We then conclude with some final remarks in Sect. 6. Appendices provide some additional technical details as well as elements that are subject of further studies. The expressions for the second flow of the Pfaff Lattice with odd and even times are provided in Appendix A; the explicit form of the coupled system (1.8) and its higher order corrections are given in Appendix B; higher order corrections to the hydrodynamic chain (1.10) are reported in Appendix C; Appendix D lists nonzero entries of the Nijenhuis tensor (1.13) for the hydrodynamic chain (1.10).

2 Symmetric Matrix Ensemble and Pfaff Lattice

In this section, we briefly review definitions and properties of the SME and its connection with the Pfaff Lattice with a focus on aspects that are relevant for the purposes of this paper [3, 5, 34]. As mentioned above, the partition function (1.1) can be reduced, up to a proportionality constant, to the integral over the eigenvalues of the form

$$\begin{aligned} \tau _{2 n}({\mathbf {t}}) := \frac{1}{(2n)!}\int _{{\mathbb {R}}^{2n}} |\Delta _{2n}(z)| \, \prod _{i=1}^{2n} \text {e}^{- \frac{z^2_i}{2} + \sum _{k\ge 1}t_k\,z_i^k}\,dz_i \,, \end{aligned}$$

where \(\tau _{2n}\) is referred to as Pfaffian \(\tau \)-function. The function \(\tau _{2n}\) is in fact the Pfaffian of the \({\mathbf {t}}\)-dependent moment matrix \(m_{2n}({\mathbf {t}})=\big (\mu _{ij}({\mathbf {t}})\big )_{0 \le i,j < 2n-1}\), i.e.

$$\begin{aligned} \tau _{2n}({\mathbf {t}}) = \text {pf}(m_{2n}({\mathbf {t}}))=\left( \text {det} \, m_{2n}({\mathbf {t}}) \right) ^{1/2} \, \end{aligned}$$
(2.1)

where entries of \(m_{2n}({\mathbf {t}})\) are constructed via the skew-symmetric inner product \(\langle \cdot \,,\,\cdot \rangle _{{\mathbf {t}}} \)

$$\begin{aligned} \mu _{ij}({\mathbf {t}}) = \langle x^i\,,\,y^j \rangle _{{\mathbf {t}}} := \int \int _{{\mathbb {R}}^2} x^i \, y^j \,\sigma (x-y)\,\text {e}^{-\frac{1}{2} (x^2 + y^2) + \sum _{k}t_k\,(x^k+y^k)}\,dx\,dy \nonumber \\ \end{aligned}$$
(2.2)

with \(\sigma (x)=\text {sign}(x)\). Noting that \(\mu _{ij}=-\mu _{ji}\), the moments matrix \(m_{2n}\) is skew-symmetric and takes the form

$$\begin{aligned} m_{2n}=\big ( \mu _{i\,j} \big )_{0\le i,j\le 2n-1}= \begin{pmatrix} 0 \quad &{} \mu _{0\,1} \quad &{} \mu _{0\,2} \quad &{} \mu _{0\,3} \quad &{} \mu _{0\,4} \quad &{} \mu _{0\,5} \quad &{} \dots \\ -\mu _{0\,1} \quad &{} 0 \quad &{} \mu _{1\,2} \quad &{} \mu _{1\,3} \quad &{} \mu _{1\,4} \quad &{} \mu _{1\,5} \quad &{} \dots \\ -\mu _{0\,2} \quad &{} -\mu _{1\,2} \quad &{} 0 \quad &{} \mu _{2\,3} \quad &{} \mu _{2\,4} \quad &{} \mu _{2\,5} \quad &{} \dots \\ -\mu _{0\,3} \quad &{} -\mu _{1\,3} \quad &{} -\mu _{2\,3} \quad &{} 0 \quad &{} \mu _{3\,4} \quad &{} \mu _{1\,5} \quad &{} \dots \\ -\mu _{0\,4} \quad &{} -\mu _{1\,4} \quad &{} -\mu _{2\,4} \quad &{} -\mu _{3\,4} \quad &{} 0 \quad &{} \mu _{4\,5} \quad &{} \dots \\ -\mu _{0\,5} \quad &{} -\mu _{1\,5} \quad &{} -\mu _{2\,5} \quad &{} -\mu _{3\,5} \quad &{} -\mu _{4\,5} \quad &{} 0 \quad &{} \dots \\ \vdots \quad &{} \vdots \quad &{} \vdots \quad &{} \vdots \quad &{} \vdots \quad &{} \vdots \quad &{} \ddots \\ \end{pmatrix}\nonumber \\ \end{aligned}$$
(2.3)

The evolution equations of matrix elements \(\mu _{ij}\) with respect to the coupling constants \(\{t_k\}_{k\in {\mathbb {N}}}\) follow from the direct differentiation from the definition (2.2) and read as

$$\begin{aligned} \frac{\partial \mu _{ij}}{\partial t_k}=\mu _{i+k,\,j}+\mu _{i,\,j+k} . \end{aligned}$$
(2.4)

Therefore, the semi-infinite moment matrix \(m_{\infty }\) satisfies the equation

$$\begin{aligned} \frac{\partial m_{\infty }}{\partial t_k} = \Lambda ^k\,m_{\infty }+m_{\infty } \, \Lambda ^k \,, \end{aligned}$$
(2.5)

where \(\Lambda \) is the shift matrix

$$\begin{aligned} \Lambda =\begin{pmatrix} 0 &{}\quad 1 &{} \quad 0 &{}\quad 0 &{}\quad \dots \\ 0 &{}\quad 0 &{}\quad 1 &{}\quad 0 &{}\quad \dots \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad 1 &{}\quad \dots \\ \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \vdots &{}\quad \ddots \\ \end{pmatrix}. \end{aligned}$$
(2.6)

As mentioned above regarding Eq.  (1.4), the matrix \(m_{\infty }\) admits the unique factorisation [1, 3, 34]

$$\begin{aligned} m_{\infty }({\mathbf {t}})= \left( Q({\mathbf {t}})^{-1} \right) \,J\,\left( Q({\mathbf {t}})^{-1}\right) ^{\top } \end{aligned}$$
(2.7)

where J is the semi-infinite skew-symmetric matrix such that \(J^2=-I\) and Q is a semi-infinite lower triangular matrix of the form

figure a

Due to the factorisation (2.7), elements of the matrix Q depend on the moments \(\mu _{ij}\) as well as the Pfaffian \(\tau \)-function \(\tau _{2n}\). From Eq.  (2.4), it follows that the evolution of the Pfaffian \(\tau \)-function with respect to the k-th time can be written as

$$\begin{aligned} \frac{\partial \tau _{2n}}{\partial t_k} = \sum _{i,j=0}^{2n-1} \frac{\partial \tau _{2n}}{\partial \mu _{i,j}} \, \frac{\partial \mu _{i,j}}{\partial t_k} = \sum _{i,j=0}^{2n-1} \frac{\partial \tau _{2n}}{\partial \mu _{i,j}} \, \left( \mu _{i+k,j} + \mu _{i,j+k} \right) . \end{aligned}$$
(2.9)

Therefore, elements of the decomposition matrix Q are expressed in terms of \(\tau _{2n}\) and suitable combinations of its derivatives with respect to the times \(\{t_k\}_{k\in {\mathbb {N}}}\) determined by the Schur’s polynomials of the differential operators \(\{\partial _{t_k}\}_{k\in {\mathbb {N}}}\) [34].

The factorisation of the moments matrix allows to define the Lax matrix

$$\begin{aligned} L({\mathbf {t}})= Q({\mathbf {t}}) \Lambda Q({\mathbf {t}})^{-1} \end{aligned}$$
(2.10)

for which the following theorem holds:

Theorem 1

([5]) The function \(\tau _{2n}\) is a \(\tau \)-function for the Pfaff Lattice, i.e. the following operator

$$\begin{aligned} L({\mathbf {t}})= & {} Q({\mathbf {t}}) \Lambda Q({\mathbf {t}})^{-1} \nonumber \\= & {} \begin{pmatrix} 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ * \quad &{}\quad \partial _{t_1} \log \tau _2 \quad &{}\quad \left( \dfrac{\tau _4 \, \tau _0}{\tau _2^2} \right) ^{1/2} \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \ddots \\ * \quad &{}\quad * \quad &{}\quad -\partial _{t_1} \log \tau _2 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \ddots \\ * \quad &{}\quad * \quad &{}\quad * \quad &{}\quad \partial _{t_1} \log \tau _4 \quad &{}\quad \left( \dfrac{\tau _6 \, \tau _2}{\tau _4^2} \right) ^{1/2} \quad &{}\quad 0 \quad &{}\quad \ddots \\ * \quad &{}\quad * \quad &{}\quad * \quad &{}\quad * \quad &{}\quad -\partial _{t_1} \log \tau _4 \quad &{}\quad 1 \quad &{}\quad \ddots \\ \vdots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \\ \end{pmatrix} \end{aligned}$$
(2.11)

satisfies the commuting equations

$$\begin{aligned} \frac{\partial L}{\partial t_k} = \left[ \,- (L^k)_{{\mathfrak {t}}}\,, L \, \right] \, \qquad k\in {\mathbb {N}}, \end{aligned}$$
(2.12)

Eq. (2.12) is referred to as the Lax equation of the Pfaff Lattice.

The star symbols \(*\) in the expression of the Lax matrix (2.11) stand for suitable differential expressions of the \(\tau \)-functions \(\tau _{2k}\) and \(A_{{\mathfrak {t}}}\) is the projection defined in (1.7). This follows from the splitting of the Lie algebra \(\text {gl}(\infty )\)

$$\begin{aligned} \text {gl}(\infty ) = {\mathfrak {t}} \oplus {\mathfrak {n}} {\left\{ \begin{array}{ll} {\mathfrak {t}} = \{\text {lower triangular matrices of the form} \, (2.8)\} \\ {\mathfrak {n}} = \text {sp}(\infty ) = \{ A \text { such that } J A^{\top } J = A\} \end{array}\right. } \end{aligned}$$
(2.13)

which yields the unique decomposition

$$\begin{aligned} A= & {} A_{{\mathfrak {t}}} + A_{{\mathfrak {n}}} \nonumber \\= & {} A_{-} - J\,(A_{+})^{\top }\,J + \frac{1}{2} \left( A_{0} - J (A_0)^{\top } J \right) + A_{+} + J\,(A_{+})^{\top }\,J \nonumber \\&+ \frac{1}{2} \left( A_{0} + J (A_0)^{\top } J \right) \, \end{aligned}$$
(2.14)

where \(A_{\pm }\) denote, respectively, the upper and lower triangular part of A, with all \(2\times 2\) diagonal blocks equal to zero, and \(A_{0}\) the block diagonal part of A with \(2\times 2\) diagonal blocks.

3 Lattice equations and initial conditions for the Pfaff hierarchy

In this section, we further investigate the structure of the Lax equation (2.12). Our main observation is that the Lax equation can be recast in the form of a two-component infinite chain.

Let us introduce the following notation for the Lax matrix (2.11)

$$\begin{aligned} L=\begin{pmatrix} 0 &{}\quad 1 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ w_1^{-1} &{}\quad v^{0}_{1} &{} \quad w_1^0 &{}\quad 0 &{} \quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ v^{-1}_{1} &{}\quad w_1^{1} &{}\quad -v^{0}_{1} &{}\quad 1 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ w_1^{-2} &{}\quad v^{1}_{1} &{}\quad w_2^{-1} &{}\quad v^{0}_{2} &{}\quad w_2^0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ v^{-2}_{1} \quad &{}\quad w_1^2 \quad &{}\quad v^{-1}_{2} \quad &{}\quad w_2^1 \quad &{}\quad -v^{0}_{2} &{}\quad 1 &{}\quad 0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ w_1^{-3} &{}\quad v^{2}_{1} &{}\quad w_2^{-2} &{}\quad v^{1}_{2} &{}\quad w_3^{-1} &{}\quad v^{0}_{3} &{}\quad w_3^0 &{}\quad 0 &{}\quad 0 &{}\quad \dots \\ v^{-3}_{1} &{}\quad w_1^3 &{}\quad v^{-2}_{2} &{}\quad w_2^{2} &{}\quad v^{-1}_{3} &{}\quad w_3^1 &{}\quad -v^{0}_{3} &{}\quad 1 &{}\quad 0 &{}\quad \dots \\ \vdots &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \\ \end{pmatrix}. \end{aligned}$$
(2.1)

The variables \(\{w^{0}_{n}\}_{n\ge 1}\) constitute the non-constant entries in the first upper diagonal (even positions) of L, and \(\{v^{0}_{n}\}_{n\ge 1}\) the entries in the main diagonal of L. In the lower triangular part, for any \(k>0\), \(\{w^{-k}_n\}_{n\ge 1}\) and \(\{w^{k}_n\}_{n\ge 1}\) occupy, respectively, odd and even positions on the \((2 k - 1)\)-th diagonal. Similarly, the variables \(\{v^{-k}_n\}_{n\ge 1}\) and \(\{v^{k}_n\}_{n\ge 1}\) occupy odd and even positions on the 2k-th diagonal. The evolution equations for the variables \(v^{k}_{n}\) and \(w^{k}_{n}\) follow from the Lax equation (2.12). For instance, the \(t_1\)-flow for the variables \(v^{\pm k}_n\) and \(w^{\pm k}_n\), respectively, is given by following equations

$$\begin{aligned} \begin{aligned} \partial _{t_1} v^{k}_n =&\,\frac{1}{2} \left( v^0_{n-1} + v^0_n - v^0_{n-k -1} - v^0_{-k+n}\right) v^{k}_n + w^{k-1}_{n} - w^0_n w^{-(k+1)}_{n+1} \\&- w^{-1}_n w^{-k}_n - w^0_{n-1}w^{-(k-1)}_{n-1}, \qquad k<-1 \\ \partial _{t_1} v^{-1}_n =&\, \frac{1}{2}\left( v^0_{n-1} - v^0_{n+1}\right) v^{-1}_n + w^{-2}_n - w^0_n - w^{-1}_n w^1_n - w^0_{n-1}w^2_{n-1} \\ \partial _{t_1} v^0_n =&\, w^0_n w^1_n\\ \partial _{t_1} v^1_n =&\, \frac{1}{2}\left( v^0_{n+1} - v^0_{n-1}\right) v^1_n - w^{-2}_n + w^0_n + w^{-1}_{n+1}w^1_n+w^0_{n+1}w^2_n \\ \partial _{t_1}v^k_n =&\, \frac{1}{2}\left( v^0_{k+n} + v^0_{k+n-1}-v^0_n - v^0_{n-1}\right) v^k_n + w^0_{n+k-1} w^{k-1}_n + w^{-1}_{n+k} w^k_n \\&+ w^0_{n+k} w^{k+1}_n - w^{-(k+1)}_n, \;\;\;\; k>1 \end{aligned} \end{aligned}$$
(2.2)

and

$$\begin{aligned} \partial _{t_1} w^{k}_n =&\,\frac{1}{2} \left( v^0_{n-k-1} + v^0_{n-k-2}+ v^0_{n} + v^0_{n-1}\right) w^{k}_n + w^{0}_{n-k-2}v^{k+2}_{n} - w^0_n v^{-(k+2)}_{n+1} \nonumber \\&+ w^{-1}_{n-k-1} v^{k+1}_n - w^{-1}_{n}v^{-(k+1)}_{n} +w^0_{n-k-1} v^{k}_n - w^0_{n-1}v^{-k}_{n-1}, \qquad k<-1\nonumber \\ \partial _{t_1} w^{-1}_n =&\,w^{0}_n v^{-1}_n - w^0_{n-1}v^1_{n-1} \\ \partial _{t_1} w^0_n =&\,\frac{1}{2}\left( v^0_{n+1} -2 v^0_n + v^0_{n-1}\right) w^0_n\nonumber \\ \partial _{t_1} w^k_n =&\, -\frac{1}{2}\left( v^0_{n+k} + v^0_{n+k-1} +v^0_n + v^0_{n-1}\right) w^k_n +v^k_n - v^{-k}_n, \qquad k > 0. \nonumber \end{aligned}$$
(2.3)

To give an idea of the increasing complexity of higher flows, the \(t_2\)-flows for both variables \(v^{\pm k}_n\) and \(w^{\pm k}_n\) are also reported in Appendix A.

We now consider the initial condition for the Lax matrix L. From (2.11), we have that the component \(w^{0}_n({\mathbf {t}})\) can be expressed in terms of \(\tau _{2n}({\mathbf {t}})\) as follows (see [34])

$$\begin{aligned} w^{0}_n({\mathbf {t}})=\bigg (\frac{\tau _{2n+2}({\mathbf {t}})\,\tau _{2n-2}({\mathbf {t}})}{\tau _{2n}^2({\mathbf {t}})}\bigg )^{1/2}\,. \end{aligned}$$
(2.4)

The function \(\tau _{2n}({\mathbf {0}})\) is given by a Selberg’s integral which can be evaluated explicitly so that

$$\begin{aligned} \tau _{2n}({\mathbf {0}})=\pi ^{n/2}\, \prod _{k=0}^{n-1} 2^{-2k} (2k)!. \end{aligned}$$
(2.5)

Therefore, Eqs. (2.4) and (2.5) imply

$$\begin{aligned} w^{0}_n({\mathbf {0}}) = 2 \sqrt{\pi } \sqrt{2n(2n-1)}. \end{aligned}$$
(2.6)

Similarly, using the expression for \(v^0_n({\mathbf {t}})\) obtained in [34], i.e.

$$\begin{aligned} v^{0}_n({\mathbf {t}}) = \partial _{t_1}\log \tau _{2n}({\mathbf {t}}). \end{aligned}$$
(2.7)

one can evaluate the initial datum \(v^0_n({\mathbf {0}})\). Hence, from the definition of \(\tau _{2n}({\mathbf {t}})\) given in (1.3), and due to the skew symmetry of the integration measure we have

$$\begin{aligned} v^{0}_n({\mathbf {0}}) = 0. \end{aligned}$$
(2.8)

In general, the variables \(v^{\pm k}_n\) are represented as suitable combinations of integrals of odd functions and therefore \(v^{\pm k}_n({\mathbf {0}}) = 0\). We conclude that the Lax matrix L evaluated at \({\mathbf {t}}= {\mathbf {0}}\) takes the following form

$$\begin{aligned} L(\mathbf{{0}}) =\begin{pmatrix} 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ w_1^{-1}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_1^0(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^{1}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ w_1^{-2}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_2^{-1}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_2^0(\mathbf{0})\quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^2(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_2^1(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ w_1^{-3}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_2^{-2}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_3^{-1}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_3^0(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^3(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_2^{2}(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad w_3^1(\mathbf{0}) \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad \dots \\ \vdots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \\ \end{pmatrix}\nonumber \\ \end{aligned}$$
(2.9)

4 The reduced even Pfaff hierarchy

In this section, we consider the SME with even power interactions specified by the partition function (1.9) and show that it provides a solution to it provides a solution to a reduction of the even Pfaff hierarchy, i.e. the commuting flows (2.12) associated with the even times \(t_{2k}\) only.

In this case, Eq. (2.1), i.e. \(\tau _{2n}({\mathbf {t}}) = \mathrm{pf}(m_{2n}({\mathbf {t}}))\), still holds with \(m_{2n} = (\mu _{ij})_{0\le i,j\le 2n-1}\) and

$$\begin{aligned} \mu _{ij}({\mathbf {t}})= & {} \langle x^i\,,\,y^j \rangle _{{\mathbf {t}}} \nonumber \\= & {} \int \int _{{\mathbb {R}}^2} x^i \, y^j \,\sigma (x-y)\,\text {e}^{\sum _{k\ge 1}t_{2k}\,(x^{2k}+y^{2k})}\text {e}^{-\frac{1}{2}(x^2 + y^2)}\,dx\,dy. \end{aligned}$$
(2.1)

Hence, the moments matrix \(m_{2n}({\mathbf {t}})\) reads as

$$\begin{aligned} m_{2n}=\big ( \mu _{i\,j} \big )_{0\le i,j\le 2n-1}= \begin{pmatrix} 0 \quad &{}\quad \mu _{0\,1} \quad &{}\quad 0 \quad &{}\quad \mu _{0\,3} \quad &{}\quad 0 \quad &{}\quad \mu _{0\,5} \quad &{}\quad \dots \\ -\mu _{0\,1} \quad &{}\quad 0 \quad &{}\quad \mu _{1\,2} \quad &{}\quad 0 \quad &{}\quad \mu _{1\,4} \quad &{} \quad 0 \quad &{} \quad \dots \\ 0 \quad &{}\quad -\mu _{1\,2} \quad &{}\quad 0 \quad &{}\quad \mu _{2\,3} \quad &{}\quad 0 \quad &{}\quad \mu _{2\,5} \quad &{}\quad \dots \\ -\mu _{0\,3} \quad &{}\quad 0 \quad &{}\quad -\mu _{2\,3} \quad &{}\quad 0 \quad &{}\quad \mu _{3\,4} \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad -\mu _{1\,4} \quad &{}\quad 0 \quad &{}\quad -\mu _{3\,4} \quad &{}\quad 0 \quad &{}\quad \mu _{4\,5} \quad &{}\quad \dots \\ -\mu _{0\,5} \quad &{}\quad 0 \quad &{}\quad -\mu _{2\,5} \quad &{}\quad 0 \quad &{}\quad -\mu _{4\,5} \quad &{}\quad 0 \quad &{}\quad \dots \\ \vdots \quad &{}\quad \vdots \quad &{}\quad \vdots \quad &{}\quad \vdots \quad &{}\quad \vdots \quad &{}\quad \vdots \quad &{}\quad \ddots \\ \end{pmatrix}.\nonumber \\ \end{aligned}$$
(2.2)

The moments (2.1) satisfy the evolution equations

$$\begin{aligned} \frac{\partial \mu _{ij}}{\partial t_{2k}}=\mu _{i+2k,j}+\mu _{i,j+2k} \end{aligned}$$
(2.3)

which imply

$$\begin{aligned} \frac{\partial m_{\infty }}{\partial t_{2k}} = \Lambda ^{2k}\,m_{\infty }+m_{\infty } \, \Lambda ^{2k} \,. \end{aligned}$$
(2.4)

We consider the reduction of the Lax equation (2.12) of the form

$$\begin{aligned} \frac{\partial L}{\partial t_{2k}} = \left[ \,- (L^{2k})_{{\mathfrak {t}}}\,, L \, \right] \,, \end{aligned}$$
(2.5)

with

$$\begin{aligned} L = \begin{pmatrix} 0 \quad &{}\quad 1 \quad &{} \quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ w_1^{-1} \quad &{}\quad 0 \quad &{}\quad w_1^0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^{1} \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{} \quad \dots \\ w_1^{-2} \quad &{}\quad 0 \quad &{}\quad w_2^{-1} \quad &{} \quad 0 \quad &{}\quad w_2^0\quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^2 \quad &{}\quad 0 \quad &{}\quad w_2^1 \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{} \quad 0 \quad &{}\quad \dots \\ w_1^{-3} \quad &{}\quad 0 \quad &{}\quad w_2^{-2} \quad &{}\quad 0 \quad &{} w_3^{-1} \quad &{}\quad 0 \quad &{}\quad w_3^0 \quad &{}\quad 0 \quad &{}\quad 0 \quad &{}\quad \dots \\ 0 \quad &{}\quad w_1^3 \quad &{}\quad 0 \quad &{}\quad w_2^{2} \quad &{}\quad 0 \quad &{}\quad w_3^1 \quad &{}\quad 0 \quad &{}\quad 1 \quad &{}\quad 0 \quad &{}\quad \dots \\ \vdots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \quad &{}\quad \ddots \\ \end{pmatrix}, \end{aligned}$$
(2.6)

that is the Lax matrix associated with the SME with even power interactions is obtained from the general one by setting the variables \(v^0_n\), \(v^{\pm k}_n\) identically equal to zero for any \({\mathbf {t}}_{k}\). In other words, the partition function gives a solution to a reduction of the even Pfaff hierarchy which preserves the zeros of the initial Lax matrix \(L({\mathbf {0}})\) given by the expression (2.9).

The first non-trivial flow of the reduced even Lax hierarchy (2.5) provides the following evolution equations for the variables \(w^k_n\)

$$\begin{aligned} \begin{aligned} \partial _{t_2} w^{k}_{n} =&~ \frac{1}{2} \left( w^{k}_{n} w^{0}_{n} w^{1}_{n} + w^{k}_{n} w^{0}_{n-k-1} w^{1}_{n-k-1} - w^{k}_{n} w^{0}_{n-1} w^{1}_{n-1} - w^{k}_{n} w^{0}_{n-k-2} w^{1}_{n-k-2} \right) \\&+w^{k+1}_{n+1} w^{0}_{n} + w^{k-1}_{n} w^{0}_{n-k-1} - w^{k-1}_{n-1} w^{0}_{n-1} - w^{k+1}_{n} w^{0}_{n-k-2}, \qquad k <- 1\\ \partial _{t_2} w^{-1}_{n}=&~ w^{0}_{n} \left( w^{-1}_{n} w^{1}_{n} + w^{-2}_{n} + w^{0}_{n} \right) - w^{0}_{n-1} \left( w^{-1}_{n} w^{1}_{n-1} + w^{-2}_{n-1} \right) - \left( w^{0}_{n-1} \right) ^{2} \\ \partial _{t_2} w_{n}^{0} =&~ \frac{1}{2} \left( w^{0}_{n+1} w^{1}_{n+1} - w^{0}_{n-1} w^{1}_{n-1} \right) w^0_n + \left( w^{-1}_{n+1} - w^{-1}_{n}\right) w^0_n \\ \partial _{t_2} w_{n}^{1} =&~ \frac{1}{2} \left( w^{0}_{n-1} w^{1}_{n-1} w^{1}_{n} - w^{0}_{n+1} w^{1}_{n} w^{1}_{n+1}\right) + w^{0}_{n+1} w^{2}_{n} - w^{0}_{n-1} w^{2}_{n-1} \\ \partial _{t_2} w^{k}_{n} =&~ \frac{1}{2} \left( w^{0}_{n-1} w^{1}_{n-1} w^{k}_{n} + w^{0}_{n+k-1} w^{1}_{n+k-1} w^{k}_{n} -w^{0}_{n} w^{1}_{n} w^{k}_{n} - w^{0}_{n+k} w^{1}_{n+k} w^{k}_{n} \right) \\&+ w^{0}_{n} w^{k-1}_{n+1} + w^{0}_{n+k} w^{k+1}_{n} - w^{0}_{n-1} w^{k+1}_{n-1} - w^{0}_{n+k-1} w^{k-1}_{n}, \qquad k> 1.\nonumber \end{aligned}\!\!\!\!\!\!\\ \end{aligned}$$
(2.7)

The derivation described above naturally compares with the case of the HME, as studied in [8], where the partition function corresponds to a particular solution of the Toda Lattice and the reduction of even power interactions provides a solution to the Volterra Lattice. The Volterra Lattice is effectively an independent integrable system as it arises from a reduction in the even Toda hierarchy and it is not compatible with the odd flows of the Toda hierarchy. Similarly, the reduction of the even Pfaff hierarchy obtained from the SME with even interactions is not compatible with the odd flows of the Pfaff hierarchy.

5 Thermodynamic limit and integrable hydrodynamic chain

We study the large n asymptotic properties of the SME via the continuum limit of the Pfaff Lattice equations. In particular, we focus on the case of even power interactions (1.9) described by the Lax matrix (2.6) which satisfies the Lax equations (2.5).

As observed above, the lattice equations for the reduced even Pfaff hierarchy (2.7) constitute an infinite chain for the variables \(w^{k}_{n}\), where \(k \in {\mathbb {Z}}\) labels the components of the chain and \(n \in {\mathbb {N}}\) labels points on the lattice. In Sect. 3, we noted that the variables \(w^{k}_{n}\) can be expressed in terms of suitable elements of the sequence of \(\tau \)-functions \(\{\tau _{2n}\}_{n\ge 1}\) and their derivatives. As \(n \rightarrow \infty \), for the variables \(w^{k}_{n}\) we have \(w^{k}_{n+1} - w^{k}_{n} = O(\varepsilon )\), with \(\varepsilon \rightarrow 0\) such that \(x = \varepsilon n\) remains finite. In the following, we derive the continuum limit equations for the chain and study the integrability at the leading order with respect to the \(\varepsilon \) expansion. We illustrate the result for the first equation of the hierarchy given by the \(t_{2}\)-flows. As mentioned in Sect. 1, our considerations extend to the \(t_4\)- and \(t_6\)-flows as well, and we conjecture they hold for any equation of the hierarchy.

Let us introduce the interpolation function \(w^{k}(x/\varepsilon )\) with \(x = \varepsilon n\) so that \(w^{k}(n) = w^{k}_{n}\), and define

$$ u^{k}(x) := w^{k}\left( \frac{x}{\varepsilon } \right) $$

with \(u^{k}(x \pm \varepsilon ) = w^{k}_{n \pm 1}\). Substituting \(u^k(x)\) into Eq. (2.7), expanding in Taylor series for \(\varepsilon \rightarrow 0\) and setting \(t = \varepsilon \, t_{2}\), at the leading order \(O(\varepsilon ^{0})\) we get the following system of PDEs

$$\begin{aligned} \begin{aligned} u^{k}_{t} =&\left( (k+2) u^{k+1} - k u^{k-1} + u^{1} u^{k} \right) u^{0}_{x} + u^{0} u^{k} u^{1}_{x} + u^{0} u^{k-1}_{x} + u^{0} u^{k+1}_{x}, \ k < 0\\ u^{0}_{t} =&~u^{0} u^{1} u^{0}_{x} + \left( u^{0} \right) ^{2} u^{1}_{x} + u^{0} u^{-1}_{x} \\ u^{1}_{t} =&\left( 2 u^{2} - \left( u^{1} \right) ^{2} \right) u^{0}_{x} - u^{0} u^{1} u^{1}_{x} + u^{0} u^{2}_{x} \\ u^{k}_{t} =&\left( (k+1) u^{k+1} - (k-1) u^{k-1} - u^{1} u^{k} \right) u^{0}_{x} - u^{0} u^{k} u^{1}_{x} + u^{0} u^{k-1}_{x} + u^{0} u^{k+1}_{x}, \ k>1\nonumber \end{aligned}\!\!\!\!\\ \end{aligned}$$
(2.1)

with the notation \(f_{t} = \partial _{t} f\), \(f_{x} = \partial _{x} f\). In particular, we note that the system (2.1) is an infinite chain of quasilinear PDEs of hydrodynamic type. In fact, the equations of the chain are of the form

$$\begin{aligned} u^{k}_{t} = a^{k}_{0} \, u^{0}_{x} +a^{k}_{1} \, u^{1}_{x} +a^{k}_{k-1} \, u^{k-1}_{x} + a^{k}_{k+1}\, u^{k+1}_{x} \end{aligned}$$
(2.2)

or equivalently

$$\begin{aligned} \mathbf{u}_{t} = A(\mathbf{u}) \mathbf{u}_{x} \end{aligned}$$
(2.3)

where \(A (\mathbf{u}) = \{ a^{k}_{j} \}_{j,k=-\infty }^{+\infty }\) is an infinite matrix such that \(a^{k}_{j} = 0\) if \(\notin \{0,1,k-1,k+1 \}\) and

$$\begin{aligned} \begin{aligned} a^{k}_{0}&= {\left\{ \begin{array}{ll} (k+2) u^{k+1} - k u^{k-1} + u^{1} u^{k} &{}if \quad k<0 \\ u^{0} u^{1} &{}if \quad k = 0 \\ (k+1) u^{k+1} - (k-1) u^{k-1} - u^{1} u^{k} &{}if \quad k \ge 1 \\ \end{array}\right. } \\ a^{k}_{1}&= {\left\{ \begin{array}{ll} u^{0} u^{k} &{}if \quad k \le 0 \\ - u^{0} u^{k} &{}if \quad k \ge 1 \end{array}\right. } \\ a^{k}_{k-1}&= {\left\{ \begin{array}{ll} u^{0} &{}if \quad k \ne 1 \\ \left( 2 u^2 - (u^1)^2\right) &{}if \quad k = 1 \end{array}\right. } \qquad \quad a^{k}_{k+1} = {\left\{ \begin{array}{ll} u^{0} &{}if \quad k \ne 0 \\ (u^0)^2 &{}if \quad k = 0 \end{array}\right. } \end{aligned} \end{aligned}$$
(2.4)

By applying the same procedure, one can construct a hierarchy of infinitely many commuting flows, each of them in the form of a hydrodynamic chain (1.11) from the thermodynamic limit of the higher flows of the hierarchy (2.5). The hydrodynamic chain (2.1) is integrable as it possesses an infinite hierarchy of commuting flows. In the following, we show that the hydrodynamic chain (2.1) is diagonalisable and integrable according to the criterion introduced in [15], namely the existence of integrable hydrodynamic reductions in an arbitrary number of components.

Following [15], the diagonalisability of the hydrodynamic chain is established by studying the Haantjes tensor

$$\begin{aligned} H^{i}_{jk} = N^{i}_{pr} \, a^{p}_{j} \, a^{r}_{k} - N^{p}_{jr} \, a^{i}_{p} \, a^{r}_{k} - N^{p}_{rk} \, a^{i}_{p} \, a^{r}_{j} + N^{p}_{jk} \, a^{i}_{r} \, a^{r}_{p} \end{aligned}$$
(2.5)

where \(N^i_{jk}\) is the Nijenhuis tensor

$$\begin{aligned} N^{i}_{jk} = a^{p}_{j} \, \partial _{u^{p}} a^{i}_{k} -a^{p}_{k} \, \partial _{u^{p}} a^{i}_{j} - a^{i}_{p} \left( \partial _{u^{j}} a^{p}_{k} -\partial _{u^{k}} a^{p}_{j} \right) . \end{aligned}$$
(2.6)

In the case of infinite matrices, both Nijenhuis and Haantjes tensors are well defined for the so-called matrices of chain class.

Definition 2

(Chain class matrices [15]) An infinite matrix \(A(\mathbf{u})\) is said to belong to the class C (chain class) if it satisfies the following two properties:

  1. (a)

    each row of \(A(\mathbf{u})\) contains finitely many nonzero elements;

  2. (b)

    each matrix element of \(A(\mathbf{u})\) depends on finitely many variables \(u^{k}\).

Bearing in mind the form of the matrix \(A(\mathbf{u})\) as specified in (2.4), we have the following

Proposition 3

Given the chain (2.1), the associated matrix \(A(\mathbf{u})\) in (2.3) belongs to the chain class.

Moreover, based on the Haantjes theorem given in [16], the following definition extends the concept of diagonalisability to the case of infinite matrices:

Definition 4

(Diagonalisable hydrodynamic chains [15]) A hydrodynamic chain from the class C is said to be diagonalisable if all components of the corresponding Haantjes tensor are zero.

We show that our chain fulfils the definition above.

Proposition 5

Given the chain (2.1), the Haantjes tensor of the associated matrix \(A(\mathbf{u})\) vanishes.

Proof

The proof proceeds by direct inspection. Observing that by definition \(N^i_{jk}\) is antisymmetric under exchange of j and k, i.e. \(N^i_{jk} = - N^i_{kj}\), a direct calculation shows that \(N^0_{jk} = 0\) for any j and k. Similarly, for \(i \ne 0\) the only nonzero elements of \(N^i_{jk}\) are

$$\begin{aligned} N^i_{0\,\pm 1}, \,N^{i}_{0\,i}, \,N^{i}_{0\,i\pm 1}, \,N^i_{1\,i \pm 1} \,N^{i}_{-1\,i\pm 1}, \,N^{i}_{-1\,+1} \end{aligned}$$
(2.7)

and their counterparts with the lower indices exchanged. The above components can be computed for a generic value of i, and their explicit expressions are listed in Appendix D. The structure of \(N^{i}_{jk}\) and \(A(\mathbf{u})\) induces constraints on the range of values the indices p and r can take in the expression of the Haantjes tensor (2.5), and consequently on potential nonzero elements. Indeed, the form of \(N^i_{jk}\), specified by the elements (2.7), implies that for any fixed i the only components of \(H^i_{jk}\) which are not trivially zero are those for \(j, k\in \{0,\pm 1,\pm 2, 3, i, i\pm 1, i\pm 2, i\pm 3\}\). Given the explicit expressions for \(a^k_j\) in (2.4) and \(N^i_{jk}\) in Appendix D, a direct calculation demonstrates that \(H^i_{jk} = 0\) for the listed values of the lower indices. This proves the statement. \(\square \)

We now study the integrability of the chain (2.3) by following the approach based on the method of hydrodynamic reductions applied to the system (2.1). We look for solutions of the form

$$\begin{aligned} u^{k} = u^{k}(R^{1},R^{2},\dots , R^{N}) \end{aligned}$$
(2.8)

for an arbitrary number N of components \(R^{i} = R^{i}(x,t)\). The functions \(\{R^i\}_{i=1}^N\) are the Riemann invariants and satisfy by definition the diagonal system

$$\begin{aligned} R^{i}_{t} = \lambda ^{i}(R^{1},\dots ,R^{N}) R^{i}_{x} \end{aligned}$$
(2.9)

where the characteristic speeds \(\lambda ^{i}\) are such that the system (2.9) possesses the semi-Hamiltonian property, that is

$$\begin{aligned} \partial _{k} \left( \frac{\partial _{j} \lambda ^{i}}{\lambda ^{j} - \lambda ^{i}} \right) = \partial _{j} \left( \frac{\partial _{k} \lambda ^{i}}{\lambda ^{k} - \lambda ^{i}} \right) , \end{aligned}$$
(2.10)

with the notation \(\partial _{i} = \partial _{R^{i}}\). The diagonal form of the system (2.9) and the condition (2.10) guarantee that Eq. (2.9) constitute a system of conservation laws [30] which is integrable via the generalised hodograph method [33]. Substituting the assumption (2.8) into the system (2.3) and using (2.9), we obtain the equations of the form

$$\begin{aligned} \lambda ^{i} \, \partial _{i} \mathbf{u} = A(\mathbf{u}) \, \partial _{i} \mathbf{u}, \qquad i = 1,2,\dots N \end{aligned}$$
(2.11)

where we used the fact that \(R^{i}_{x}\) for \(i=1,\dots ,N\) are independent. We observe that, due to the specific sparse structure of the matrix \(A(\mathbf{u})\), the components of the eigenvectors \(\partial _{i}{\mathbf{u}}\) can be parametrised in terms of the components \(\partial _{i} u^{0}\) and \(\partial _{i} u^{1}\).

Let us consider, for example, the equations for \(\partial _{i} u^{-2}\), \(\partial _{i} u^{-1}\), \(\partial _{i} u^{2}\) and \(\partial _{i} u^{3}\):

$$\begin{aligned} \begin{aligned} \partial _{i}{u^{-2}} =&\, \frac{1}{(u^{0})^{2}} \left( (\lambda ^{i})^{2} - u^{0} u^{1} \lambda ^{i} - u^{0} (2 u^{0} + u^{-2} + u^{-1} u^{1}) \right) \partial _{i} u^{0} - \left( \lambda ^{i} + u^{-1} \right) \partial _{i} u^{1} \\ \partial _{i}{u^{-1}} =&\, \left( \frac{\lambda ^{i}}{u^{0}} - u^{1} \right) \partial _{i} u^{0} - u^{0} \partial _{i} u^{1} \\ \partial _{i}{u^{2}} =&\, \frac{1}{u^{0}} \left( (u^{1})^{2} - 2 u^{2} \right) \partial _{i} u^{0} + \frac{1}{u^{0}} \left( \lambda ^{i} + u^{0} u^{1} \right) \partial _{i} u^{1} \\ \partial _{i}{u^{3}} =&\, \frac{1}{(u^{0})^{2}} \left( \left( (u^{1})^{2} - 2 u^{2} \right) \lambda ^{i} + u^{0} \left( u^{1} (1 + u^{2}) - 3 u^{3} \right) \right) \partial _{i} u^{0} \\&+ \frac{1}{(u^{0})^{2}} \left( (\lambda ^{i})^{2} + u^{0} u^{1} \lambda ^{i} + (u^{0})^{2} (u^{2} -1) \right) \partial _{i} u^{1}\,, \qquad i = 1,\dots , N. \end{aligned} \end{aligned}$$

The compatibility conditions

$$\begin{aligned} \partial _{j}\partial _{i} u^{-2} = \partial _{i}\partial _{j} u^{-2} \qquad \partial _{j}\partial _{i} u^{-1} = \partial _{i}\partial _{j} u^{-1} \qquad \partial _{j}\partial _{i} u^{2} = \partial _{i}\partial _{j} u^{2} \qquad \partial _{j}\partial _{i} u^{3} = \partial _{i}\partial _{j} u^{3} \end{aligned}$$

lead to a so-called Gibbons–Tsarev system. For our chain, this takes the form

$$\begin{aligned} \begin{aligned} \partial _{j} \lambda ^{i}&= \frac{4 (u^{0})^{2} - \lambda ^{i} \lambda ^{j}}{u^{0} (\lambda ^{i} - \lambda ^{j})} \partial _{j} u^{0} \\ \partial _{i} \lambda ^{j}&= \frac{4 (u^{0})^{2} - \lambda ^{i} \lambda ^{j}}{u^{0} (\lambda ^{j} - \lambda ^{i})} \partial _{i} u^{0}\\ \partial _{i} \partial _{j} u^{0}&= \frac{(\lambda ^{i})^{2} + (\lambda ^{j})^{2} - 8 (u^{0})^{2}}{u^{0} (\lambda ^{i} - \lambda ^{j})^{2}} \partial _{i} u^{0} \partial _{j} u^{0} \\ \partial _{i} \partial _{j} u^{1}&= - \frac{(\lambda ^{j} - 2 \lambda ^{i}) \lambda ^{j} + 4 (u^{0})^{2}}{u^{0} (\lambda ^{i} - \lambda ^{j})^{2}} \partial _{i} u^{0} \partial _{j} u^{1} - \frac{(\lambda ^{i} - 2 \lambda ^{j}) \lambda ^{i} + 4 (u^{0})^{2}}{u^{0} (\lambda ^{i} - \lambda ^{j})^{2}} \partial _{j} u^{0} \partial _{i} u^{1}. \end{aligned} \end{aligned}$$
(2.12)

A direct calculation shows that the system of Eq. (2.12) is in involution, i.e. compatibility conditions of the form

$$\begin{aligned} \partial _{k} \partial _{j} \lambda ^{i} = \partial _{j} \partial _{k} \lambda ^{i} \qquad \partial _{k} \partial _{i} \partial _{j} u^{0} =\partial _{i} \partial _{k} \partial _{j} u^{0} \qquad \partial _{k} \partial _{i} \partial _{j} u^{1} =\partial _{i} \partial _{k} \partial _{j} u^{1} \end{aligned}$$

are satisfied modulo the equations (2.12) for all permutations of the derivatives with respect to \(R^{i}\), \(R^{j}\), \(R^{k}\). A first classification of Gibbons–Tsarev systems has been provided by Odesskii and Sokolov [24, 25]. We note that, at the best of our knowledge, the system (2.12) has not appeared before in the literature and it is not included in the class considered in [24, 25].

The compatibility of the Gibbons–Tsarev system (2.12) guarantees that for any solution of the Riemann invariants system (2.1) it is possible to construct a solution of the hydrodynamic chain. This property was proposed in [15] as definition of integrability of a hydrodynamic chain:

Definition 6

(Integrable hydrodynamic chains [15]) A hydrodynamic chain of class C is integrable if it admits N-phase solutions of the form (2.8) for arbitrary N.

Therefore, the above calculations prove the following

Theorem 7

The hydrodynamic chain (2.1) is integrable in the sense of the hydrodynamic reductions.

6 Concluding remarks

Extensive studies of Random Matrix Ensembles and their connection with the theory of integrable systems (see, e.g. [34] and references therein) show that the order parameters, defined as derivatives of the partition function, and their suitable combinations appear as entries of the Lax matrix and the associated Lax equation. For example, in the case of HME, one has the Lax equations for the Toda Lattice. The matrix ensemble of interest is specified by a particular solution of the hierarchy, obtained from a suitable initial condition. Such initial condition is evaluated by considering the partition function and its derivatives in the case where all coupling constants \(t_k\) vanish. Similarly, in the case of SMEs the underlying integrable system is constituted by the Pfaff Lattice and the equations of its hierarchy. These equations specify the behaviour of the order parameters, namely the entries of the Lax matrix, as functions of the coupling constants.

For even power interactions, the thermodynamic limit of the HME is constituted by an order parameter that evolves according to a scalar integrable hierarchy (the Hopf hierarchy) [8]. On the other hand, in the \(n\rightarrow \infty \) limit, the even SME is specified by infinitely many order parameters that satisfy an integrable hydrodynamic chain. This result follows from the key observation that the components of the reduced even Pfaff Lattice can be rearranged in the form of a chain of equations, where the state of each site is given by a vector of infinitely many components. From this point of view, the SME reveals a higher level of complexity compared to the HME due to the existence of integrable reductions in any number of components and associated critical scenarios.

It is indeed well known that, for generic initial conditions, solutions of systems of hydrodynamic type break down in finite time, namely, in the context of SME, for finite values of the coupling constants \(t_{k}\). In the case of the HME, the critical behaviour of the order parameter at the leading order is described by the Whitney cusp and, as observed in [8], finite size corrections resolve the singularity via the onset of a modulated highly oscillating quasi-periodic wave, known as dispersive shock. The dispersive shock characterises a new type of phase transition where asymptotic stable states are connected by an intermediate state where order parameters develop fast oscillations induced by the dispersive nature of finite size corrections. The case of SME presents a similar scenario with a potentially richer variety of behaviours due to the higher number of components. Further studies in this direction will entail the detailed analysis of the solution with the specific initial condition induced by the partition function (1.9) calculated at \({\mathbf {t}} ={\mathbf {0}}\).

The application of the Haantjes tensor test and the method of hydrodynamic reductions allow to prove the integrability of the hydrodynamic chain by considering the first nontrivial flow only. Integrability implies the existence of infinitely many commuting flows which describe the evolution of the order parameters in the space of coupling constants.

We finally note that above considerations are concerned with a direct comparison between HME and SME when restricted to even power interactions. It is important to note that, with the given scaling, the hydrodynamic chain arises in the case of even power interactions only. As mentioned earlier and further specified in Appendix B, the first flow (B.1) associated with \(t_{1}\) does not lead to an infinite chain of quasilinear PDEs. The system (B.1) and its relation with the large n scaling properties of the initial condition will be analysed in detail in a separate work. The previously unseen connection between matrix ensembles and hierarchies of hydrodynamic chains discussed in this paper, together with the aforementioned results for the HME discussed in [8], suggests that the study of random matrix models may lead to the discovery of new interesting integrable hydrodynamic PDEs. The study of the PDEs so obtained arises as a general framework and a new methodology to classify and describe asymptotic properties of complex systems.