Thermodynamic Limit and Dispersive Regularisation in Matrix Models

TL;DR

This paper investigates phase transitions in Hermitian matrix models, revealing a dispersive regularisation mechanism that resolves singularities and leads to complex, chaotic behaviors, with implications for understanding nonlinear matrix dynamics.

Contribution

It introduces a novel dispersive regularisation framework for phase transitions in matrix models, linking them to Toda hierarchy reductions and dispersive shocks.

Findings

Dispersive regularisation resolves singularities at critical points.

Chaotic behaviors emerge from dispersive shocks in matrix models.

The mechanism applies to nonlinearities of arbitrary order.

Abstract

We show that Hermitian matrix models support the occurrence of a new type of phase transition characterised by dispersive regularisation of the order parameter near the critical point. Using the identification of the partition function with a solution of a reduction of the Toda hierarchy, known as Volterra system, we argue that the singularity is resolved by the onset of a multi-dimensional dispersive shock of the order parameter in the space of coupling constants. This analysis explains the origin and mechanism leading to the emergence of chaotic behaviours observed in M^6 matrix models and extends its validity to even nonlinearity of arbitrary order.

Equations

$$Z_n(\mathbf{t})={\int }_{{\mathcal{H}}_n}{\mathrm{e}}^{H(M)}dM,$$

$$H(M)=\mathrm{Tr}\left(-M^2/2+\sum _{j=1}^{\infty }{t}_{2j}{M}^{2j}\right)$$

$${\tau }_n(\mathbf{t})=\frac{1}{n!}{\int }_{{\mathbb{R}}^n}{\Delta }_n(\mathbf{\lambda }{)}^2\prod _{i=1}^n\left({\mathrm{e}}^{H({\lambda }_i)}d{\lambda }_i\right)$$

$$\frac{\partial L}{\partial t_{2k}}=\left[\frac{1}{2}{\left(L^{2k}\right)}_s,L\right]k=1,2,\dots .$$

$$L=\left(\begin{array}{ccccc}0& b_1& 0& 0& \dots \\ b_1& 0& b_2& 0& \dots \\ 0& b_2& 0& b_3& \dots \\ ⋮& \ddots & \ddots & \ddots & \ddots \end{array}\right)$$

$$\frac{\partial b_{n}}{\partial t_{2k}}=\frac{b_{n}}{2}\biggl{(}b_{n+1}(L^{2k-1})_{n+1,n+2}-b_{n-1}(L^{2k-1})_{n-1,n}\biggr{)}.$$

$$\frac{\partial B_n}{\partial t_{2k}}=B_n(V_{n+1}^{(2k)}-V_{n-1}^{(2k)})k=1,2,\dots$$

$$V_n^{(2)}$$

$$V_n^{(4)}$$

$$V_n^{(6)}$$

$$n=B_n-\sum _{j=1}^{\infty }2j{t}_{2j}{V}_n^{(2j)}.$$

$$\left[L,P\right]=1$$

$$P=-\frac{1}{2}{L}_s+\sum _{k\ge 1}k{t}_{2k}{\left(L\right)}_s^{2k-1}$$

$$\frac{n}{N}=u_n-\sum _{j=1}^{\infty }2j{T}_{2j}{W}_n^{2j}$$

$${\Omega }_{ϵ}=0$$

$$\begin{array}{rl}{\Omega }_{ϵ}:=& -x+(1-2T_2)u-12T_4{u}^2-60T_6{u}^3\\ & +{ϵ}^2\left(p_1{u}_{xx}+p_2{u}_x^2\right)\\ & +{ϵ}^4\left(q_1{u}_{xxxx}+q_2{u}_x{u}_{xxx}+q_3{u}_{xx}^2\right)+O({ϵ}^6),\end{array}$$

$$p_1=-4T_{4}u-60T_6{u}^2$$

$$q_1=-\frac{T_4}{3}u-11T_6{u}^2$$

$$q_3=-\frac{33T_6}{2}u.$$

$${\Omega }_0:=-x+(1-2T_2)u-12T_4{u}^2-60T_6{u}^3=0.$$

$$u_{T_{2k}}=\sum _{n=0}^{\infty }{ϵ}^n{g}_n^{(k)}\left(u;u_x,\dots ,{\partial }_x^{n}u\right)$$

$$u_{T_2}=2u\left[\frac{1}{ϵ}\sinh \left(ϵ{\partial }_x\right)\right]u$$

$$u_{T_2}=2u{u}_x.$$

$$u_{T_{2k}}=C_k{u}^k{u}_x$$

$$f_0[u]=-xu+\frac{1}{2}\left(1-2T_2\right)u^2-4T_4{u}^3-15T_6{u}^4.$$

Full text

I Introduction

Random matrix models, originally introduced as an attempt to model the complex structure of heavy nuclei energy spectra, have since become a universal paradigm for modelling complex phenomena. They naturally arise in connection with different areas of mathematics and physics, from quantum field theory to the theory of completely integrable dynamical systems wigner ; dysonI ; mehta .

For instance, in the context of quantum field theory, Hermitian matrix models can be introduced as a discrete approximation of 2D quantum gravity. Based on this result, a celebrated conjecture of Witten witten , proven by Kontsevich konts , establishes the connection between the free energy of 2D quantum gravity and a particular solution of the Korteweg-de Vries (KdV) equation. The KdV equation, the prototypical example of soliton equation ablowitz , first appeared in water wave dynamics to model small amplitude elevation waves in shallow water. Along with solitons, another important class of solutions of the KdV equation is constituted by Dispersive Shock Waves (DSW) whitham ; gurevich-pitaevskiy . DSWs occur as a universal regularisation mechanism of singularities in dispersive hydrodynamics, and effectively explain the emergence of a variety of complex behaviours in hydrodynamic systems whitham ; gurevich-pitaevskiy ; el ; eh ; ehs ; congy . The remarkable mathematical structure of completely integrable systems, of which the KdV equation is possibly the most celebrated example, such as the existence of an infinite number of conservation laws, allows for an effective and detailed description of dispersive shock waves (see e.g. gurevich-pitaevskiy ; el ; eh ; ehs ). The above mentioned Witten’s conjecture established an unexpected connection between quantum field theory and dispersive hydrodynamics via the identification of correlators of 2D gravity with conserved densities of a particular solution of the KdV equation. Thereafter, similar correspondence between different matrix models on Hermitian, Unitary, Symmetric and Symplectic ensembles and completely integrable dynamical systems have been discovered (see e.g. gerasimov ; adler ; adler2 ; okounkov ; takasaki and references therein). The relevance of these connections lies in the fact that methods developed for solving completely integrable dynamical systems are indeed effective for the study of the associated matrix models and unveil further intriguing mathematical structure and elegance. Moreover, extensive studies of properties of partition functions for matrix models led to remarkable connections between the theory of integrable systems, statistical mechanics, quantum field theory, algebraic and enumerative geometry witten ; adler ; okounkov ; takasaki ; dubrovin ; dubrovin2 . In this paper we investigate Hermitian matrix models and observe a novel type of phase transition, resulting from a dispersive regularisation mechanism of the order parameter in the space of coupling constants and the consequent onset of a dispersive shock. We provide the physical conditions and scaling regime under which this phenomenon occurs. In particular, we find that, unlike classical mean field models, multivalued solutions at the leading order of the asymptotic expansion do not necessarily correspond to a phase transition. For the sake of simplicity, we focus on a case of Hermitian Matrix Model with even nonlinear interaction terms and and its formulation in terms of the 1D Toda hierarchy. However, our considerations can be extended to other matrix ensembles. The Toda Lattice is an important example of completely integrable nonlinear dynamical system which, in the continuum limit, contains as a particular case several examples of soliton equations, including the KdV equation. We also note that asymptotic properties of partition functions in the thermodynamic limit of one-matrix models with even and odd nonlinearity and their relation with the Toda Lattice have been previously considered in ercolani1 ; ercolani2 . A key point in our analysis is that the sequence of partition functions $Z_{n}$ for the one-matrix model of $n\times n$ matrices can be expressed in terms of a particular solution of a suitable restriction of the Toda Lattice equations, known as Volterra Lattice (or discrete KdV equation). The complete integrability of the Volterra Lattice system implies the existence of infinitely many conservation laws. The Volterra Lattice together with its set of symmetries associated to the conservation laws constitutes the Volterra hierarchy. We show that the identification of the Volterra hierarchy with the matrix model is based on a one to one correspondence between coupling constants and equations of the hierarchy. Each equation of the hierarchy provides the evolution of the order parameter of the theory as a function of the associated coupling constant, that is identified with the time variable of the chosen equation. The partition function $Z_{n}$ is therefore specified by the state of the $n-$th point of the lattice for the corresponding values of the coupling constants (i.e. time parameters of the hierarchy). Most importantly, the dynamics on the lattice is uniquely specified by the initial conditions that are given in terms of the partition function of the free model, i.e. where all coupling constants vanish. In this respect, the model is relatively simpler than the case of 2D gravity studied in witten where the initial condition is specified by additional symmetries that are compatible with the hierarchy, namely the Virasoro constraints dubrovin2 ; adler .

We thus exploit the relation between the matrix model under consideration and the Volterra Lattice in order to explore the phase diagram of the model. In his pioneering work jurkiewicz2 , Jurkiewicz observed that a natural order parameter can be introduced by using orthogonal polynomial decompositions and combinatorial considerations ParisiEtAl . Such order parameter develops, in the thermodynamic limit, a singularity that is regularised by oscillations revealing an apparent underlying chaotic behaviour of the system jurkiewicz2 ; senechal . Rigorous proof of the occurrence of asymptotic oscillations of the partition function has been found in ken ; grava . In this work we explain the oscillatory behaviour of the order parameter as the occurrence of a dispersive shock. Using the fact that the thermodynamic limit ($n\rightarrow\infty$) of the random matrix model is equivalent to the continuum limit of the Volterra Lattice, given by a system of partial differential equations of hydrodynamic type deift ; kodama ; michor , we show that the order parameter evolves in the space of the coupling constants as a shock wave solution of the associated hydrodynamic system. The chaotic phase is therefore interpreted as a dispersive shock propagating through the chain in the continuum/thermodynamic limit. The intriguing complexity of its phase diagram can hence be explained in the context of dispersive hydrodynamics. The physical constraint on the signature of the order parameter determines whether a catastrophe evolves into a dispersive shock.

Considerations above outline the following general scenario: when a thermodynamic system undergoes a phase transition, some specific quantities, the order parameters, develop singularities. In the context of conservation laws of hydrodynamic type, when a singularity (hydrodynamic catastrophe) occurs, viscosity and dispersion underpin two different mechanisms of regularisation of such singularity. In presence of small viscosity the solution develops a sharp but smooth wavefront whitham . If small viscosity is replaced by small dispersion, when the wavefront approaches the point of gradient catastrophe the dispersion induces initially small oscillations that further evolve into a dispersive shock gurevich-pitaevskiy ; el ; ehs . In classical mean field fluid and spin models phase transitions are associated to classical shocks of order parameters in the space of thermodynamic parameters moro1 ; barra ; genovese . In this work we show that the chaotic behaviour observed in jurkiewicz2 is indeed a phase transition where the order parameter develops a singularity that is resolved via dispersion rather than viscosity as in classical spin models. The observation that phase transitions in matrix models are explicitly connected to the occurrence of a dispersive shock in the order parameter paves the way to the application of current methods of dispersive hydrodynamics in the context of quantum field theories.

II Hermitian Matrix Model

We study the model defined by the partition function

[TABLE]

where $M$ are Hermitian matrices of order $n$,

[TABLE]

is the Hamiltonian, with $\mathbf{t}=\{t_{2j}\}_{j\geq 1}$ the coupling constants, and $dM$ is the Haar measure in the space of Hermitian Matrices ${{\mathcal{H}}_{n}}$. Based on a classical result by Weyl weyl , the partition function (1) is proportional to an integral over the eigenvalues of the matrix $M$, that is $Z_{n}(\mathbf{t})=c_{n}\tau_{n}(\mathbf{t})$ where $c_{n}$ is a constant and

[TABLE]

where $\Delta_{n}(\boldsymbol{\lambda})=\prod_{1\leq i<j\leq n}(\lambda_{i}-\lambda_{j})$ is the Vandermonde determinant. A theorem by Adler and van Moerbeke adler implies that the quantity (2) can be interpreted as a tau-function of the Toda hierarchy restricted to the even flows

[TABLE]

with $L$ the tridiagonal symmetric Lax matrix of the form

[TABLE]

where $b_{i}=\sqrt{\tau_{i+1}\tau_{i-1}/\tau_{i}^{2}}$ and $\left(L^{2k}\right)_{s}$ denotes the skew-symmetric part of the matrix $L^{2k}$ (see e.g. adler ). The solution of interest is specified by the initial conditions $b_{i}(\mathbf{0})=\sqrt{n}$ obtained via a direct calculation of Gaussian integrals for the quantities $\tau_{n}(\mathbf{0})=(2\pi)^{n/2}\prod_{j=1}^{n}j!/n!\;$. We note that the Lax matrix of the type (4), originally considered in kac , and more recently in ercolani1 , corresponds to a reduction of the even Toda hierarchy known as Volterra hierarchy. The description of the Volterra hierarchy in terms of the matrix resolvent, its $\tau$ structure and application to ribbon graphs and Hodge integrals has been recently studied in Yang . Incidentally, we also mention that the model with odd nonlinearities is different from the present case and its relation with the Toda hierarchy has been considered in ercolani2 .

Writing the equations of the hierarchy (3) in terms of lattice variables $b_{n}$ we have

[TABLE]

Introducing the notation $B_{n}=b_{n}^{2}$ and $V^{(2k)}_{n}=b_{n}(L^{2k-1})_{n,n+1}$ and multiplying both sides by $b_{n}$, the above equation (5) reads as follows

[TABLE]

One can simply prove by induction that $V^{(2k)}_{n}$ is a linear combination of products of the variable $B_{n}$ evaluated at different points on the lattice. For instance, for the first three non-trivial flows we have

[TABLE]

In the following we refer to $B_{n}$ as the order parameter of the theory.

Based on a result in bonora , we have proven that the required solution to the above reduction of the even Toda hierarchy is given by the recursive formula (string equation)

[TABLE]

Indeed, Eq. (7) follows from the the string equation for the Toda lattice

[TABLE]

where

[TABLE]

by considering its restriction to even times only bonora . Eq. (7) allows to evaluate the order parameter of the $M^{2q}$ model for arbitrary $q$ and generalises the formula obtained by Jurkiewicz for $q=3$ jurkiewicz ; jurkiewicz2 .

III Thermodynamic limit

We analyse the Matrix Model in the large $n$ (thermodynamic) regime via the continuum limit of the solution (7) of the reduced Toda hierarchy. Introducing the scale given by a suitable large integer $N$ and the rescaled variables $u_{n}=B_{n}/N$, $T_{2k}=N^{k-1}t_{2k}$, Eq. (7) reads as follows

[TABLE]

where $W_{n}^{2j}=V_{n}^{(2j)}/N^{j}$. We then define the interpolating function $u(x)$ such that $u_{n}=u(x)$ for $x=n/N$ and $u_{n\pm 1}=u({x\pm\epsilon})$ with the notation $\epsilon=1/N$. Expanding in Taylor series for small $\epsilon$, we obtain an ODE as a formal series

[TABLE]

where $\Omega_{\epsilon}$ has the following form

[TABLE]

Here for simplicity we have fixed $T_{2k}=0$ for $k>3$, with $p_{i}$ and $q_{i}$ as follows.

[TABLE]

At the leading order we have the polynomial equation in $u$ of the form

[TABLE]

In order to obtain further insights on the behaviour of the solution of the recurrence equation (9) in the thermodynamics limit, it is interesting to study the continuum limit of the Volterra hierarchy (6). Proceeding as above, i.e. introducing the interpolating function $u(x)$ in (6) and expanding $u(x\pm\epsilon)$ in Taylor series, one obtains a compatible hierarchy of dispersive partial differential equations of the form

[TABLE]

where functions $g^{(k)}_{n}$ are differential polynomials of $u$. For instance, the first member of the hierarchy (for $k=1$), which gives the flow with respect to $T_{2}$, takes the following compact form

[TABLE]

where the operator stays for the formal McLaurin expansion of $\sinh$. In the thermodynamic limit, i.e. $\epsilon\to 0$, equation  (14) gives, at the leading order, the Hopf equation

[TABLE]

Similarly, higher flows in $T_{2k}$ lead to higher members of the so-called Hopf hierarchy

[TABLE]

for suitable constants $C_{k}$, the solution of which, for the assigned initial condition, is implicitly given by the polynomial equation (12). The effect of the corrections in the parameter $\epsilon$ will be discussed in the following section.

For illustrative purposes, we consider the case $t_{2k}=T_{2k}=0$ for $k>3$, i.e. Eq. (12) is a cubic equation for the order parameter $u$. In this case, Eq. (12) provides the condition for extremising the free energy functional $F=\int_{0}^{\beta}f_{0}(u)\;dx$, for some $\beta>0$, of density

[TABLE]

In particular, local minima and maxima depend on the signature of the discriminant $\Delta(x,T_{2},T_{4},T_{6})$ of the cubic equation (12). If $\Delta>0$ the free energy has two local minima and one local maximum, if $\Delta<0$ the free energy presents one local minimum only. The set in the space of parameters such that $\Delta=0$ corresponds to the critical set where a phase transition occurs. For example, in Fig. 1 we plot the set $\Delta=0$ in the $x$-$T_{6}$ plane for a given choice of $T_{2}$ and $T_{4}$ (notice that in order to ensure convergence of the integral (1), we have $T_{6}<0$). The convex sector corresponds to the region where the equation of state (12) admits three real solutions that correspond to the stationary points of the free energy density. Fig. 2 shows the free energy density as a function of $u$ for two different values of $T_{6}$, with fixed $T_{2}$, $T_{4}$ and $x$. For these values of the parameters the discriminant of the cubic is positive, and the free energy has three extremal points. We see that for $T_{6}=-0.0051$ the difference of the values of the free energy density at its local minima is particularly pronounced compared with the case $T_{6}=-0.0067$.

The above scenario is compatible with the well known fact that the generic solution of the Hopf hierarchy (16) develops singular behaviour for finite value of the “time” variables $T_{2k}$ whitham . In the next section, we study these singularities in relation to the occurrence of phase transitions. A phase transition is associated to the occurrence of a dispersive shock induced by dispersive corrections to the Hopf hierachy. Eq. (10) obtained as the continuum limit of (9) provides quasi-trivial deformations of the Hopf hierarchy and the behaviour near the singularity that are universally described by a solution of the fourth order analogue of the Painlevé I equation dubrovin3a ; dubrovin3 ; dgkm .

IV Dispersive regularisation

In the thermodynamic limit, the evolution of the order parameter in the space of coupling constants is governed by Eq. (16) provided derivatives of $u$ are bounded. In the vicinity of the gradient catastrophe, dispersive corrections in the equation (13) induce fast oscillations responsible for a rich and interesting phenomenon known as dispersive shock wave gurevich-pitaevskiy ; eh .

In this section we illustrate the general phenomenology by considering the particular case $T_{2k}$ $=$ $t_{2k}$ $=0$ for all $k>3$ so that $T_{2}$, $T_{4}$ and $T_{6}$ are the only non zero coupling constants. This choice allows for a simple but sufficiently general analysis demonstrating as chaotic behaviours observed in jurkiewicz correspond to a type of phase transition comprised by a dispersive shock of the order parameter. The shock occurs as a dispersive regularisation mechanism of a particular solution of the hierarchy (6) in the continuum limit.

In Fig. 3(b) we compare the order parameter $u(x)$ obtained as solution of the recurrence equation (9) and the limit equation (12). Values $T_{2}$, $T_{4}$ and $T_{6}$ are chosen in such a way that the solution of the cubic equation (12) is single valued. Fig. 3(a) shows that the two solutions fully overlap for sufficiently small $\epsilon$, but, as shown in Fig. 3(b), a relevant deviation is observed in the vicinity of the point of gradient catastrophe of the solution to Eq. (12).

Figs. 4(a) and 4(b) show a comparison between the cubic solution (12) and the exact solution (9) for different values of $T_{6}$ within the convex region shown in Fig. 1, where the solution of (12) is multivalued. Both figures demonstrate the onset of a DSW. This behaviour is intriguing as, unlike classical statistical mechanical systems, e.g. magnetic and fluid models baxter , the order parameter $u(x)$ develops oscillations in the form of a dispersive shock in conjunction with the existence of additional stationary points for the free energy such as unstable and metastable states. The emergence of such oscillations is therefore explained as a result of higher order corrections in the string equation (10) of which the cubic equation (12) is the leading order approximation. Based on the results in dubrovin3a ; claeys ; dubrovin3 ; dgkm , the mechanism of formation of such oscillations in the vicinity of the of the critical points of the solution to the equation (12) is universal and it is given by a particular real analytic solution of the second member of the Painlevé I hierarchy, known as Painlevé I2 equation.

V Analysis of scenarios

As discussed above, a dispersive shock in the order parameter develops for values of the couplings such that the density of free energy admits multiple stationary points. Hence, the signature of the discriminant of the cubic in Eq. (12) determines a necessary condition for the occurrence of a dispersive shock and therefore a phase transition. We now focus on the different scenarios in regimes where the discriminant $\Delta$ is positive and the Eq. (12) admits three real and distinct solutions. In particular, different cases need to be considered depending on whether the coefficients of the cubic equation (12) are negative or positive. Necessarily, in order to ensure convergence of the integral (1), it is $T_{6}<0$. Hence, we have four distinct cases, depending on the signs of the coefficient $1-2T_{2}$ and $-12T_{4}$ in Eq. (12).

Scenario 1: $T_{2}<1/2$ and $T_{4}>0$.

This choice corresponds to the case analysed in jurkiewicz ; jurkiewicz2 , hence it allows for a direct comparison. As by definition $u(x)\geq 0$, only non negative branches of $u(x)$ correspond to admissible states of the system. In fact, in both Figs. 4(a) and 4(b), the three branches of the cubic which correspond to stationary points of the free energy are positive.

**Scenario 2: $T_{2}>1/2$ and $T_{4}<0$. **

In this case, as illustrated for example in Fig. 5, the solution to Eq. (12), shown in Fig. 5(a), it is three-valued but one root, associated to a local minimum, is negative and therefore does not correspond to a state of the system. However, two concurrent states, although of different nature, one stable and one unstable, underlie a dispersive shock, shown in Fig. 5(b). Notice that for $x>0$ the solution to Eq. (12) has one non negative branch only. Nonetheless, $u(x)$ develops a dispersive shock profile at positive $x$, although this is originated by a catastrophe located at $x<0$. The solution to Eq. (12) is multivalued with two non negative branches for a small interval of negative values of $x$.

Scenario 3: $T_{2}>1/2$ and $T_{4}>0$. Similarly to Case 2, $u(x)$ has only one positive solution for $x>0$ (see Fig. 6(a)) and as shown in Fig. 6(b), a dispersive shock arises in correspondence of two non negative roots for negative $x$.

Scenario 4: $T_{2}<1/2$ and $T_{4}<0$. The last scenario is given by the case, illustrated in Fig. 7(a), where the cubic solution is multivalued with one positive branch for all values of $x$ and therefore only one solution corresponds to a state that is accessible by the system. Interestingly, as shown in Fig. 7(b), the solution of the recurrence equation (9) overlaps with the cubic solution and no oscillations occur. This suggests that the dispersive regularisation in the form of a dispersive shock of the order parameter is related to the existence of accessible (meta-)stable/unstable states. This phenomenon is rather unexpected as, from the point of view of the governing dispersive PDE, one would expect that for a generic initial condition multivaluedness would be replaced by a dispersive shock originated by the gradient catastrophe.

VI Conclusions

We have analysed the Hermitian matrix model with even degree of nonlinearity for which the order parameter of the model is obtained as a particular solution of the Volterra hierarchy. The critical properties of the model are obtained in terms of a solution to the continuum limit of Volterra hierarchy. We have shown that the chaotic behaviour observed in previous literature can be described as a dispersive shock propagating in the space of thermodynamic parameters. Also, the profile of the order parameter, specifically the form of the envelope, appears to be highly sensitive to the choice of the parameters $T_{2k}$. For instance, Fig.s 4(a), 4(b), 5(b) and 6(b) show the onset of dispersive shocks whose envelope displays very distinctive features.

A further detailed study of this intriguing behaviour will involve the construction of the asymptotic genus expansion of the solution (9) and the Whitham modulation theory for solutions of related integrable hierarchy. We finally anticipate that the rich phenomenology described in this paper reflects the fact that the dispersive shock given by the solution (9) is an intrinsic multidimensional object arising from the simultaneous solution of equations of the hierarchy (6) in the continuum limit. In fact, the solution to the Volterra hierarchy (6) satisfies the modified Kadomtsev-Petviashvili hierarchy similarly to how the solution of the Toda hierarchy satisfies the Kadomtsev-Petviashvili hierarchy. Hence, the description of the dispersive shock solution entails the development of the Whihtam modulation theory of a $2+1$-dimensional integrable dispersive equation. This further development is in progress and will be presented elsewhere BBMP .

Acknowledgements. This work is dedicated to the memory of Boris Dubrovin (1950-2019) whose magnificent scientific foresight and generosity have inspired the original ideas of this research. Authors are grateful to The Leverhulme Trust RPG-2017-228 Integrable Dressed Networks for supporting this research project.