UNIVERSAL SPECTRAL SHOCKS IN RANDOM MATRIX THEORY — LESSONS FOR QCD ∗

Following Dyson, we treat the eigenvalues of a random matrix as a sys-tem of particles undergoing random walks. The dynamics of large matrices is then well-described by ﬂuid dynamical equations. In particular, the inviscid Burgers’ equation is ubiquitous and controls the behavior of the spectral density of large matrices. The solutions to this equation exhibit shocks that we interpret as the edges of the spectrum of eigenvalues. Going beyond the large N limit, we show that the average characteristic polynomial (or the average of the inverse characteristic polynomial) obeys equations that are equivalent to a viscid Burgers’ equation, or equivalently a diﬀusion equation, with 1 /N playing the role of the viscosity and encoding the entire ﬁnite N eﬀects. This approach allows us to recover in an elementary way many results concerning the universal behavior of random matrix theories and to look at QCD spectral features from a new perspective.

(Received January 31, 2015) Following Dyson, we treat the eigenvalues of a random matrix as a system of particles undergoing random walks. The dynamics of large matrices is then well-described by fluid dynamical equations. In particular, the inviscid Burgers' equation is ubiquitous and controls the behavior of the spectral density of large matrices. The solutions to this equation exhibit shocks that we interpret as the edges of the spectrum of eigenvalues. Going beyond the large N limit, we show that the average characteristic polynomial (or the average of the inverse characteristic polynomial) obeys equations that are equivalent to a viscid Burgers' equation, or equivalently a diffusion equation, with 1/N playing the role of the viscosity and encoding the entire finite N effects. This approach allows us to recover in an elementary way many results concerning the universal behavior of random matrix theories and to look at QCD spectral features from a new perspective.

Introduction
The physics motivations of (some of) the authors of this contribution have their roots in the study of various aspects of Quantum Chromodynamics (QCD). A particularly inspiring example is the case of Yang-Mills theory in two dimensions and the order-disorder phase transition first identified by Durhuus and Olesen, or the chiral symmetry breaking and its restoration at finite temperature. There were also attempts to apply the Random Matrix Theory to QCD evolution equations at high energy, but so far these attempts have not been successful. In all cases considered, some aspects of the physics are captured by the Random Matrix Theory. This is so, in particular, in situations where, as a parameter is varied, qualitative changes of behavior present a universal behavior. The parameter may be the area of a Wilson loop, the volume of the system, the number of colors, the rapidity, etc. We shall generically refer to this parameter as a "time", and, indeed, our main effort will be to arrive at dynamical equations that describe the evolution of the system. These equations will follow Dyson's prescription [1], that: The x i should be interpreted as positions of particles in Brownian motion. This means that the particles do not have well-defined velocities, nor do they possess inertia. Instead, they feel frictional forces resisting their motion.

Burgers' equation
Much of the work to be discussed in this note was inspired by the ubiquitous emergence of Burgers' equation in the Random Matrix Theory. We shall, therefore, start by recalling some properties of Burgers' equation [2]. This will give us acquaintance with concepts that will be used later.
Burgers' equation describes the evolution of the velocity field u(x, t) of a fluid. In one dimension, it reads where ν is the viscosity. This equation differs from the Navier-Stokes equation by the absence of pressure forces (in the Navier-Stokes equation, a term ∇P/ρ, with P the pressure and ρ the density, would be present on the righthand side of Eq. (1)). The non-linear term in Eq. (1) can be viewed as of the purely "kinematical" origin. When ν = 0, the equation is referred to as the inviscid Burgers equation. It is often used as a one-dimensional model of turbulence. Its solution can exhibit shocks, and this will play an important role in our discussion. The inviscid Burgers equation can be solved by the method of characteristics: We look for solutions of the form of u(x, t) = u(x(t), t), such that u is constant, du = 0, along the characteristic lines The solution of the equation can then be written as u(x, t) = u 0 (ξ(x, t)), or equivalently as the solution of the implicit equation u(x, t) = u 0 (x−tu(x, t)). A shock starts to develop when the characteristics cross each other, which happens for the smallest value t * of t for which the condition is satisfied. A generic initial velocity field for which this happens is u 0 (ξ) = aξ + bξ 3 , with a < 0 and b > 0. Then, it is easily shown that t * = −1/a. As we shall see, we can learn a lot about the general (universal) structure of the solution by analyzing its behavior in the vicinity of the shocks. In order to solve the viscid Burgers equation, it is convenient to perform a so-called Cole-Hopf transform with φ(x, t) obeying a diffusion equation The solution can be then written as a convolution of the heat kernel and the initial condition u 0 (x) When ν → 0, one can evaluate the integral using the saddle point method.
The saddle point equation is nothing but the characteristic equation, as one can easily verify: we recover the inviscid Burgers equation.
Most of these elements (characteristics and shocks, Cole-Hopf transform, heat kernel, etc.) will be present in the foregoing analysis of specific random matrix models. In particular, we shall identify the shocks in the Burgers equation with the edges of the spectrum of eigenvalues, and the viscosity will be seen to be simply related to the size N of the matrices, ν = 1/2N . However, while we found the fluid analogy inspiring, a word of caution is in order: we should emphasize that the Burgers equations that we shall meet will involve complex valued fields and coordinates, and the viscosity may turn out to be sometimes negative.

Yang-Mills theory in d = 2
The first example that we shall discuss is that of the Yang-Mills theory in d = 2 dimensions. This is not the simplest example, but it is the one for which the main observations on which our work is based were made [3]. The fundamental object to consider in this context is the average of the Wilson loop along a (simple) curve C where P is a path ordering operator, and the average (denoted by the angular brake · ) is taken over the two-dimensional Yang-Mills field configurations. The gauge field A µ is taken in the fundamental representation so that W is an N c ×N c matrix, with N c the number of colors. To within a normalization, W depends only on the area enclosed by the loop. The matrix W is unitary (for each realization of the gauge field), with its eigenvalues of the form of λ = e iθ , living on the unit circle. The density of these eigenvalues, ρ(θ), evolves as a function of the size of the loop. For small loops (which probe short distance, perturbative physics), the spectrum covers a small fraction of the unit circle around θ = 0: the spectrum exhibits a gap. As the area increases, the gap eventually closes, with the eigenvalues occupying uniformly the unit circle in the limit of large areas. Remarkably, the region of the crossover where the closing of the gap takes place, becomes infinitely thin as N c , the number of colors (the size of the matrix) becomes infinite, suggesting a phase transition at infinite N c , for a particular critical loop area A = A * [4]. This behavior can be reproduced by a simple random matrix model, in which the building of Wilson loops of increasing areas is mimicked by multiplying random unitary matrices [5] where H i is a random Hermitian matrix drawn from a Gaussian distribution with second cumulant 1 N TrH 2 = m 2 . The square root dependence on time reflects the underlying diffusion process. This model exhibits a phase transition for a critical area A * ∼ 1/m 2 .
It turns out that this phenomenon can be described by a Burgers equation. This equation has been derived in a number of independent ways, none of them very direct (see references in [3]). It reads In this equation, A plays the role of time, and the angular variable α the role of the position (on the unit circle). The function F (which plays the role of the velocity field) is related to the resolvent G(z) The density ρ(θ) can be obtained from the imaginary part of F (z = e iα ). This relation between the resolvent G and the function F that satisfies the Burgers equation is complicated here by the structure of the unitary group and its compact support. We shall later treat examples where it is the resolvent itself that satisfies the (inviscid) Burgers equation. The main issue here is to observe that F is a complex function of a complex variable. By solving this equation using the method of (complex) characteristics (see [6]), one finds that for the trivial initial condition corresponding to the unit matrix, i.e. G 0 (z) = 1/(z −1), or F (α) = (1/2) cot α/2, shocks develop. The characteristics are straight lines α = ξ + AF 0 (ξ), where F 0 (α) = F (A = 0, α), and the locations of the shocks are given by the solutions of dα dξ ξc = 1 + AF 0 (ξ c ) = 0. By expanding F 0 (ξ) in the vicinity of ξ c , one gets This equation allows us to invert the characteristic equation and get ξ(α, A), and hence the solution in the vicinity of the singularity (recall that This is the familiar square root singularity corresponding to the edges of the spectrum in the gapped phase. At the closure of the gap, when the two edges of the spectrum on the unit circle meet at α = π, F 0 (ξ c ) = 0, and the cubic term is dominant, leading to a cubic root singular behavior in the spectral density. This point corresponds to the Durhuus-Olesen transition. In the vicinity of this point, the spectral density is of the form of ρ(α) ∼ (π − α) 1/3 . This completes the analysis of the spectral density in the large N limit.
In order to go beyond the infinite N limit, and capture the effect of finite N corrections of the spectrum of eigenvalues, it is convenient to consider the average characteristic polynomial This object was shown to be convenient in order to obtain spectral information from numerical simulations [7]. In the large N limit, the average characteristic polynomial is simply related to the resolvent. We have indeed (with λ denoting an eigenvalue of W ) and in the large N limit, we have approximately ln det(z−W ) ≈ ln det(z− W ) . The function F that satisfies the Burgers equation is related to Q N (z) by (see Eq. (10)) For the unitary group, the average characteristic polynomial can be calculated via a character expansion, and an integral representation has been given by Neuberger [8,9] with Expression (15) has the form of a convolution of the heat kernel and some initial condition (it is easily verified that this is indeed the initial condition that was mentioned earlier). Thus, q N satisfies a diffusion equation with a diffusion constant equal to 1/2N , while −(1/N )∂ ln q N /∂y satisfies the corresponding viscid Burgers equation. In terms of the function F introduced above, and for z = e iα , this equation reads Note that in this case, the viscosity is negative. One should be careful, however, with the interpretation of the viscosity term in this complex setting, since the sign of the viscosity term depends on the direction of the complex plane in which one is looking (it is positive, for instance, if one choose z = −e y , with y real). One may argue that the negative sign is indeed what is expected in order to amplify the spectral oscillations of the eigenvalue density that are observed in the vicinity of the shocks, while the positive viscosity would lead naturally to damping. The viscosity term in Eq. (17) allows us to take into account the finite N corrections in a compact way, and to zoom into the structure of the spectrum in the close vicinity of the shocks. As emphasized already, the viscid Burgers equation controls the behavior, not of the spectral density itself, but of the average characteristic polynomial (or a function closely related to it), that evolves smoothly towards the resolvent as N → ∞. As is well-known, to capture the universal behavior of the spectrum close to its edges, one needs to introduce special scalings that take into account the way the spectral density at large N drops at the edge. The density generically drops as |α − α c | η , so that an interval that contains a fixed number of eigenvalues scales as N δ , with δ = 1/(1+η). For the typical values η = 1/2 and η = 1/3, this yields δ = 2/3 and δ = 3/4, respectively. These exponents indeed control the universal behavior of the average characteristic polynomial and, in particular, the double scaling limit that holds near the Durhuus-Olesen transition. We refer to Ref. [7] for a thorough analysis. We shall present an example of the similar double scaling limit (so-called Pearcey universality class) in the (simpler) context of random Hermitian matrices in the next sections.

Random walk of Hermitian matrices
The analysis of the simpler case of random Hermitian matrices [10] will reveal why the emergence of the inviscid Burgers equation is quite natural in the study of the spectral density of matrices of large sizes N . In the next section, we shall argue that this equation naturally extends to the viscid Burgers equation for the average characteristic polynomial, with 1/2N playing the role of the viscosity.
Following Dyson, we regard a random matrix as the result of independent random walks undergone by each independent matrix element. As is wellknown, the change of variables from the matrix elements to the eigenvalues introduces a Coulomb-like repulsion between the eigenvalues, so that the random walks of the eigenvalues x i obey Equivalently, the joint probability distribution P (x 1 , . . . , x N , t) obeys the Smoluchowski-Fokker-Planck equation Although the solution of this equation, for the trivial initial condition corresponding to a null matrix, is easy to obtain we shall be interested in simpler objects, such as reduced densities and, in particular, the eigenvalue densitỹ Knowing the equation satisfied by P (x 1 , . . . , x N , t), it is easy to establish the equation satisfied by the eigenvalue density in which the two particle density appears:ρ(x, . This equation is the first equation of an infinite hierarchy of equations, somewhat analogous to the BBGKY hierarchy of equations for the n-point functions in statistical mechanics, or the Schwinger-Dyson equations in quantum field theory. In order to solve such equation, some truncation is needed. We shall base the present truncation on the large N limit. In this limit, one expects the connected part of the two particle density, in the expressionρ(x, y) =ρ(x)ρ(y) +ρ c (x, y) to be subleading (of the order of 1/N ) with respect to the disconnected part, which is simply the product of single eigenvalue densities. It is then convenient to redefinẽ and to also redefine the time, setting τ = N t, to take into account the fact that the rate for "collisions" among eigenvalues (the last term in Eq. (22)) is N times larger than the rate of diffusion. We then obtain where the terms of the order of 1/N have been moved to the right-hand side of the equation, and can be ignored at this stage (we return to finite N corrections in the next section). On the left-hand side, one recognizes a non-linear term that is reminiscent of that present in the Burgers equation, except for the denominator that reflects the repulsion among the eigenvalues. But it is easy to get rid of this denominator, at the price of introducing the complex resolvent of the random matrix H where ρ is the eigenvalue density of H. Then, simple manipulations (involving taking the Hilbert transform of the equation) yield the following equation This equation is a complex Burgers equation. It is the analog of Eq. (9) of the unitary case, however here it is the resolvent itself that enters the equation (instead of the function F of the unitary case). This (small) difference is to be attributed to the specificity of the unitary group. At this point, we note that the characteristic equation (in the complex domain) takes the form of z = ξ + τ G 0 (ξ). For the initial condition corresponding to a vanishing matrix, G 0 (z) = 1/z. Then, the characteristic equation yields an implicit equation for G, G = 1/(z − τ G) which defines G(z, τ ) as an algebraic curve, with possible square root singularities. These singularities, we associate with shocks and they correspond to the edges of the spectrum (in the large N limit).

Beyond the large N limit with the average characteristic polynomial
It appears that the fluid analogy that guided us in the previous section to the inviscid Burgers equation cannot be pursued in a simple way. This is because in Eq. (24) we do not know how to treat the connected 2-point function in order to truncate the hierarchy in a consistent fashion.
Our first attempt to study the finite N corrections was to use other well-known techniques in Random Matrix Theory, namely the technique of orthogonal polynomials, from which we can reconstruct the spectral information. For the Hermitian matrices, these polynomials are Hermite polynomials. By imposing that these polynomials be orthogonal with respect to the distribution (20), one finds that these polynomials differ from the usual ones by a simple scaling. They are given explicitly by One can verify, by a direct calculation, that these polynomials satisfy the differential equation or, by taking the Cole-Hopf transform f k (z, τ ) ≡ 2ν s ∂ z ln π k (z, τ ), the viscid Burgers equation This is now an exact equation, with all finite N effects captured by the viscosity term. The knowledge of the orthogonal polynomials allows us to reconstruct the spectral density (and other correlation functions) using standard techniques of Random Matrix Theory. However, we would like to proceed in a different way. Let us observe first that the (monic) polynomial with n = N is the average characteristic polynomial, det(z − H(τ )) = π N (z, τ ) which, as we have already emphasized, coincides with the resolvent in the large N limit. We have indeed and Tr ln (z − H(τ )) = ln det (z − H(τ )) ≈ ln det (z − H(τ )) , where the last, approximate, identity holds in the large N limit. The viscid Burgers equation (29) constitutes, therefore, a natural extension of the inviscid Burgers equation for the resolvent, albeit for a slightly different object. At this point, before we proceed further, we can show how we can use the viscid Burgers equation in order to analyze the behavior of the average characteristic polynomial near a shock. To that aim, recall that we need to look at interval in the eigenvalues spectrum that scales as N δ , with δ = 1/(1 + η). As for the correction to the large N density, it is of the order of N γ with γ = δ − 1. Focusing, for the sake of illustration, on the case η = 1/2, corresponding to δ = 2/3 and γ = −1/3, one is then led to look for solutions of the viscid Burgers equation in the form A simple analysis then allows us to recover the well-known Airy behavior in the vicinity of the edges of the spectrum. The method of orthogonal polynomial is powerful, but it has limitations. For instance, the polynomials depend on the initial condition (here trivial), and there are cases where such polynomials are unknown. However, the existence of the diffusion and Burgers equations seems to be generic. We have already seen that these equations can be derived in the unitary case by using a character expansion of the determinant (or its inverse). In other cases, we can derive the equations by expressing the determinant as Gaussian integrals over Grassmann or complex variables. This has been achieved in a number of cases. Let us just focus on the Hermitian case, which is the simplest. In this case, the equations are simply where Q N denotes the average characteristic polynomial, and θ N the average of the inverse characteristic polynomial. One can then easily construct the solutions of these equations, for arbitrary initial conditions, as convolutions of the heat kernel with the initial conditions. We get and This contour can be deformed to the real axis and a half circle enclosing the pole at 0 from above (Imz > 0) or below (Imz < 0). These integral representations provide an alternative method to study the universal behavior near the shocks to that just mentioned base on the Burgers equation. Here, the analysis will follow the saddle points and their merging (recall that the saddle point equation coincides with the equation for the characteristics of the Burgers equation). For detailed analysis of large N limit of above integral representations, we refer to [11].

Spontaneous breakdown of the chiral symmetry in QCD
The fluid analogy outlined in the previous sections helps also to understand another challenging problem of strong interactions -i.e. the spontaneous breakdown of the chiral symmetry in QCD. The order parameter, known as a quark condensate Σ ≡ | qq | (i.e. the expectation value of the density of quark-antiquark pairs in the vacuum state of QCD) is directly related to the average spectral density ρ(λ) of the Euclidean-Dirac operator near the vanishing eigenvalue, by the so-called Banks-Casher relation [12] where V 4 = L 4 represents a Euclidean, four-dimensional volume. In order to get the non-zero value of the quark condensate on the l.h.s. of the Banks Casher formula, in the limit when the volume of the Euclidean space-time tends to infinity, a dramatic accumulation of small eigenvalues has to take place in the vicinity of zero. We argue, that this sudden increase in the density is achieved by the formation of a spectral shock wave at zero eigenvalue [13]. In order to demonstrate this phenomenon, we have to tune the random matrix model in such a way that it includes chiral properties, i.e. non-zero eigenvalues come in pairs (−λ, λ). The simplest choice reads where the entries of K, an M ×N (M ≥ N ) matrix, evolve in time t according to a Brownian walk [13,14]. The fundamental object of our studies is, again, the characteristic polynomial Q N where the averaging represents the Brownian motion, and for simplicity, we have put M = N . Using similar methods as in the Hermitian case, we arrive at the equation i.e. of the diffusion type. The next step is the corresponding Cole-Hopf transformation (f N = 1 N ∂ w ln Q N (w, t)) and a rescaling of time, τ = N t. Then, one gets from (39) therefore, again a Burgers type equation. In the large N limit, we recover the inviscid Burgers equation known the from previous sections, since the r.h.s. of (40) vanishes and f N becomes the resolvent (Green's function) G(z, τ ). Let us study now complex characteristics ξ corresponding to the case of non-trivial initial conditions, i.e. where the spectrum at τ = 0 is represented by a chirally symmetric pair (−a, a). In the physical context, a may e.g. represent the lowest Matsubara frequency πT , where T is the temperature. First, we notice that the singularities (shock waves) may belong to three distinct classes, corresponding to the case when τ is smaller, equal or larger comparing to some critical value τ c = a 2 , respectively. Using the parametrization δ from the section devoted to Yang-Mills studies, we arrive at δ = 2/3, δ = 3/4 or δ = 1, respectively. The microscopic properties of the spectrum can be now rederived by, first, setting second, plugging-in these quantities into, exact for any N , equation (40), and finally, performing the large N limit. In the first case (τ < τ c ), corresponding to the scenario when the spectrum forms two symmetric bumps with respect to zero, averaged spectral determinant of the Dirac takes the shape of the Airy function -chiral symmetry is unbroken. In the third case (τ > τ c ), where gap at zero is closed, averaged spectral determinant takes the shape of the Bessel function. Temperatures are still low enough, so the chiral symmetry is broken, and the spectral density within the zero eigenvalues is completely determined by the so-called Bessel kernel [15]. The most interesting is the middle case, when the gap closes (opens) at the critical value τ = τ c , i.e. the chiral symmetry is just to be broken (restored). In this case, in the large N limit, all three saddle points corresponding to the integral representation of the characteristic polynomial merge, yielding, after proper rescalings [14]: We note the qualitative similarity to the integral representation of the Pearcey function, encountered at the strong-weak coupling transition for the Wilson loop. This is not unexpected, since in both cases δ = 3/4, corresponding to the closure of the gap. However, quantitatively, the behavior is different. This is caused by the fact, that in the chiral case, the diffusion operator on the r.h.s. of (40) represents the radial part of two-dimensional Laplace operator, comparing to one-dimensional Laplace operator in (17). The two-dimensional pattern of the diffusion equation origins form the chiral symmetry -the spectrum of (37) is invariant under the transformation K → Ke iφ , where φ represents azimuthal angle. Historically, Bessel functions were derived to describe the cylindrically symmetric propagation of the heat, hence the omnipresence of the Bessel functions in the description of universal properties of the chiral models should not be puzzling. The chirally-symmetric analog of the Pearcey function is sometimes called the spun cusp or the Bessoid. The Pearcey function has a direct analogy in the diffractive optics, corresponding to the so-called cusp catastrophe [16]. The spun cusp (Bessoid) can be also, in principle, measured in certain dichroic crystals [17]. In Fig. 1, we compare the "interference" patterns of both functions, plotting the modulus of the Bessoid function versus the modulus of the Pearcey function.

Conclusions
In these notes, we concentrated on links between the diffusive behavior of certain random matrix models and QCD. In particular, we stressed that two most spectacular phenomena of strong interactions -the weak-strong coupling transition and the spontaneous breakdown of the chiral symmetry -require qualitatively similar massive rearrangement of the eigenvalues of pertinent operators (here Wilson loop and Euclidean-Dirac operator, respectively). Such rearrangements belong to the universality classes of certain random matrix models. In the case of the spontaneous breakdown of the chiral symmetry, agreement of spectral properties of QCD with chiral Random Matrix Theory was confirmed in an impressive way in several lattice studies, for various realizations of fermions, various values of number of flavors, topological charges, finite masses and temperatures, etc. [18]. As far as we know, the Bessoid-like kernel was not yet measured, despite it is known how to construct the microscopic spectral densities in such cases [19][20][21]. The Pearcey character of the Durhuus-Olesen phase transition was convincingly confirmed in large N c Yang-Mills lattice studies [22]. The joint simulation of the microscopic Dirac spectrum and the Wilson loop spectrum at the closure of the gap was also not done. We believe that such simulation could shed more light on the mutual relation between both phase transitions in QCD. Last but not least, we express the belief that the presented here "hydrodynamic" way of approaching the spectral properties of QCD may also turn out beneficial in still mysterious domains of strong interactions, alike finite density, strong CP violating angle or strongly correlated quark-gluon plasma.