Unveiling the significance of eigenvectors in diffusing non-hermitian matrices by identifying the underlying Burgers dynamics

Following our recent letter, we study in detail an entry-wise diffusion of non-hermitian complex matrices. We obtain an exact partial differential equation (valid for any matrix size $N$ and arbitrary initial conditions) for evolution of the averaged extended characteristic polynomial. The logarithm of this polynomial has an interpretation of a potential which generates a Burgers dynamics in quaternionic space. The dynamics of the ensemble in the large $N$ is completely determined by the coevolution of the spectral density and a certain eigenvector correlation function. This coevolution is best visible in an electrostatic potential of a quaternionic argument built of two complex variables, the first of which governs standard spectral properties while the second unravels the hidden dynamics of eigenvector correlation function. We obtain general large $N$ formulas for both spectral density and 1-point eigenvector correlation function valid for any initial conditions. We exemplify our studies by solving three examples, and we verify the analytic form of our solutions with numerical simulations.


Introduction
The concept of matrices filled with entries subject to the diffusion process was first introduced by Dyson [2] and applied in the context of both Gaussian Unitary Ensemble (GUE) and Circular Unitary Ensemble (CUE). The arising Coulomb gas analogy had a major impact on understanding of random matrices [3]. Today the marriage of stochastic processes and random matrices still brings several new insights. Examples include the study of determinantal processes [4,5,6], Loewner diffusion [7] or the fluctuations of non-intersecting interfaces in thermal equilibrium [8].
Recently, several authors [9,10] have approached the diffusion in the GUE from a new perspective. They found a viscid complex Burgers equation for the logarithmic derivative of the averaged characteristic polynomial f N (z, τ ) ≡ 1 N ∂ z ln U N (z, τ ) (associated with a Hermitian matrix filled with entries performing Brownian motion in the complex space): where τ is the diffusion time and z is a complex variable. The role of viscosity is played by the inverse size of the matrix N . In the N → ∞ limit, f N (z, τ ) becomes the Green's function G(z, τ ) and the partial differential equation becomes inviscid: ∂ τ G(z, τ ) + G(z, τ )∂ z G(z, τ ) = 0.
A solution of the latter equation (by the method of characteristics) requires an introduction of shocks, which turn out to coincide with the edges of the spectra. This phenomenon leads to a novel interpretation of known matrix results, since microscopic universal behavior of the spectra emerges as an expansion around the shock wave of the viscous equation. Nontrivial initial conditions give rise to shock collisions which are equivalent to the merging of the spectrum boundaries. This way, not only Airy but also the Pearcey functions are captured in the same formalism. A similar Burgers equation was also obtained for the Wishart ensemble and chiral GUE yielding a universal scaling associated with the Bessoid function [11]. The equivalent phenomenon appears also at the level of CUE diffusion, providing new insight for order-disorder transition of Wilson loops in Yang-Mills theory [12].
Recently, this program got extended to the realm of matrices with complex eigenvalues [1]. The success of such an extension is a priori surprising, since hermitian and non-hermitian matrix models seem to be hardly comparable. In the former, the hermiticity condition confines the eigenvalues to the real axis. In the latter, there are no constraints, and the eigenvalues spread over the whole complex plane. Furthermore, non-hermitian models develop discontinuities at the spectral density boundaries, a feature observed even for the well-known Ginibre Ensemble [13], for which the spectral density is given by: whereas e.g. in the GUE case, the Wigner semicircle is continuous across the spectral edge, and only the derivatives are discontinuous. Other differences arise when one considers the evolution of the diffusing matrices. In the hermitian case, the evolution is determined by the initial eigenvalues only, whereas in the non-hermitian models, the information on initial eigenvectors additionally affects the shape of the spectral density. This paper is a continuation and an extension of the ideas of non-hermitian diffusion announced briefly in [1].
Our approach is rooted in the standard electrostatic analogy, but requires a novel setting (we call it the quaternionic method), which can be viewed as an extension of the standard Dysonian strategy applied originally to the hermitian (or unitary) matrix diffusion case. The main objects of our interest are the spectral density (obtained from the resolvent) and a certain one-point eigenvector correlation function. We stress that the aforementioned eigenvector correlator is crucial for understanding the diffusion process of nonhermitian matrices. We also point out, why the importance of this correlator was disregarded in majority of the studies of non-hermitian random matrix models.
The basic object of our studies is an averaged "extended" characteristic polynomial (AECP). An extension follows from an introduction of two pairs of complex variables (compared to one complex variable in standard treatments). Surprisingly, AECP obeys a certain partial differential equation akin to the diffusion equation, for arbitrary size of the matrix and for arbitrary initial conditions, and is exactly integrable. The diffusion happens in the auxiliary plane "perpendicular" to the complex plane where the eigenvalues reside. In the large matrix size limit, the logarithm of the AECP can be viewed as an electrostatic potential and its derivatives with respect to the two complex variables yield a pair of coupled Burgers-like equations for the non-hermitian Green's function and eigenvector correlation function. We would like to mention, that in the standard electrostatic analogy [14,15,16] the "second variable" is treated as an infinitesimal regulator only. This is the reason why the dynamics, as a function of this variable, remained hidden, and the complementary information on the eigenvector correlator co-evolving with the spectra was absent.
To illustrate our findings, we consider a couple of examples of initial conditions. In most of them we demonstrate the explicit solutions of Burgers equations, obtaining formulas for the spectral density, eigenvector correlators and the electrostatic potential. We note that an inspection of the Burgers-like equation identifies the shock line with the non-holomorphic sector of the spectral density. Moreover, we show the insensitivity of the shock formation to the initial condition chosen. This hints to a lack of truly new universality classes in this type of models, which is furthermore corroborated in the study of the universal behavior in the vicinity of the spectral collision in one of the examples.
The paper is organized as follows. In Section 2 we define Dysonian non-hermitian diffusion. In Section 3 we briefly review the electrostatic analogy and the quaternionic method. We proceed in Section 4 by deriving the partial differential equation for the AECP and presenting its integral representation. In Section 5 we derive a pair of coupled Burgers equations for the on-and off-diagonal parts of the quaternionic Green's function (in the large N limit) thus making a link to the quaternionic method. Subsequently we solve them with the method of complex characteristics. Finally, we obtain an implicit solution to the equation for the potential in terms of the Hopf-Lax formula and derive large N formulas for the spectral density, eigenvector correlation functions and the boundary of the spectrum valid for an arbitrary initial matrix. Section 6 is devoted to the examples of a) Ginibre, b) spiric and c) 1-band non-normality matrices. We apply previously described methods to these cases, depict the characteristics picture and obtain the large N limit spectral density and eigenvector correlators. We also comment on critical behavior of the EACP. Section 7 is devoted to a curious observation by Osada [17], which actually has triggered our interest in the diffusion of the Ginibre ensemble. We provide also an explanation of the Osada's observation. Section 8 summarizes the paper and outlines some possibilities of further investigations of the observed patterns.
Three appendices hide technicalities: Appendix A demonstrates the derivation of the key diffusion equation for an AECP, Appendix B clarifies the link between Ginibre and Wishart/chiral ensembles, and Appendix C determines weights and normalizations needed for establishing the universal scaling at the shock.

Dysonian non-hermitian diffusion
where B x ij and B y kl are independent Wiener processes. We restrict ourselves to deterministic initial conditions, that is we assume that each element of the matrix has a given initial value x ij = (x 0 ) ij and y ij = (y 0 ) ij for τ = 0. This condition can be concisely written as X = X 0 , for τ = 0. Clearly, this model is a straightforward extension of the Dyson random walk [2] to the realm of non-hermitian random matrices.
The joint probability density for matrix elements evolves according to the 2N 2 -dimensional diffusion equation with the initial condition P (X, 0) = δ(X − X 0 ). The probability measure for random matrices at time τ is defined by dµ τ (X) = D[X]P (X, τ ) where D[X] = ij dx ij dy ij , and the statistical averages by The hermitian version of the model, discussed by Dyson, reduces to a model of evolution of eigenvalues. In that case eigenvectors can be integrated out. What makes the non-hermitian extension interesting is that in addition to eigenvalues one has to control also the evolution of eigenvectors. We present a systematic method to do so.

Electrostatic analogy and quaternions
In this section we briefly recall the method to calculate eigenvalue distribution of random matrices in the limit N → ∞. The method is based on "electrostatic" analogy [14,15,16]. One defines a quantity which can be interpreted in the limit w → 0 as an electrostatic potential of a cloud of N electric charges interacting on the z-complex plane. The corresponding electric field is Identifying the real and imaginary part of as vector components G = (E x − iE y )/2 one can rewrite the last equation in the vector notation as E = (E x , E y ) = ∇ z Φ. The minus sign in front of E y and the scale factor 1/2 in the relation of G to electric field E is a matter of convention. We are interested in the eigenvalue distribution where z i 's are the eigenvalues of X. The limiting eigenvalue density can be calculated from the Gauss law This relation follows from a standard representation of the complex Dirac delta function πδ (2) The expression in the brackets on the r.h.s. of (8) can be cast into the standard form of resolvent (z − X) −1 at the price of introducing 2N × 2N matrices in place of the original N × N ones. The resolvent is a 2 × 2 matrix where the block-trace is defined as We refer to G(z, w, τ ) (12) as to generalized Green's function or quaternionic resolvent [18,19]. Note, that we use the representation of the quaternion in terms of Pauli matrices, i.e. Q = q 0 1 2 + i 3 j=1 σ j q j , so z = q 0 + iq 1 and −w = q 2 + iq 3 . The diagonal element of the quaternionic resolvent G 11 is equal to G(z, w, τ ) (8). The extended Green's function G(z, w) is an advantageous object since one can apply geometric series expansion to (Q − X ) −1 which has a diagrammatic interpretation. This leads to a closed set of Dyson-Schwinger equations enumerating underlying planar Feynman diagrams. From these equations one can derive an exact form of the Green's function in the limit N → ∞, as well as matrix-valued addition and multiplication laws [20]. We mention that the quaternionic extension is equivalent to another approach known under the name of hermitization method [21,22,23], in which the diagonal and off-diagonal blocks of matrices Q and X are flipped before the block-trace operation.
Having determined the quaternionic resolvent G(z, w, τ ) one can determine the potential Φ(z, w, τ ) or vice versa, since the two objects are related by a simple relation: As follows from (10), the eigenvalue density can be derived from the potential Φ(z, w, τ ) using the Poisson equation It turns out that the potential Φ(z, w, τ ) encodes also information about the correlations of eigenvectors [24], being the special case of the Bell-Steinberger matrix [25,26,27]. One defines the correlation function as 1 with O αβ = L α |L β R α |R β where |L α (|R α ) are the left (right) eigenvectors of matrix X. It can be shown [28] that, in the N → ∞ limit, this correlation function is related to the off-diagonal elements of the resolvent as O(z, τ ) = − 1 π G 11 (z, 0, τ )G1 1 (z, 0, τ ). Applying (13) we have where V (z, w, τ ) = ∂ w Φ(z, w, τ ) is the velocity field, which plays the same role in the w-complex plane as the electric field G(z, w, τ ) = ∂ z Φ(z, w, τ ) in the z-complex plane. It is a vector field. If we parametrize positions on the w-complex plane as w = a + ib and V = (V a − iV b )/2, then V = ∇ w Φ. The term "velocity" is related to the underlying Burgers dynamics to be discussed later. To summarize, the limiting eigenvalue density and the eigenvector correlation function can be calculated from the electrostatic potential using eqs. (14) and (16), respectively, and taking the limit w → 0 which project quaternions to the z-plane. It remains to show how to calculate the electrostatic potential Φ(z, w, τ ) (7). The standard method is based on enumeration of planar diagrams as mentioned above. The main object in this method is the Green's function G.
Here we propose an alternative approach which is based on the diffusion equation in the quaternionic (2 + 2)-dimensional space in the direction perpendicular to the z-complex plane. The primary object in this calculation is an extended characteristic polynomial defined in the next section.

Averaged extended characteristic polynomial
In order to calculate the potential (7) in the limit N → ∞ we rewrite it as and define an associated object -an effective potential where 1 Note that we introduced an additional 1/N factor as compared to the definition given in [24] to obtain a limiting density.
We refer to D(z, w, τ ) as to averaged extended characteristic polynomial (AECP). For the Gaussian process (4) the determinant D(z, w, τ ) self-averages for N → ∞ and Φ can be replaced by φ in this limit. The advantage of using the latter is that the averaged extended characteristic polynomial D(z, w, τ ) (19), which appears in the definition of φ = 1 N log D obeys a simple diffusion equation with respect to the variable w as shown in Appendix A. Note that, from the point of view of this equation, z is a dummy parameter. The z-dependence appears only in the initial condition which is completely determined by the initial matrix X 0 . In other words, for each z we have an independent diffusion in the perpendicular w-complex plane. The problem is therefore exactly integrable, and the solution of the diffusion equation (20) reads where D 0 (z, w ) = D(z, w , 0). The solution can be equivalently written as where φ 0 (z, w ) = 1 N log D 0 (z, w ). In the limit N → ∞ the last equation assumes the form of the Hopf-Lax formula [29] φ(z, w, τ ) = max This equation describes the evolution of the electrostatic potential in the limit N → ∞ for the given initial configuration φ 0 (z, w), so this equation solves our original problem. A few remarks are in order. The characteristic polynomial D(z, w, τ ) and the potential φ(z, w, τ ) (18) depend on w only through the norm |w| 2 . The diffusion preserves the spherical symmetry of these quantities in the w-complex plane, so it is convenient to rewrite these equations in the radial part r of w = re iα skipping the dependence on the phase α. In particular (22) takes the form and (24) simplifies to The second remark is on the role of the parameter w. Originally it was introduced as a regulator to the expression for the potential (7) and eventually sent to zero. Here we promote w to a full complex variable and analyse the dynamics of the model on the entire w-complex plane. This approach allows one to trace not only eigenvalues but also eigenvectors of the random matrix X and to break the symmetry between matrices having identical eigenvalues but different eigenvectors. A complex valued matrix can be Schur decomposed X = U (Λ + T )U † where U is a unitary matrix, Λ is a diagonal matrix containing the complex eigenvalues, and T is a strictly upper-triangular matrix encoding information about eigenvectors. Two different matrices X 1 , X 2 having the same eigenvalues Λ but different eigenvectors have different T 1 and T 2 . The averaged extended characteristic polynomial (19) for these matrices differs D 1 (z, w) = D 2 (z, w) when |w| = 0. From this difference one can read off information about eigenvectors.
Moreover, the dynamics of the model in the quaternionic space has a beautiful physical interpretation in terms of the Burgers dynamics. The behavior of the model on the z-complex plane is a shadow of this dynamics. In particular, the support of the eigenvalue density ρ(z, τ ) coincides with the location of shocks of the quaternionic Burgers dynamics in the full quaternionic (z, w)-space.

Burgers dynamics
In the previous section we have found a solution to the diffusion equation (20). Here for completeness we relate the diffusion equation to Burgers dynamics. Using the definition (18) it is easy to see that the effective potential φ = φ(z, w, τ ) and its gradient v(z, w, τ ) = ∂ w φ(z, w, τ ) fulfill the following differential equations and respectively. These equations describe Burgers dynamics on the w-complex plane for a two dimensional The coefficients of the Laplacian term can be identified as a hydrodynamic viscosity parameter ν = 1 N . Equation (28) follows from (20) by an inverse Cole-Hopf transformation [30]. One can also write an equation for the z-gradient, g(z, w, τ ) = ∂ z φ(z, w, τ ): The two gradients are related to each other as ∂ w g = ∂ z v. The effective potential φ and the gradients reproduce the electrostatic potential and the quaternionic Green's function in the limit N → ∞ Let us now discuss the inviscid limit N → ∞. The effective potential (27) obeys the equation which after applying rotational symmetry (in the variable w) simplifies to an equation for the radial variable r = |w|. The solution is given by the Hopf-Lax formula (26) which in our case is equivalent to with r * being the location r * = r of the maximum (26) given by the usual extremum condition To complete the scheme, the maximizing parameter r * has to be calculated from (34) and the result r * = r * (z, r, τ ) has to be inserted to (33). We solve this equation Eliminating r * from this set of equation we obtain the effective potential φ(z, r, τ ) for an arbitrary initial matrix X 0 . From it we derive the limiting density ρ(z, τ ) (14) and the eigenvector correlations O(z, τ ) (16) For the initial condition of the form φ 0 (z, r) = 1 N Tr log M(z, r) we arrive, after differentiation and some algebraic manipulations, at where M = M(z, r * ). In the final formulas we set r = 0 to project the results to the z-complex plane. The equation for r * (35) simplifies for r = 0 to Equations (38), (39) are valid inside the boundary given by which corresponds to r * = 0. Outside this boundary O(z, τ ) = 0 and ρ(z, τ ) = 0. For a normal initial matrix X 0 , the second term in the spectral density (38) drops out since [M −1 , z − X 0 ] = 0. By inspecting the boundary equation (41) one finds a surprising connection to the so-called pseudospectrum [31], a mathematical concept generalizing the notion of the eigenvalue spectrum. Pseudospectrum of a matrix A is defined as an subset σ A of the complex plane z such that where the symbol || · || is some arbitrary matrix norm and is the parameter of the pseudospectrum. In the → 0 limit, one recovers the standard eigenvalues as poles of the resolvent (z − A) −1 . For the initial matrix A = X 0 , the boundary of the pseudospectrum subset σ X0 is exactly the eigenvalue boundary (41) with 2 = τ N and a Frobenius norm ||X|| F = √ TrX † X. From this simple identification we conclude that a diffusion model with an initial matrix X 0 is also a probabilistic realization of the pseudospectrum for the same matrix.
We finish this section by discussing an equation for the gradient v = ∂ w φ in the inviscid limit. The equation is equivalent to the one for the potential φ that we discussed above but in some situations the equation for the gradient is more handy to use. The inviscid version of the (28) reads It is an inviscid Burgers equation in 2 + 1 dimensions. A general solution to this equation for smooth (differentiable) initial conditions can be deduced from the Hopf-Lax formula (24) for the effective potential φ for N → ∞. The maximum in (24) is achieved for w = w * fulfilling the standard extremum condition which is equivalent to which is given by an implicit equation for v = v(z, w, τ ) depending only on the initial condition v 0 (z, w). The parameter w * can be viewed as a labeling parameter for the family of the characteristic lines. These lines start to cross when the labeling fails to be bijective i.e. dw dw * = 0. This singularity condition defines a caustic surface -a boundary across which an ambiguity of the solutions arises. The development of multivalued solutions is an unwanted feature of the inviscid Burgers equation and is circumvented by constructing shock lines along which one cuts the characteristics rendering the solution unique.
The two-dimensional Burger's evolution (43) can be simplified in our case due to the rotational symmetry. We are interested only in the solutions which depend on the modulus r = |w|. In this case the vector velocity field is a central field v =w r ν, with ν = |v|, and the vector equation (43) reduces to a scalar equation for the modulus of the velocity field with a solution given by From this solution we can reconstruct the full 2d-solution: v(z, w, τ ) =w r ν(z, r, τ ) with r = |w|.

Examples
In this section we discuss three examples: (1) the canonical Ginibre evolution for which the initial matrix is X 0 = 0, (2) the spiric evolution for which X 0 = diag(−a, −a..., a, a...) with an equal number of ±a and (3) an evolution initiated from a non-normal matrix: (X 0 ) ij = αδ i,j−1 . This matrix has eigenvalues equal zero as the initial matrix in the first example but it is not a normal matrix.
The first example serves as a proof-of-concept. We solve the Burgers equations by the method of characteristics to obtain the spectral density, the eigenvector correlator and the potential function in the large N limit. The discussion is accompanied by Appendices B,C where we consider finite w results and a relation of the averaged extended characteristic polynomial to the two-point kernel of the underlying determinantal process.
The second example illustrates an evolution of the eigenvalue density initiated from two disconnected eigenvalue "islands" which grow in the course of time to collide at some critical time. We discuss a novel universality arising in the vicinity of the collision.
The third example demonstrates the dependence of the evolution of the eigenvalue distribution on the initial information which goes beyond the eigenvalues themselves.

Ginibre evolution
The evolution is initiated from the matrix X 0 = 0. The evolution of the eigenvalue density is shown in Fig. 1. The spectral density forms a circular eigenvalue "island" expanding in time. For each time τ the ensemble of matrices in this evolution is equivalent to a Ginibre Ensemble with a time-rescaled dispersion. The initial condition for the determinant (21) is D 0 (z, w) = |z| 2 + |w| 2 N , for the effective potential (18) φ 0 (z, w) = log(|z| 2 + |w| 2 ), for the velocity v 0 (z, w) =w/(|z| 2 + |w| 2 ), and for its modulus ν 0 (z, r) = r/(|z| 2 + r 2 ), where r = |w|, respectively. We solve the inviscid 2+1 Burgers equation (43)  place in the w complex plane however the position on z-plane defines the initial condition. The z-variable acts as a spectator on the eigenvalue plane "observing" the evolution in the perpendicular w-direction. We identify the cone-like caustic surface (left plot) whose apex is located at r = 0 and at critical time τ c = |z| 2 . This surface is the boundary along which the characteristics start to cross, making the Burgers solution multi-valued. This ambiguity develops for τ > τ c as the expanding eigenvalue boundary "swallows" the spectator at z 0 . We find the position of the shock line as a locus distanced equally from the caustic surface at each given time (the Rankine-Hugoniot condition). Since our problem is radially symmetric, the shock is positioned exactly at r = 0 and starts from the critical time τ c . Therefore although the dynamics takes place in the whole r space, the shock is always confined to the "physical" r = 0 region. Moreover if we confine it to r = 0, the spectator at z stays on the shock line at every time τ > τ c . We conclude that inside the bulk of the spectrum (i.e. the non-holomorphic sector), the observer is constantly on the shock line.
Having identified the positions of the shock we can write down a solution of (48) for the reduced 1+1 Burgers equation for our initial conditions ν 0 (z, r) = r/(|z| 2 + r 2 ). It reads It is an implicit algebraic equation for ν = ν(z, w, τ ). It can be rewritten as a cubic equation. Solutions for different τ is plotted in Fig. 3. Rather than showing the solution for the modulus ν we show a cross section of the vector field which has two symmetric branches ±ν. The right plot shows the solution inside the caustic surface with an unphysical branch depicted as dashed lines. For r = 0, that is on the z-plane, we have The boundary |z| 2 = τ is found by the sewing condition of zero and non-zero solutions. From ν we readily obtain the eigenvector correlation function (15) in the large N limit We can also find the second gradient g = ∂ z φ of the effective potential. To this end we can use the equation (29) which in the limit N → ∞ simplifies to ∂ τ g = ∂ w |v| 2 . For our initial condition g 0 (z, r) =z/(|z| 2 + r 2 ) it gives g(z, τ ) = 1/z for τ < |z| 2 z/τ for τ > |z| 2 .
Now we can use the Gauss law (10) to obtain the limiting eigenvalue density It is given by a uniform distribution on a disc of radius √ τ . For τ = 1, results (51,53) reproduce results of [24,13], respectively. We could alternatively find the same formulas by directly applying (38) and (39).
The maximizer r * can be found by solving the constraint (40): We find φ(z, τ ) = ln |z| 2 for τ < |z| 2 ln τ + |z| 2 τ − 1 for τ > |z| 2 . (55) The last formula was obtained for τ = 1 in [32]. So far we have discussed the limit N → ∞, in which as the effective potential φ and its gradients are equal to the electrostatic potential and the quaternionic Green's function (30). What about finite N ? For finite N the order of calculating the average and the logarithm in (17) and (18) matters , so φ is at best only an approximation of Φ for large but finite N . So clearly our method is not able to get insight into the finite N corrections. Surprisingly as we show below, the behavior of the diffusion kernel D(z, w, τ ) on the z complex plane, that is for w = 0 reveals the same type of finite size effects as known from exact calculations of the eigenvalue density for Ginibre ensemble [24]. For our initial conditions the diffusion kernel (25) is On the z-complex plane that is for r = 0 this integral significantly simplifies and it can be for large N calculated using the saddle point method. There are three saddle points When τ approaches |z| 2 from above, the two points r ± approach each other along the real axis to coalesce for τ = |z| 2 at r = 0. When τ becomes smaller than |z| 2 the two points move on the imaginary axis. Near the critical value τ = |z| 2 it is convenient to introduce a rescaled parameters and use them in the calculations. One finds that for large N the diffusion kernel behaves as where erfc is the complementary error function. This result holds not only on the z-complex plane but also sufficiently close to the z-complex plane that is for r approaching zero as N −3/4 or faster: r = O(N −3/4 ). Indeed in this case the argument of the function I 0 is a finite number and this function behaves as a constant for large N . The error function behavior has the same form as the finite size expression for the eigenvalue density function. This indicates that there may be some deeper relation between the diffusion kernel (averaged extended characteristic polynomial) and the two-point kernel known from the considerations of the underlying determinantal process for Ginibre matrices. We work out this analogy at heuristic level in Appendices B,C.

Spiric case
We now consider diffusion initiated from a diagonal matrix X 0 = diag(a, ..., a, −a, ..., −a) with the same number of a and −a eigenvalues. In the course of time two initial eigenvalue "islands", initially concentrated around ±a, expand to collide at some critical time τ c as shown in Fig. 4. The initial condition corresponds to D 0 (z, r) = (r 2 + |z − a| 2 ) N 2 (r 2 + |z + a| 2 ) N 2 ,   or if we write it for the modulus of the velocity field (47) ν 0 (z, r) = 1 2 r r 2 + |z − a| 2 + 1 2 r r 2 + |z + a| 2 . (60) In Fig. 5 we show charecteristics for three different values of the observer position z, including z = 0 where the collision takes place. As we can see the dynamics in the r plane behaves qualitatively in the same way as for the Ginibre case. In particular , characteristics form the same type of caustic surfaces with a shock line present after some critical time τ c = |a 2 −z 2 | 2 |a| 2 +|z| 2 . This critical time can be found by sewing the solutions ν or by finding the position of the caustic surface apex as a function of z and τ . In Fig. 6 we show a comparison of the caustic surfaces for Ginibre and spiric evolution for τ = 0.5. For this inititial condition the solution (48) of the inviscid Burgers equation takes the form 2ν = r + τ ν |z + a| 2 + (r + τ ν) 2 + r + τ ν |z − a| 2 + (r + τ ν) 2 . (61) For r = 0 the solution can be written explicitly where S a = τ 2 + 4Z 2 a , Z a =za + zā. The symbol S stands for the interior of the spiric section defined by the contour Generally spiric section is a curve obtained by intersection of a torus and a plane parallel to its rotational symmetry axis. We plotted spiric curves as contours around the scatter plots for eigenvalue densities in Fig.  4. The eigenvector correlation function inside the spiric section S is We can also calculate the diagonal element of the Green's function using the equation (29) that in the limit N → ∞ simplifies to ∂ τ g = ∂ z ν 2 . For the initial condition as g 0 = 1 2 z+ā r 2 +|z+a| 2 +z −ā r 2 +|z−a| 2 the solution reads with a constant c(z, a) =ā 2Za obtained by the sewing condition along the spiric section S. We use the Gauss law to obtain the spectral density as The same results can be obtained from the calculations of the effective potential φ using the equations (38), (39) and (40) for r * : Outside S we get and inside S φ in (z, τ ) = 1 2 ln We recall that S a = τ 2 + 4Z 2 a , Z a =za + zā, as defined after (62). In Fig. 7 we show a comparison of theoretical predictions for the limiting eigenvalue density (66) and the eigenvector correlation function (64) with Monte-Carlo simulations. The agreement is very good.
We complete the discussion of the spiric case by deriving a finite N formula for the diffusion kernel D(z, r, τ ) We are interested in the behavior close to the origin z = 0, r = 0 for τ close to the collision time. Without loss of generality we perform calculations for a = 1. In this case the collision time is τ = 1. We set r = 0 and zoom into the vicinity of the critical region where the saddle points merge. After expanding the solution (70) we obtain an asymptotic formula which is plotted in Fig. 8. Again, the result holds not only for r = 0 but more generally for r = O(N −3/4 ).

Non-normal Ginibre case
We consider now diffusion initiated from a matrix (X 0 ) ij = αδ i,j−1 . The matrix X 0 has all eigenvalues equal to 0. The initial eigenvalue distribution coincides with the one for the Ginibre case but afterwards the eigenvalue density obeys completely different evolution. Three snapshots of this evolution are shown in Fig. 9. The initial distribution concentrated initially at zero instantaneously expands to a circle of radius |α|, which corresponds to the pseudospectrum of the matrix. In the course of evolution the support of the density takes the form of the growing annulus. After a finite time τ = |α| 2 the inner radius of the annulus shrinks to zero and eigenvalues fill up a full disk.
Let us show this by direct calculations. The matrix M = (z − X 0 )(z − X 0 ) † + r 2 (35) has for our choice  of X 0 a tridiagonal form with a = |z| 2 + r 2 + |α| 2 , b = −zα and d = |z| 2 + r 2 . The determinant of this matrix can be calculated explicitly: , where ∆ = a 2 − 4|b| 2 and a ± = 1 2 (a ± ∆). The initial effective potential for large N is where a = |z| 2 + r 2 + |α| 2 . We have neglected 1/N terms which disappear in the limit N → ∞. The value r * which maximizes the expression in the Hopf-Lax formula (26) can be calculated from the equation (34). For r = 0 we get where we used the notation T α = τ 2 + 4|α| 2 |z| 2 . The boundary of the annulus A is given by the radii We see that at the beginning of the evolution, τ ≈ 0, the annulus is infinitely narrow forming a onedimensional circular pseudospectrum. On the other hand for τ → |α| 2 the inner radius tends to zero and the annulus becomes a disk. Inserting r * to (26) we obtain the potential for the three regions of the annulus. The spectral density is obtained from the Poisson equation (36) by differentiation the effective potential twice on the support of the annulus and zero otherwise. Using (39) we obtain the eigenvector correlation function to the support of the annulus In Fig. 10 we compared the prediction of the two formulas given above with numerical simulations. The agreement is good. Deviations are observed only close to the boundaries and can be attributed to finite-size effects. It is easy to check that all the expressions transform to the corresponding expressions for the Ginibre evolution when α → 0. Let us shortly discuss finite N effects for the averaged extended characteristic polynomial. We are looking for an universal behaviour near the origin for τ → 0. Around this space-time point an instantaneous transition from N eigenvalues positioned at the origin to the ring of radius |α| happens. In our example the characteristic polynomial is given explicitly as We consider the following scaling around the origin Inserting these formulas to the equation above we find an asymptotic function for large N and for |α| = xN −1/6 : We see that it reveals a non-perturbative divergent character near x = 0. For the normal Ginibre case we do not find the second term and therefore we identify it as a consequence of the non-normal initial condition.
Concluding this section, we note that the formula (78) was obtained in [22] for an initial matrix with N initial equidistant eigenvalues lying on a circle of radius |α|. This result is not surprising since the matrix X C 0 is unitarily equivalent to a circulant matrix which in turn differs from the non-normal initial matrix considered in this paper by just one element whose effect can be neglected in the large N limit. We also mention that finite N spectral densities and large deviations have been recently studied in [33].

Spectral stochastic equations
An interesting issue arises, if the stochastic dynamics of the elements of non-hermitian random matrices corresponds to some stochastic equation for the spectra of the underlying matrix. Such phenomenon holds in the case of Gaussian Unitary Ensembles (GUE), as well as in the case of Circular Unitary Ensembles (CUE), as shown in the seminal paper by Dyson [2]. In the simpler case of GUE, the corresponding Langevin equation for real eigenvalues reads where dB i reflects the Brownian dynamics and the second, "drift" term comes from the Van der Monde determinant. We neglect the optional Ornstein-Uhlenbeck term −λ i dt on the r.h.s. of (86) since it only contributes to freezing the diffusing end-points of the spectra in the stationary limit. The corresponding Smoluchowski-Fokker-Planck equation written for the resolvent G(z, τ ) of the spectral density ρ(λ, t), takes (after rescaling the time τ = N t and in the large N limit) the form of complex, inviscid Burgers equation, i.e. the eq. (2) [9,34]. Equations for non-hermitian Gaussian ensemble obtained in this work exhibit structural similarity to the hermitian "Burgulence", hence a question arises, if some stochastic dynamics for complex eigenvalues does exist as well. Recently, Osada [17] had examined a stochastic two-dimensional Coulomb system defined as where dB (2) i is a two-dimensional (complex) Brownian walk. In particular, he has shown that such interacting Brownian motion equipped with trivial initial condition (all λ i put to zero at t = 0) leads to limiting distribution of λ i representing a uniform disk, therefore resembling the Ginibre Ensemble spectrum.
From the point of view of our analysis this result is curious, since we have proven that the dynamics of eigenvalues is intimately connected to the dynamics of eigenvectors, which seem to decouple completely from "hypothetical" eigenvalues of Ginibre ensemble in (87). We suggest the resolution of this puzzle.
It is well known, that in so-called Normal Random Matrix model [35] with axially symmetric potentials all correlation functions can be expressed in terms of holomorphic functions of a single variable. Moreover, the exact integrability of such models can be linked to (2 + 1)-dimensional Burgers equations. Quite remarkably, in the case of the potential V (z,z) = |z| 2 , the correlations of the Normal Random Matrix model are identical to the correlations of the Ginibre Ensemble [36]. We therefore conjecture, that the stochastic equation (87) corresponds rather to Gaussian Normal Random Matrix model than to the Ginibre Ensemble.
Because normal matrices are diagonalizable by a single unitary transformation, alike the hermitian matrices, the eigenvectors decouple from the eigenvalues which explains the lack of these degrees of freedom in (87) and its formal similarity to hermitian stochastic equation (86). Based on this and the significance of eigenvectors presented in this paper we find it highly probable that the Ginibre-like dynamics of (87) is accidental and proper stochastic equation governing the dynamics of non-hermitian random matrices is not known.
The speculative links between our approach for generic complex matrices and non-Gaussian Normal Random Matrix models represent a challenge, which we plan to address in the forthcoming publications.

Conclusions
We have shown that a consistent description of non-hermitian Gaussian ensemble requires the knowledge of the detailed dynamics of co-evolving eigenvalues and eigenvectors. Moreover, the dynamics of eigenvectors seems to play a superior role (at least in the N → ∞ limit) and leads directly to the inference of the spectral properties. This is a dramatically different scenario as compared to the standard random matrix models, where the statistical properties of eigenvalues are of primary importance, and the properties of eigenvectors are basically trivial due to the their decoupling from the spectra. We have considered examples of Ginibre, spiric section and non-normal Ginibre where the formulas for spectral density, the 1-point eigenvector correlation function, electrostatic potential and universal functions were obtained and positively crosschecked with numerical simulations. By studying the dynamics of characteristics we anticipated the novel universality obtained for the spiric case as an error function type. We obtained compact formulas for both spectral density and eigenvector correlator for which the latter unraveled a promising determinantal structure. The diffusion equation, as an equation exact for finite N , was used mainly to obtain the universal behavior.
We conjecture that the hidden dynamics of eigenvectors discovered by us and described for the Gausssian non-hermitian ensemble, is a general feature of all non-hermitian random matrix models, and has to appear systematically in 1/N expansion.
Our formalism could be exploited to expand the area of application of non-Hermitian random matrix ensembles within problems of growth, charged droplets in quantum Hall effect and gauge theory/geometry relations in string theory beyond the subclass of complex matrices represented by normal matrices.
One of the challenges is an explanation, why, despite being so different, the Smoluchowski-Fokker-Planck equations for hermitian and non-hermitian random matrix models exhibit structural similarity to simple models of turbulence, where so-called Burgers equation plays the vital role, establishing the flow of the spectral density of eigenvalues in the case of the hermitian or unitary ensembles and the flow of certain eigenvector correlator in the case of non-hermitian ensembles.
We believe that our findings will contribute to understanding of several puzzles of non-hermitian dynamics, alike extreme sensitivity of spectra of non-hermitian systems to perturbations [24,27]. We also hope that the quaternion extension used in our paper may help to understand better the mathematical subtleties of the measure of non-hermitian operators [37].

Acknowledgments
MAN, JG, ZB, WT and PW are supported by the Grant DEC-2011/02/A/ST1/00119 of the National Centre of Science. MAN thanks Neil O'Connell for pointing out the reference [17] and appreciates discussions with Neil O'Connell, Boris Khoruzhenko, Andu Nica, Dima Savin, Roger Tribe and Oleg Zaboronski. The authors are also grateful to Charles Bordenave for pointing the relevant reference [33] on mathematical subtleties of normal Gaussian models. JG thanks prof. Thomas Guhr for the warm hospitality at the University Duisburg-Essen where part of this work was done.
Appendix A. Derivation of the diffusion equation (20) In this appendix we demonstrate that the averaged extended characteristic polynomial (19) obeys the diffusion equation (20). The determinant in (19) can be expressed with help of Grassmann variables as with the object T G given by where η,η, ξ andξ are Grassmann variables and X ij = x ij + iy ij . With the help of heat equation (5) for the joint probability density function P (X, τ ) we obtain where we integrated by parts twice and we used On the other hand we have We see that the expressions on the right hand side of (A.3) and (A.7) are identical and thus we have 1 N ∂ ww D(z, w, τ ) = ∂ τ D(z, w, τ ) (20). As a side remark we note that this calculation can be almost verbatim repeated for the averaged extended "inverse" characteristic polynomial The majority of results discussed in this paper are confined to the "physical" region where r → 0. One can however keep the parameter r nonzero. In the Coulomb gas interpretation, this deformation introduces a complex nonlinear interaction between the eigenvalues of unknown interpretation.
We consider the Ginibre case. The expression in the Hopf-Lax equation (26) φ(z, r, τ ) = max is maximized by r * obeying a cubic equation The corresponding spectral density ρ(z, r, τ ) and the eigenvector gradient ν(z, r, τ ) are shown in Fig. B.11. There are no critical points for r = 0 and instead we see a smooth crossover between two phases. This is in accordance with the fact that the shock is present only on the r = 0 plane. Consider the formula for the Applying the Newton binomial formula to |z| 2 + r 2 N an using the following integral representation of Laguerre polynomials Appendix C. The kernel structure We argue that in the case of Ginibre matrices, the characteristic polynomial D in the r → 0 limit has essentially the same information as the microscopic kernel of the underlying determinantal process. This stems from the observation made in [39] where the author connects the n-point correlation function of a Ginibre matrix model to a random matrix QCD partition function with n quarks with appropriately decreased matrix size: det(K N (z i , z j )) i,j=1...n = c(z 1 , ..., z n ) where det(K N (z i , z j )) i,j=1...n = n i=1 Trδ 2 (z i − X) X N is the correlation function averaged over N dimensional matrix and c denotes a known z dependent proportionality factor.
We study closer the n = 1 case because for this parameter the prefactor c is the weight function w(z) and the averaged determinants are exactly the characteristic polynomial D in the r → 0 limit. Therefore, in this particular case, we obtain a on-diagonal kernel formula K N (z, z) = C N w(z)D (N −1) (z, r = 0, τ ), (C.2) with some numerical constant C N . The off-diagonal kernel is not so easily obtainable by the above formula however we make an educated guess based on the symmetry of arguments and the result of Akemann and Vernizzi [38]. We write the full kernel as To complete the argument demonstrating that the full information resides in the characteristic polynomial D we will find an a priori unknown weight function w(z) present in formula (C.3) from the D alone. This is done by using decomposition in the biorthogonal basis where X i , Y i are biorthogonal polynomials with respect to an unknown measure W d 2 zW (z)X i (z)Ȳ j (z) = g ij . (C.7) The expansion terms c ij and the bi-orthogonality matrix g ij are inverses g ij = c −1 ij . We proceed by finding a basis X i , Y j in which c kl are diagonal. Then we infer the formula for W by considering its moments (C.7). We start with a formula (C.5) which is already in a diagonal form with The moments are therefore given by By assuming the radial symmetry and setting |z| 2 = p we obtain The kernel is therefore equal to (C.14)