Burgers-like equation for spontaneous breakdown of the chiral symmetry in QCD

We link the spontaneous breakdown of chiral symmetry in Euclidean QCD to the collision of spectral shock waves in the vicinity of zero eigenvalue of Dirac operator. The mechanism, originating from complex Burger's-like equation for viscid, pressureless, one-dimensional flow of eigenvalues, is similar to recently observed weak-strong coupling phase transition in large $N_c$ Yang-Mills theory. The spectral viscosity is proportional to the inverse of the size of the random matrix that replaces the Dirac operator in the universal (ergodic) regime. We obtain the exact scaling function and critical exponents of the chiral phase transition for the averaged characteristic polynomial for $N_c \ge3$ QCD. We reinterpret our results in terms of known properties of chiral random matrix models and lattice data.


Introduction
It was a recently emphasized that the Burgers equation could be used to understand universal features of the weak to strong-coupling transition in two-dimensional Yang-Mills theory with a large number of colors N c [1,2,3]. This transition, first studied by Durhuus and Olesen [4], can be pictured, in the language of the Burgers equation, as resulting from the collision of two spectral shock waves at the closure of the gap. The emergence of the Burgers equation in such kind of problems, and related questions in random matrix theory, seems generic, and can be simply understood by exploiting Dyson's original idea of matrix random walks. Then, after a proper rescaling of time that separates the fast motion of the eigenvalues caused by their mutual repulsion, from the diffusion that results from the random walks of the matrix elements, one can easily show, by using standard tools of statistical mechanics, that the average resolvent obeys indeed a Burgers equation (or its simple generalizations) in the limit of large size of the matrices. Here we shall obtain the Burgers equation from an exact equation, akin to a diffusion equation, satisfied by the average characteristic polynomial. The average resolvent and average characteristic polynomial are simply related in the large N limit. However, the equation for the characteristic polynomial is exact for any matrix size, while no such equation seems to exist for the average resolvent.
Returning to two-dimensional Yang-Mills theory, we note that the average characteristic determinant for a Wilson loop of a given area was indeed shown to obey an exact diffusion equation for any number of colors N c (which indeed yields a Burgers equation in the infinite number of colors limit). For finite number of colors, the flow of the eigenvalues is viscid, with a negative spectral viscosity v s = 1 2N c . It is the negative sign of the viscosity, or of the diffusion constant, that causes the rapid spectral oscillations at the closure of the gap. This picture was confirmed by extensive numerical studies of full Yang-Mills theory in three dimensions and recently, in four dimensions [3,5].
In this letter, we demonstrate that a similar mechanism, governed by similar viscous spectral Burgers-like equation, is responsible for the universal spectral oscillations of the spectrum of the Dirac operator in QCD that accompany the spontaneous breakdown of the chiral symmetry [6,7]. In the case of the Durhuus-Olesen transition, where the role of the time is played by the area of the Wilson loop, an explicit construction of a random matrix model was provided by Janik and Wieczorek [8], matrices attached to large loops being formed by multiplying random matrices representing small loops. Here, as already mentioned, following Dyson [9], we add a fictitious time (somewhat analogous to Schwinger's proper time) in order to describe the diffusion of the random matrices in 4 dimensional Euclidean space.

Random matrix theory and QCD
A cornerstone for the microscopic understanding of the spontaneous breakdown of chiral symmetry is the Banks-Casher [10] relation where the quark condensate < qq > is an order parameter for chiral symmetry, ρ(0) is the averaged (over the gauge field configurations) level density of the Euclidean Dirac operator near the vanishing eigenvalue, and V 4 ≡ L 4 is the Euclidean volume. This relation shows that chiral symmetry breaking requires a strong accumulation of eigenvalues near zero, i.e. a level spacing ∆ ∼ 1/L 4 , much larger than the level spacing ∆ ∼ 1/L of a free system [11]. This accumulation of eigenvalues leads to universal properties, that are well captured by random matrix theory: For eigenvalues smaller than a characteristic energy scale, referred to as the Thouless scale E T h , the fluctuations of eigenvalues are described by chiral random matrix models respecting the global symmetries of the Dirac Hamiltonian. In QCD, the condition E QCD T h /∆ = F 2 π L 2 ≫ 1, where F π is the pion decay constant, determines the regime of applicability of random matrix theory [12,13].
In Euclidean QCD, all four Dirac matrices can be chosen to be anti-hermitian, hence the spectrum of the massless Dirac operator D ≡ iγ µ (∂ µ − igA µ ) is purely imaginary. The partition function, for a fixed topological sector, reads where the averaging is done with respect to gluonic configurations of a given topological charge ν, ±iλ k are the eigenvalues of D, and m f is the mass of a quark with flavor f . Due to the chiral symmetry non-zero eigenvalues of D come in pairs, and the number of fermionic zero modes is related to the topological charge. In the chiral Gaussian random matrix model (hereafter χGUE), corresponding to QCD with N c ≥ 3, the role of the massless Dirac operator is played by a random matrix W(= −iD) of the following form: Here K is a rectangular M × N (M > N) matrix with complex entries, K i j ≡ x i j + iy i j , where x i j and y i j are drawn from a Gaussian distribution. Note that W is hermitian, so that its eigenvalues κ i 's are real. The block-diagonal structure of (3) reflects the chiral symmetry of the Dirac operator: W anticommutes with the analogue of the Dirac matrix γ 5 , defined here as γ 5 = diag(1 N , −1 M ). This implies in particular that the eigenvalues come in pairs of opposite values, (κ, −κ). By construction, W has in addition ν ≡ M − N zero eigenvalues. These mimic the zero modes of quarks propagating in gauge fields of non trivial topology. General spectral properties of the random matrices can be obtained form correlation functions containing both products and ratios of the characteristic polynomial [14] where < ... > denote averaging with respect to the χGUE measure, and the characteristic polynomial is Z(w) = det(w − W). Directly related to (4) is the resolvent, R(z) = ∂ w C(z; w)| w=z , whose imaginary part yields the spectral density. In the vicinity of zero, the microscopic (unfolded) resolvent predicts for QCD [15] where b = N f +|ν|, with N f the number of quark flavors, and x = zV 4 Σ. The appearance of Bessel functions I b , K b in the correlation functions, such as in Eq. (5), is generic. They encode the universal behaviors that show up already in the simplest objects, the average characteristic polynomial and/or the average of the inverse of the characteristic polynomial. As the first result of the paper, we shall obtain an exact differential equation for the averaged characteristic polynomial. This equation is akin to a diffusion equation, from which a Burgers equation can be derived from a simple transformation.

Burgers equation and the characteristic polynomial for Dirac operator
We assume now that the entries of the matrix K follow independent random walks. Let us denote by P(X, Y, t) the joint probability that the entries of K take the values The random walks of the matrix elements translates into the following diffusion equation for P(X, Y, t): We shall be interested in this paper with the time evolution of the averaged characteristic polynomial where n = N + M, w is an arbitrary complex number, and dX = i j dx i j , and similarly for dY. In order to get the equation obeyed by Q ν n (w, t), we consider first the equation satisfied by the average characteristic polynomial of the associated N × N Wishart where To derive the equation for M, we take a time derivative of the expression above. This acts on P(X, Y, t), which, using Eq. (6), we transform into derivatives with respect to x i j and y i j . Then, we integrate by parts, and use standard Grassmann calculus to obtain (after a somewhat tedious but straightforward calculation) the following differential equation This equation is valid for any N and M, and arbitrary initial conditions. Note that for the trivial initial condition K i j (t = 0) = 0, its solution is given by time dependent associated Laguerre polynomial [16]. From the equation for M ν N (z, t), Eq. (9) above, one easily obtains the equation for It will be also useful to consider the equation for the Cole-Hopf transform of Q ν n (w, t), t)). This object identifies with the average resolvent in the large n limit. After a rescaling of the time, τ = Mt, one gets from Eq. (10) where we have separated on the left hand side the terms that survive the large n limit, and on the right hand side the terms that are explicitly suppressed by powers of 1/n. Note the crucial role played by the rescaling of time in arriving at this equation. The motivation behind this rescaling is that the diffusion associated with the random walks is taking place over a time scale that is larger, typically by a factor n, than the time scale corresponding to the local rearrangements of the eigenvalues due to their mutual repulsion [9,16]. After rescaling, the diffusion terms are dwarfed by a factor 1/n, and the large n dynamics is dominated by repulsion. The last term, of order 1/n 2 finds it origin in the kinematical zero modes present when ν 0.

Large n limit
We consider now the limit n → ∞, with ν constant. We set g(w, τ) = lim n→∞ f ν n (w, τ). Eq. (11) reduces then to the inviscid Burgers equation [1,17,18], independent of ν: It can be solved using complex characteristics. We choose the system to be initially in a chiral symmetric state, with eigenvalues localized at ±a.The characteristic lines are given by w = ξ + τg 0 (ξ), with The solution g(w, τ) is constant along the characteristics, meaning g(w, τ) = g 0 ξ(w, τ) . By eliminating ξ, one obtains an implicit equation for g: This equation can be solved by elementary means, and well known results recovered. In fact, the change of variables an equation for a time-independent resolvent G st that has been obtained in this context using different techniques [19,20,21]. In previous studies, the parameter d was introduced as the "deterministic" part of the chiral matrix (with K in Eq. (3) replaced by K + d, d being fixed and K random), in order to control the approach to the chiral transition. The aforementioned change of variables renders transparent the dynamics captured by the Burgers equation: For small time, i.e. τ < a 2 , the spectral density remains localized in humps centered around the values ±a. As time reaches the critical value τ c = a 2 (corresponding to d = 1 in the static approach), the two domains of the spectrum merge at the origin, which we picture as the collision of two spectral shock waves. A finite "condensate" then develops, and chiral symmetry is spontaneously broken.
We now recover these features by studying the singularities of the characteristics, and the behavior of the solution in the vicinity of these singularities. This brief discussion will pave the way for the scaling analysis to be performed in the next section. Singularities appear when characteristics start to cross (appearance of a spectral shock wave). This occurs for values of ξ that obey the equation, that is The values of ξ c correspond to the edges of the spectrum when τ < τ c = a 2 . At the critical time τ c = a 2 , where the two humps of the spectrum start to merge, the equation for ξ c admits a double solution at ξ c = 0, which splits into two purely imaginary and opposite solutions when τ > a 2 .
In order to study the behavior of g at the edge of the spectrum, we expand g 0 in the vicinity of a singular point When τ < a 2 , g ′ 0 (ξ c ) = −1/τ and, using g 0 = (w − ξ)/τ, we get One can then easily invert the relation between w and ξ, and get (for the rightmost edge) so that which exhibits the familiar square root behavior of the spectrum near its (right) edge. For τ = a 2 , a similar analysis taking into account that g 0 (ξ c ) = 0 = g ′′ 0 (ξ c ), so that the cubic term must be kept, yields (for w c = 0) For time τ > a 2 , we have g(w = 0, τ) = g 0 (ξ c (w = 0, τ)) = − √ a 2 −τ τ , which is imaginary and hence directly proportional to the spectral density ρ(0).
The behavior of g(w) at the edge of the spectrum determines the average eigenvalue spacing in the limit of large matrices. A singularity ∼ |w − w c | α yields a level spacing ∼ n −δ with δ = 1/(1 + α). We have therefore δ = 2/3 for τ < a 2 , δ = 3/4 for τ = a 2 and δ = 1 for τ > a 2 . We shall exploit these properties in the next section.

Critical properties of the characteristic polynomial
In this section we carry out a scaling analysis of the average characteristic polynomial Q(w, τ), or its Cole-Hopf transform f (w, τ), in the vicinity of the singular points. To that aim, we set with γ = 1 − δ, and s and χ(s, τ) remain finite as n → ∞.

Bessel universality
For τ > a 2 , there is no singular behavior (δ = 1). So, we set w = n −1 s and f ν n = χ. In the large n limit (at ν constant) we obtain, following the same manipulations as above, the following partial differential equation: which integrates to χ 2 + ∂ s χ + χ s − ν 2 s 2 + u(τ) = 0. Then, setting χ(s, τ) = ∂ s ln φ(s, τ) we obtain: whose solution is The determination of the arbitrary function u(τ) proceeds as in the previous case, by matching the asymptotic χ(s, τ) ∼ −i √ u(τ) with the large N solution. This gives √ u(τ) = √ τ − a 2 /τ. We recover the scaling of the ratio of spectral densities discussed for instance in [22,23,24].

Bessoid (axially symmetric Pearcey) universality
Finally, we move to the case of τ = τ c = a 2 , which will lead to the second new result of this paper. As we have shown above, at τ = a 2 , the two pre-shocks collide. Before the collision, these pre-shocks are accompanied by oscillations of the Airy type. Our purpose now is to describe the modification of the pattern of oscillations for τ close to τ c . To perform this analysis, it is most convenient to start from the diffusion equation obeyed by the characteristic polynomial, i.e., Eq. (10) which, after changing t into τ = Mt, we rewrite as This equation is to be solved with the initial condition Q ν n (w, τ = 0) = w ν w 2 − a 2 N (where the function of w, (w 2 − a 2 ) N is defined with a cut between −a and a). It can be verified by a direct calculation (that does not require the explicit calculation of the integral below), that where C = (−1) ν+1/2 2M, is a solution with the proper initial condition. (For ν = 0, this solution agrees with a known solution [25].) The y-integral runs over a half-line that starts at the origin and goes to infinity, making a constant angle φ y = arg(y) with the real axis, with −π ≤ arg(y) < π. For the integral to be convergent, we require π 4 < |φ y | < 3π 4 . The modified Bessel function I ν (x) had the following asymptotic expansion lim |x|→∞ I ν (x) ≃ 1 √ 2πx e x , valid for | arg(x)| < π 2 (see [26]; here x = 2Myw τ ). This is useful in particular to verify the initial behavior. Indeed, as τ → 0, one may estimate the integral using the saddle point method. The saddle point equation yields y + w = 0, which fixes in particular arg(y) = arg(−w). This new condition for φ y , together with the convergence condition noted above, are easily seen to be compatible with the condition of validity of the asymptotic expansion of the Bessel function. In turn, these conditions limit the allowed arguments of w to π 4 < | arg(w)| < 3π 4 . The integral representation (30) of the characteristic polynomial allows us to study the vicinity of the critical point. We note that the saddle point equation reads (in the large n limit) Identifying y = −ξ, we recognize the equation for the characteristic lines. This indicates how the large n dynamics is coded in this integral. We shall focus more specifically at the critical point, w = 0, y = 0, τ = a 2 . In this regime, we may expand ln(a 2 − y 2 ) ≈ ln(a 2 ) − y 2 a 2 − y 4 2a 4 , and obtain with C ′ = C(−1) N . To capture the critical behavior also as a function of time, we set τ = a 2 + θ. Then Eq. (32) becomes This expression suggests the following change of variables that will ensure a smooth large n limit: We then define where the limit is taken with ν constant. Finally This exact scaling function for the characteristic polynomial is the second important result of this paper. It has a form very similar to the Pearcey function, which gives the asymptotic behavior of the characteristic polynomial in the case when gap closes for GUE [27] or in the case of unitary diffusion on the circle [5] P(q, t) = ∞ −∞ dy exp(−y 4 − ty 2 + qy).
Here q is the rescaled angle representing the position of the eigenvalue on the unitary circle and t parameterizes the fluctuations around the critical area. The critical indices that determine the scaling with n in Eq. (34) are identical to those in the Pearcey integral, but the form of the integral is different. The reason is the chiral symmetry, which imposes an additional polar symmetry of the spectrum, trading the exponential function of q in the Pearcey integral (37) for a Bessel function in Eq. (36).

Conclusions
In this letter we have obtained an exact differential equation for the average characteristic polynomial of a chiral random matrix, and its Cole-Hopf transform. For the latter, the equation takes the form of a generalized viscid Burgers equation, where the viscosity is proportional to the inverse of the size of the matrix, but with a negative sign. This allowed us to provide a complete description of the full critical behavior of averaged characteristic polynomials in chiral QCD, based on a single equation. In particular, we considered the case of chiral Gaussian Unitary Ensembles and we have identified the exact universal scaling function (Bessoid B ν (q, θ)) in the vicinity of the chiral critical point for the average characteristic polynomial. We did not analyze in this letter the properties of the average of the inverse characteristic polynomial, but we have checked that if fulfills similar equations, albeit with initial conditions singular at w = 0, alike in the cases of unitary and GUE diffusions.
Since F π scales like √ N c , more and more eigenvalues of the Dirac operator fall into the universal window when the number of colors tends to infinity, the volume of the lattice being kept finite. This suggests, that the present study is relevant also for analyzing the spontaneous breakdown of chiral symmetry at finite volume and large N c QCD, which was observed, and explained by Neuberger and Narayanan [28]. For small lattice sizes, chiral symmetry is unbroken, while at some critical scale L c a condensate is formed. The same authors [29] have also observed, that the N c -dependence of the level spacing closest to zero goes from 1/N c in the broken chiral symmetry phase to 1/N 2/3 c in the symmetric (gapped) phase. At L = L c , the critical scaling changes to 1/N 3/4 behavior and the condensate vanishes at the critical size L c as √ L − L c . These results are in agreement with our analysis.
The critical universal scaling function (Bessoid) for large N c Dirac operator resembles closely the critical universal scaling function (Pearcey's cuspoid) for the weak to strong coupling transition in Yang-Mills theory at large N c . It would be interesting to study numerically both transitions simultaneously (at least in some simple model like [28]) to see the interplay between the cuspoid and the Bessoid, or, in other words, the relation between the critical size L c for chiral symmetry breakdown and the critical area ∼ L 2 c for the weak to strong coupling transition in Yang-Mills theory. Of particular interest would be to measure on the lattice the microscopic spectral density exactly at the point of the transition. As far as we know, the microscopic spectral density at criticality was constructed explicitly only for b = N f + |ν| = 0 case [30], and was never checked by the lattice simulation. Taking into account the still ongoing discussion on the nature of chiral phase transition, its relation to confinement and Anderson localization [31,32], lattice verification of analytic predictions for microscopic densities at the critical point may be a powerful tool to shed more light on this aspect of strong interactions.