The unperturbed GUE random matrix, corresponding to a i ≡ 0, is unitary invariant in the sense that there is equality in distribution G d = U GU * for any fixed N × N unitary matrix U . This implies that the distribution of the eigenvalues of H is symmetric in the perturbation parameters a 1 , . . . , a N . The overall goal of the paper is to explore probabilistic consequences of this symmetry property.

Together with the eigenvalues λ N = (λ N N ≤ . . . ≤ λ N 1 ), λ N i ∈ R, of the full matrix H = [h i j ] N i, j=1 , one can also consider its corners process, 1 that is, the interlacing collection of eigenvalues of the principal corners [h i j ] k i, j=1 of H for all k = 1, 2, . . . , N . (See Fig. 1 for an illustration.) The distribution of the corners process of H is not symmetric in the parameters a i . Moreover, assuming that the a i 's are all distinct, there are N ! different probability distributions on interlacing collections of eigenvalues at N levels.

In this paper we present explicit couplings between these N ! distributions, by showing that each nearest neighbour transposition a k ↔ a k+1 , k = 1, . . . , N -1, of the parameters is equivalent in distribution to a rather simple Markov swap operator S a k -a k+1 k . This swap operator randomly changes the entries λ k i on the k-th level of the array given the two adjacent levels λ k-1 , λ k+1 , while leaving all other entries intact. If a k > a k+1 , S a k -a k+1 k is realized as an independent collection of instantaneous exponential type jumps of each λ k i to the left:foot_1

where E i a k -a k+1 's are independent exponential random variables with parameter a ka k+1 (and mean 1/(a ka k+1 )). Here by agreement, λ k-1 k = -∞. In particular, these left jumps are constrained by the interlacing. For a k < a k+1 , the same jumps are performed to the right in a symmetric way:

where, by agreement, λ k-1 0 = +∞. Finally, if a k = a k+1 , then S a k -a k+1 k is the identity operation.

Theorem 1.1 (Follows from Theorem 4.4 below) Assume that a k = a k+1 . Then the action of the Markov operator S a k -a k+1 k (with left jumps for a k > a k+1 , and right jumps otherwise) turns the corners distribution of G + diag(a 1 , . . . , a k , a k+1 , . . . , a N ) into the one of G + diag(a 1 , . . . , a k+1 , a k , . . . , a N ), where G is the N × N GUE random matrix.

We establish this theorem by relying on a perturbed Gibbs structure of the corners distribution of the matrix H . Namely, it is well-known that in the unperturbed case, the conditional distribution of the eigenvalues λ k i , 1 ≤ i ≤ k ≤ N -1, given λ N , is uniform on the polytope defined by all the interlacing inequalities (known as the Gelfand-Tsetlin polytope). In the perturbed case, the Gibbs structure should be deformed in a certain way by means of the parameters a i (see Sect. 3.1). The coupling follows by considering the conditional distribution of λ k given two adjacent levels λ k±1 , which reduces to a collection of independent exponential random variables confined to the corresponding intervals. Producing a suitable Markov swap operator for a single such variable (see Proposition 4.1 below), we arrive at the result of Theorem 1.1. does not change the distribution of G + diag(a 1 , . . . , a N ). However, this composition is not an identity transformation: two random jumps return a particle to the original location with probability 0.

The perturbed GUE corners distributions are compatible for various N , and so one can define the corresponding perturbed GUE corners distribution on infinite interlacing arrays. It depends on an infinite parameter sequence a = {a i } i∈Z ≥1 . One particular interesting case is when the perturbation parameters form an arithmetic progression a i = -(i -1)α, where α > 0. Swapping a 1 with a 2 , then a 1 with a 3 , and so on all the way to infinity leads to an additive shift in the perturbation parameters, which is equivalent in distribution to a global shift: Theorem 1.3 (Theorem 5.2 below) The action of a sequence of left exponential jumps (where the parameter at level k is taken to be kα), from level 1 up to infinity, is equivalent in distribution 3 to shifting all the elements of the interlacing array by α to the left.

Let α > 0, and let let X k (t), k = 1, 2, . . ., be reflected Brownian motions constructed as follows. First, X 1 (t) is the standard driftless Brownian motion started from 0. Inductively, let X k (t), k = 2, 3, . . ., be a new independent Brownian motion with drift -(k -1)α, and reflected down off of X k-1 (t) by means of subtracting local time when

is the standard Brownian motion, and

s) is the continuous non-decreasing process which increases only at times when X 1 (s) = X 2 (s) (in other words, it is twice the semimartingale local time of X 1 -X 2 at zero. We refer to [5,19] for further details on the reflection mechanism, and for an explanation on how to start all these reflected processes from zero (which formally results in infinitely many collisions in finite time). Almost surely we have X 1 (t) ≥ X 2 (t) ≥ X 3 (t) . . . for all t.

Fix t and define

where E kα , k = 1, 2, . . ., are independent exponential random variables with parameters kα (and mean 1/(kα)).

Theorem 1. 4 For each fixed t, we have equality of joint distributions

In particular,

) is a normal random variable with mean (-αt) and variance t. To the best of our knowledge, even this result for two processes (one a usual Brownian motion, and one reflected off it) is new.

Theorem 1.4 follows from Theorem 1.3 combined with the connection between the reflected drifted Brownian motions and the perturbed GUE corners process due to [5]. We recall this connection in detail in Sect. 2.3 below, and prove Theorem 1.4 in the end of Sect. 5.

As stated, Theorem 1.4 assumes that the time t is fixed. Indeed, naively taking independent exponential shifts at different times t would make the functions t → X k (t) discontinuous. It is interesting to see whether a stochastic process version of Theorem 1.4 holds:

Presumably, if such a Markov operator on processes exists, then its construction could be accomplished using the sticky Brownian motion, 4 as exponential random variables arise in the study of this process, e.g., see Theorem 1 in [18]. It seems plausible that the difference process t → X 1 (t) -X 1 (t) ≥ 0 itself could be distributed as the sticky Brownian motion, as the single-time distributions coincide thanks to the results of [18] and [10,Proposition 14]. However, it is less clear how to extend this idea to all differences t → X k (t) -X k (t) ≥ 0.

The results of this paper might be viewed as a random matrix limit of the ones from the recent work [17]. There, similar Markov swap operators were considered on discrete interlacing arrays as in Fig. 1. A combination of these swap operators together with a certain Poissontype limit (cf. Sect. 1.5 below) has lead to a Markov chain on distributions of TASEP (totally asymmetric simple exclusion process) which decreases the time parameter. The shifting result for reflected drifted Brownian motions (Theorem 1.4) may be viewed as a certain analogue of the TASEP reversal property. In the Brownian case, instead of decreasing the time, the exponential jumps lead to a deterministic shift.

It should be pointed out that even though the discrete stochastic systems in [17] converge to the reflected Brownian motions [8] (and [5] in the drifted case), here we do not rely on this convergence or the results of [17]. Instead we obtain the results independently using basic mechanisms related to the (perturbed) Gibbs property.

In the arithmetic progression setting a i = -(i -1)α with α > 0, when α 0, the perturbed GUE corners process of H = G + diag(0, -α, -2α, . . .) becomes the usual GUE corners process, and the system of reflected Brownian motions {X k (t)} k∈Z ≥1 becomes driftless. It would be very interesting to see whether the Markov operators considered in the present paper have meaningful limits as α 0. However, this limit presents certain immediate issues which we discuss now.

For simplicity, consider the Brownian motion setup. Fix t > 0 and suppress this parameter in the notation. As α → 0, the Markov operator X k → X k (1.1) turns into the (deterministic) identity operator X k → X k . Indeed, this is because Prob(E kα > x) = e -kαx ∼ 1 -αkx for all k and x, and so the minimum in (1.1) is equal to X k -X k+1 with probability of order 1 -O(α). Arguing similarly to the discrete case considered in [17,Sect. 6], one can apply the map (1.1) a large number τ/α of times, where τ ∈ R >0 is the scaled time.

Taking a Poisson-type limit should lead to a continuous time Markov process (with τ as the new time parameter) under which X k has an exponential clock of rate k(X k -X k+1 ), and when the clock rings, X k instantaneously jumps into

This jumping mechanism is known as the Hammersley process [1,9]. However, applying this continuous time jumping process to the whole system {X k } k∈Z ≥1 is problematic, as it leads to infinitely many jumps in finite time due to the growing jump rates k(X k -X k+1 ) as k → ∞. Moreover, under this hypothetical process X k would depend on all X j for j > k, and so one cannot simply restrict the dynamics to finitely many particles where it would make sense.

On the other hand, by Theorem 1.4, the hypothetical continuous time dynamics should be equivalent in distribution to a deterministic shift of the (driftless) reflected Brownian motions by -αt τ/α ∼ -tτ . It is reasonable to expect that such a deterministic shift of infinitely many X k 's cannot be achieved only by finitely many jumps in finite time. To summarize, Open problem 1.6 Do there exist well-defined α 0 limits of the Markov operators acting on the GUE corners process perturbed by an arithmetic progression a i = -(i -1)α or on the reflected drifted Brownian motions? These hypothetical limits should act on (much more studied) unperturbed GUE corners process and driftless reflected Brownian motions.

This section is preliminary. We recall the perturbed GUE corners process [5] (also called the GUE corners process with external source [3]), and its connection to reflected Brownian motions with drifts. The original, unperturbed GUE corners process is due to [11,12], and it was linked to driftless reflected Brownian motions in [19].

Take a time parameter t > 0 and an infinite sequence of parameters a = (a 1 , a 2 , . . .), a i ∈ R.

Unless otherwise indicated, we assume that the parameters a i are pairwise distinct. Consider a random matrix H = t 1/2 • G + t • diag (a) of infinite size with entries:

Here g kk are independent real standard normal random variables, and g kl are independent complex standard normal random variables (that is, their real and imaginary parts are independent real normal random variables each with mean 0 and variance 1 2 ). The matrix H is Hermitian.

For each m ∈ Z ≥1 , take the m × m principal corner

, be its eigenvalues. At adjacent levels, the eigenvalues interlace (notation λ m ≺ λ m+1 ):

We call the joint distribution of all {λ k j } 1≤ j≤k<∞ the perturbed GUE corners process.

A standard application of the Harish-Chandra-Itsykson-Zuber integral shows that the joint eigenvalue density of {λ N i } N i=1 at a fixed level N is given by

where the normalizing constant does not depend on a 1 , . . . , a N . Here and throughout the paper we use the notation

for the Vandermonde determinant. Observe from (2.2) that the distribution of {λ N j } N i=1 depends on the parameters a i in a symmetric way. This should indeed be the case, since the distribution of the eigenvalues of the N × N matrix t 1/2 G N ×N + t diag(a 1 , . . . , a N ) does not depend on the order of the a i 's due to the unitary invariance of G N ×N . The main goal of this paper is to explore this distributional symmetry from a Markov operator point of view. For this, we will need the joint distribution of eigenvalues of all corners: Proposition 2.1 ([5, Proposition 2.3]) The joint density of the eigenvalues {λ k j } 1≤ j≤k≤N at the first N levels, where N ∈ Z ≥1 is arbitrary, has the following form:

3) where we use the notation

and the normalizing constant does not depend on a 1 , . . . , a N .

Fix a perturbation sequence a = {a i } i∈Z ≥1 . Consider a family of interacting Brownian motions B k j , 1 ≤ j ≤ k < ∞, such that: • All processes start from zero.

• The processes B k j , j = 1, . . . , k have the same drift a k . • The evolution of each B k j does not depend on any of the B l i 's with l > i. • The processes interact only through their local times. That is, when the processes are sufficiently far apart, each B k j behaves as an independent Brownian motion with drift a k .

foot_4 and reflects off both B k-1 j and B k-1 j-1 . Therefore, at each time t, the random variables {B k j (t)} 1≤ j≤k<∞ almost surely form an interlacing array as in Fig. 1.

We refer to [5,Sect. 4] (and [19] in the driftless case) for details on the reflection mechanism.

At each fixed time moment t ∈ R ≥0 , we have equality of joint distributions of two infinite interlacing arrays:

, where the right-hand side is the perturbed GUE corners process with the same time parameter t and perturbation sequence a.

In this section we place the perturbed GUE corners process into a wider family of Gibbs measures on interlacing arrays.

A measure on infinite interlacing arrays {λ k j } 1≤ j≤k<∞ (satisfying inequalities (2.1) between any two consecutive levels) is called a-Gibbs if for each N and any fixed configuration λ N at level N , the density of the conditional distribution of all the lower entries of the array has the form

(if some of the λ N i 's are equal, the density would have delta components and formula (3.1) should be understood in a limiting sense). Here and below 1 B stands for the indicator of an event B. Proposition 2.1 implies that the perturbed GUE corners process is an example of an a-Gibbs measure.

The fact that the density (3.1) integrates to 1 in λ 1 , . . . , λ N -1 can be checked by induction on N . Remark 3.2 When a i ≡ a are all equal to each other, the a-Gibbs property becomes the usual Gibbs property, with (3.1) replaced by the uniform conditioning provided that the configurations λ 1 , . . . , λ N -1 , λ N interlace. A classification of uniform Gibbs measures on interlacing arrays is due to [16]. In fact, performing a suitable exponential change of variables, one can see that when a is an arithmetic progression, the space of a-Gibbs measures is essentially the same as in the uniform case. This is somewhat parallel to how the two-sided q-Gelfand-Tsetlin graph degenerates to the "graph of spectra" [7,15].

To each a-Gibbs measure we can associate a family of a-harmonic functions as follows:

Density(λ N ),

where Density(λ N ) is the marginal density of λ N . The term "harmonic function" comes from the Vershik-Kerov theory of boundaries of branching graphs, cf. [13]. Harmonicity means that the functions satisfy a version of a mean value theorem associated to a directed graph Laplacian. In the context of random matrices the discrete graph is replaced by a suitable continuous analogue, and the graph Laplacian becomes an integral operator. In other words, since the a-harmonic functions ϕ N for different N 's come from the same Gibbs measure, they must be consistent in the following sense:

3) should be viewed as a version of the mean value theorem, as discussed before Lemma 3.3.

The claim follows by writing down the joint distribution of λ 1 , . . . , λ N through ϕ N and the conditional distribution (3.1), and then integrating out λ 1 , . . . , λ N -2 (this produces the factor V(a 1 , . . . ,

) and λ N to get the marginal density of λ N -1 . The resulting marginal density is expressed through ϕ N -1 via (3.2), which yields the result.

For an a-Gibbs measure, let each ϕ k depend on a 1 , . . . , a k in a symmetric way. Then the distribution of λ k , where k ∈ Z ≥1 is fixed, depends on the parameters a 1 , . . . , a k in a symmetric way, too.

Proof An immediate consequence of (3.2).

Proof Follows from the Kolmogorov extension theorem.

Let us emphasize that the results of this subsection (Lemmas 3.3 and 3.4 and Proposition 3.5) are valid not only for the perturbed GUE corners process (which, as we see next, is an example of an a-Gibbs measure), but hold in the full generality for any a-Gibbs measure.

One readily sees that for the perturbed GUE corners process we have the following harmonic functions:

(3.4) where the constant is the same as in (2.3) and does not depend on the a j 's. One readily checks that the a-Gibbs property (Lemma 3.3) for the perturbed GUE corners process is equivalent to the well-known integral identity for the Vandermonde determinants:

where the constant does not depend on the a j 's. The shift by a N t in the exponents in both sides by changing the variables in the integral and renaming the λ N -1 i 's, can also be removed (or replaced with any other shift bt) since the Vandermonde is translation invariant.

In particular, (2.3) together with Lemma 3.4 implies the symmetry (as in this lemma) of the perturbed GUE corners distribution with respect to the parameters a i .

In this section we explore the Gibbs property and prove Theorem 1.1 on Markov swap operators.

The next observation plays a key role: Proposition 4.1 Take real numbers c < d and α > 0. Let X be distributed as E α (c, d), and E α ∈ (0, +∞) be an independent usual exponential random variable with parameter α (i.e., with density αe -α y , y > 0). Then the random variable

is distributed as E -α (c, d).

Proof We have for the conditional distribution of Y given X = x:

The distribution of Y has an atom at y = x (coming from the event E α > Xc in (4.1)) and an absolutely continuous part on (0, x). The overall density of Y in the variable y is obtained from the following integral:

-α e -dαe -cα e -α y , which completes the proof.

We will view the operation of passing from X to Y as in (4.1) as a one-step Markov transition operator. One can think that the "particle" X ∈ (c, d) jumps left into the new location Y , by means of the new exponential random variable E α . Note that the new location Y depends only on X and not on the right endpoint d of the interval. We call this Markov transition operator the elementary swap operator and denote it by S α . This operator acts on distributions (in our case, densities) as Density Y = Density X S α .

The swap operator S α is analogous to the jump operator L α in the discrete situation considered in [17,Sect. 4]. Let us make a number of remarks. Remark 4.2 (1) When α = 0, the swap operator S α should be understood as the identity map, which is evident from (4.2). ( 2) For α = -β < 0, algebraic manipulations in the proof of Proposition 4.1 make sense, but the new random variable Y obtained by applying S -β to X ∼ E -β (c, d) does not admit a probabilistic interpretation as in (4.1). (3) Instead of applying S -β to E -β (c, d), let us consider the operator which moves X to the right symmetrically to how S α moves the "particle" X to the left. That is, this new operator acts as

, where E β is an independent exponential random variable. One can show similarly to Proposition 4.

All our results for Markov operators built from the left jumps S α have straightforward analogues for these right jumping operators, and so we will only focus on the left jumps in the paper.

Let us fix a perturbation sequence a, and let {λ m j } 1≤ j≤m<∞ be a random interlacing array distributed according to some a-Gibbs measure (for example, the perturbed GUE corners process with an arbitrary time parameter t ≥ 0).

Next, fix a level k ∈ Z ≥1 , and consider the conditional distribution of λ k given the two adjacent levels λ k-1 , λ k+1 (if k = 1, the conditioning is only on λ 2 ). From (3.1) one readily sees that this conditional distribution takes the form

where we have denoted α := a ka k+1 . Equivalently, we can describe distribution (4.3) as follows.

Proposition 4. 3 The conditional distribution of λ k given λ k-1 and λ k+1 is such that each λ k i , i = 1, . . . , k, is an independent random variable distributed as

where α = a ka k+1 . (For i = k we set λ k-1 k = -∞, and for i = 1 we set λ k-1 0 = +∞, but both ends of the interval in (4.4) are always finite.) Proof Readily follows from (4.3).

Assume that α = a ka k+1 > 0, and take an array {λ m j } 1≤ j≤m<∞ as above. Let us define a new random interlacing array {ν m j } 1≤ j≤m<+∞ for which ν m j = λ m j for all m = k, j = 1, . . . , m, and such that

where E 1 α , . . . , E k α are independent usual exponential random variables with parameter α. Note that almost surely we have ν k i ≤ λ k i , i = 1, . . . , k. In other words, in (4.5) we independently apply the elementary swap operator S α to each λ k i which is confined to the corresponding interval as in Proposition 4.3. Denote this combination of the swap operators applied at level k by S α k . As in Remark 4.2, the Markov operator S α k makes sense only for α > 0. Let τ k denote the elementary transposition (k, k + 1). For a perturbation sequence a, let τ k a = (a 1 , . . . , a k-1 , a k+1 , a k , . . .) be the permuted sequence. Theorem 4.4 (Theorem 1.1 in Introduction) Take an a-Gibbs measure for which each harmonic function ϕ N depends on the parameters a 1 , . . . , a N in a symmetric way. If a k > a k+1 , then the action of the Markov operator S α k (with α = a ka k+1 ) on this a-Gibbs measure results in a τ k a-Gibbs measure which corresponds to harmonic functions modified as follows:

Proof Since the action of S α k does not change levels j = k (and hence distributions of these levels), we clearly have ϕ j = ϕ j for j = k.

Thus, it remains to show that under S α k the a-Gibbs property becomes τ k a-Gibbs. This can be seen by representing the conditional distributions as

The left-hand side depends on a 1 , . . . , a k+1 in a symmetric way. One can readily check that Prob(λ 1 , . . . , λ k-1 | λ k+1 ) depends on the parameters a k , a k+1 in a symmetric way, too. Indeed, this conditional distribution corresponds to integrating (3.1) (with N = k + 1) over λ k . The non-exponential prefactor in (3.1) is already symmetric, and for the exponential part we have

where we used the normalizing constant for the confined exponential distribution, and α = a ka k+1 . Swapping the parameters as a k ↔ a k+1 brings exp -α(|λ k+1 | + |λ k-1 |) from the exponential factor in front of the product in (4.8). This factor compensates the product of the expressions exp

over i = 1, . . . , k, coming out of the product in (4.8) after the same swap. Thus, (4.8) is symmetric under a k ↔ a k+1 .

The action of S α k affects only the part Prob(λ k | λ k-1 , λ k-1 ) in the right-hand side of (4.7). Before the action of S α k , each λ k i was distributed as E α on the corresponding interval (see Proposition 4.3). By Proposition 4.1, after the action of S α k , these random variables turn into the E -α 's, which corresponds to a τ k a-Gibbs structure. Combining this with the symmetries in (4.7) described above and using Lemma 3.4 and Proposition 3.5, we arrive at the claim.

In particular, for a k > a k+1 , the perturbed GUE corners process coming from the random matrix

, turns into the corners process for the random matrix

where T k is the permutation matrix of τ k = (k, k+1). In other words, applying the exponential jumps S a k -a k+1 k on the level of eigenvalues is equivalent in distribution to the change of basis e k ↔ e k+1 in the space corresponding to the random matrix.

Acting on a-Gibbs measures with a = (0, -α, -2α, . . .) (5.1), S α 1 interchanges 0 with -α, then S 2α 2 interchanges 0 (which is now the new a 2 ) with -2α, and so on. After infinitely many swaps, the parameter 0 disappears, and one expects that the resulting distribution would be a-Gibbs with a = (-α, -2α, -3α, . . .). For the special choice of the perturbed GUE corners process (5.2), the action of S α is, moreover, equivalent in distribution to a global shift. We establish the following result: where I is the infinite identity matrix.

Proof Informally, one can think that (5.3) follows by a sequential application of the change of basis (4.9) under a single-level action S kα k . The shift by αt is precisely the difference between t • diag(a) before and after the modification of a. We will now prove this claim more formally, using Theorem 4.4 on how Gibbs measures change under swap operators.

Take the harmonic functions ϕ N = ϕ pertGUE(a;t) N as in (3.4) with a = (0, -α, -2α, . . .). The action of each S kα k changes only the k-th function ϕ k as in (4.6) and leaves all other functions intact. Therefore, the action of the whole S α replaces {ϕ k } by the family

Here we took a k = 0 because this is precisely the perturbation parameter that is being swapped with a k+1 = -kα under the action of S kα k . The integral in the right-hand side of (5.4) can be computed using (3.5) (with a N = 0 in that formula), and we obtain

k+1 j=1 e -t(( j-1)α) 2 /2 = C 0 ϕ k (λ k ) e -tk 2 α 2 /2 .

Here both const and C 0 are some constants which are independent of α. Sequentially applying Theorem 4.4, we see that the new the harmonic functions {ϕ k } satisfy Gibbs property with the sequence a = (-α, -2α, -3α, . . .), and hence (by Proposition 3.5) correspond to a Gibbs measure with shifted parameters. Let us now identify this particular Gibbs measure. The modified density of λ k after the application of S α reads, by (3.2),

V(-α, -2α, . . . , -kα)

V(-α, -2α, . . . , -kα)

V(-α, -2α, . . . , -kα) where Density(•) is the original density before applying S α . In the last step, the two Vandermondes are equal by their translation invariance, and the ratio of the determinants is e -α|λ k | (indeed, factor out e -λ k j from each j-th column of the determinant in the numerator). Now, using (2.2) we have C 0 e -α|λ k |-tk 2 α 2 /2 Density(λ k ) = C 0 const × e -α|λ k |-tk 2 α 2 /2 det exp -(λ k i + t( j -1)α) 2 2t k i, j=1

V(λ k 1 , . . . , λ k k ) V(0, -α, . . . , -(k -1)α) .

Here const is the normalizing constant in (2.2) which is independent of α. Observe that in the exponents inside the determinant we have

(2 j -1).

Factoring out the last two terms from each j-th column, we get a factor which precisely cancels with e -α|λ k |-tk 2 α 2 /2 . Therefore, we see that

Normalizing, this implies that C 0 = 1. Thus, we see that applying S α is indeed equivalent to the global shift by αt to the left, as desired.

We can now establish the shifting property for the reflected Brownian motions: