A Multi-Layer Extension of the Stochastic Heat Equation

Motivated by recent developments on solvable directed polymer models, we define a ‘multi-layer’ extension of the stochastic heat equation involving non-intersecting Brownian motions. By developing a connection with Darboux transformations and the two-dimensional Toda equations, we conjecture a Markovian evolution in time for this multi-layer process. As a first step in this direction, we establish an analogue of the Karlin-McGregor formula for the stochastic heat equation and use it to prove a special case of this conjecture.


Introduction and Summary
We consider the fundamental solution to the stochastic heat equation in one dimension with initial condition u(0, x, y) = δ(x − y), where W is a standard L 2 (R) cylindrical Brownian motion and the term u dW is interpreted as an Itô integral. Alternately, we can write this as whereẆ denotes space-time white noise associated with the Brownian motion W . This is a distribution-valued Gaussian field on [0, ∞) × R with covariance function For more background see, for example [7,24,41,53]. We will write stochastic integrals Then the solution u(t, x, y) to (1) is given by the chaos expansion u(t, x, y) = p(t, x, y) where k (t) = {0 < t 1 < · · · < t k < t} and For each t > 0 and x, y ∈ R, the expansion (3) is convergent in L 2 (W ). It satisfies (1) in the sense that it satisfies the integral equation u(t, x, y) = p(t, x, y) + For more details see, for example [41]. The chaos series representation for the solution can be viewed as an expansion of the Feynman-Kac formula. To emphasize this interpretation, the solution (3) is often written as a generalised Wiener functional where the expectation is with respect to a Brownian bridge (X s , s ≤ t), which starts at x at time 0 and ends at y at time t [24]. The factor p(t, x, y) arises because of this use of the bridge in the representation, corresponding to the δ-function initial condition, rather than the more usual free Brownian motion. The Wick exponential, see Sect. 3, is denoted by exp . It is important to note that the integral t 0 dW (s, X s ), representing the integral of the white noiseẆ along the path s → (s, X s ), is not well-defined, and equation (5) is purely formal.
This solution arises as a scaling limit of partition functions associated with lattice directed polymers, in the 'intermediate disorder' regime [1,31]. In fact, as suggested by (5), it can be interpreted directly as a partition function for the continuum random polymer [2]. On the other hand, it has been known for some time that h = log u arises as the scaling limit of the height profile of the weakly asymmetric simple exclusion process, at least for equilibrium initial conditions [8], recently extended to include the present initial condition in [3]. With this 'surface growth' interpretation, h is understood to be the physically relevant solution (also known as the Cole-Hopf solution) to the KPZ (Kardar-Parisi-Zhang) equation [28] with 'narrow wedge' initial condition.
In a remarkable recent development the exact distribution of the random variable u(t, x, y) has been determined. This has been acheived via two distinct approaches. One of these [3,[43][44][45][46] uses the asymmetric simple exclusion process approximation together with recent work by Tracy and Widom [49][50][51][52] in which exact formulas have been obtained for that process using an approach based on the Bethe ansatz. Another approach [11,[16][17][18] is based on replicas, where the moments of u(t, x, y) are related to the attractive δ-Bose gas and computed via the Bethe ansatz. These developments indicate that there is an underlying integrable structure behind the KPZ and stochastic heat equations that is not yet fully understood.
On the other hand, it has recently been found that there are exactly solvable discrete (or semi-discrete) directed polymer models [15,32,[34][35][36]47,48], yielding yet another approach. For these models, there is a direct connection to integrable systems (specifically to the quantum Toda lattice), which one might hope to understand in the continuum scaling limit. We will describe here one of the main results of the paper [34], which provided the motivation for the present work.
Define an 'up/right path' in R × Z to be an increasing path that either proceeds to the right or jumps up by one unit. For each sequence 0 < t 1 < · · · < t N −1 < t we can associate an up/right path π from (0, 1) to (t, N ) which has jumps between the points (t i , i) and (t i , i + 1), for i = 1, . . . , N − 1, and is continuous otherwise. Let B(t) = (B 1 (t), . . . , B N (t)), t ≥ 0, be a standard Brownian motion in R N and define and the integral is with respect to Lebesgue measure on the Euclidean set {(t 1 , . . . , t N −1 ) ∈ R N −1 : 0 < t 1 < · · · < t N −1 < t} (6) of all such paths. This is the partition function for the model, which was introduced in [35]. In [34] an explicit integral formula is obtained for the Laplace transform of the distribution of Z N (t), via the following 'multi-layer' construction. For n = 1, 2, . . . , N , define Z N n (t) = e E(π 1 )+···+E(π n ) dπ 1 . . . dπ n , where the integral is with respect to Lebesgue measure on the set of n-tuples of nonintersecting (disjoint) up/right paths with respective initial points (0, 1), . . . , (0, n) and respective end points (t, N − n + 1), . . . , (t, N ). Here, the notion of Lebesgue measure arises by identifying this set of paths as a suitable subset of R n(N −n) , as in (6). Define . The relevance of this construction is analogous to the role of the RSK correspondence in the study of last passage percolation and longest increasing subsequence problems; in this setting it is based on a geometric variant of the RSK correspondence. The main result in [34] is the following.
, is a diffusion process in R N with infinitesimal generator given by where ψ 0 is a (particular) ground state eigenfunction of the quantum Toda lattice Hamiltonian The function ψ 0 is a (class-one) G L(n, R)-Whittaker function. The diffusion process with generator L is the analogue of Dyson's Brownian motion in this setting. The law of the logarithmic partition function X N 1 (t) = log Z N (t), which can now be seen as the analogue of the top line (or largest eigenvalue) in the Dyson process, is determined as a corollary. Similar results have been obtained in [15,36] for a lattice directed polymer model with log-gamma weights, which was introduced by Seppäläinen [47]. In that setting, the eigenfunctions of the quantum Toda lattice (also known as Whittaker functions) continue to play a central role, as does the geometric lifting of the RSK correspondence which was introduced and studied in the papers [30,33].
Motivated by these developments, in this paper we introduce continuum versions of the partition functions Z N n (t), which we expect will play an important role in our understanding of the integrable structure that appears to lie behind the KPZ and stochastic heat equations. The continuum partition functions are defined as follows. For n = 1, 2, . . . , t ≥ 0 and x, y ∈ R, define where R (n) k is the k-point correlation function for a collection of n non-intersecting Brownian bridges which all start at x at time 0 and all end at y at time t, as defined in Sect. 2 below. Note that Z 1 = u is the solution of the stochastic heat equation defined by (3) above. To explain the above definition we note that, just as in (5), we can formally write (9) as where (X 1 s , . . . , X n s , 0 ≤ s ≤ t) denote the trajectories of n non-intersecting Brownian bridges which all start at x at time 0 and all end at y at time t. These should be compared with the partition functions (7). The first main result of this paper is that the continuum partition functions Z n (t, x, y) are well-defined. Theorem 1.2. The series (9) is convergent in L 2 (W ).
The proof will given in Sect. 4. Define u 1 = u and, for n ≥ 2, u n = Z n /Z n−1 . From Theorem 1.1 we expect that, for each fixed t > 0, the process u t (x) = {u n (t, 0, x), n ∈ N}, indexed by x ∈ R, is a diffusion process in R N which is a scaling limit at the edge of the diffusion with generator (8), just as the multi-layer Airy process introduced in [38] is a scaling limit at the edge of Dyson's Brownian motion. Note that it is not at all clear a priori that this process should have the Markov property: in fact, this is quite a remarkable property and can be regarded as the analogue, in this setting, of Pitman's celebrated '2M − X ' theorem. Moreover, for large t, the process u t (x), x ∈ R should rescale (after taking logarithms) to the multi-layer Airy process. At present, we only know this to be the case for the one-dimensional distributions of the first layer: it has been shown in the papers [3,46] that the distribution of log[u(t, 0, x)/ p(t, 0, x)] (which is independent of x) converges in a suitable scaling limit to the Tracy-Widom distribution. This result on the first layer has been tentatively extended to the finitedimensional distributions by Prolhac and Spohn [39,40]. As such, it is natural to regard the process u t (x), x ∈ R as the analogue of the multi-layer Airy process in the setting of the KPZ and stochastic heat equations. Similarly, for fixed t > 0, the random sequence {u n (t, 0, 0), n ∈ N} can be regarded as the analogue of the Airy point process. There are many things to understand in these directions, especially the law of the process (u t (x), x ∈ R). For further recent progress in this direction, see [10,14].
It is also natural to ask about the evolution of u t as t varies: after all, this is an extension of the process (u(t, 0, ·), t > 0), which evolves according to the stochastic heat equation and (more or less as an immediate consequence) has the Markov property. This Markov property can also be seen by means of the Feynman-Kac representation (5), by splitting the exponential into two factors: Using the fact that a Brownian bridge over the time interval [0, t + u], conditioned on its value at the intermediate time t, splits into two independent bridges, one obtains the flow property of the stochastic heat equation: where u(u, z, y; t) denotes the solution the stochastic heat equation driven by the shifted white noise W (t + ·, ·). The Markov property follows immediately. Turning to the multilayer process we see such an argument fails dramatically: the definition of u t involves non-intersecting Brownian bridges with the same starting and ending points but when we condition on intermediate values we obtain bridges with distinct starting and end points and consequently quantities that depend on the white noise in a more general manner. In view of this there is really no reason to expect the extended process (u t , t > 0) to have the Markov property. Nevertheless, we will present a number of results which strongly indicate that it does indeed, somewhat remarkably, have the Markov property. In fact, we believe that for each n, the process {(u 1 (t, 0, ·), . . . , u n (t, 0, ·)), t ≥ 0} is Markov and give a proof of this claim in the case n = 2 (see Corollary 6.4).
In order to develop a better understanding of the continuum partition functions Z n (t, x, y), we consider their analogues when the space-time white-noise potential W (t, x) is replaced by a smooth time-varying potential φ(t, x). For example, we can take φ to be in the Schwartz space E of rapidly decreasing smooth (C ∞ ) functions on R + × R. For each n = 1, 2, . . . , t > 0 and x, y ∈ R, define where once again (X 1 s , . . . , X n s , 0 ≤ s ≤ t) denote the trajectories of n non-intersecting Brownian bridges which all start at x at time 0 and all end at y at time t.
In the following, we drop the superscript φ and write Z n = Z φ n . Now, as above, we define u = u 1 = Z 1 and, for n ≥ 2, u n = Z n /Z n−1 . Note that, by Feynman-Kac, u(t, x, y) satisfies the heat equation with initial condition u(0, x, y) = δ(x − y). In this setting we will show, using a straightforward generalisation of the Karlin-McGregor formula (see Propositions 3.1 and 3.2 below), that, for n ≥ 2, where c n,t = t n(n−1)/2 ( n−1 j=1 j!) −1 . But this determinant is the Wronskian associated with the solutions u, ∂ x u, . . . , ∂ n−1 x u of the heat Eq. (12). It follows (see for example [4]) that the functions u n are in fact Darboux transformations and satisfy the coupled system of heat equations with initial conditions u n (0, x, y) = δ(x − y). These equations are not immediately meaningful if we replace the smooth potential φ by space-time white noise. However they do suggest that, for each n, the multi-layer process (in the white noise setting) has a Markov evolution. In order to pursue this, there are two possible directions one could take. The first, and most obvious one, is to try to make sense of these evolution equations with white-noise replacing φ. To this end, it is helpful to consider the following change of variables which, as it happens, also reveal some deeper structure as we shall see. Set τ 0 = 1 and, for n ≥ 1, define , a n = τ n−1 τ n+1 τ 2 n .
We will show (see Propositon 3.7) that the evolution of the a n is given by ∂ t a n = 1 2 ∂ 2 y a n + ∂ y [a n ∂ y log u n ].
It seems to be the case that these evolution equations will make sense as stochastic partial differential equations in the white-noise setting [21]. We note the following connection to the 2D Toda equations. We will show (see Lemma 3.6) that a n = ∂ xy log τ n , for each n. Thus, if we define, for n ≥ 1, then a n = e q n+1 −q n and the q n satisfy the 2D Toda equations ∂ xy q n = e q n+1 −q n − e q n −q n−1 , n ≥ 1, with the convention that q 0 = +∞. In this notation, from (15), the time-evolution of the a n is given by ∂ t a n = 1 2 ∂ 2 y a n + ∂ y [a n ∂ y q n ].
A second, and quite different, approach one can take in order to understand (and prove) the Markov property of the multi-layer process starts from the discussion we had above regarding flow property and the Feynman-Kac representation for the stochastic heat equation. This suggests considering a natural extension of the partition functions Z n to collections of non-intersecting Brownian paths that start at distinct points and end at distinct points. In the first instance we continue to replace the space time white noise by the smooth potential φ. Set and denote the interior of n by • n . For each t > 0 and for x = (x 1 , . . . , x n ) and y = (y 1 , . . . , y n ) in • n , define where (X 1 s , . . . , X n s , 0 ≤ s ≤ t) denote the trajectories of a collection of nonintersecting Brownian bridges which start at positions x = (x 1 , . . . , x n ) and end at positions y = (y 1 , . . . , y n ) at time t, and p * n (t, x, y) is the transition density of a Brownian motion in n killed when it first hits the boundary, given by the Karlin-McGregor formula [29], In fact, K n satisfies the following generalisation of the Karlin-McGregor formula, which we record here as it may be of independent interest.
Then we find that (see Theorem 3.4) for t > 0, x ∈ R and y ∈ • n , In the case φ = 0, this identity reduces to the fact that the volume of GT (y) is proportional to (y). Now, if the analogous formula holds in the white-noise setting, with K n (t, x, y) defined by the chaos series (44) below, then this implies the Markov property of the multi-layer process This is because for each x ∈ R and for each n, the process (M n (t, x1, ·), t ≥ 0) clearly has the Markov property, in fact a version of the flow property holds, see Corollary 6.2. One important ingredient in this argument is the Karlin-McGregor formula (18). A priori, its not clear whether or not this will hold in the white-noise setting. In Sect. 5 we show that it does. Then in Sect. 6 we argue that, modulo technical considerations, this indicates that the analogue of the integral formula (21) should hold in the white-noise setting. We give a proof in the case n = 2, just to demonstrate that it can be done. The main technical issue concerns the continuity of M n (t, x, y) at the boundary of n × n . A proof of the existence of an almost surely continuous extension to the boundary based on Kolmogorov's criterion would be long and technical and will not be included in the present work. Here we satisfy ourselves with a continuous extension in L 2 , which allows us to establish the analogue of the integral formula (21), and hence the Markov property of the multi-layer process, in the special case n = 2. We remark that an interesting consequence of our proof of the L 2 -continuity property is that the ratio of two solutions to the stochastic heat equation is in H 1 . In fact, such ratios have recently been shown (in a slightly different setting, for smooth initial data and periodic boundary conditions) by Hairer [21] to be in C 3/2− . The index 3/2 is consistent with our expectation that the continuum partition functions Z n (t, x, y) are locally Brownian in the space variable y.
We conclude this section with some remarks on the RSK interpretation. As remarked above, the multi-layer construction presented in this paper is based on a geometric lifting of the RSK correspondence, so it is natural to consider such an interpretation in the continuum setting. The analogue of RSK in the context of smooth potentials is the mapping In the language of RSK, is the Q-tableau, and their common 'shape' is the sequence {u n (t, 0, 0), n ≥ 1}. We note the following symmetry, which corresponds to a well known symmetry property of the RSK correspondence.
Similarly, in the white noise setting, we define where W [0,t] denotes the restriction of W to [0, t]×R. As explained above, we expect that, , which are conditionally independent given their starting position {u n (t, 0, 0), n ≥ 1}. This would be the analogue, in this setting, of Pitman's '2M − X ' theorem. Since the first version of the present paper appeared, substantial progress has been made in this direction by Corwin and Hammond [14], where a natural candidate for this infinitedimensional diffusion process has been constructed.
The outline of the paper is as follows. In the next section we provide some background on non-intersecting Brownian motions and their bridges. Following this we study the analogue of the partition functions when the space-time white noise is replaced by a smooth time-varying potential. In this setting we establish a connection with Darboux transformations of solutions to the heat equation, which give rise to evolution equations for the multi-layer process of partition functions. These equations are not directly meaningful in the white noise setting, but suggest that the multi-layer process has a Markovian evolution. We also give proofs of the integral formula (21) above, the evolution Eq. (15) and remark on the connection with the 2D Toda equations. In Sect. 4 we present the proof of Theorem 1.2 on the existence of the continuum partition functions in the whitenoise setting. In Sect. 5, we show that the Karlin-McGregor formula (5.3) holds in the white-noise setting. In Sect. 6 we consider the evolution of the multi-layer process in the white-noise setting and give a proof of the Markov property for n = 2.

Non-intersecting Brownian Motions
Non-intersecting Brownian motions play a large role in this paper, and we record here definitions and facts concerning them that will be useful to us.
Recall that and denote the interior of n by • n . Standard n-dimensional Brownian motion killed on exiting • n has transition densities given by the Karlin-McGregor forumla [29], Dyson Brownian motion is obtained as a Doob-h transform of this killed process by the harmonic function (x) = i< j (x i − x j ). It has transition densities with respect to the measure (y) 2 dy given by Lemma 2.1. For each t > 0, the transition density q n (t, x, y) extends continuously to a uniformly bounded strictly positive function on n × n .
Proof. By the Harish-Chandra/Itzykson-Zuber formula [22,25], we can write where 1/c n = n−1 j=1 j!, X and Y are diagonal matrices with entries given by the vectors x and y, and the integral is with respect to normalised Haar measure on the group of n × n unitary matrices. By bounded convergence, the RHS defines a continuous function on n × n , and is bounded by (2π) −n/2 t −n 2 /2 c n . The strict positivity follows from the strict positivity of the integrand.
The semigroup property thus also extends to the boundary, by continuity and dominated convergence. Moreover, it follows from the representation (25) that q n (t, x, y) (y) 2 dy defines a probability measure on n for every t > 0 and x ∈ n . Consequently, Dyson Brownian motion can be started from any point on the boundary of n . In fact, as was shown by Cépa and Lépingle [12], it almost surely never subsequently returns to the boundary.
Let (H t , t ≥ 0) be a standard Brownian motion in the space of n × n Hermitian matrices. Then the vector of ordered real-valued eigenvalues of H t evolves as Dyson Brownian motion. This can be verified by deriving the transition density for the eigenvalues using the Harish-Chandra/Itzykson-Zuber formula. Alternatively, following Dyson's original approach, applying Itô's formula shows that, denoting the vector of eigenvalues by (X 1 t , X 2 t , . . . , X n t ), the following stochastic differential equations are satisfied.
where β i , i = 1, 2, . . . , n are a collection of independent standard one-dimensional Brownian motions. Note that these equations hold even if the intitial value H 0 of the Hermitian Brownian motion has repeated eigenvalues, in which case the Dyson Brownian motion is starting from the boundary of n . This can been seen by the following argument. Applying Itô's formula from some strictly positive time onwards ( recalling the process of eigenvalues does not visit the boundary) we obtain for t ≥ , where β i, are Brownian motions. Now for < the increments of β i, − +· and β i, agree, and by virtue of this consistency there exist Brownian motions β i starting from 0 such that β i, t = β i +t − β i for all > 0. Now returning to (27), writing it using β i , and letting tend down to 0 gives (26) as desired even in the case when the Dyson Brownian motion starts from the boundary.
We can construct bridges for Dyson Brownian motion using the standard Markovian framework, see for example Proposition 1 of [19]. Specifically given points x and y belonging to n , we define the bridge from x at time 0 ending at y at time t, to be a process (X s , 0 ≤ s ≤ t) whose law over [0, s], for any s < t, is absolutely continuous with respect to that of Dyson Brownian motion starting from x, with a density Note that this is well-defined as the denominator is strictly positive, by Lemma 2.1 above. We will also refer to this bridge, somewhat informally, as a collection of nonintersecting Brownian bridges. In the special case x = x1, y = y1, where x, y ∈ R and we denote by 1 ∈ R n the vector with all coordinates equal to 1, the process is also often referred to as a watermelon. The correlation function (9) is defined to be the sum over i 1 , i 2 , . . . , i k of the (continuous) probability densities of (X i 1 t 1 , . . . , X i k t k ) with respect to Lebesgue measure evaluated at (x 1 , x 2 , . . . , x k ). Correlation functions for non-intersecting Brownian bridges with arbitrary starting and ending positions, which appear in Sect. 5 below, are defined analogously.

Darboux Transformations
In this section we replace the white noise potentialẆ by a smooth potential φ, which we assume for convenience to be in the Schwartz space E of rapidly decreasing smooth (C ∞ ) functions on R + × R.
For each n = 1, 2, . . . , t > 0 and x, y ∈ R, define where (X 1 s , . . . , X n s ), 0 ≤ s ≤ t, denote the trajectories of n non-intersecting Brownian bridges which all start at x at time 0 and all end at y at time t. On one hand, these are the analogues of the partition functions Z n introduced in the previous section with the white noiseẆ replaced by a smooth potential φ. On the other, they are directly related to the Z n by the formula where and exp (W (φ)) is the Wick exponential of W (φ) defined by In other words, as a function of φ, Z φ n (t, x, y) is the S-transform of the white noise functional Z n (t, x, y) (see, for example [24]). To see that (30) holds, on the RHS replace Z n (t, x, y) by the series (9) and exp (W (φ)) by its Wiener chaos expansion; computing the expectation of the product of these two series, and making use of the orthogonality of multiple integrals of different order, we obtain For the remainder of this section we will only consider the case of smooth potential φ. For notational convenience we will drop the superscript and simply write Z n (t, x, y) = Z φ n (t, x, y). By the Feynman-Kac formula, u := Z 1 satisfies the heat equation with initial condition u(0, x, y) = δ(x − y).
where c n,t = t n(n−1)/2 c n and 1/c n = n−1 j=1 j!. We will prove this via a generalisation of the Karlin-McGregor formula.
For each t > 0 and for x = (x 1 , . . . , x n ) and y = (y 1 , . . . , y n ) in • n , define where (X 1 s , . . . , X n s , 0 ≤ s ≤ t), denote the trajectories of a collection of nonintersecting Brownian bridges which start at positions x 1 , . . . , x n and end at positions y 1 , . . . , y n at time t, and p * n (t, x, y) is the transition density of a Brownian motion in n killed when it first hits the boundary, given by the Karlin-McGregor formula (23).
Proof. According to the Feynman-Kac formula, K n satisfies the equation with Dirichlet boundary conditions on ∂ n and initial condition K n (0, x, y) = i δ(x i − y i ). Moreover it is the unique solution to this initial-boundary value problem which vanishes as |y| → ∞ uniformly for t in compact intervals. This follows from a variant of the maximum principle (see, for example [20, Chapter 2, Theorem 2]), which applies in this setting since φ is bounded and continuous.
On the other hand, det[Z (t, x i , y j )] n i, j=1 satisfies the same initial-boundary value problem and vanishes as |y| → ∞ uniformly for t in compact intervals. So the identity follows by uniqueness.
We remark that the same argument can be applied to more general expressions than (33) in which the potential (s, x) → φ(s, x i ) is replaced by a bounded continuous potential ψ(t, x) which is a symmetric function of the coordinates of x; we will make use of this fact in the proof of Theorem 5.3 below. Z n (t, a, b)  p(t, a, b) , x, y) .

Proof of Proposition 3.1. It is immediate from the definitions that
.
Define u n (t, x, y) recursively by Z n = u 1 u 2 · · · u n .
The Eq. (36) follow from Proposition 3.1 together with known properties of Darboux transformations of solutions to one-dimensional heat equations with time-varying potentials, see for example [4,5]. For completeness, we will include a direct proof of Proposition 3.3 just after the statement of Proposition 3.7 below. The coupled heat equations of Proposition 3.3 are not immediately meaningful if we replace the smooth potential φ by space-time white noise. However they do suggest that the multi-layer process (in the white noise setting) t, x, ·), . . . , Z n (t, x, ·), t ≥ 0) has the Markov property. In the following, we introduce a natural extension of the Z n which will play an important role in our understanding of the Markov property when we return to the white noise setting.
This extends continuously to the boundary of n × n ; by Proposition 3.2 and [54, Rather surprisingly, we will now show that the apparently richer object M n (t, x1, ·) is, for a fixed x ∈ R and t > 0, given as a function of (Z 1 (t, x, ·), . . . , Z n (t, x, ·)).
as required.
Lemma 3.6. Let g(x, y) be a smooth function and define W 0 = 1, W 1 = g and, for . Suppose that W n is strictly positive for all n ≥ 1 and define T n = W n−1 W n+1 /W 2 n . Then the following identities hold: proving the first identity. Essentially the same argument shows that, for any k ≥ 1, This implies that In particular, and so on, yielding the second identity. Now, by Lemma 3.5, Applying Lemma 3.5 again, using . Now apply Lemma 3.5 again, using T 2 = ∂ y (∂ y (∂ 2 x g/g)/T 1 ), and so on, to obtain the third identity.
Note that, by Lemma 3.6, a n = ∂ xy log τ n .
Thus, if we define, for n ≥ 1, then a n = e q n+1 −q n and, by (41), the functions q n satisfy the 2D Toda equations ∂ xy q n = e q n+1 −q n − e q n −q n−1 , n ≥ 1, with the convention that q 0 (t, x, y) = +∞ for all x, y ∈ R.
The time-evolution of the a n is given by the following proposition.

Proof (of Propositions 3.3 and 3.7)
. First we note that the initial condition u n (0, x, y) = δ(x − y) follows immediately from the definition (29) of Z n . We will verify the equations (36) and (42) simultaneously by induction over n. Write h n = log u n and note that (36) is equivalent to For n = 1, the Eq. (36) holds by hypothesis. Thus (t, y).
It follows that a 1 = ∂ xy h 1 satisfies as required. Now assume the induction hypothesis (with two parts): ∂ t u n = 1 2 ∂ 2 y u n + φ(t, y) + ∂ 2 y log Z n−1 / p n−1 u n and ∂ t a n = 1 2 ∂ 2 y a n + ∂ y [a n ∂ y h n ].
Note that u n+1 = (t/n)a n u n and Thus, ∂ t u n+1 = 1 n a n u n + t n [a n ∂ t u n + u n ∂ t a n ] = 1 t u n+1 + t n a n 1 2 ∂ 2 y u n + φ(t, y) + ∂ 2 y log Z n−1 / p n−1 u n + t n u n 1 2 ∂ 2 y a n + ∂ y [a n ∂ y h n ] as required. Note that this implies By (41) we can write a n+1 = a n + ∂ xy h n+1 . Thus, using the second part of the induction hypothesis again, ∂ t a n+1 = 1 2 ∂ 2 y a n+1 + ∂ y [a n ∂ y h n ] + ∂ y [(a n+1 − a n )∂ y h n+1 ] + ∂ 2 y a n .

Proof of Theorem 1.2
We return now to the white noise setting, denoting by u(t, x, y) the solution to the stochastic heat equation (1) with initial condition u(0, x, y) = δ(x − y). In this section we will show that for each n ≥ 2, Z n (t, x, y) defined by (9) is convergent in L 2 (W ) or, equivalently, The first step is to show that this is equivalent to Ee L < ∞, where L is the total intersection local time between two independent copies of the system of n non-intersecting Brownian bridges. Let (X i s , 0 ≤ s ≤ t, i = 1, . . . , n) be a collection of nonintersecting bridges which all start at x at time 0 and all end at y at time t, and let to be the semimartingale local time process at 0 of (X i − Y j )/2, as defined for example in [42,Chapter VI]. The total intersection local time is defined by L = L t = n i, j=1 L i j t .
Lemma 4.1. In the above notation, Ee L is given by (43).
Proof. We show, by induction on k, that the kth term of (43) is equal to E[L k t ]/k!. First recall that R (n) 1 ((t 1 , x 1 )) is the sum over i of the marginal probability densities for each X i t 1 evaluated at x 1 . Consequently (R (n) 1 ((t 1 , x 1 )) 2 can be expressed as the sum over all i and j of the joint density of (X i t 1 , Y j t 1 ) evaluated at (x 1 , x 1 ). Integrating over x 1 and t 1 gives an expression which agrees with that for E[L] given by the occupation time formula (see [42,Chapter 6]).
Similarly R ((t 1 , x 1 ), . . . , (t k , x k ))) 2 can be expressed as the sum over all i 1 , . . . , i k and j 1 , . . . j k of the joint density of This is integrated over all x i and t i to give an expression for the kth term of (43). On the other hand we can derive the same expression for E[L k t ] by writing L k t = k t 0 L k−1 s d L s , and evaluating the expectation of this using Proposition 3 and Lemma 1 of [19], together with inductive hypothesis.
Next will show that, in fact, all exponential moments of L t are finite. First note that L t = A + B, where A is the intersection local time on the time interval [0, t/2] and B is the remainder. Thus, by Cauchy-Schwartz, it suffices to show that A and B each have finite exponential moments of all orders. Now, on the time interval [0, t/2], the joint law of the bridges X = (X 1 , X 2 , . . . X n ) and Y = (Y 1 , Y 2 , . . . , Y n ) is absolutely continuous to the law of two independent copies of Dyson Brownian motion with Radon-Nikodym density a product of two factors each given by (28) with s = t/2, x = x1, y = y1. By Lemma 2.1, this Radon-Nikodym density is a bounded random variable. A similar statement holds on [t/2, t] after time reversal. It therefore suffices to show that, for two independent Dyson Brownian motions the total intersection local time has finite exponential moments of all orders. This is established as a special case of the following lemma which controls the intersection local time for arbitrary starting points of the Dyson Brownian motions. Proof. The processes X and Y satisfy a system of SDEs Recall that L i j is the local time process of (X i −Y j )/2 at zero. Consequently Tanaka's formula, see [42,Chapter 6], states that, where Thus, it suffices to show that each of the random variables has finite exponential moments of all orders, each of which can be bounded arbitrarily close to 1, uniformly in u and v, by choosing t small. The third of the above random variables is the absolute value of a Gaussian random variable with mean zero and variance 2t and so the desired property holds straightforwardly. (To see this, note that the stochastic integral is a continuous martingale with quadratic variation process 2t, hence is a Brownian motion by Lévy's characterisation theorem [42, Chapter IV, Theorem 3.6].) To control the exponential moments of the first and second of the above random variables, we recall that Dyson's Brownian motion arises as the process of eigenvalues of a Brownian motion (H t , t ≥ 0) in the space of n × n Hermitian matrices. Theñ H t = H t − H 0 defines a Hermitian Brownian motion starting from the zero matrix, and applying Weyl's eigenvalue inequalities to the sum H 0 +H t we deduce that where σ t is the spectral radius ofH t . Since the Brownian motion (H t , t ≥ 0) can be taken to not depend on u, and the spectral radius ofH t has exponential moments that approach 1 as t tends 0, this gives the desired control for |X i t − u i |. The same argument applies of course to |Y j t − v j |. The fourth and fifth random variables are essentially the same, so it remains to show that, for each i, has finite exponential moments which can be made arbitrarily close to 1. We will prove this by induction over i. For i < j, define First we note that has finite exponential moments which may be bounded as desired. Now, since ξ 1 = ξ 12 + · · · + ξ 1n and each term is non-negative, this implies that ξ 1 j ≤ ξ 1 and hence that ξ 1 j has exponential moments satisfying the same bound for each j = 2, . . . , n. Now Thus ξ 2 and ξ 23 , . . . , ξ 2n all have exponential moments of all orders which may be bounded as desired. Similarly, the fourth term is handled similarly, and so on.
We note that, by the stationarity of space-time white noise, the law of Z n (t, x, y)/ p(t, x, y) n does not depend on x, y and the above proposition implies that C t := E|Z n (t, x, y)/ p(t, x, y) n | 2 < ∞.

Karlin-McGregor Type Formula in the White Noise Setting
For n = 1, 2, . . . and x, y ∈ • n , define where R Proof. We need to show that E[e L 1 A ] < ∞ where L is the total intersection local time between two independent sets of n independent Brownian bridges started at positions x and ending at positions y at time t, and A is the event that each set is non-intersecting. Note that P(A) > 0 since x, y ∈ • n . So it suffices to show that Ee L < ∞. By considering pairwise intersection local times and applying Hölder's inequality one obtains Ee L ≤ Ee n 2 √ t/2R where R is the local time at zero of a standard Brownian bridge on [0, 1], which has the Rayleigh distribution P(R > r ) = e −r 2 /2 , r > 0 (see, for example [37]).
In fact, Proposition 4.2 yields the following stronger statement.

Proposition 5.2.
For each t > 0, there is a constant d t < ∞ such that for all x, y ∈ • n , E|K n (t, x, y)| 2 ≤ d t p * n (t, x, y) (x) (y).
Proof. As in the proof of Theorem 1.2, let L be the total intersection local time between two independent copies of the system of n non-intersecting Brownian bridges started at positions x and ending at positions y at time t. Then (cf. Lemma 4.1) E|K n (t, x, y)| 2 = p * n (t, x, y) 2 Ee L .
whereÊ denotes the expectation with respect to a Dyson Brownian motion started at x. By Lemma 2.1 and Proposition 4.2 this is bounded by d t /q n (t, x, y) where d t is a constant independent of x, y. The second term is treated similarly, and the statement of the proposition follows.
Before proceeding to the Karlin-McGregor formula, which is the main result of this section, we recall the approach of Bertini and Cancrini to the stochastic heat equation. In [7], these authors make sense of the formal Feynman-Kac representation (5) by means of smoothing the white noise. For κ > 0 introduce the mollified white noise W κ defined by where δ κ (·) is the centered Gaussian density of variance 1/κ. Then the analogue of (5) is then meaningful and defines random variables u κ (t, x, y).
Moreover, for each (t, x, y), and any p ≥ 1, Proof. Recall that E denotes Schwartz space of rapidly decreasing smooth (C ∞ ) functions on R + × R. Let φ ∈ E, multiply both sides by exp W (φ) and take expectations. The left-hand side becomes K φ n (t, x, y), defined earlier by the Feynman-Kac expression (33). The right-hand side becomes which will we now argue is also given by (33), and since φ ∈ E is arbitrary the statement of the theorem will follow. Consider the quantity Replacing each u κ (t, x i , y j ) by its Feynman-Kac representation, using Fubini, and integrating over W we obtain where for each permutation σ , B σ is a bridge of standard n-dimensional Brownian motion starting from x and ending at σ y = (y σ (1) , . . . , y σ (n) ), and ψ κ is given by Because ψ κ is invariant under permutations of the coordinates, as remarked earlier the argument for the Karlin-McGregor formula of Proposition 3.2 allows this to be rewriten as where X = (X 1 s , . . . , X n s , 0 ≤ s ≤ t) denotes a collection of non-intersecting Brownian bridges which start at positions x = (x 1 , . . . , x n ) and end at positions y = (y 1 , . . . , y n ) at time t, and p * n (t, x, y) is the transition density of a Brownian motion in n killed when it first hits the boundary.
We now let κ tend to infinity. On the one hand, from their definitions and (45), we have C κ n (t, x, y) → C n (t, x, y). On the other hand, since φ κ converges uniformly to φ, and X does not visit the boundary of n , the right-hand side of equation (46) converges to that of (33). This is justified by an application of the Dominated Convergence Theorem using the fact that sup x∈R l x t , where l x t denotes the local time of a one-dimensional Brownian motion at level x, has finite exponential moments.
It is known [7] that for each x the solution to the stochastic equation u(t, x, y) admits a version that is almost surely continuous in t and y and moreover is strictly positive. It follows from the above theorem that for each x ∈ R n , K n (s, x, z) admits a version that is almost surely continuous in t and z. Define K n (t, z, y; s) via the chaos expansion (44) but with the shifted white noiseẆ (s + ·, ·). For each y ∈ R n , this similarly admits a version that is almost surely continuous in t and z. In the following we assume that we are using these versions. We conclude this section with the following.
Proof. In the above notation, we first claim that for each x, y ∈ • n and κ > 0, it holds almost surely that To see this, we use the Feynman-Kac representation, from [7, (2.17)], where the expectation is with respect to a Brownian bridge b starting at x and ending at y at time t, and F κ s (b), 0 ≤ s ≤ t is an almost surely continuous, strictly positive, multiplicative functional of b. It follows, by a standard path-switching argument (see, e.g., [27,Sect. 1.2]), that where the expectation is with respect to a collection of n non-intersecting Brownian bridges (X 1 , . . . , X n ) started at positions x and ending at positions y at time t. Indeed, multiplying both sides of (47) by 1 A (y), where A is a measurable subset of n , and integrating with respect to y ∈ n , gives where B is a standard Brownian motion in R n started at x and T is the first exit time of B from n . To prove (47) it suffices to show that (48) holds for every measurable A ⊂ n . Let us write and note that (48) is equivalent to for each i = 1, . . . , n − 1. Fix i and definẽ where s i denotes the adjacent transposition (i, i + 1). By the strong Markov property,B has the same law as B. Moreover, since F κ t is a multiplicative functional, we also have (B) = (B). Hence, and it follows that as required. The result now follows from (47), letting κ → ∞.

On the Evolution of the Z n in the White Noise Setting
In this section we discuss the analogue of Theorem 3.4 in the white noise setting, and the implication that (Z 1 (t, x, ·), . . . , Z n (t, x, ·)), t ≥ 0 is Markov.
We expect, but will not prove here, that for each t > 0, has a version which almost surely extends continuously to a strictly positive function on n × n . In particular, for each t > 0, almost surely, extends continuously to n × n this representation naturally defines the extension of M n (t, x, y) to n × n . Our task is show continuity in L 2 (W ). For this it is enough to show that is continuous. Now, as in the proof of Theorem 1.2 this expectation is equal to where L is the total intersection local time of two independent sets of non-intersecting bridges, X and X say, starting at positions x = (x 1 , . . . where the kernel p(·, ·) can be written as a product of transition densities for nonintersecting Brownian motions, and is thus seen to be continuous. From this it follows by a dominated convergence argument that depends continuously on (x, x , y, y ) also. To deduce the continuity of we must show that the difference can be made uniformly small for z within compact sets by choosing δ small enough. Applying the Cauchy-Schwartz inequality this amounts to showing that can be made uniformly close to 1. This follows from Proposition 4.2, noting that the joint law of X and X over the time interval [0, δ] is absolutely continuous to that of a pair of independent Dyson Brownian motions, with a density, specified by equation (28) that, by virtue of Lemma 2.1, is bounded uniformly for (x, x , y, y ) belonging to compact sets.
We can now establish the identity (50) in the case n = 2. Theorem 6.3. For x ∈ R and y ∈ • 2 , M 2 (t, x1, y) u(t, x, y 1 )u(t, x, y 2 ) = 1 y 1 − y 2 Proof. It is known [7] that the solution to the stochastic equation u(t, x, y) admits a version that is almost surely continuous in t and y and moreover is strictly positive. We assume in the following that we are using this version. In particular, having fixed t, x and y 1 > y 2 we let A (x) be the event {u(t, x, z) > for all z ∈ [y 2 , y 1 + 1]}. Then as ↓ 0 we have P(A (x)) ↑ 1. By Theorem 5.3, for x, y ∈ • 2 , [u(t, x 1 , y 1 )u(t, x 2 , y 2 ) − u(t, x 1 , y 2 )u(t, x 2 , y 1 )]. Hence, Writing y = (z + h, z) where h > 0 this becomes Integrating this equation with respect to z over the interval [y 2 , y 1 ] we obtain y 1 y 2 M 2 (t, x, (z + h, z)) u(t, x 2 , z + h)u(t, x 2 , z) dz Now let h tend to zero. By the continuity of u(t, x 1 , ·) and u(t, x 2 , ·) the RHS converges almost surely to We want to identify the limit of the LHS. Consider By virtue of the uniform continuity in L 2 of the mappings (z 1 , z 2 ) → M 2 (t, x, (z 1 , z 2 )) and z → u(t, x 1 , z) these integrals tend to zero as h ↓ 0, and consequently E tends to 0 in probability. Thus we have proven Next let x = (x + h, x) and let h ↓ 0. The LHS of (56) can be rewritten as (y 1 − y 2 ) M 2 (t, (x + h, x), (y 1 , y 2 )) u(t, x, y 1 )u(t, x, y 2 ) ; as h ↓ 0 this converges in probability to (y 1 − y 2 ) M 2 (t, x1, (y 1 , y 2 )) u(t, x, y 1 )u(t, x, y 2 ) .
On the other hand, if we consider We remark that the identity (56) shows that the ratio of two solutions to the stochastic heat equation is in H 1 ; in fact, such ratios have recently been shown (in a slightly different setting) by Hairer [21] to be in C 3/2− . Corollary 6.4. For each x ∈ R, the process (Z 1 (t, x, ·), Z 2 (t, x, ·)), t ≥ 0 has the Markov property.
Proof. Fix times 0 ≤ s < t. Suppose that F = F(Z 1 (t, x, ·), Z 2 (t, x, ·)) is a bounded random variable (depending on the random fields at a finite number of points). We wish to show that the conditional expectation given the white noise W [0,s] of this random variable is measurable with respect to the random fields (Z 1 (s, x, ·), Z 2 (s, x, ·)). To see this, note firstly that by Lemma 6.1 the same random variable F is a function of the fields (M 1 (t, x1, ·), M 2 (t, x1, ·)). Now the flow property for M n obtained in Corollary 6.2 together with the independence of M(t, ·, ·; s) from W [0,s] implies the conditional expectation E F|W [0,s] has a version G of the form G = G (M 1 (s, x1, ·), M 2 (s, x1, ·)). We have M 1 (s, x1, ·) is proportional Z 1 (t, x, ·), and more profoundly, by the preceeding Theorem, M 2 (s, x1, ·)) can be expressed in terms of Z 1 (t, x, ·) and Z 2 (t, x, ·). Thus we see that G is of the required form.
It would be interesting to understand the evolution of the multi-layer process in terms of a system of stochastic partial differential equations. Motivated by the evolution equations obtained in Sect. 2 in the case of a smooth potential, it is natural to consider and to try to make sense of the system of equations ∂ t a n = 1 2 ∂ 2 y a n + ∂ y [a n ∂ y log u n ], where Z n = u 1 . . . u n . For recent progress in this direction, in the case of smooth initial data and periodic boundary conditions, see [21].