Solution of the Kolmogorov equation for TASEP
Mihai Nica, Jeremy Quastel, Daniel Remenik
TL;DR
This paper offers a straightforward proof that the TASEP transition probabilities satisfy the Kolmogorov backward equation, also extending the method to the PushASEP particle system, thereby clarifying the mathematical structure of these stochastic processes.
Contribution
It provides an elementary proof confirming the TASEP transition probabilities solve the Kolmogorov backward equation and extends this approach to PushASEP.
Findings
Confirmed TASEP transition probabilities solve the Kolmogorov backward equation
Extended the proof technique to PushASEP system
Simplified understanding of TASEP and PushASEP dynamics
Abstract
We provide a direct and elementary proof that the formula obtained in [MQR17] for the TASEP transition probabilities for general (one-sided) initial data solves the Kolmogorov backward equation. The same method yields the solution for the related PushASEP particle system.
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Solution of the Kolmogorov equation for TASEP
Mihai Nica
Department of Mathematics
University of Toronto
40 St. George Street
Toronto, Ontario
Canada M5S 2E4
,
Jeremy Quastel
Department of Mathematics
University of Toronto
40 St. George Street
Toronto, Ontario
Canada M5S 2E4
and
Daniel Remenik
Departamento de Ingeniería Matemática and Centro de Modelamiento Matemático (UMI-CNRS 2807)
Universidad de Chile
Av. Beauchef 851, Torre Norte, Piso 5
Santiago
Chile
Abstract.
We provide a direct and elementary proof that the formula obtained in [MQR17] for the TASEP transition probabilities for general (one-sided) initial data solves the Kolmogorov backward equation. The same method yields the solution for the related PushASEP particle system.
1. Introduction
The totally asymmetric simple exclusion process (TASEP) consists of particles on the lattice performing totally asymmetric nearest neighbour random walks with exclusion: Each particle independently attempts jumps to the neighbouring site to the right at rate , the jump being allowed only if that site is unoccupied. We will always consider initial conditions in which there is a rightmost particle; this always remains so, and the positions of the particles can be denoted . The dynamics of the first particles is independent of the rest, so the infinite system clearly makes sense. If we let
[TABLE]
denote the label of the rightmost particle which sits to the left of, or at, at time , then the TASEP height function associated to can be defined, for , by
[TABLE]
which fixes . The height function itself is a simple random walk path with if there is a particle at at time and if there is no particle at at time . The dynamics of the height function is that local max’s become local min’s at rate ; i.e. if then at rate , the rest of the height function remaining unchanged. The rate of decrease is
[TABLE]
where is the indicator function of a local max. Hence the TASEP height function can be seen as a simple discretization of the Kardar-Parisi-Zhang (KPZ) equation,
[TABLE]
Because of its amenability to computations, TASEP has become the most popular model in the KPZ universality class.
-particle TASEP was solved by [Sch97] using the coordinate Bethe ansatz. The transition probabilities are given by a determinant,
[TABLE]
with
[TABLE]
where is any positively oriented simple loop which includes and . A direct proof of this formula is not difficult and can be obtained in a couple of pages (see, for example [MQ17]).
On the other hand, one is generally less interested at a later time in the exact positions of the particles, or the exact height function, but in the joint distribution of the height function at a finite number of points , where is fixed, and is large, or infinite. Or, what amounts to the same thing, one would like to compute, for some sequence of indices and any vector , the joint probability
[TABLE]
where the subscript in denotes the initial condition. In Schütz’s formula, this would involve a sum over the positions of the other particles. The resulting formula is not useful, and, in particular, not conducive to the limit.
This was overcome by [Sas05, BFPS07], who were able to rewrite the right hand side of (1.4) as a (signed) determinantal point process on a space of Gelfand-Tsetlin patterns. This allowed them to employ the Eynard-Mehta technology [EM98] to conclude that (1.6) can be written as the Fredholm determinant which (in principle) could be obtained from the operators and acting on given by
[TABLE]
where and, for ,
[TABLE]
where is any positively oriented simple loop including the pole at but not the one at . The functions , , are defined implicitly by: I. The biorthogonality relation ; II. is a polynomial of degree at most in for each .
Except for a few very special choices of initial data, the solution for the was not discovered until [MQR17]. Without discussing their exact form, we state the final result after performing some manipulations to get a nice formula. Let be the real analytic extension of , let be the backwards discrete difference operator, and observe that is invertible (with inverse ). Define
[TABLE]
the first one being defined for all and the second one for . Here is the semigroup of a Poisson process with jumps to the left at rate , which we may think of as an integral operator on with kernel ; this formula is actually valid for all and it defines the whole group of operators (in particular is invertible, with inverse ). Define also, for ,
[TABLE]
where is defined to be the hitting time of the strict epigraph of the ‘curve’ \big{(}X_{0}(k+1)\big{)}_{k=0,\dotsc,n-1} by a discrete time random walk with transition probabilities . Then, as proven in [MQR17], the kernel from (1.8) can be expressed as
[TABLE]
Now let denote the generator of TASEP. It acts on bounded cylinder (i.e. depending on finitely many coordinates) functions , where , as follows:
[TABLE]
where we take the convention in the indicator . Since is a cylinder function, the sum has only finitely many non-zero terms. If is a Markov process with generator , its transition probabilities satisfy the Kolmogorov backward equation, with final condition . The generator acts on the variable. In our case, we have a time homogeneous process and a generating family of sets and , so the backward equation for reads,
[TABLE]
with initial condition
[TABLE]
The solution of the Kolmogorov backward equation for a continuous time Markov chain on a countable state space is unique under the condition that the supremum over all states of the rate of leaving that state is finite (see, for example, Thm. 2.26 and Cor. 2.34 of [Lig10], where it is written for the point-to-point probabilities ; being countable, the probability measure is therefore prescribed, and satisfies uniquely the backward equation in , with initial condition ). Our state space is not countable; however, in TASEP (and PushASEP), the evolution of particles is unaffected by particles , . Hence the backward equation (1.15), (1.16) actually takes place on the countable set . Therefore we have
Theorem 1.1**.**
The unique solution of the Kolmogorov backward equation (1.15), (1.16) for TASEP is given explicitly by the Fredholm determinant of the kernel introduced in (1.8):
[TABLE]
However, a detailed proof along the historical lines sketched above would run about 30 pages. A natural and important question is whether one could just prove directly that the determinant satisfies the Kolmogorov equation. The purpose of this article is to provide such a proof, for TASEP and its variant PushASEP. We comment that the direct proof is relatively short, but far from obvious. It is not known if it can be obtained directly from the biorthogonal representation (1.8) without knowledge of the special form of the . It is also worth noting that there are examples (such as discrete time TASEP, see [MQR+]) with explicit for which the present short proof does not work.
TASEP has a rich history, some of which can be found in the introduction of [Sch97]. Since that paper there has been significant progress in exact solvability, with much interest stemming from TASEP’s role as one of the fundamental models in the KPZ universality class. There are many results; we mention here only a few besides the ones which fit directly into our story above. The one-point distribution of TASEP with step initial data (, ) was first solved in [Joh00] by exploiting a Toeplitz structure which is only available for that choice of initial condition; the main goal of that paper was to prove the now-classic fact that the fluctuations of the position of an appropriately chosen particle converge, as time goes to infinity, to the Tracy-Widom GUE distribution [TW94]. This was later extended in [Joh03] to the multi-point distributions of TASEP, which were shown there to converge to the Airy2 process [PS02]. Much work was devoted during the last fifteen years to extending this type of results to two other choices of initial data: periodic (, , which was the subject of the papers [Sas05, BFPS07] discussed above) and stationary ( corresponding to placing particles on according to a product measure) [IS04, BFS08, BFP10]. In the past ten years there has been a huge effort to extend some of these results to the partially asymmetric simple exclusion process (see, for example, [TW08, IS11, BCS14]). These still depend on very particular initial data. In the case of TASEP on a ring there is also a huge literature with some recent breakthroughs: [BL18], where asymptotics is done for a type of step initial data, and [Pro13, Pro14], which actually computes the entire spectrum.
2. Kolmogorov equation for TASEP
In this section we provide a short self-contained proof of Thm. 1.1. The strategy is to compute the two sides of (1.15) with given by (1.17) and given by (1.13), and check that they are equal.
We begin with the right hand side of (1.15). We will first consider the effect of moving a single particle, and then we will sum over all particles to obtain the effect of the generator . Fix a particle label and consider the original initial condition and the initial condition where particle is moved to the right by one, i.e. (for now it does not even matter whether or not this particle can actually be moved without violating the strict particle order condition). We will compare the kernel with these two initial conditions; for notational convenience we put tildes on top of all objects when they depend on the modified initial condition (e.g. refers to the kernel started from the initial condition ). By (1.13) we have
[TABLE]
while by (1.12) we have
[TABLE]
Recall that means the hitting time of the strict epigraph of \big{(}X_{0}(\ell+1)\big{)}_{\ell\geq 0}; is the same but with the modified . It is clear from the definitions that unless and (see Fig. 1 for a visual explanation), and thus the above expectation equals
[TABLE]
which by the Markov property we may write as with
[TABLE]
where the subscript in expectation on the right hand side indicates that the random walk is started now at at time and the superscript in indicates that the hitting time is now restricted to times in (we will use the same notation for other hitting times below). This means that
[TABLE]
Crucially, this is a rank one kernel (acting on ), and then the same is true of the extended kernel when thought of as acting on the extended space , so using the fact that for any rank one operator, we deduce that
[TABLE]
(here as in (2.6)). This difference corresponds to the bracket in the -th term of the sum (1.14) defining the action of on our function . The conclusion then is that if is the set of labels of particles which can be moved (i.e. those so that ), then equals the sum over of the above difference, i.e.
[TABLE]
Note that, as above, the sum only has finitely many non-zero terms, because for all (as is clear, for instance, from (2.5) and (2.6)). Finally we note if (i.e. if ), then since the random walk has increments , we have and thus . Thus summing over has no effect, and we have then
[TABLE]
Now we turn to the left hand side of (1.15). In general, if a kernel depends smoothly on a parameter one has . Thus
[TABLE]
and then in view of (2.9) all we need in order to deduce (1.15) is to check that
[TABLE]
We have from (1.13) that
[TABLE]
From (1.10) we get
[TABLE]
(we have used here that commutes with and since they all have Toeplitz kernels). This gives
[TABLE]
with
[TABLE]
Shifting the curve defined by up by one, the first term on the second line can be expressed as \mathbb{E}_{B_{0}=z_{1}}\big{[}\bar{\mathcal{S}}_{-t,n-\hat{\tau}}(B_{\hat{\tau}}-1,z_{2})\mathbf{1}_{\hat{\tau}<n}\big{]} with the hitting time of the strict epigraph of \big{(}X_{0}(m+1)+1\big{)}_{m\geq 0}, and thus
[TABLE]
Note that , so the difference inside the brackets vanishes for , while if and we have , so the difference vanishes again. We deduce that
[TABLE]
with as in (2.4) and given by the bracket in the middle line (see also the comment after (2.5)).
Shifting now the curve defining back to in the last expectation we get
[TABLE]
We have then
[TABLE]
Comparing this formula and (2.6), and since for by its definition (2.5), all that remains in order to prove (2.10) is that
[TABLE]
But from (1.11) we have , so applying on the right in (2.17) and using the definitions (2.5) and (2.15) of and we see that (2.17) is equivalent to
[TABLE]
for . By definition of , if we multiply the equation by then both sides are polynomials in , so it is enough to prove the equality for . But it is easy to check that
[TABLE]
(see [MQR17, Eqn. (LABEL:fp-eq:QQext)]), so for such , since is decreasing in then is also smaller than , and in (2.18), and we deduce that all the ’s in the identity can be replaced by . As a consequence, what we need to prove is that
[TABLE]
or simply
[TABLE]
which is easy to see since the walk takes Geom steps to the left: in fact, since and and may only differ on the event that the walk hits at time , then for , decomposing according to the first step,
[TABLE]
the case is even simpler. This yields the desired equality (2.17).
3. Initial condition
Now we check that the initial condition (1.16) is satisfied. It is simpler to deal first with the one-point case (), so we start there.
3.1. One-point case
Let
[TABLE]
We need to prove that
[TABLE]
Remark 3.1**.**
In the one point case, the initial condition (3.1) can be checked rather easily from the biorthogonal ensemble representation (1.8), as mentioned by one of the referees; the proof using this biorthogonalization is given below. However, since our goal in this paper is to prove (1.17) directly from the representation (1.13) (and thus avoid having to use the result of [MQR17] to connect this formula back to (1.8)), we provide the direct (but relatively more complicated) proof from (1.13) in this subsection.
From (2.6) we have with in the definition (1.9). Let be given by and , so that the left hand side of (3.1) equals
[TABLE]
where we have used the cyclic property of the Fredholm determinant. Note that on the determinant on the right hand side is finite dimensional, with the matrix inside being given by
[TABLE]
Assume first that . Then for all , which implies from (1.9) that for all and all such , so as desired. Assume otherwise that . From (1.9) it is straightforward to check that and thus the first column of is given by . But the biorthogonalization condition implies that . We deduce that the first column of equals zero, so vanishes as desired.
We turn now to the proof of (3.1) directly from (1.13). Consider first the case , which implies that in (3.1) we need to evaluate only for . By definition of the hitting time in (1.12) we have for , so we may write
[TABLE]
Now
[TABLE]
so, in particular, for . If and we have , so and the first term on the right hand side of (3.2) vanishes. For the second term we consider \bar{\mathcal{S}}^{\operatorname{\text{{epi}}}(X_{0})}_{0,n}(y,z_{2})=\sum_{k=0}^{n-1}\mathbb{E}_{B_{0}=y}\big{[}\bar{\mathcal{S}}_{0,n-k}(B_{k},z_{2})\mathbf{1}_{\tau=k}\big{]} with and note that implies , which cannot happen unless the starting location satisfies (because the walk takes downward steps of size at least ). But for such (and ) we have , and thus again we get . We have proved that
[TABLE]
and then the left hand side of (3.1) equals for , as desired.
Next we turn to the case , and focus on the column corresponding to the index in , that is \big{(}K^{(n)}_{0}(z,X_{0}(n))\big{)}_{z\leq a}. We claim that, in this case,
[TABLE]
If that is the case then the matrix has a column which is identically 0, and thus its determinant vanishes as required.
In order to prove (3.5), observe first that the argument we used to prove (3.4) actually shows that for all . Consider next the case . For the first term in (3.2), a simple residue computation using the contour integral formulas for and which follow from (1.10) and (1.11) with gives (remembering )
[TABLE]
Since we want to show that , then from (3.6) together with (3.2) we see that we need to show that
[TABLE]
We consider the two possibilities for separately. Assume first that . This means that the first particles are packed one next to the other, i.e. , and thus, since the walk starts below and takes downward steps of size at least , it cannot hit the strict epigraph of the curve by time , so as desired. The other possibility is that . Using we may express the kernel as \sum_{k=0}^{n-1}\mathbb{E}_{B_{0}=y}\big{[}\bar{Q}^{(n-k)}(B_{k},X_{0}(n))\mathbf{1}_{\tau=k}\big{]} and then note that, since inside the expectation is the location where the strict epigraph is hit, , so by (2.19) we may replace by there to get . But implies , which in turn implies , so the last sum equals simply , and then
[TABLE]
And the last sum can be extended to all ’s because if then, under our assumption and using (3.3), the first factor in the summand vanishes. This shows that, in this case, as desired.
Up to here we have proved (3.5) for . Consider finally the case . The exact same argument used in the last paragraph to prove (3.7) also shows
[TABLE]
In particular, in the case we get using (3.2) that which equals [math] because in fact (see [MQR17, Eqn. (LABEL:fp-eq:QinvQext)]). It remains to prove that the first term in (3.2) vanishes when , but this is straightforward, because in this case we have (using the contour integral formulas (1.10) and (1.11) as before)
[TABLE]
since .
3.2. Multi-point case
We let and think of as an operator-valued matrix with entries . The initial condition (1.16) will follow from the arguments in the one-point case and the identity
[TABLE]
which is (LABEL:fp-eq:checkQKK) in [MQR17]. Consider first the case for each . Here (3.4) gives , and then in the sum only terms with survive. On the other hand, we have for (any ) so for . This shows that , and since this holds for all , we get with a strictly upper triangular operator-valued matrix, which yields as desired.
Suppose next that for some , so we need to prove . We may assume without loss of generality that is the largest particle label among for which this inequality holds. We will focus on the -th column of , and more precisely on the subcolumn of this column of kernels corresponding to the index , that is v^{(j)}\coloneqq\big{(}v^{(j)}_{\ell}(z)\big{)}_{z\in\mathbb{Z},\,i=1,\dotsc,m} with . In this case (3.5) gives us for all , and then , so that using (3.8) we get
[TABLE]
The case is direct and yields . Otherwise, if , then necessarily by our choice of , and hence for we have so as before . The conclusion is that . But then has a column (namely ) which is identically [math], and thus as required.
4. PushASEP
We consider now PushASEP, a generalization of TASEP and the Toom model [DLSS91] introduced in [BF08]. Again we have particles on the lattice with positions . Each particle attempts to jump one step to the right at rate , with jumps being permitted only if the neighboring site is empty; this is the same TASEP dynamics considered above (except run at rate ). On the other hand, each particle has another (independent) exponential clock running at rate , and when it rings the particle jumps to the nearest vacant site on its left. Relabeling the particles in order to keep the ordering after such a jump, we may think of the effect of a particle jumping left as pushing all its left neighbors one step to the left. TASEP and the Toom model (or PushTASEP) are recovered from PushASEP by setting and , respectively.
The generator of this process can be written similary to the one for TASEP: letting denote the length of the block of nearest neighbor occupied sites lying to the left of (and including) particle , meaning that for and ,
[TABLE]
where again we take . As for the TASEP case we define
[TABLE]
again we have that satisfies the Kolmogorov backward equation
[TABLE]
Theorem 4.1**.**
The unique solution of (4.3), (4.4) is given by
[TABLE]
with
[TABLE]
where the operators and are now given by
[TABLE]
here and is defined similarly to after (1.11).
This result, which is new, can be proved along the same lines as the solution of TASEP given in [MQR17]; such a proof will appear in [MQR+] (where convergence to the KPZ fixed point will also be proved). Our goal here, as for the Kolmorogov equation for TASEP in Thm. 1.1, is to prove directly that the Fredholm determinant provides a solution to the Kolmogorov backward equation (4.3), (4.4).
Note first that for the formulas for and given in Thm. 4.1 are the same as those for TASEP, so the initial condition (4.4) follows from the TASEP proof. Next we claim that in order to prove that the Fredholm determinant solves (4.3) it is enough to consider the case of PushTASEP, i.e. , since the TASEP part of the dynamics was already handled in Sec. 2. To see this it is convenient to define with
[TABLE]
, , and defined analogously using , so that (4.3) will follow (after setting ) if we prove that . But since , it is actually enough to prove that
[TABLE]
The first of these is a version of the Kolmogorov equation for TASEP (1.1) where in the contour integrals in (1.10) and (1.11) there are respectively additional factors and inside the integrands, both independent of , and it is not hard to see that the proof in Sec. 2 works without any difficulty after this modification. We will prove below that (4.3) holds in the PushTASEP case , and inspecting the proof it is easy again to seen that this implies the second equation in (4.10). We conclude that it suffices to check the case , and we can also take as it is just a scaling factor.
Given this reduction, one may naturally wonder why the proof works for PushASEP and not ASEP. It just turns out that when one tries to adapt Schütz’s formula [Sch97] to include jumps to the left (for example, as done in [BF08]) the resulting formula ends up solving the push dynamics instead of the leftward TASEP dynamics.
Having made these comments, we begin the proof. Proceeding as for TASEP, let and consider the initial condition and the one where particle is moved one step to the left along with its neighboring block. Write for the length of this block. As for TASEP, we need to compute
[TABLE]
where now the tildes refer to the initial condition modified in this new way. We clearly have unless and . But since the particles for lie in a block, implies (because the walk always jumps down by at least one step). This means that the above expectation equals
[TABLE]
with
[TABLE]
(note that does not depend on but does, through ). Since this is a rank one kernel, the same arguments as in the TASEP case give
[TABLE]
where is defined exactly as in (2.6) but using the new versions of and (note that here acts on all particles, since all particles are allowed to move in the PushTASEP dynamics).
Now we consider the left hand side. As for TASEP, all we need to check in order to finish the proof is that
[TABLE]
Now we have
[TABLE]
and using (4.6) and (4.8) we have \tfrac{\mathrm{d}}{\mathrm{d}t}(\mathcal{S}_{-t,-n})^{*}=2(\mathcal{S}_{-t,-n}\big{)}^{*}\nabla^{+} and \tfrac{\mathrm{d}}{\mathrm{d}t}(\bar{\mathcal{S}}_{-t,n})^{*}=-2\nabla^{+}(\mathcal{S}_{-t,n}\big{)}^{*}. Therefore
[TABLE]
with
[TABLE]
The first term on the second line can be expressed as \mathbb{E}_{B_{0}=z_{1}}\big{[}\bar{\mathcal{S}}_{-t,n-\hat{\tau}}(B_{\hat{\tau}}+1,z_{2})\mathbf{1}_{\hat{\tau}<n}\big{]} with the hitting time of the strict epigraph of \big{(}X_{0}(m+1)-1\big{)}_{m\geq 0}, and thus
[TABLE]
Now we have , so the difference inside the brackets vanishes for , and it also vanishes when and . Thus equals
[TABLE]
with
[TABLE]
and
[TABLE]
We have then
[TABLE]
and so it remains to prove that for all and for (using as for TASEP that for ).
For fixed let denote the hitting time of the initial condition introduced above where the block starting at is moved one step to the left, so that . We obviously have , and thus and implies , so that . This gives
[TABLE]
Next we prove the identity for . Proceeding as for TASEP, we apply on the right of this identity to see that it is equivalent to
[TABLE]
Following again the TASEP argument, the two sides are of the form times a polynomial in , so it is enough to prove the equality for , in which case the identity becomes
[TABLE]
which is easy to prove again since the walk takes Geom steps to the left and .
Acknowledgements. MN and JQ were supported by the Natural Sciences and Engineering Research Council of Canada. DR was supported by Programa Iniciativa Científica Milenio grant number NC120062 through Nucleus Millenium Stochastic Models of Complex and Disordered Systems, by Conicyt Basal-CMM Proyecto/Grant PAI AFB-170001, and by Fondecyt Grant 1160174.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[BFP 10] Jinho Baik, Patrik L. Ferrari and Sandrine Péché “Limit process of stationary TASEP near the characteristic line” In Comm. Pure Appl. Math. 63.8 Wiley Subscription Services, Inc., A Wiley Company, 2010, pp. 1017–1070 DOI: 10.1002/cpa.20316 · doi ↗
- 2[BL 18] Jinho Baik and Zhipeng Liu “Fluctuations of TASEP on a ring in relaxation time scale” In Comm. Pure Appl. Math. 71.4 , 2018, pp. 747–813 DOI: 10.1002/cpa.21702 · doi ↗
- 3[BCS 14] Alexei Borodin, Ivan Corwin and Tomohiro Sasamoto “From duality to determinants for q 𝑞 q -TASEP and ASEP” In Ann. Probab. 42.6 , 2014, pp. 2314–2382 DOI: 10.1214/13-AOP 868 · doi ↗
- 4[BF 08] Alexei Borodin and Patrik L. Ferrari “Large time asymptotics of growth models on space-like paths. I. Push ASEP” In Electron. J. Probab. 13 , 2008, pp. no. 50, 1380–1418 DOI: 10.1214/EJP.v 13-541 · doi ↗
- 5[BFPS 07] Alexei Borodin, Patrik L. Ferrari, Michael Prähofer and Tomohiro Sasamoto “Fluctuation properties of the TASEP with periodic initial configuration” In J. Stat. Phys. 129.5-6 , 2007, pp. 1055–1080 DOI: 10.1007/s 10955-007-9383-0 · doi ↗
- 6[BFS 08] Alexei Borodin, Patrik L. Ferrari and Tomohiro Sasamoto “Transition between Airy 1 and Airy 2 processes and TASEP fluctuations” In Comm. Pure Appl. Math. 61.11 , 2008, pp. 1603–1629
- 7[DLSS 91] B. Derrida, J. L. Lebowitz, E. R. Speer and H. Spohn “Dynamics of an anchored Toom interface” In J. Phys. A 24.20 , 1991, pp. 4805–4834 URL: http://stacks.iop.org/0305-4470/24/4805
- 8[EM 98] Bertrand Eynard and Madan Lal Mehta “Matrices coupled in a chain. I. Eigenvalue correlations” In J. Phys. A 31.19 , 1998, pp. 4449–4456 DOI: 10.1088/0305-4470/31/19/010 · doi ↗