A Determinantal structure of the correlation functions

# Anisotropic growth of random surfaces in 2+1 dimensions

## Abstract

We construct a family of stochastic growth models in dimensions, that belong to the anisotropic KPZ class. Appropriate projections of these models yield dimensional growth models in the KPZ class and random tiling models. We show that correlation functions associated to our models have determinantal structure, and we study large time asymptotics for one of the models.

The main asymptotic results are: (1) The growing surface has a limit shape that consists of facets interpolated by a curved piece. (2) The one-point fluctuations of the height function in the curved part are asymptotically normal with variance of order for time . (3) There is a map of the -dimensional space-time to the upper half-plane such that on space-like submanifolds the multi-point fluctuations of the height function are asymptotically equal to those of the pullback of the Gaussian free (massless) field on .

## 1 Introduction

In recent years there has been a lot of progress in understanding large time fluctuations of driven interacting particle systems on the one-dimensional lattice, see e.g. [4, 2, 1, 47, 30, 29, 38, 45, 25, 48, 46, 8, 9, 11, 7, 10]. Evolution of such systems is commonly interpreted as random growth of a one-dimensional interface, and if one views the time as an extra variable, the evolution produces a random surface (see e.g. Figure 4.5 in [43] for a nice illustration). In a different direction, substantial progress have also been achieved in studying the asymptotics of random surfaces arising from dimers on planar bipartite graphs, see the review [32] and references therein.

Although random surfaces of these two kinds were shown to share certain asymptotic properties (also common to random matrix models), no direct connection between them was known. One goal of this paper is to establish such a connection.

We construct a class of two-dimensional random growth models (that is, the principal object is a randomly growing surface, embedded in the four-dimensional space-time). In two different projections these models yield random surfaces of the two kinds mentioned above (one reduces the spatial dimension by one, the second projection is fixing time). We partially compute the correlation functions of an associated (three-dimensional) random point process and show that they have determinantal form that is typical for determinantal point processes.

For one specific growth model we compute the correlation kernel explicitly, and use it to establish Gaussian fluctuations of the growing random surface. We then determine the covariance structure.

Let us describe our results in more detail.

### 1.1 A two-dimensional growth model

Consider a continuous time Markov chain on the state space of interlacing variables

 S(n)={{xmk}k=1,…,mm=1,…,n⊂Zn(n+1)2|xmk−1

can be interpreted as the position of particle with label , but we will also refer to a given particle as . As initial condition, we consider the fully-packed one, namely at time moment we have for all , see Figure 1.1.

The particles evolve according to the following dynamics. Each of the particles has an independent exponential clock of rate one, and when the -clock rings the particle attempts to jump to the right by one. If at that moment then the jump is blocked. If that is not the case, we find the largest such that , and all particles in this string jump to the right by one. For any denote by the resulting measure on at time moment .

Informally speaking, the particles with smaller upper indices are heavier than those with larger upper indices, so that the heavier particles block and push the lighter ones in order for the interlacing conditions to be preserved. This anisotropy is essential, see more details in Section 1.4.

Let us illustrate the dynamics using Figure 1.2, which shows a possible configuration of particles obtained from our initial condition. If in this state of the system the -clock rings, then particle does not move, because it is blocked by particle . If it is the -clock that rings, then particle moves to the right by one unit, but to keep the interlacing property satisfied, also particles and move by one unit at the same time. This aspect of the dynamics is called “pushing”.

Observe that for , and the definition of the evolution implies that is a marginal of for any . Thus, we can think of ’s as marginals of the measure on . In other words, are measures on the space of infinite point configurations .

Before stating the main results, it is interesting to notice that the Markov chain has different interpretations. Also, some projections of the Markov chain to subsets of are still Markov chains.

1. The evolution of is the one-dimensional Poisson process of rate one.

2. The row evolves as a Markov chain on known as the Totally Asymmetric Simple Exclusion Process (TASEP), and the initial condition is commonly referred to as step initial condition. In this case, particle jumps to its right with unit rate, provided the arrival site is empty (exclusion constraint).

3. The row also evolves as a Markov chain on that is sometimes called “long range TASEP”; it was also called PushASEP in [7]. It is convenient to view as particle locations in . Then, when the -clock rings, the particle jumps to its right and pushes by one unit the (maybe empty) block of particles sitting next to it. If one disregards the particle labeling, one can think of particles as independently jumping to the next free site on their right with unit rate.

4. For our initial condition, the evolution of each row , , is also a Markov chain. It was called Charlier process in [37] because of its relation to the classical orthogonal Charlier polynomials. It can be defined as Doob -transform for independent rate one Poisson processes with the harmonic function equal to the Vandermonde determinant.

5. Infinite point configurations can be viewed as Gelfand-Tsetlin schemes. Then is the “Fourier transform” of a suitable irreducible character of the infinite-dimensional unitary group , see [13]. Interestingly enough, increasing corresponds to a deterministic flow on the space of irreducible characters of .

6. Elements of can also be viewed as lozenge tiling of a sector in the plane. To see that one surrounds each particle location by a rhombus of one type and draws edges through locations where there are no particles, see Figure 1.2. Our initial condition corresponds to a perfectly regular tiling, see Figure 1.1.

7. The random tiling defined by is the limit of the uniformly distributed lozenge tilings of hexagons with side lengths , when so that , and we observe the hexagon tiling at finite distances from the corner between sides of lengths and .

8. Finally, Figure 1.2 has a clear three-dimensional connotation. Given the random configuration at time moment , define the random height function

 h:(Z+12)×Z>0×R≥0→Z≥0, (1.2) h(x,n,t)=#{k∈{1,…,n}|xnk(t)>x}.

In terms of the tiling on Figure 1.2, the height function is defined at the vertices of rhombi, and it counts the number of particles to the right from a given vertex. (This definition differs by a simple linear function of from the standard definition of the height function for lozenge tilings, see e.g. [33, 32].) The initial condition corresponds to starting with perfectly flat facets.

Thus, our Markov chain can be viewed as a random growth model of the surface given by the height function. In terms of the step surface of Figure 1.2, the evolution consists of removing all columns of -dimensions that could be removed, independently with exponential waiting times of rate one. For example, if jumps to its right, then three consecutive cubes (associated to ) are removed. Clearly, in this dynamics the directions and do not play symmetric roles. Indeed, this model belongs to the anisotropic KPZ class of stochastic growth models, see Section 1.4.

### 1.2 Determinantal formula, limit shape and one-point fluctuations

The first result about the Markov chain that we prove is the (partial) determinantal structure of the correlation functions. Introduce the notation

 (n1,t1)≺(n2,t2)iffn1≤n2,t1≥t2, and (n1,t1)≠(n2,t2). (1.3)
###### Theorem 1.1.

For any , pick triples

 ϰj=(xj,nj,tj)∈Z×Z>0×R≥0

such that

 t1≤t2≤⋯≤tN,n1≥n2≥⋯≥nN. (1.4)

Then

 \mathbbmP{For each j=1,…,N there exists a% kj,1≤kj≤nj such that xnjkj(tj)=xj}=det[K(ϰi,ϰj)]Ni,j=1, (1.5)

where

 K(x1,n1,t1;x2,n2,t2)=−12πi∮Γ0dwwx2−x1+1e(t1−t2)/w(1−w)n2−n1\mathbbm1[(n1,t1)≺(n2,t2)]+1(2πi)2∮Γ0dw∮Γ1dzet1/wet2/z(1−w)n1(1−z)n2wx1zx2+11w−z, (1.6)

the contours , are simple positively oriented closed paths that include the poles and , respectively, and no other poles (hence, they are disjoint).

This result is proved at the end of Section 2.8. The above kernel has in fact already appeared in [7] in connection with PushASEP. The determinantal structure makes it possible to study the asymptotics. On a macroscopic scale (large time limit and hydrodynamic scaling) the model has a limit shape, which we now describe, see Figure 1.3. Since we look at heights at different times, we cannot use time as a large parameter. Instead, we introduce a large parameter and consider space and time coordinates that are comparable to . The limit shape consists of three facets interpolated by a curved piece. To describe it, consider the set

 D={(ν,η,τ)∈R3>0|(√η−√τ)2<ν<(√η+√τ)2}. (1.7)

It is exactly the set of triples for which there exists a nondegenerate triangle with side lengths . Denote by the angles of this triangle that are opposite to the corresponding sides (see Figure 3.1 too).

Our second result concerns the limit shape and the Gaussian fluctuations in the curved region, living on a scale.

###### Theorem 1.2.

For any we have the moment convergence of random variables

 Missing or unrecognized delimiter for \bigr (1.8)

with .

We also give an explicit formula for the limit shape:

 limL→∞\mathbbmEh([(ν−η)L]+12,[ηL],τL)L=:{h}(ν,η,τ)=1π(−νπη+η(π−πν)+τsinπνsinπηsinπτ). (1.9)

Theorem 1.2 describes the limit shape h of our growing surface, and the domain describes the points where this limit shape is curved. The logarithmic fluctuations is essentially a consequence of the local asymptotic behavior being governed by the discrete sine kernel (this local behavior occurs also in tiling models [31, 24, 42]). Using the connection with the Charlier ensembles, see above, the formula (1.9) for the limit shape can be read off the formulas of [5].

Using Theorem 1.1 it is not hard to verify (see Proposition 3.2 below) that near every point of the limit shape in the curved region, at any fixed time moment the random lozenge tiling approaches the unique translation invariant measure on lozenge tilings of the plane with prescribed slope (see [35, 32, 16] and references therein for discussions of these measures). The slope is exactly the slope of the tangent plane to the limit shape, given by

 ∂{h}∂ν=−πηπ,∂{h}∂η=1−πνπ. (1.10)

This implies in particular, that are the asymptotic proportions of lozenges of three different types in the neighborhood of the point of the limit shape. One also computes the growth velocity (see (1.12) for the definition of )

 ∂{h}∂τ=1πsinπνsinπηsinπτ=Im(Ω(ν,η,τ))π. (1.11)

Since the right-hand side depends only on the slope of the tangent plane, this suggest that it should be possible to extend the definition of our surface evolution to the random surfaces distributed according to measures ; these measures have to remain invariant under evolution, and the speed of the height growth should be given by the right-hand side of (1.11). This is an interesting open problem that we do not address in this paper.

### 1.3 Complex structure and multipoint fluctuations

To describe the correlations of the interface, we first need to introduce a complex structure. Set and define the map by

 |Ω(ν,η,τ)|=√η/τ,|1−Ω(ν,η,τ)|=√ν/τ. (1.12)

Observe that and . The preimage of any is a ray in that consists of triples with constant ratios . Denote this ray by . One sees that ’s are also the level sets of the slope of the tangent plane to the limit shape. Since for any , the height function grows linearly in time along each . Note also that the map satisfies

 (1−Ω)∂Ω∂ν=Ω∂Ω∂η=−∂Ω∂τ, (1.13)

and the first of these relations is the complex Burgers equation, cf. [34].

From Theorem 1.2 one might think that to get non-trivial correlations we need to consider . However, this is not true and the division by is not needed. To state the precise result, denote by

 G(z,w)=−12πln∣∣∣z−wz−¯w∣∣∣ (1.14)

the Green function of the Laplace operator on with Dirichlet boundary conditions.

###### Theorem 1.3.

For any , let be any distinct triples such that

 τ1≤τ2≤⋯≤τN,η1≥η2≥⋯≥ηN. (1.15)

Denote

 HL(ν,η,τ):=√π(h([(ν−η)L]+12,[ηL],τL)−\mathbbmEh([(ν−η)L]+12,[ηL],τL)), (1.16)

and . Then

 limL→∞\mathbbmE(HL(ϰ1)⋯HL(ϰN))=⎧⎪ ⎪⎨⎪ ⎪⎩∑σ∈FNN/2∏j=1G(Ωσ(2j−1),Ωσ(2j)),N%iseven,0,N is odd, (1.17)

where the summation is taken over all fixed point free involutions on .

The result of the theorem means that as , is a Gaussian process with covariance given by , i.e., it has correlation of the Gaussian Free Field on . We can make this statement more precise. Indeed, in addition to Theorem 1.3, a simple consequence of Theorem 1.2 gives (see Lemma 5.4),

 \mathbbmE(HL(ϰ1)⋯HL(ϰN))=O(Lϵ),L→∞, (1.18)

for any and any . This bounds the moments of for infinitesimally close points . A small extension of Theorem 1.3 together with this estimate immediately implies that on suitable surfaces in , the random function converges to the -pullback of the Gaussian free field on , see Theorem 5.6 and Theorem 5.8 in Section 5.5 for more details.

###### Conjecture 1.4.

The statement of Theorem 1.3 holds without the assumption (1.15), provided that -images of all the triples are pairwise distinct.

Theorem 1.3 and Conjecture 1.4 indicate that the fluctuations of the height function along the rays vary slower than in any other space-time direction. This statement can be rephrased more generally: the height function has smaller fluctuations along the curves where the slope of the limit shape remains constant. We have been able to find evidence for such a claim in one-dimensional random growth models as well [23].

### 1.4 Universality class

In the terminology of physics literature, see e.g. [3], our Markov chain falls into the class of local growth models with relaxation and lateral growth, described by the Kardar-Parisi-Zhang (KPZ) equation

 ∂th=Δh+Q(∂xh,∂yh)+white noise, (1.19)

where is a quadratic form. Relations (1.10) and (1.11) imply that for our growth model the determinant of the Hessian of , viewed as a function of the slope, is strictly negative, which means that the form in our case has signature . In such a situation the equation (1.19) is called anisotropic KPZ or AKPZ equation.

An example of such system is growth of vicinal surfaces, which are naturally anisotropic because the tilt direction of the surface is special. Using non-rigorous renormalization group analysis based on one-loop expansion, Wolf [52] predicted that large time fluctuations (the roughness) of the growth models described by AKPZ equation should be similar to those of linear models described by the Edwards-Wilkinson equation (heat equation with random term)

 ∂th=Δh+white noise. (1.20)

Our results can be viewed as the first rigorous analysis of a non-equilibrium growth model in the AKPZ class. (Some results, like logarithmic fluctuations, for an AKPZ model in a steady state were obtained in [44]. Some numerical numerical results are described in [28, 36, 27]). Indeed, Wolf’s prediction correctly identifies the logarithmic behavior of height fluctuations. However, it does not (at least explicitly) predict the appearance of the Gaussian free field, and in particular the complete structure (map ) of the fluctuations described in the previous section.

On the other hand, universality considerations imply that analogs of Theorems 1.2 and 1.3, as well as possibly Conjecture 1.4, should hold in any AKPZ growth model.

### 1.5 More general growth models

It turns out that the determinantal structure of the correlations functions stated in Theorem 1.1 holds for a much more general class of two-dimensional growth models. In the first part of the paper we develop an algebraic formalism needed to show that. At least three examples where this formalism applies, other than the Markov chain considered above, are worth mentioning.

1. In the Markov chain considered above one can make the particle jump rates depend on the upper index in an arbitrary way. One can also allow the particles jump both right and left, with ratio of left and right jump rates possibly changing in time [7].

2. The shuffling algorithm for domino tilings of Aztec diamonds introduced in [21] also fits into our formalism. The corresponding discrete time Markov chain is described in Section 2 below, and its equivalence to domino shuffling is established in the recent paper [39].

3. A shuffling algorithm for lozenge tilings of the hexagon (also known as boxed plane partitions) has been constructed in [12] using the formalism developed in this paper, see [12] for details.

Our original Markov chain is a suitable degeneration of each of these examples.

We expect our asymptotic methods to be applicable to many other two-dimensional growth models produced by the general formalism, and we plan to return to this discussion in a later publication.

### 1.6 Other connections

We have so far discussed the global asymptotic behavior of our growing surface, and its bulk properties (measures ), but have not discussed the edge asymptotics. As was mentioned above, rows and can be viewed as one-dimensional growth models on their own, and their asymptotic behavior was studied in [7] using essentially the same Theorem 1.1. This is exactly the edge behavior of our two-dimensional growth model.

Of course, the successive projections to and then to a fixed (large) time commute. In the first ordering, this can be seen as the large time interface associated to the TASEP. In the second ordering, it corresponds to considering a tiling problem of a large region and focusing on the border of the facet.

Interestingly enough, an analog of Theorem 1.1 remains useful for the edge computations even in the cases when the measure on the space is no longer positive (but its projections to and remain positive). These computations lead to the asymptotic results of [48, 9, 8, 7, 11, 10] for one-dimensional growth models with more general types of initial conditions.

Another natural asymptotic question that was not discussed is the limiting behavior of when but remains fixed. After proper normalization, in the limit one obtains the Markov chain investigated in [51].

Two of the four one-dimensional growth models constructed in [20] (namely, “Bernoulli with blocking” and “Bernoulli with pushing”) are projections to and of one of our two-dimensional growth models, see Section 2 below. It remains unclear however, how to interpret the other two models of [20] in a similar fashion.

Finally, let us mention that our proof of Theorem 1.1 is based on the argument of [18] and [50], the proof of Theorem 1.3 uses several ideas from [33], and the algebraic formalism for two-dimensional growth models employs a crucial idea of constructing bivariate Markov chains out of commuting univariate ones from [19].

Outline. The rest of the paper is organized as follows. It has essentially two main parts. The first part is Section 2. It contains the construction of the Markov chains, with the final result being the determinantal structure and the associated kernel (Theorem 2.25). Its continuous time analogue is Corollary 2.26, whose further specialization to particle-independent jump rate leads to Theorem 1.1. The second main part concerns the limit results for the continuous time model that we analyze. We start by collecting various geometric identities in Section 3. We also shortly discuss why our model is in the AKPZ class. In Section 4 we first give a shifted version of the kernel, whose asymptotic analysis is the content of Section 6. These results then allow us to prove Theorem 1.2 in Section 4 and Theorem 1.3 in Section 5.

Acknowledgments. The authors are very grateful to P. Diaconis, E. Rains, and H. Spohn for numerous illuminating discussions. The first named author (A. B.) was partially supported by the NSF grant DMS-0707163.

## 2 Two dimensional dynamics

All the constructions below are based on the following basic idea. Consider two Markov operators and on state spaces and , and a Markov link that intertwines and , that is . Then one can construct Markov chains on (subsets of) that in some sense has both and as their projections. There is more than one way to realize this idea, and in this paper we discuss two variants.

In one of them the image of under the Markov operator is determined by sequential update: One first chooses according to , and then one chooses so that the needed projection properties are satisfied. A characteristic feature of the construction is that and are independent, given and . This bivariate Markov chain is denoted ; its construction is borrowed from [19].

In the second variant, the images and are independent, given , and we say that they are obtained by parallel update. The distribution of is still , independently of what is. This Markov chain is denoted for the operator that plays an important role.

By induction, one constructs multivariate Markov chains out of finitely many univariate ones and links that intertwine them. Again, we use two variants of the construction — with sequential and parallel updates.

The key property that makes these constructions useful is the following: If the chains , , and , are -Doob transforms of some (simpler) Markov chains, and the harmonic functions used are consistent, then the transition probabilities of the multivariate Markov chains do not depend on . Thus, participating multivariate Markov chains may be fairly complex, while the transition probabilities of the univariate Markov chains remain simple.

Below we first explain the abstract construction of , , and their multivariate extensions. Then we exhibit a class of examples that are of interest to us. Finally, we show how the knowledge of certain averages (correlation functions) for the univariate Markov chains allows one to compute similar averages for the multivariate chains.

### 2.1 Bivariate Markov chains

Let and be discrete sets, and let and be stochastic matrices on these sets:

 ∑y∈SP(x,y)=1,x∈S;∑y∗∈S∗P∗(x∗,y∗)=1,x∗∈S∗. (2.1)

Assume that there exists a third stochastic matrix such that for any and

 ∑x∈SΛ(x∗,x)P(x,y)=∑y∗∈S∗P∗(x∗,y∗)Λ(y∗,y). (2.2)

Let us denote the above quantity by . In matrix notation

 Δ=ΛP=P∗Λ. (2.3)

Set

 SΛ={(x∗,x)∈S∗×S|Λ(x∗,x)>0}, SΔ={(x∗,x)∈S∗×S|Δ(x∗,x)>0}.

Define bivariate Markov chains on and by their corresponding transition probabilities

 PΛ((x∗,x),(y∗,y))=⎧⎪⎨⎪⎩P(x,y)P∗(x∗,y∗)Λ(y∗,y)Δ(x∗,y),Δ(x∗,y)>0,0,otherwise, (2.4)
 PΔ((x∗,x),(y∗,y))=P(x,y)P∗(x∗,y∗)Λ(y∗,x)Δ(x∗,x). (2.5)

It is immediately verified that both matrices and are stochastic.

The chain was introduced by Diaconis-Fill in [19], and we are using the notation of that paper.

One could think of and as follows.

For , starting from we first choose according to the transition matrix , and then choose using , which is the conditional distribution of the middle point in the successive application of and provided that we start at and finish at .

For , starting from we independently choose according to and according to , which is the conditional distribution of the middle point in the successive application of and provided that we start at and finish at .

###### Lemma 2.1.

For any , we have

 ∑y∗∈S∗:(y∗,y)∈SΛPΛ((x∗,x),(y∗,y)) =P(x,y), (2.6) ∑y∗∈S∗:(y∗,y)∈SΔPΔ((x∗,x),(y∗,y)) =P(x,y),

and for any ,

 ∑x∈S:(x∗,x)∈SΛΛ(x∗,x)PΛ((x∗,x),(y∗,y)) =P∗(x∗,y∗)Λ(y∗,y), (2.7) ∑x∈S:(x∗,x)∈SΔΔ(x∗,x)PΔ((x∗,x),(y∗,y)) =P∗(x∗,y∗)Δ(y∗,y).

Proof of Lemma 2.1. Straightforward computation using the relation . ∎

###### Proposition 2.2.

Let be a probability measure on . Consider the evolution of the measure on under the Markov chain and denote by the result after steps. Then for any the joint distribution of

 (x∗(0),x∗(1),…,x∗(k),x(k),x(k+1),…,x(k+l)) (2.8)

coincides with the stochastic evolution of under transition matrices

 (P∗,…,P∗k,Λ,P,…,Pl). (2.9)

Exactly the same statement holds for the Markov chain and the initial condition with replaced by in the above sequence of matrices.

Proof of Proposition 2.2. Successive application of the first relations of Lemma 2.1 to evaluate the sums over , and of the second relations to evaluate the sums over . ∎

Note that Proposition 2.2 also implies that the joint distribution of and has the form , where is the result of -fold application of to .

The above constructions can be generalized to the nonautonomous situation.

Assume that we have a time variable , and our state spaces as well as transition matrices depend on , which we will indicate as follows:

 S(t),S∗(t),P(x,y|t),P∗(x∗,y∗|t),Λ(x∗,x|t),P(t),P∗(t),Λ(t). (2.10)

The commutation relation 1.3 is replaced by or

 Δ(x∗,y|t):=∑x∈S(t)Λ(x∗,x|t)P(x,y|t)=∑y∗∈S∗(t+1)P∗(x∗,y∗|t)Λ(y∗,y|t+1). (2.11)

Further, we set

 SΛ(t) ={(x∗,x)∈S∗(t)×S(t)|Λ(x∗,x|t)>0}, (2.12) SΔ(t) ={(x∗,x)∈S∗(t)×S(t+1)|Δ(x∗,x|t)>0},

and

 PΛ((x∗,x),(y∗,y)|t)=⎧⎪⎨⎪⎩P(x,y|t)P∗(x∗,y∗|t)Λ(y∗,y|t+1)Δ(x∗,y|t),Δ(x∗,y|t)>0,0,otherwise, (2.13)
 PΔ((x∗,x),(y∗,y)|t)=P(x,y|t+1)P∗(x∗,y∗|t)Λ(y∗,x|t+1)Δ(x∗,x|t). (2.14)

The nonautonomous generalization of Proposition 2.2 is proved in exactly the same way as Proposition 2.2. Let us state it.

###### Proposition 2.3.

Fix , and let be a probability measure on . Consider the evolution of the measure on under the Markov chain , and denote by the result after steps. Then for any the joint distribution of

 (x∗(t0),x∗(t0+1),…,x∗(t0+k),x(t0+k),x(t0+k+1),…,x(t0+k+l)) (2.15)

coincides with the stochastic evolution of under transition matrices

 P∗(t0),…,P∗(t0+k−1),Λ(t0+k),P(t0+k),…,P(t0+k+l−1) (2.16)

(for only remains in this string).

A similar statement holds for the Markov chain and the initial condition : For any the joint distribution of

 (x∗(t0),x∗(t0+1),…,x∗(t0+k),x(t0+k+1),x(t0+k+2),…,x(t0+k+l+1)) (2.17)

coincides with the stochastic evolution of under transition matrices

 P∗(t0),…,P∗(t0+k−1),Δ(t0+k),P(t0+k+1),…,P(t0+k+l). (2.18)
###### Remark 2.4.

Observe that there is a difference in the sequences of times used in (2.8) and (2.17). The reason is that for nonautonomous , the state space at time is a subset of , and we denote its elements as . In the autonomous case, an element of the state space at time was denoted as .

### 2.2 Multivariate Markov chains

We now aim at generalizing the constructions of Section 2.1 to more than two state spaces.

Let be discrete sets, be stochastic matrices defining Markov chains on them, and let be stochastic links between these sets:

 Pk:Sk×Sk→[0,1], ∑y∈SkPk(x,y)=1,x∈Sk,k=1,…,n; (2.19) Λkk−1:Sk×Sk−1→[0,1], ∑y∈Sk−1Λkk−1(x,y)=1,x∈Sk,k=2,…,n.

Assume that these matrices satisfy the commutation relations

 Δkk−1:=Λkk−1Pk−1=PkΛkk−1,k=2,…,n. (2.20)

The state spaces for our multivariate Markov chains are defined as follows

 S(n)Λ ={(x1,…,xn)∈S1×⋯×Sn|n∏k=2Λkk−1(xk,xk−1)≠0}, (2.21) S(n)Δ ={(x1,…,xn)∈S1×⋯×Sn|n∏k=2Δkk−1(xk,xk−1)≠0}.

The transition probabilities for the Markov chains and are defined as (we use the notation , )

 P(n)Λ(Xn,Yn)=⎧⎪ ⎪⎨⎪ ⎪⎩P1(x1,y1)n∏k=2Pk(xk,yk)Λkk−1(yk,yk−1)Δkk−1(xk,yk−1),n∏k=2Δkk−1(xk,yk−1)>0,0,otherwise, (2.22)
 P(n)Δ(Xn,Yn)=P(x1,y1)n∏k=2Pk(xk,yk)Λkk−1(yk,xk−1)Δkk−1(xk,xk−1). (2.23)

One way to think of and is as follows. For , starting from , we first choose according to the transition matrix , then choose using