Tristan Gautié

Jean-Philippe Bouchaud

Pierre Le Doussal

Matrix Kesten Recursion, Inverse-Wishart Ensemble and Fermions in a Morse Potential

Keywords: numbers: 05.30.Fk, 02.10.Yn, 02.50.-r, 05.40.-a

The random variable 1 + z1 + z1z2 + . . . appears in many contexts and was shown by Kesten to exhibit a heavy tail distribution. We consider natural extensions of this variable and its associated recursion to N × N matrices either real symmetric β = 1 or complex Hermitian β = 2. In the continuum limit of this recursion, we show that the matrix distribution converges to the inverse-Wishart ensemble of random matrices. The full dynamics is solved using a mapping to N fermions in a Morse potential, which are non-interacting for β = 2. At finite N the distribution of eigenvalues exhibits heavy tails, generalizing Kesten's results in the scalar case. The density of fermions in this potential is studied for large N , and the power-law tail of the eigenvalue distribution is related to the properties of the so-called determinantal Bessel process which describes the hard edge universality of random matrices. For the discrete matrix recursion, using free probability in the large N limit, we obtain a self-consistent equation for the stationary distribution. The relation of our results to recent works of Rider and Valkó, Grabsch and Texier, as well as Ossipov, is discussed.

I. INTRODUCTION

The dynamical generation of random variables with heavy tails is a fascinating and ubiquitous phenomenon. In physics it occurs, among many others, in problems such as diffusion in random media [START_REF] Bouchaud | Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications[END_REF], directed polymers in random media [START_REF] Derrida | Polymers on disordered trees, spin glasses, and traveling waves[END_REF], Anderson localization in an external field [START_REF] Delyon | From Power-Localized to Extended States in a Class of One-Dimensional Disordered Systems[END_REF][START_REF] Prigodin | One-dimensional disordered system in an electric field[END_REF][START_REF] Crosnier De Bellaistre | Expansion of a quantum wave packet in a one-dimensional disordered potential in the presence of a uniform bias force[END_REF], growing networks [START_REF] Dorogovtsev | Evolution of networks: From biological nets to the Internet and WWW[END_REF], avalanches in driven elastic systems [START_REF] Alessandro | Domain-wall dynamics and Barkhausen effect in metallic ferromagnetic materials I[END_REF] or in self-organized critical systems such as sandpiles [START_REF] Bak | Self-organized criticality: an explanation of 1/f noise[END_REF]. Power-law tails are also ubiquitous in economics and finance [START_REF] Bouchaud | Power laws in economics and finance: some ideas from physics[END_REF][START_REF] Gabaix | Power laws in economics and finance[END_REF]: wealth distribution in models of economic agents [START_REF] Bouchaud | Wealth condensation in a simple model of economy[END_REF], returns of financial markets [START_REF] Gabaix | Power laws in economics and finance[END_REF], city and firm sizes [START_REF] Axtell | Zipf Distribution of U.S. Firm Sizes[END_REF][START_REF] Gabaix | Zipf 's law for cities: an explanation[END_REF][START_REF] Saichev | Theory of Zipf 's law and beyond[END_REF], etc. One of the simplest example giving rise to power-law tails is a stochastic recursion equation studied in a seminal paper by Kesten [START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF], which we will call here the Kesten recursion. It was studied in a scalar and a vector form. It occurs naturally in its scalar form in models of random walks in random environments [START_REF] Solomon | Random walks in a random environment[END_REF][START_REF] Kesten | Spitzer A limit law for random walk in a random environment[END_REF] and in its vector form to describe the mean population in multi-specie branching processes in random environments. Its continuum limit appears in the problem of Arrhenius diffusion of a particle in a one-dimensional random force field, where it predicts the generation of broad distributions of waiting times, which can be seen as effective traps for the particle [START_REF] Bouchaud | Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications[END_REF][START_REF] Bouchaud | Classical diffusion of a particle in a one-dimensional random force field[END_REF].

In this paper our main aim is to introduce and study a matrix realization of the Kesten recursion, which we call the matrix Kesten recursion. There are several ways to generalize the scalar recursion (N = 1) to the case of a N × N matrix. Here we will study two such generalizations, one with real symmetric matrices (denoted β = 1) and the other with Hermitian matrices (β = 2). We will obtain two sets of exact results. First we will provide a solution for the continuum limit of this recursion for arbitrary size N of the matrix. It shows that the marginal distribution of the few largest eigenvalues of this matrix exhibit a power-law distribution with a similar exponent as for the scalar case. These eigenvalues however also exhibit strong correlations, as is typical in random matrix theory [START_REF] Mehta | Random matrices[END_REF][START_REF] Forrester | Log-Gases and Random Matrices[END_REF][START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Allez | Eigenvector dynamics: General theory and some applications[END_REF]. We will also provide some exact results in the limit of a large size of the matrix both for the discrete recursion and its continuum limit. These results unveil a connection to a classical random matrix ensemble, the inverse-Wishart ensemble, well studied in various other contexts [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF][START_REF] Allez | Invariant β-Wishart ensembles, crossover densities and asymptotic corrections to the Marcenko-Pastur law[END_REF].

Here we also unveil a remarkable connection between the continuum limit of the matrix Kesten recursion and the quantum mechanics of N fermions in a Morse potential. The eigenvalues of the matrix are related to the positions of the fermions. The large time limit is obtained by considering the ground state of these fermions which can be obtained exactly in both cases. For Hermitian matrices β = 2 the fermions are non-interacting and the full dynamics at finite time can be obtained in a simple way. Such fermion problems are important in the context of cold atoms in traps, and the present results thus add the Morse potential to the list of potentials in one dimension for which the quantum correlations can be computed using their connection to random matrices [START_REF] Mehta | Random matrices[END_REF][START_REF] Forrester | Log-Gases and Random Matrices[END_REF][START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF][START_REF] Lacroix-A-Chez-Toine | Statistics of fermions in a d-dimensional box near a hard wall[END_REF][START_REF] Lacroix-A-Chez-Toine | Non-interacting fermions in hard-edge potentials[END_REF][START_REF] Nadal | Nonintersecting Brownian interfaces and Wishart random matrices[END_REF][START_REF] Cunden | Free fermions and the classical compact groups[END_REF].

Our study also has connections with previous works on matrix Brownian functionals by Rider and Valkó in probability theory and by Grabsch and Texier in the study of mesoscopic quantum transport in disordered wires. Our model, in its continuum version, leads to a different matrix dynamics compared to the ones considered in the aforementioned works. Remarkably however, we show that the evolution of the eigenvalues can be made to coincide by choosing the appropriate discretization (Itô versus Stratonovich). Hence we find that two different models with two different noise terms and two different discretization prescriptions lead to the exact same dynamics for the eigenvalues.

Let us now introduce the main objects of interest and summarize the known results.

A. Overview of the Kesten recursion

The scalar Kesten recursion is the following system of equations on the sequence (Z n ) n 0 :

Z 0 = 0 Z n = ξ n (1 + Z n-1 ) (1)
where (ξ n ) n 0 is a sequence of independent, identically distributed positive random variables. The solution of this recursion is formally:

Z n = n-1 j=0 ξ n • • • ξ n-j = n-1 j=0 n k=n-j ξ k . (2)
Note that for a given integer n, the iid hypothesis gives, by permuting the ξ k variables as (0, 1, • • •, n) → (n, • • •, 1, 0), the following equality in law:

Z n law = n-1 j=0 j k=0 ξ k (3)
Harry Kesten studied the limit as n → ∞ of such recursions in [START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF][START_REF] Kesten | Spitzer A limit law for random walk in a random environment[END_REF][START_REF] Kesten | Convergence in distribution of products of random matrices[END_REF]. Under certain conditions, depending on the sign of E[log ξ], the sum in (3) may either grow unboundedly with n, if E[log ξ] 0, or, if E[log ξ] < 0, converge to a positive random variable denoted Z ∞ . Kesten showed that in that case the tail of the distribution of Z ∞ exhibits a power-law behaviour with exponent 1 + ν such that

E [ξ ν ] = 1 . (4)
It is easy to see why this condition arises. From the definition (1) one must have that in law Z ∞ law = ξ(1 + Z ∞ ), where in the r.h.s ξ and Z ∞ are taken as independent random variables. Denoting p(ξ) the probability density function (PDF) of ξ, we see that the PDF P (Z) of Z = Z ∞ must obey the integral equation

P (Z) = dZ P (Z ) dξp(ξ)δ(Z -ξ(1 + Z )) = dξ ξ p(ξ)P ( Z ξ -1) (5)
If one assumes that P (Z) ∼ Z -(1+ν) at large Z, the relation [START_REF] Prigodin | One-dimensional disordered system in an electric field[END_REF] follows. This and related integral equations admit a variety of possible behaviors which were studied in details in physics [START_REF] De Calan | On the distribution of a random variable occurring in 1D disordered systems[END_REF] and in mathematics [START_REF] Goldie | Implicit renewal theory and tails of solutions of random equations[END_REF]. Large deviations and rate of convergence were also studied in [START_REF] Buraczewski | Large deviations for solutions to stochastic recursion equations under Kesten's condition[END_REF][START_REF] Buraczewski | On the rate of convergence in the Kesten renewal theorem[END_REF].

The Kesten recursion has a close connection to the problem of a Random Walk in a Random Environment (RWRE) [START_REF] Solomon | Random walks in a random environment[END_REF][START_REF] Kesten | Spitzer A limit law for random walk in a random environment[END_REF][START_REF] De Calan | On the distribution of a random variable occurring in 1D disordered systems[END_REF][START_REF] Sinai | The limiting behavior of a one-dimensional random walk in a random medium[END_REF][START_REF] Derrida | Classical diffusion on a random chain[END_REF][START_REF] Derrida | Velocity and diffusion constant of a periodic one-dimensional hopping model[END_REF][START_REF] Kesten | The limit distribution of Sinai's random walk in random environment[END_REF][START_REF] Bouchaud | The relaxation-time spectrum of diffusion in a one-dimensional random medium: an exactly solvable case[END_REF][START_REF] Doussal | First-passage time for random walks in random environments[END_REF][START_REF] Comtet | Exact results on Sinai's diffusion[END_REF][START_REF] Le Doussal | Random walkers in one-dimensional random environments: exact renormalization group analysis[END_REF]. The RWRE setting is the study of a random walk where the hopping probabilities are random variables that constitute a random environment, that the random walk explores. This problem has been widely studied on the lattice, and in the continuum. Let us explain how the Kesten recursion occurs from RWRE in one dimension, following an argument by Solomon [START_REF] Solomon | Random walks in a random environment[END_REF]. We denote by α n = P n,n+1 and β n = P n,n-1 = 1 -α n the transition probabilities of a random walk on the Z lattice, from site n to its two neighbours. In the random environment setting, these are random variables. We suppose the sequence (α n , β n ) n∈Z to be fixed, such that a random environment has been chosen beforehand. Let µ n be the mean first passage time from site n to site n + 1 for a random walk in this environment: the first passage time is 1 with probability α n , if the first step from site n is taken to the right, or (1 + µ n-1 + µ n ) in mean, with probability β n , if the first step is taken to the left:

µ n = α n + β n (1 + µ n-1 + µ n ) (6)
Defining ξ n = βn αn , a prominent random variable to understand the behaviour of the RWRE, and Z n = 1 2 (µ n -1), we obtain directly the Kesten recursion Z n = ξ n (1 + Z n-1 ), with the correct initialization Z 0 = 0 if we choose α 0 = 1 such that the random walk is constrained to the positive sites. The mean first passage time µ n from site n to site n + 1 is thus distributed as the Kesten random variable. The quantity -E[log ξ] thus representes the effective bias of the environment. If it is positive, the random walker typically moves to the right, and µ ∞ is a finite random variable, which can be large with a power-law distribution corresponding to all possible explorations backward from site n before reaching site n + 1.

The Kesten recursion also appears in a simplified model of Directed Polymers in Random Media [START_REF] Derrida | Polymers on disordered trees, spin glasses, and traveling waves[END_REF], where the polymer lives on a complete graph [START_REF] Bouchaud | Wealth condensation in a simple model of economy[END_REF]. The time evolution of the partition function of polymers ending on site a = 1, • • • , N at time n is given by:

z a n+1 = η a n (1 -ε)z a n + ε N N b=1 z b n , (7)
where N is the total number of sites and ε the hopping rate of the polymer. On the complete graph, the Laplacian takes a mean-field form that allows one to simplify the analysis considerably. Indeed, introducing z n = b z b n /N , one can then show that for large N , the rescaled partition function Z a n := (1 -ε)z a n /εz n obeys a Kesten recursion:

Z a n+1 = ξ a n (1 + Z a n ), (8)
with ξ n := (1 -ε)η n /E[η], see [START_REF] Bouchaud | Wealth condensation in a simple model of economy[END_REF] for details. The Kesten recursion is again useful in the study of the random-field Ising chain at low temperature, as shown by Derrida and Hilhorst in [START_REF] Derrida | Singular behaviour of certain infinite products of random 2×2 matrices[END_REF]. In this setting, the free energy of the chain involves infinite products of random 2 x 2 matrices and can be expressed as F (ε) = lim N →∞ (1/N ) log (a N /a 0 ), where the sequence of random variable ( ai+1 ai ) i satisfies a recursion, which can be mapped to the Kesten recursion in the low temperature limit. A similar connection to quantum models, such as random Dirac operators in 1D, was studied in e.g. [START_REF] Comtet | On the distribution of the Wigner time delay in one-dimensional disordered systems[END_REF][START_REF] Texier | Universality of the Wigner time delay distribution for one-dimensional random potentials[END_REF][START_REF] Steiner | Statistical properties of a localization-delocalization transition in one dimension[END_REF].

In [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF], Dufresne exhibits the Kesten recursion in an investment management problem: consider that a savings account is credited every year by a unit currency. Define the random variable ξ n = 1 + R n , with R n the rate of return during the nth year. The accumulated value in the account S n is then given as a function of the preceding value S n-1 , through the Kesten recursion:

S n = ξ n (1 + S n-1 ) . (9)
Finally, these recursions also appeared early on in biosciences [START_REF] Paulson | Limit laws of a sequence determined by a random difference equation governing a one-compartment system[END_REF][START_REF] Cavalli-Sforza | Models for cultural inheritance, I. Group mean and within group variation[END_REF].

B. Scalar, vector and matrix models

The discrete recursion defined in (1) was studied for scalars and vectors by Kesten in [START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF], in the slightly more general form Z n = ξ n Z n-1 + q n . We will be interested throughout this work only in the special case where q n = ξ n . In the scalar case, as mentioned above, the variable Z ∞ is in the domain of attraction of a stable law and there exists an exponent ν such that the following limit exists and is nonzero:

0 < lim t→∞ t ν P (Z ∞ > t) < ∞ (10)
Furthermore, ν is characterized by (4) (Th. 5 of [START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF]). De Calan et al. give in [START_REF] De Calan | On the distribution of a random variable occurring in 1D disordered systems[END_REF] an argument for this characterization by studying poles of the Mellin transforms associated with the relevant distributions in this problem.

The vector case for this discrete recursion was also studied by Kesten, promoting Z n and q n to d-dimensional vectors and ξ n to a d × d matrix. It is shown, under suitable hypothesis, that in this case also there exists ν such that for any unit d-dimensional vector x the following limit exists and is nonzero:

0 < lim t→∞ t ν P x T . Z ∞ > t < ∞ (11)
In this case however ν is not characterized by a simple moment equation.

A continuous version of the scalar recursion was studied in physics as a Langevin equation describing the diffusion of a particle in a one-dimensional white noise Gaussian random field [START_REF] Bouchaud | Classical diffusion of a particle in a one-dimensional random force field[END_REF] (see section 4.3 there), the interpretation of the Kesten variable being the mean sojourn time of a particle in some region of space. It was studied in mathematics at about the same time by Dufresne in [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF][START_REF] Dufresne | On the integral of geometric Brownian motion[END_REF]. In these works the exact distribution was obtained as an inverse Gamma law (see below). Note that, by contrast, in the discrete case, apart from special cases [START_REF] De Calan | On the distribution of a random variable occurring in 1D disordered systems[END_REF], only the power-law behaviour is accessible. One way to obtain a meaningful continuous version is to consider the variable U t = Z n= t dt dt, i.e. the Kesten variable multiplied by the vanishing time interval. In the context of RWRE it corresponds to going from number of steps to continuous time duration. Writing the random variable ξ in exponential form as ξ n = e ηn , the reordered equivalent form (3) gives for the variable U at time t = ndt :

U t law = n-1 j=0 e j k=0 η k dt dt→0 -→ t 0 e As ds ( 12
)
where the process A s is the continuous limit of the random walk j k=0 η k . In order to ensure a well-defined limit for the random walk, the iid η random variables must scale as:

η = γdt + α √ dtX, X ∼ N (0, 1) (13)
where X follows the standard centered normal law. In that case, the limit process is a Brownian motion with a drift A t = γt + αW t , where we denote throughout this work W t as the (real scalar) standard Brownian motion. The infinite-time limit of the U variable can then be written as the following Brownian exponential functional:

U ∞ law = ∞ 0 e γt+αWt dt ( 14
)
where γ = E [log ξ] /dt and α 2 = V [log ξ] /dt have well defined limits. Such objects have been widely studied in the mathematical literature [START_REF] Matsumoto | Exponential functionals of Brownian motion, I: Probability laws at fixed time[END_REF][START_REF] Matsumoto | Exponential functionals of Brownian motion, II: Some related diffusion processes[END_REF]. For some applications in physics see e.g. [START_REF] Comtet | Exponential functionals of Brownian motion and disordered systems[END_REF][START_REF] Comtet | Functionals of Brownian motion, localization and metric graphs[END_REF].

It is also possible to derive the stochastic evolution of U from the random evolution of Z. Note that the scaling for η shows that one must choose the i.i.d random variables ξ of the following form in the continuous limit:

ξ = 1 + (γ + α 2 2 )dt + α √ dtX, X ∼ N (0, 1) (15)
The random recursion

Z n = ξ n (1 + Z n-1 ) becomes U t+dt = (1 + (γ + α 2
2 )dt + αdW t )(dt + U t ) which, expanded up to order dt using dW 2 t = dt, gives the following Itô stochastic differential equation (SDE) for U t in the continuous limit, denoting m = γ + α 2 2 and σ = α √ 2

1 :

dU t = (1 + m U t ) dt + √ 2σ U t dW t (Itô) (16)
It is possible to obtain the exact stationary large time limit, e.g. using the associated Fokker Planck equation, i.e. the PDF of the variable U ∞ which exists whenever γ < 0 (or equivalently m < σ 2 ). One finds [START_REF] Bouchaud | Wealth condensation in a simple model of economy[END_REF][START_REF] Bouchaud | Classical diffusion of a particle in a one-dimensional random force field[END_REF][START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF] that the inverse of U ∞ follows a gamma law

U -1 ∞ ∼ Γ 1 - m σ 2 , σ 2 (17)
such that U ∞ itself follows an inverse-gamma law with the following PDF denoted by P

P (U ∞ ) = σ -2+2 m σ 2 Γ 1 -m σ 2 U -2+ m σ 2 ∞ e -1 σ 2 U∞ . ( 18
)
with the heavy tail exponent 1

+ ν = 2 -m σ 2 .
The aim of the present work is to propose a random matrix extension of the Kesten recursion, and to study it both in the discrete and continuous settings. In the continuous setting, the main result will be the convergence to the inverse-Wishart distribution which is the natural matrix extension of the inverse-gamma law, see Appendix C for a precise definition. Before introducing the model, let us mention results from relevant random matrix models. Until now only continuum models have been considered, while here we extend the discrete Kesten recursion to the matrix world. Rider and Valkó have extended some identities related to the Dufresne inverse-gamma law [START_REF] Bouchaud | Classical diffusion of a particle in a one-dimensional random force field[END_REF] to matrices in [START_REF] Rider | Matrix Dufresne identities[END_REF]. In particular they show that, defining the matrix geometric Brownian motion M t such that dM t = M t dH t + 1 2 -µ M t dt, where dH t are a matrix of N 2 independent standard Brownian motions (see Section IV A), the random matrix defined by the following integral

∞ 0 M s M T s ds (19)
follows an inverse-Wishart law. This model is a direct matrix extension of the integral expression of the continuous Kesten problem in Eq. ( 14) in which the integrand is a scalar geometric Brownian motion with drift. We note that Rider and Valkó further extend the geometric M -X Lévy and 2M -X Pitman theorems, by considering the

M -1 t t 0 M s M T s ds M -T t and M -1 t t 0 M s M T
s ds matrix processes. The matrix Dufresne property and the 2M -X theorems were also covered and extended by O'Connell in [START_REF] O'connell | Interacting diffusions on positive definite matrices[END_REF].

Another important matrix model in relation to this paper is the one introduced by Grabsch and Texier in [START_REF] Grabsch | Topological phase transitions in the 1D multichannel Dirac equation with random mass and a random matrix model[END_REF] to model the stochastic evolution of the Riccati matrix associated to a random-mass Dirac equation on the half-line, as a model for a multichannel disordered wire. Furthermore, the same authors have recently recovered in [START_REF] Grabsch | Wigner-Smith matrix, exponential functional of the matrix Brownian motion and matrix Dufresne identity[END_REF] a similar stochastic evolution for a different model, namely the Wigner-Smith time delay matrix of a disordered multichannel model obeying a Schrödinger equation with random potential (see also [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF]). In this last setting, their model has very close connections with the integral representation of Rider and Valkó. In fact, the general problem of the Wigner-Smith time delay in chaotic cavities has close connections with the present work, and power law tails have also been obtained in that problem [START_REF] Texier | Wigner time-delay distribution in chaotic cavities and freezing transition[END_REF][START_REF] Fyodorov | Statistics of resonance poles, phase shifts and time delays in quantum chaotic scattering: Random matrix approach for systems with broken time-reversal invariance[END_REF]. For a general review on the Wigner-Smith time delay see [START_REF] Texier | Wigner time delay and related concepts: Application to transport in coherent conductors[END_REF].

A precise explanation of the difference between the matrix stochastic model introduced in this work with the ones from Rider-Valkó and Grabsch-Texier will be adressed in section IV.

C. Definition of the matrix Kesten recursion

For N ∈ N * , let us define the Kesten recursion for a sequence of symmetric positive semidefinite N × N matrices (Z n ) n 0 , with (ξ n ) n 0 a sequence of independent identically distributed (iid) symmetric positive semidefinite N × N matrices and I the N × N identity matrix, as:

Z n+1 = I + Z n ξ n I + Z n (20)
We have introduced a real parameter > 0 in the model, which will be given a value of order 1 in the discrete setting, and taken to 0 in the continuous limit. We notice that the symmetric positive semidefinite property of the sequence ξ implies by recursion the same property for the sequence (Z n ) n 0 , as we denote for each n by √ I + Z n the principal square-root of I + Z n . The properties of ξ n and Z n can be summed up by stating that they are symmetric matrices with only nonnegative eigenvalues. One expects that for large n the Kesten matrix Z n tends to a random matrix Z ∞ and one of the aim of this paper is to study the properties of the distribution of Z ∞ .

A generalization where all matrices are Hermitian positive semidefinite will be considered as well. We define the usual Dyson index to be β = 1 in the real symmetric case and β = 2 in the complex Hermitian case.

The initial motivation for studying such a matrix generalization of the Kesten problem arose from an attempt to define a matrix analogue of the Directed Polymer problem, for which the Kesten recursion is a mean-field version (as we explained in the previous section).

D. Outline

The outline of this work is as follows. In section II, we study the continuous limit of this model for a fixed arbitrary integer N . We obtain the stationary distribution for the Kesten matrix in section II A and provide a mapping to a quantum problem of interacting fermions in the Morse potential in section II B. The case β = 2 corresponding to noninteracting fermions is studied in details in section II C both for finite N and in the large N limit. In addition to their relation to the matrix Kesten recursion, the results in that section are of interest for the physics of trapped fermions. The section III is devoted to the large N limit, where free probability results can be obtained for the discrete model, as well as for the continuous limit. We characterize the stationary random matrix Z ∞ by its S-transform, which obeys the self-consistent equation (103). In the continuum limit a resolvent analysis is performed which shows agreement with the finite N results. In section IV we discuss the links with the Rider-Valkó and Grabsch-Texier models. In particular we show in section IV C that the flow of the eigenvalues coincide provided appropriate correspondence is performed between the discretization prescriptions of the stochastic equations. Finally the six appendices contain additional details and information.

II. FINITE N MATRIX KESTEN RECURSION, IN CONTINUOUS TIME

In this section we study the continuous limit of the Kesten recursion for N × N matrices, with N a fixed arbitrary integer. We consider here only Gaussian noise and restrict to strictly positive definite matrices ξ n and Z n (including the initial condition Z 0 ). Remembering the picture laid out in section I B for the scalar case, we choose the following structure for the matrices ξ n which ensures their positive definite nature

ξ n = exp γ I + σ √ B n (21)
where is the parameter introduced in [START_REF] Forrester | Log-Gases and Random Matrices[END_REF].

The definition of the iid noise matrices B n will be different in the real and complex cases. If β = 1, the matrices B n are independent and sampled from the Gaussian Orthogonal Ensemble (GOE), i.e. B n = (H n + H T n )/ √ 2 where the entries of H n are N 2 independent real variables following N (0, 1). If β = 2, the matrices B n are independent and sampled from the Gaussian Unitary Ensemble (GUE), i.e. B n = (H n + H T n )/2 + i( Hn -HT n )/2 where H n and Hn are independent samples with the same distribution as above. The linear combinations are chosen such that the off-diagonal entries of B n are real (if β = 1) or complex (if β = 2) Gaussian variables with variance 1 in both cases. The PDF of the matrix B n is rotationally invariant and proportional to exp -β TrB n B † n /4 . With this choice the spectrum of B n / √ N converges in density at large N to a semi-circle of support [-2, 2] for any β. In summary:

B n = H n + H T n √ 2β + i Hn -HT n 2 δ β=2 (22)
Using the above definitions for B n , the small-expansion of the matrix exponential definition for ξ n gives to order :

ξ n = (1 + m )I + σ √ B n (23)
where parameters m = γ + 1 2 σ 2 (N + δ β=1 ) and σ respectively tune the mean and the noise in ξ n . Indeed for the sake of the continuum limit obtained below we can replace B 2 n → I(N + δ β=1 ). Note that for N = 1, (23) with = dt identifies with [START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF]. With the definition [START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF], the Kesten recursion for Z, Eq. ( 20), becomes:

Z n+1 = (1 + m )( I + Z n ) + σ √ I + Z n B n I + Z n (24)
To first order in , the square-root of the sum of two commuting matrices can be expanded as:

I + Z n = Z n I + Z -1 n = Z n I + 2 Z -1 n + O( 2 ) = Z n + 2 Z -1 n + O( 2 ) ( 25
)
since we have restricted to the case where Z n has only strictly positive eigenvalues Finally, the evolution of Z reads:

Z n+1 = Z n + I + m Z n + σ √ Z n B n Z n + O( 3/2 ) (26)
In the → 0 limit, we define a continuous time matrix process U t = Z t/ , defining time as t = n . Notice that the insertion of in the definition of the recursion avoids the need for a rescaling between U and Z. We obtain the time evolution for U t as the following Itô matrix SDE:

dU t = (I + mU t ) dt + σ U t dB t U t (Itô) (27)
with B t = (H t + H T t )/ √ 2β + i( Ht -HT t )/2δ β=2 where the entries of H t and Ht are 2N 2 independent standard Brownian motions. We remark that taking N = 1 in the real case β = 1 recovers exactly the continuous Kesten evolution [START_REF] Solomon | Random walks in a random environment[END_REF] as expected, because the diagonal terms in the matrix B t are √ 2 times a standard Brownian. Fixing N = 1 in the complex case β = 2 requires a rescaling of σ by √ 2 to fit the scalar case because a diagonal term of B t is then a standard Brownian motion.

One method to study the matrix stochastic equation ( 27) is to write the evolution of the eigenvalues of U t . The flow that we study is such that U t is non-degenerate (see below) and its eigenvalues, denoted {λ i (t)} 1 i N , are distinct and positive. One can then use non-degenerate perturbation theory, as detailed in Appendix A. The evolution of the set {λ i (t)} 1 i N is obtained as the joint stochastic equations, for all 1 i N , see (A13),

dλ i = (1 + mλ i ) dt + σ 2 1 j N j =i λ i λ j λ i -λ j dt + 2 β σλ i dW i (Itô) (28)
where the (W i (t)) 1 i N are N independent standard Brownian motions. Note that the evolution of the eigenvalues is obtained with no a priori knowledge of the evolution of the eigenvectors. As compared to the case N = 1, which recovers the scalar continuous Kesten evolution, we note the interaction between the eigenvalues in the form λ i λ j /λ i -λ j for each pair. Note the similarity with the usual Dyson Brownian motion (DBM) [START_REF] Dyson | On the Brownian-motion model for the eigenvalues of a random matrix[END_REF]. However the difference is that here, because of the multiplicative noise, the repulsion depends also multiplicatively on the eigenvalues. From the stochastic equations [START_REF] Nadal | Nonintersecting Brownian interfaces and Wishart random matrices[END_REF] one can obtain the evolution of the joint probability distribution function (JPDF) P ( λ, t) of the set of eigenvalues {λ i } 1 i N denoted by the vector λ. One finds that the Fokker-Planck equation associated with the eigenvalue evolution of the Itô process U t is:

∂P ( λ, t) ∂t = 1 2 N i=1 ∂ 2 ∂λ 2 i 2 β σλ i 2 P ( λ, t) - N i=1 ∂ ∂λ i         1 + mλ i + σ 2 1 j N j =i λ i λ j λ i -λ j     P ( λ, t)     (29)
where the initial condition is specified by the eigenvalues of U 0 (which we keep arbitrary at this stage).

A. Large-time limit: stationary solution

In this section we obtain the stationary solution of the [START_REF] Cunden | Free fermions and the classical compact groups[END_REF]. The equation ( 29) can be written as ∂P ( λ,t) ∂t = i ∂ λi J i , which defines the components of the current (or flux) J i . The stationary solution P stat ( λ) is such that i ∂ λi J i = 0, and here it will be found as the normalizable solution of the zero-flux condition J i = 0 for all i, namely 0

= J i = ∂ ∂λ i 1 β (σλ i ) 2 P stat -     1 + mλ i + σ 2 1 j N j =i λ i λ j λ i -λ j     P stat (30)
Inspired by the DBM we will look for a solution of the form

P stat ( λ) = F ( λ) ∆( λ) β (31)
in terms of ∆( λ) = 1 i<j N (λ j -λ i ) the determinant of the Vandermonde matrix (λ j-1 i ) 1 i,j N . Using the following identities:

∂ λi log |∆( λ)| = j =i 1 λ i -λ j ; j =i λ 2 i λ i -λ j = (N -1)λ i + j =i λ i λ j λ i -λ j , (32)
where the second one is used when computing the derivative in the first term in [START_REF] Kesten | Convergence in distribution of products of random matrices[END_REF], one finds that F ( λ) must obey the set of equations for 1 i N

σ 2 β λ 2 i ∂ λi F = (1 + λ i (m - 2σ 2 β -σ 2 (N -1)))F ( 33
)
where the interaction term has cancelled. Integrating this equation one finally obtains the stationary measure for any N as

P stat ( λ) = K ∆( λ) β N i=1 λ β( m σ 2 -N +1)-2 i e -β σ 2 N i=1 1 λ i (34)
with K a normalization constant. We note that for N = 1 the Vandermonde factor is absent and setting β = 1 one recovers exactly the inverse gamma law for λ 1 ≡ U ∞ , i.e. Eq. ( 18), β = 1 (and also for β = 2 setting in [START_REF] Buraczewski | On the rate of convergence in the Kesten renewal theorem[END_REF] σ → σ √ 2). One recognizes in [START_REF] Buraczewski | On the rate of convergence in the Kesten renewal theorem[END_REF] the joint PDF of the eigenvalues of the inverse-Wishart ensemble. We recall the (white) inverse-Wishart matrix distribution over N × N real (β = 1) or Hermitian (β = 2) positive matrices, with parameters T, C

P iw (M ) = Q T,N,C,β det(M ) -β 2 (T +N -1)-1 e -C β 2 tr M -1 (35)
and its joint eigenvalue distribution:

P iw ( λ) = K T,N,C,β ∆( λ) β N i=1 λ -β 2 (T +N -1)-1 i e -C β 2 N i=1 1 λ i (36)
where the definition and normalization constants are detailed in Appendix C. We denote q = N T . We thus see in (34) that P stat is the Inverse-Wishart eigenvalue joint distribution with the following correspondence of parameters:

                     T = N -2m σ 2 + 2 β -1 1 q = T N = 1 -1 N 2m σ 2 -2 β + 1 C = 2 σ 2 K = K T,N,C,β = 1 N ! σ -N (β(N-1-2m σ 2 )+2) Γ( β 2 ) N π β 2 N (N -1) Γ N,β ( β 2 N ) Γ N,β ( β 2 (N -2m σ 2 -1)+1) (37)
The eigenvector-basis rotation invariance of dB t in [START_REF] Lacroix-A-Chez-Toine | Non-interacting fermions in hard-edge potentials[END_REF] then allows us to conclude that the full matrix distribution P stat (U ) is a stationary distribution for this model, with:

P stat (U ) = Q det(U ) β( m σ 2 -N +1)-2 e -β σ 2 tr U -1 (38)
where

Q = Q T,N,C,β = σ -N ( β(N -1-2m σ 2 )+2 ) Γ N,β ( β 2 (N -2m σ 2 -1)+1
) . We note that this distribution is well-defined only when the arguments of the Γ functions in the denominator of K T,N,C,β are positive in order to avoid the poles of Γ at negative integers: 1 -

β m σ 2 + β 2 (j -1) > 0 , ∀j ∈ [1, N ].
As a consequence, the normalizability condition for the stationary distribution does not depend on N :

m < σ 2 β ( 39
)
It is indeed similar to the N = 1 case obtained in continuous scalar Kesten evolution ( 16), as can be seen in the real case β = 1.

Another way to interpret this constraint is that the Inverse-Wishart distribution is well-defined only if the corresponding Wishart samples are almost surely invertible, i.e. they have full-rank. Since Wishart samples are obtained as X † X where X is a T × N matrix with independent Gaussian entries, they have full-rank only if T > N -1. Again, this yields [START_REF] Bouchaud | The relaxation-time spectrum of diffusion in a one-dimensional random medium: an exactly solvable case[END_REF].

In order to shed light on the stationary solution, we detail the marginal density for the largest eigenvalue λ 1 . In the limit where λ 1 is large, for a fixed value of N , its distribution reads

P (λ 1 ) λ1→∞ ≈ N (2/σ 2 ) 1-βm σ 2 Γ(1 + β 2 ) Γ(1 + β 2 (N -2m σ 2 )) Γ(1 + β 2 N ) Γ(1 -βm σ 2 ) Γ(1 -βm σ 2 + β 2 ) λ -2+β m σ 2 1 , (40)
see Appendix C for the derivation. We see that the tail exponent of the maximal eigenvalue distribution is independent of N and coincides with the N = 1 scalar Kesten result [START_REF] Bouchaud | Classical diffusion of a particle in a one-dimensional random force field[END_REF] as can be seen simply by injecting β = 1. We will show below (see section II C 2 b) that, in the limit of large N and fixed σ, this tail is related to the universal statistics at the hard-edge of the Wishart random matrix ensemble, which is given by the Bessel determinantal point process.

Note here some connections with the statistics of the Wigner-Smith time delay studied in the context of chaotic cavities [START_REF] Brouwer | Quantum mechanical time-delay matrix in chaotic scattering[END_REF][START_REF] Brouwer | Distribution of the quantum mechanical time-delay matrix for a chaotic cavity[END_REF][START_REF] Cunden | Large-N expansion for the time-delay matrix of ballistic chaotic cavities[END_REF]. In Ref. [START_REF] Texier | Wigner time-delay distribution in chaotic cavities and freezing transition[END_REF], the proper time delay τ was studied. It has the same distribution as the following sum over the eigenvalues λ i studied here, i.e. τ = 1 N N i=1 λ i . The parameters used in [START_REF] Texier | Wigner time-delay distribution in chaotic cavities and freezing transition[END_REF] correspond here to σ 2 = 2 and m = -N . With these parameters we see that our result (40) exhibits the same power law as the one conjectured in [START_REF] Fyodorov | Statistics of resonance poles, phase shifts and time delays in quantum chaotic scattering: Random matrix approach for systems with broken time-reversal invariance[END_REF] and obtained in [START_REF] Texier | Wigner time-delay distribution in chaotic cavities and freezing transition[END_REF] in the large N limit from a Coulomb gas calculation. This suggests that the sum is dominated by the largest eigenvalue. We will return to the large N limit below.

B. Finite-time solution: Morse-Sutherland quantum mapping

In the case of the usual Dyson Brownian motion, there exists a mapping onto a quantum system, the Calogero-Sutherland model [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF]. In the present case, as we now show, there also exists such a mapping, albeit onto a different quantum model. To elicit this mapping we first perform a change of variables by defining

λ i = e xi . (41)
From Itô rules, the stochastic evolution of the x variables has additive noise following:

dx i = F i dt + 2 β σdW i F i := e -xi + m - σ 2 β + σ 2 1 j N j =i 1 e xi-xj -1 , (42)
with an initial condition x i = x 0 ∀i ∈ [1, N ] at t = 0 with x 0 a large negative constant, and F i the force felt by particle x i . Let us denote P ( x, t) the probability distribution at time t for the set of x variables, obtained from P ( λ, t) through the change of variables. The Fokker-Planck equation verified by P ( x, t) is:

∂ P ∂t = N i=1 σ 2 β ∂ 2 P ∂x 2 i - ∂ ∂x i (F i P ) = -ĤF P P (43)
where we have introduced the generator ĤF P . One can map this generator to a Schrödinger problem in the fashion of [START_REF] Gautié | Non-crossing Brownian paths and Dyson Brownian motion under a moving boundary[END_REF] as follows. Denoting Pstat ( x) the distribution obtained from P stat ( λ) under this change of variables, we define ψ( x, t) such that:

P ( x, t) = Pstat ( x) 1 2 ψ( x, t) (44)
Injecting P into (43) yields the following Schrödinger equation on ψ( x, t):

∂ψ ∂t = N i=1 σ 2 β ∂ 2 ψ ∂x 2 i -β F 2 i 4σ 2 + 1 2 ∂F i ∂x i ψ = - σ 2 β ( Ĥ -E 0 )ψ ( 45
)
where the constant E 0 is given in [START_REF] Dufresne | On the integral of geometric Brownian motion[END_REF]. The stochastic system is hereby mapped to a quantum system with the following Hamiltonian:

Ĥ = i     - ∂ 2 ∂x 2 i + V (x i ) + 1 j N j =i V int (x i , x j )     (46)
with a 1-particle potential V of the Morse form Ae -2x -Be -x [START_REF] Morse | Diatomic molecules according to the wave mechanics. II. Vibrational levels[END_REF] and with Sutherland sinh -2 interaction potential V int [START_REF] Sutherland | Beautiful models: 70 years of exactly solved quantum many-body problems[END_REF]:

V (x i ) = β 2 e -2x i 4σ 4 -β σ 2 β 2 (N -1 -m σ 2 ) + 1 e -x i V int (x i , x j ) = β(β-2) 16 1 sinh x i -x j 2 2 (47)
We thus see that the stochastic generator ĤF P and the quantum hamiltonian Ĥ are related by the operator transformation

σ 2 β ( Ĥ -E 0 ) = Pstat ( x) -1/2 ĤF P Pstat ( x) 1/2 (48)
The eigenenergies E of ĤF P and E of Ĥ are also related through

E = σ 2 β (E -E 0 ) (49)
The lowest Fokker-Planck eigenvalue E is of course equal to 0, the eigenvalue of the stationary state, such that the constant E 0 introduced in ( 45) is indeed the ground-state energy of Ĥ and is computed from the correspondence to be

E 0 = -N (βm -σ 2 ) 2 4σ 4 + N (N -1) 4 ( β 2 σ 2 m -β β + 2 2 ) -N (N -1)(N -2) β 2 12 ( 50
)
It turns out that the ground state of the quantum problem is N ! degenerate [START_REF] Calogero | Solution of a three-body problem in one dimension[END_REF]. Indeed one can choose the sign of the wavefunction for each ordering of the x i . This can be seen from the common feature in Random Matrix Theory that the eigenvalues of matrix processes do not intersect. It can be verified in particular in [START_REF] Le Doussal | Random walkers in one-dimensional random environments: exact renormalization group analysis[END_REF]: fixing two indices i and j, let us assume that x i > x j at some time t; the evolution of the difference follows, when x i -x j → 0, with W a standard Brownian Motion:

d(x i -x j ) = (e -xi -e -xj )   1 + βσ 2 k =i,j e x k (e xi -e x k )(e xj -e x k )   dt + σ 2 e xi-xj + 1 e xi-xj -1 dt + 2 σ √ β d W (51) ∼ 2σ 2 x i -x j dt + 2 σ √ β d W , (x i -x j → 0) (52)
Scaling the time variable by an adequate factor, we see that this equation becomes the SDE 2XdX = βdt + 2Xd W , i.e. a (β + 1)-Bessel process describing the norm of a Brownian motion in β + 1 dimensions. For any β 1 this process never hits zero as a consequence of the recurrence properties of Brownian motion, and we can conclude that x i = x j for all times. As a consequence the wave-function ψ can be chosen to have the fermionic symmetry (i.e antisymmetric in all its arguments). Hence we are led to study a a system of N fermions in the Hamiltonian described above.

C. The case β = 2: non-interacting fermions in the Morse potential

In the special case β = 2, the mapping is simplified as the interaction term vanishes, leading to a system of N non-interacting spinless fermions in the Morse potential. From Eq. ( 41) the position of the fermions x i are related to the eigenvalues λ i of the Kesten matrix via x i = log λ i . Using the known spectrum of the single-particle model, we detail in section II C 1 the ground-state wavefunction for the system of N fermions. It is constructed as usual as a Slater determinant where the N single-particle states of energy below the Fermi energy F are occupied. This provides an interpretation of the stationary measure of the model studied in this paper, which is related to the inverse-Wishart ensemble. It also unveils a correspondence between the inverse-Wishart random matrix ensemble and non-interacting fermions in a Morse potential.

The present study thus adds the Morse potential to the list of potentials in one dimension which are related to random matrices [START_REF] Mehta | Random matrices[END_REF][START_REF] Forrester | Log-Gases and Random Matrices[END_REF][START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF][START_REF] Lacroix-A-Chez-Toine | Statistics of fermions in a d-dimensional box near a hard wall[END_REF][START_REF] Lacroix-A-Chez-Toine | Non-interacting fermions in hard-edge potentials[END_REF][START_REF] Nadal | Nonintersecting Brownian interfaces and Wishart random matrices[END_REF][START_REF] Cunden | Free fermions and the classical compact groups[END_REF]. It includes the harmonic oscillator (related to the GUE), the inverse square potential (related to the Wishart-Laguerre ensemble) and the Jacobi box potential (related to the Jacobi ensemble). Note however that a distinct feature of the Morse potential is that it can only accommodate a finite number of bound states, N b . We say that the Morse potential is full when all the bound states are occupied by the N fermions. We study in II C 2 this correspondence in the large N limit. There are two interesting scaling regimes depending how the parameter σ is chosen in that limit. If σ is scaled as O(1/ √ N ) the support of the density of the Fermi gas is a finite interval. If σ remains of O(1) in that limit, the Morse potential is almost full and the Fermi gas extends very far to the right. In that regime, we find that the fluctuations of the largest eigenvalues of the Kesten matrix, equivalently of the rightmost fermions in the Morse potential, are described by the so-called Bessel determinantal point process up to a change of variable.

Finally we write an exact formula for the full dynamics of the model in section II C 3. It is interesting to note that the integral functional of the geometric Brownian motion, which is equivalent to the continuous scalar Kesten recursion as we recalled in Section I, has a prominent interest in the field of quantitative finance [START_REF] Geman | Bessel processes, Asian options, and perpetuities[END_REF]. In the specific setting of the Asian option, the quantum dynamics in the Morse potential has proved useful to compute the fair price of the financial product in [START_REF] Zhang | Morse Potential, Contour Integrals, and Asian Options[END_REF]. Our present study can thus be viewed as an N -particle generalization of this analytical correspondence, although there is no clear financial applications of this generalization.

Energy spectrum and ground state of N fermions

The Morse potential is described by the following single-particle Hamiltonian (in units where the mass is set to 1/2 and = 1)

Ĥ1 = - ∂ 2 ∂x 2 + V (x) , V (x) = g 2 (e -2(x-x0) -2e -(x-x0) ) ( 53
)
where the parameters g and x 0 denote respectively the strength of the potential and the position of the minimum. The value at the minimum V (x 0 ) = -g 2 is negative and the potential tends to zero at x = +∞, see Fig 1 . For the Kesten matrix model studied in this paper the values of these parameters are fixed as

g = N - m σ 2 and e -x0 = σ 2 N - m σ 2 , (54)
and depend explicitly on N , the number of fermions, and this has important consequences discussed below. However the study of the fermion problem that we perform below is valid for general g and x 0 . It is shown in Appendix D that this potential exhibits a finite number

N b = g + 1 2
-of bound states, the largest integer strictly smaller than g + 1 2 , as well as a continuous branch of diffusion states with positive energies. The bound states in the discrete branch of the spectrum have the following eigenenergies k and eigenfunctions Ψ k with

0 k < g -1 2 : k = -g -k - 1 2 2 = -N - m σ 2 -k - 1 2 2 (55) ψ k (x) = N k 2 σ 2 e -x N -m σ 2 -k-1 2 e -1 σ 2 e -x L (2N -2m σ 2 -2k-1) k 2 σ 2 e -x (56)
where

L (α) k
is the generalized Laguerre polynomial and

N k = k!(2g-2k-1) Γ(2g-k) 1 2 .
The ground-state wavefunction of the N -fermion model, denoted 0, is then simply the Slater determinant obtained from the N lowest 1-particle eigenstates k = 0, 1,

• • • , N -1: Ψ 0 ( x) = 1 √ N ! det 1 i,j N (ψ i-1 (x j )) (57)
E 0 = N -1 k=0 k (58)
It can be verified as expected that E 0 is exactly the ground-state energy obtained by inserting β = 2 in [START_REF] Dufresne | On the integral of geometric Brownian motion[END_REF]. It can also be verified that the stationary solution of the matrix Kesten problem corresponds to the ground state through Pstat ( x) = Ψ 0 ( x) 2 . Indeed, we can develop the modulus square expression:

Ψ 0 ( x) 2 = 1 N ! N k=1 N 2 k-1 2 σ 2 e -x k 2(N -m σ 2 -1 2 ) e -2 σ 2 e -x k det 1 i,j N ( 2 σ 2 e -xj ) -i+1 L (2N -2m σ 2 -2i+1) i-1 2 σ 2 e -xj 2 (59)
The polynomials

X i-1 L (2N -2m σ 2 -2i+1) i-1
X -1 forming a family of increasing degree i -1 in X, the determinant in this expression is of the Vandermonde form. More precisely, the leading term in these polynomials have coefficient

2N -2m σ 2 -i i-1 , the constant term of L (2N -2m σ 2 -2i+1) i-1
, where the binomial is extended to non-integers by a standard analytic continuation. This gives, extracting constants from the Vandermonde:

Ψ 0 ( x) 2 = σ 2 /2 N (N -1) N ! N k=1 N 2 k-1 2N -2m σ 2 -k k -1 2 2 σ 2 e -x k 2N -2 m σ 2 -1 e -2 σ 2 e -x k ∆ ({e x }) 2 (60)
where the set {e xi } 1 i N is denoted here as {e x } instead of a vector, for ease of notation. Rearranging the constant terms, we obtain in this case β = 2:

Ψ 0 ( x) 2 = K ∆ ({e x }) β N k=1 (e x k ) β( m σ 2 -N +1)-1 e -β σ 2 e -x k , β = 2 ( 61
)
where K is the constant given earlier in [START_REF] Derrida | Velocity and diffusion constant of a periodic one-dimensional hopping model[END_REF]. Changing variables to λ = e x then gives exactly the stationary distribution obtained earlier in [START_REF] Buraczewski | On the rate of convergence in the Kesten renewal theorem[END_REF], for the case β = 2:

P stat ({λ}) = K ∆ λ 2 N k=1 λ -2N +2 m σ 2 k e -2 σ 2 1 λ k (62)
This shows that the modulus square of the ground-state wavefunction is indeed the stationary distribution on the x variables:

Pstat ( x) = Ψ 0 ( x) 2 .
As we mentioned above the Morse potential has a finite number of bound states

N b = g + 1 2
-. Hence for a generic value of g it can only accommodate N ≤ N b fermions in their ground state. Since for the matrix Kesten problem g = N -m σ 2 depends itself on N , we must check for consistency and for the existence of the solution. The ground-state wavefunction, as a Slater determinant of quantized states, exists if all states are well-defined i.e. N -1 < g - 1 2 , or equivalently m < σ 2 2 . We see that, in this case β = 2, this condition agrees perfectly with the condition obtained in the stochastic study for the normalizability of the stationary solution [START_REF] Bouchaud | The relaxation-time spectrum of diffusion in a one-dimensional random medium: an exactly solvable case[END_REF]. In the ground state, we further note that the Fermi energy F , the last occupied energy level, is independent of the number of particles:

F = N -1 = - 1 2 - m σ 2 2 ( 63
)
which holds only for m < σ 2 2 , and is strictly negative. Note that when -

1 2 m σ 2 < 1 2 the potential is full with N = N b , while for m σ 2 < -1 2 one has N < N b . The region of parameters m σ 2
2 corresponds in the fermion problem to a positive Fermi energy, in which case there are continuum eigenstates extending to infinity, and no finite-N ground state. It can be given a sense for finite N at finite time (see below) and in the Kesten problem it corresponds to runaway behavior of the Kesten matrix variable to infinity, with no stationary state.

The fermionic setting suggests a determinantal rewriting of the stationary solution, as in [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF][START_REF] Gautié | Constrained non-crossing Brownian motions, fermions and the Ferrari-Spohn distribution[END_REF]:

Pstat ( x) = 1 N ! det 1 i,j N (ψ i-1 (x j )) 2 = 1 N ! det 1 i,j N N -1 k=0 ψ k (x i )ψ k (x j ) = 1 N ! det 1 i,j N K N (x i , x j ) ( 64
)
where the kernel K N is:

K N (x, y) = (2/σ 2 ) 2(N -m σ 2 -1 2 ) N -1 k=0 σ 4k N 2 k 4 k e -(x+y)(N -m σ 2 -k-1 2 ) e -1 σ 2 (e -x +e -y ) L (2N -2m σ 2 -2k-1) k ( 2 σ 2 e -x )L (2N -2m σ 2 -2k-1) k ( 2 σ 2 e -y ) (65)
By the orthonormalization of the eigenfunctions ψ * k (x)ψ k (x)dx = δ k,k , this kernel satisfies the reproducibility property, which can be written as:

K N (x, z)K N (z, y)dz = K N (x, y) (66)
This property implies that the joint distribution Pstat describes a determinantal point process (DPP) with kernel K N . From standard results on DPP's [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF][START_REF] Johansson | Random matrices and determinantal processes[END_REF][START_REF] Borodin | Determinantal point processes. The Oxford HandBook of Random Matrix Theory[END_REF], this implies that all n-point correlations can also be written as determinants in terms of K N . In particular, the 1-particle density in the ground state P 1 (x) is given as

P 1 (x) := 1 N E N k=1 δ(x -x k ) = 1 N K N (x, x) (67)
= (2/σ 2 ) 2(N -m σ 2 -1 2 ) N N -1 k=0 σ 4k N 2 k 4 k e -2x(N -m σ 2 -k-1 2 ) e -2 σ 2 e -x L (2N -2m σ 2 -2k-1) k ( 2 σ 2 e -x ) 2 (68)
where in the definition the average is taken over the ground-state JPDF. This formula has been used to plot the fermion density in Figs. 3 and4. The above results were obtained for β = 2. In the case β = 1, which corresponds to fermions in the Morse potential in presence of Sutherland like interactions given in [START_REF] Dufresne | The distribution of a perpetuity, with applications to risk theory and pension funding[END_REF], the ground-state wavefunction is also obtained from (61) substituting β = 1.

Large N behavior

When N is large, the density of the Fermi gas in the potential V (x) (53) is well described in the bulk by the semi-classical limit, also called local density approximation (LDA) [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF], P 1 (x) P bulk (x). The LDA density P bulk is supported on an interval [x -, x + ] and vanishes at the edges:

P bulk (x) = 1 πN E eff -V (x) = ge x0 πN (e -x -e -x+ )(e -x--e -x ) (69)
where the edges x ± verify V (x ± ) = E eff . Here E eff is an effective Fermi energy defined such that the normalization condition 1 =

x +

x-P bulk (x)dx holds. One obtains the relations e -x-+ e -x+ = 2e -x0 and E eff = -g 2 e -(x++x-)+2x0 which leads to e -x± = e -x0 (1 ∓ 1 + E eff g 2 ), hence the positions of the edges are given by

x ± = x 0 -ln 1 ∓ 1 + E eff g 2 . ( 70
)
Now we note that one can rewrite (69) upon the change of variable Λ = 2ge -(x-x0) and using the above result for x ± as

P bulk (x)dx = g - √ -E eff N ρ MP (Λ)dΛ , Λ ± = 2ge -(x∓-x0) (71)
where

ρ MP (Λ) = 2 πΛ( Λ + -Λ -) 2 (Λ -Λ -)(Λ + -Λ) (72)
is the Marcenko-Pastur distribution, which is normalized to unity

Λ + Λ-dλρ MP (Λ) = 1, see Appendix C.
Hence the normalization condition for P bulk (x) immediately implies that the the effective Fermi energy is related to the number of fermions as

E eff = -(g -N ) 2 (73)
The appearance of the Marcenko-Pastur distribution originates from the fact that the variables Λ i = 2ge -(xi-x0) obey, in the fermionic ground state, the statistics of the eigenvalues of a Wishart matrix. This fact will be used again below. Inserting the result for E eff in Eq. ( 70) we find that the positions of the edges become

x ± = x 0 -ln 1 ∓ N g 2 - N g . ( 74
)
for N < g. One sees that as N/g → 1 -, the variable x + -x 0 diverges which means that, relative to the minimum of the well, the right edge moves to +∞. Hence one recovers the condition on the number of particles discussed above, which in the large N limit reads N < g. Thus to study the large N limit we need to choose the strength of the Morse potential g to scale as N .

We now apply these results to our Kesten matrix problem and use the parameter values g = N -m σ 2 and e -x0 = σ 2 g given in [START_REF] Comtet | Functionals of Brownian motion, localization and metric graphs[END_REF]. In that case the effective Fermi energy [START_REF] Gautié | Constrained non-crossing Brownian motions, fermions and the Ferrari-Spohn distribution[END_REF] becomes

E eff = - m 2 σ 4 (75)
At this point we should make an important distinction: depending on the scaling of σ with N , the asymptotic behavior of the system will be drastically different. ) and e -x0 = σ2 (1 -m σ2 ). Because of the exact condition m < σ 2 2 ∼ 1/N , here we must choose m < 0 which ensures that N < g is satisfied. We see that the effective Fermi energy (75) coincides in this small σ limit with the last occupied level [START_REF] Dyson | On the Brownian-motion model for the eigenvalues of a random matrix[END_REF]. In that regime we see from [START_REF] Johansson | Random matrices and determinantal processes[END_REF] that the asymptotic support of the fermion density is [x -, x + ] such that the positions of the two edges

x ± = x 0 -ln 1 ∓ √ 1-2 m σ2 1-m σ2
remain fixed and of order O(1), see Fig. 2. In the bulk, the density is given by the LDA in (69):

� � � � � -� � � � � � � � � � � -
P bulk (x) = 1 πσ 2 (e -x -e -x+ )(e -x--e -x ) (76)
At these points, the potential has the following derivative:

V (x ± ) = 2N 2 1 -2 m σ2 ± 1 - m σ2 -1 -2 m σ2 (77)
Near the two edges, denoting w ± = V (x ± ) -1/3 = O N -2/3 the length scale which characterizes the width of the edge regions, the density is correctly described by the one obtained from the Airy kernel such that [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF]:

P edge ± (x) 1 N |w ± | F 1 x -x ± w ± , F 1 (z) = Ai (z) 2 -z[Ai(z)] 2 (78)
The agreement of the 1-particle density with the approximations in the bulk and at the edges, in the case of the σ ∼ 1/ √ N scaling, is illustrated in Fig. 3. b. Constant σ. In the case where σ is a constant of order 1, the position of the minimum of the potential is asymptotically x 0 = -ln σ 2 N +O( 1N ). g scales as N , such that the left edge is at a fixed distance of the minimum x -= x 0ln 2+O( 1N 2 ), while the right edge is at diverging distance from the minimum

x + = x 0 + 2 ln σ 2 N -ln m 2 2 + O( 1 N
). In the translated frame where the position of the minimum is fixed, the Fermi gas only has one edge to the left, while extending infinitely to the right for N = +∞.

In the vicinity of the left edge, the Airy kernel characterizes the behavior of the Fermi gas. With V (x -) = -4g 2 -4N 2 and w = V (x -) -1/3 = O(N -2/3 ), the density is then described as above by:

P left edge (x) 1 N |w| F 1 x -x - w , F 1 (z) = Ai (z) 2 -z[Ai(z)] 2 (79)
In order to obtain the behavior of the Fermi gas to the far right, we will use again the mapping to the Wishart matrix ensemble. As mentioned above, it is easily seen that the ground-state JPDF (61) can be mapped, under the change of variables Λ = 2ge -(x-x0) = 2 σ 2 e -x , to the following Wishart eigenvalue JPDF of (fixed) parameter -2m σ 2 :

P Wishart ( Λ) = K ∆ Λ 2 N k=1 Λ -2m σ 2 k e -Λ k (80)
Under this mapping, the large x behavior can be accessed by studying the left edge of the Wishart model Λ 0. It is well-known that this hard-edge behavior in the Wishart ensemble is characterized by the Bessel kernel [START_REF] Dean | Noninteracting fermions in a trap and random matrix theory[END_REF][START_REF] Tracy | Level spacing distributions and the Bessel kernel[END_REF], such that the kernel K W of the DPP of the Λ i 's is, in the large N limit for Λ = O(1/N ):

K W (Λ 1 , Λ 2 ) 4N K Be,-2m σ 2 (4N Λ 1 , 4N Λ 2 ), K Be,ν (u, v) = √ v J ν ( √ v)J ν ( √ u) - √ u J ν ( √ u)J ν ( √ v) 2(u -v) (81)
where J ν is the Bessel function of order ν. This is equivalent to the fact that in the large N limit the eigenvalues near the edges can be written as Λ i bi

4N
where the b i 's form the Bessel DPP of Kernel K Be,ν . The Bessel DPP is defined as a limit process for

N = +∞ with b ∈ [0, +∞[. Its mean density ρ(b) = E[ +∞ i=1 δ(b -b i )] decays at large b as ∼ 1/ √ b.
Many of its properties are known, such as the PDF of the minimum b min = min i b i [START_REF] Tracy | Level spacing distributions and the Bessel kernel[END_REF]. We can now translate these properties to describe the fermions which are very far to the right: their positions can be written as

e -xi = σ 2 8N b i , x i = ln N -ln σ 2 b i /8 , x i = x + + ln 4m 2 σ 4 b i (82)
We deduce the fermion kernel K N (65) in the large x, y regime, ie x -x + = O(1) and y -x + = O(1) :

K N (x, y) 8N σ 2 e -(x+y)/2 K Be,-2m σ 2 8N σ 2 e -x , 8N σ 2 e -y (83)
and in particular the density:

P far right (x) 8 σ 2 e -x K Be,-2m σ 2 8N σ 2 e -x , 8N σ 2 e -x (84)
= 2 σ 2 e -x J -2m σ 2 2|m| σ 2 e -x-x + 2 2 -J -2m σ 2 +1 2|m| σ 2 e -x-x + 2 J -2m σ 2 -1 2|m| σ 2 e -x-x + 2 (85)
The fluctuations of the position of the rightmost fermion x max = x + + ln 4m 2 σ 4 bmin are of order 1 and can be obtained from the PDF of b min :

P bmin (s) = σ(s) s e -s 0 σ(x)
x dx [START_REF] Blaizot | Universal shocks in the Wishart random-matrix ensemble[END_REF] where σ satisfies a Painlevé III equation (Eq. (1.21) in [START_REF] Tracy | Level spacing distributions and the Bessel kernel[END_REF]). The PDF of b min vanishes at s = 0 as s -2m σ 2 . Taking into account the prefactor (Eq. (1.22) in [START_REF] Tracy | Level spacing distributions and the Bessel kernel[END_REF]), we deduce the tail of the distribution of the rightmost fermion as:

P max (x) x→∞ (m 2 /σ 4 ) 1-2m σ 2 Γ(1 -2m σ 2 )Γ(2 -2m σ 2 ) e -(1-2m σ 2 )(x-x+) (87)
Finally, we note that changing variables to λ = e x and injecting the value of x + given at the beginning of the paragraph, we recover exactly the tail of the marginal distribution of the largest eigenvalue λ 1 in (40) for the present case β = 2, with matching prefactor in the large N limit:

P (λ 1 ) λ1→∞ (2N/σ 2 ) 1-2m σ 2 Γ(1 -2m σ 2 )Γ(2 -2m σ 2 ) λ -2+ 2m σ 2 1 ( 88
)
The agreement of the 1-particle density with the approximations in the bulk and at the edges, in the case of constant σ, is illustrated in Fig. 4.

Note again that in the interacting case β = 1, the ground-state wavefunction, (61) with β = 1, can also be analyzed in the large N limit. It corresponds to an inverse Wishart matrix model with β = 1. For the case σ = O(1), the left edge is now described, upon the same change of variable, by the universality of the Gaussian orthogonal ensemble. The fermions on the far right, corresponding to the largest eigenvalues of the continuous matrix Kesten recurrence, are now described, upon this change of variable, by the β = 1 Bessel point process (which is not determinantal, see [START_REF] Ramirez | Diffusion at the random matrix hard edge[END_REF] for its description). 84). The inset shows the region near the right edge for N = 50, with σ = 2, m = -3 as above. For these parameters, we have x+ 3.81. We see that the exact density is very closely approximated by P far right (x) in that region.

Finite time results

The mapping of the matrix Kesten recursion onto the quantum problem allows us to study its time evolution. In the case β = 2 one can obtain an exact formula for this time evolution.

Taking into account the bound states and the continuous branch of the spectrum, the Euclidean propagator, or Green's function in imaginary time, of a 1-particle system in the Morse potential (53) is, see App. D:

G(x, x , t) = (g-1 2 ) - k=0 ψ k (x)ψ k (x )e -(g-k-1 2 ) 2 t + +∞ 0 dp 2gπ 2 p sinh(2πp)|Γ(ip -g + 1 2 )| 2 e x+x 2 -x0 W g,ip (2ge -(x-x0) )W g,ip (2ge -(x -x0) )e -p 2 t (89)
in terms of the Whittaker W function. From completeness, this propagator satisfies G(x, x , 0) = δ(x -x ). The propagator of the N -particle system of non-interacting fermions is then simply:

G N ( x, x , t) = x| e -t Ĥ | x = det 1 i,j N G(x i , x j , t) (90)
which satisfies G N ( x, x , 0) = δ( x -x ). Assuming a known initial position for the fermions x 0 , the N -particle wavefunction at a finite time t is ψ( x, t) = G N ( x, x 0 , t). From this quantum result, we obtain using Eq. ( 48) the finite time solution of the stochastic equation (42):

P ( x, t) = x|e -tH F P |x 0 = Pstat ( x) Pstat ( x 0 ) 1 2 G N ( x, x 0 , t)e E0t , P ( x, 0) = δ( x -x 0 ) (91)
Upon change of variable this gives directly the following finite-time solution for the stochastic evolution of the λ i = e xi variables:

P ( λ, t) = P ( x = ln λ, t) N k=1 1 λ i ( 92
)
where the vector ln λ denotes the set {ln λ i } i i N .

III. DISCRETE RECURSION IN THE LARGE N LIMIT

In the large N limit, we study the discrete matrix Kesten recursion [START_REF] Forrester | Log-Gases and Random Matrices[END_REF] with tools from free probability theory and through the analysis of the resolvent evolution. We treat first the case of the discrete recursion, and then study the continuum version. In the continuum case, the results obtained by these two methods are compared with the N → ∞ limit of the finite N results obtained in the previous section (noting that here we scale the parameter σ = O(1/ √ N ) in that limit). We note that the Dyson index β does not have the same significance in the free case, since all the information is contained in the eigenvalue density: results apply to both β = 1, 2 cases. In the last subsection, we complete the study with the analysis of the expected characteristic polynomial for arbitrary N .

A. Free probability approach

The property of freeness is a generalization of the concept of independence to non-commuting random variables, such as random matrices [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF][START_REF] Voiculescu | Symmetries of some reduced free product C*-algebras[END_REF][START_REF] Voiculescu | Limit laws for random matrices and free products[END_REF][START_REF] Voiculescu | Free probability theory: random matrices and von Neumann Algebras[END_REF][START_REF] Mingo | Free Probability and Random Matrices[END_REF][START_REF] Tulino | Verdú Random Matrix Theory and Wireless Communications[END_REF][START_REF] Novak | Three lectures on free probability[END_REF]. The field of free probability theory has seen the development of numerous results and tools relevant to the study of free variables, see Appendix B for a brief summary of the main transforms used in this work and [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF] for a precise introduction to freeness and detailed consequences for random matrices. The general result that we apply in this section to establish a connection to free probability is the following: two large symmetric matrices whose eigenbasis are randomly rotated with respect to one another are free, see [START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF].

Discrete matrix recursion

With help from free probability tools, applicable in the large N limit, we study the discrete matrix recursion of Eq. ( 20):

Z n+1 = I + Z n ξ n I + Z n (93)
We make two assumptions on the distribution of ξ: i) its support is of order 1 in the large N limit, ii) it is invariant under rotations. By the rotational invariance of ξ, I + Z n and ξ n in the Kesten recursion (20) become free random matrices in the N → ∞ limit. The model under study then defines Z n+1 as the free (symmetrized) multiplication of I + Z n and ξ n . Free multiplication is endowed with a very powerful property [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF]: the S-transform -defined below and in App. B -of the product matrix Z n+1 factorizes to the product of the S-transforms of I + Z n and ξ n :

S Zn+1 = S I+Zn × S ξ (94)
where the dependence on n in S ξ has been dropped as the ξ n matrices are iid. In the following, we use this property to write a recursion on S Zn , in order to characterize the evolution of the spectrum through the evolution of the S-transform. See App. B for more details on the free probability transforms used in this section. We introduce Y n = I + Z n for ease of notation. Let us express the S-transform of this variable in terms of S Zn .

For a matrix M , the eigenvalue density is defined as ρ := E 1 N λ∈Sp(M ) δ λ . For u ∈ R, the spectrum densities of Z n and Y n are related as:

ρ Yn (u) = ρ Zn (u -) (95)
For a matrix M , the Stieltjes transform is defined as g M (z) := ρ M (u) z-u du. For z ∈ C away from the support of ρ Yn , the Stieltjes transforms of Y n and Z n are related as:

g Yn (z) = ρ Yn (u) z -u du = ρ Zn (u) z --u du = g Zn (z -) (96)
The T -transform is defined from the Stieltjes transform as T (z) := zg(z) -1. As a consequence of (96):

T Yn (z) = z g Yn (z) -1 = z g Zn (z -) -1 = T Zn (z -) + g Zn (z -) (97)
We then have the following relationship between the T -transforms of Z n and Y n :

T Yn (z) = z T Zn (z -) + z - (98)
Defining finally the S-transform as S(ω) := ω+1 ωT -1 (ω) , we have S(T (z)) = T (z)+1 zT (z) . Applying this to Y n and injecting (98) yields:

S Yn (T Yn (z)) = T Zn (z -) + 1 z T Zn (z -) + (99)
such that:

1 S Yn (T Yn (z)) = (z -)T Zn (z -) T Zn (z -) + 1 + = 1 S Zn (T Zn (z -)) + (100)
Applying this to z = T -1 Yn (z) and applying [START_REF] Dong | Factorization methods in quantum mechanics[END_REF], we finally obtain the desired relation between the S-transforms of Z n and Y n :

1 S Yn (ω) = 1 S Zn ( ω(1 -S Yn (ω)) ) + ( 101
)
Applying the factorization of the S-transforms in the free product S Zn+1 = S Yn × S ξn , we obtain the main result of this section in the form of the following recursion for S Zn , the S-transform of the matrix evolving under Kesten recursion:

S ξ (ω) S Zn+1 (ω) = 1 S Zn ω(1 - S Z n+1 (ω) S ξ (ω) ) + (102)
where we recall that S ξ is the S-transform of the noise matrix that enters the recursion.

If we assume the existence of a stationary solution Z ∞ to the Kesten recursion, then the S-transform of Z ∞ must verify the following self-consistent relation:

S ξ (ω) S ∞ (ω) = 1 S ∞ ω(1 -S∞(ω) S ξ (ω) )
+ (103)

Continuum limit

We now consider again the continuum limit of the Kesten recursion which is obtained by taking the noise matrix as ξ n = (1 + m )I + σ √ B n , in the same way as in ( 23), and then taking the limit → 0. Here, we study the large N limit with a different method from the previous section II C 2, and we use the recursion relation for the S-transform which we derived above in Eq. ( 102).

We study the case where σ = √ N σ remains fixed as N → +∞, such that the spectrum of σB n has semi-circle density supported on [-2σ, 2σ], for both β = 1, 2 as defined in [START_REF] Allez | Eigenvector dynamics: General theory and some applications[END_REF]. As discussed in Sec II C 2, we need to assume m < 0. Note that we also assume 1 + m > 2σ

√ , i.e. < c = ( σ- √ σ2 -m m
) 2 so that ξ n is positive definite as required. It is shown in Appendix B that the S-transform corresponding to the resulting distribution for ξ is then (B27):

S ξ (ω) = -1 -m + (1 + m ) 2 + 4σ 2 ω 2σ 2 ω (104)
Let us define the following -expansion for S Zn and S Zn+1 :

S Zn (ω) = S (0) n (ω) + S (1) n (ω) + 2 S (2) n (ω) + • • • (105)
Inserting these expressions in the S-transform recursion relation (102) yields the following relations at zeroth and first orders in :

S (0) n+1 (ω) = S (0) n (ω) (106)

S

(1)

n+1 (ω) -S (1) n (ω) = -S (0) n (ω) m + σ2 ω + S (0) n (ω) + ωS (0) n (ω) (107)
As a consequence, the following evolution equation for S Zn holds at first order in :

S Zn+1 (ω) -S Zn (ω) = -S Zn (ω) m + σ2 ω + S Zn (ω) + ωS Zn (ω) (108)
In the continuous limit, we denote as in section II the continuous-time matrix process U t = Z t/ , defining time as t = n . The S-transform of U t is then a function of two variables ω and t:

∂ ∂t S(ω, t) = -S(ω, t) m + σ2 ω + S(ω, t) + ω ∂ ∂ω S(ω, t) (109)
The stationary solution is the solution of 0 = m + σ2 ω + S(ω) + ω ∂ ∂ω S(ω), namely:

S ∞ (ω) = -m - σ2 2 ω = |m| 1 - σ2 2|m| ω (110)
From the scaling property S aM (ω) = 1 a S M (ω) detailed in Appendix B, S ∞ is the transform of:

U ∞ law = 1 |m| W ( 111
)
where W follows an inverse-Wishart matrix probability distribution with parameter κ = |m| σ2 , see Appendix C. Indeed, for such a distribution, the S-transform is linear and equal to S

W (ω) = 1 -ω 2κ = 1 - σ2 2|m| ω [23]
. For completeness, we note that W is distributed as the inverse of a Wishart matrix with parameter q given by 1 q = 2κ + 1 = 1 + 2|m| σ2 , up to a scaling factor of 1 -q, see App. C for details. The eigenvalue density, or spectral distribution, of U ∞ has thus a finite support [λ + , λ -] and is given by:

ρ U∞ (λ) = |m| πλ 2 σ2 (λ + -λ) (λ -λ -), λ ∈ [λ -, λ + ] , with λ ± = 1 |m| + σ2 m 2 1 ± 2 |m| σ2 + 1 (112)
where we recall that |m| = -m. One can check that this agrees, upon the change of variables λ = e x , with the density (69). This equivalence shows a perfect matching between the stationary solution in the continuous limit for the discrete-time/large-N setting and the one obtained at finite N in Eq. [START_REF] Kesten | The limit distribution of Sinai's random walk in random environment[END_REF] as N becomes large. Indeed, as N → ∞, all the information about the matrix distribution is contained in the eigenvalue density.

The matching with the previous situation can also be proved directly from the finite-N stationary distribution as follows: let us denote X a matrix sampled from the stationary distribution in Eq. [START_REF] Kesten | The limit distribution of Sinai's random walk in random environment[END_REF]. We see by a rescaling that X law = 2V /σ 2 where V is distributed proportionally to det(V ) β( m σ 2 -N +1)-2 exp -β tr V -1 /2 . The factor in the exponential is now correct to apply the convergence of the eigenvalue density of V to the inverse Marcenko-Pastur density for this inverse-Wishart matrix, see App. C. The parameter q = N/T , and thus the related parameter κ = ( 1 q -1)/2, are asymptotically, from (37) by injecting σ = σ √ N : 1 q = 1 + 2|m| σ2 and κ = |m| σ2 . We have then:

X law = 2 σ 2 V = 2q σ2 T V = 1 |m| 2|m| σ2 1 + 2|m| σ2 T V = 1 |m| (1 -q)T V law → 1 |m| W ( 113
)
where W is an inverse-Wishart matrix with parameter κ = |m| σ2 . The last equation expresses that the rescaled matrix (1 -q)T V converges in law to the inverse-Wishart distribution, since its eigenvalue density converges to the inverse Marcenko-Pastur density as shown in (C9). This is in perfect agreement with (111), where the continuous-time and large-size limits are taken in the other order.

It is interesting to study the convergence of (109) to the stationary solution. Let us define the deviation from stationarity as ∆S(ω, t) = S(ω, t) -S ∞ (ω). One can linearize the evolution equation (109) around the stationary solution. One finds that the solution of this linearized equation reads

∆S(ω, t) = e mt 1 + σ2 ω 2m (1 -e mt ) ∆S 0, ωe mt 1 + σ2 ω 2m (1 -e mt ) (114)
Since here m < 0 we see that ∆S(ω, t) vanishes for any ω at large time, such that the convergence occurs. At finite time t for large ω, we see that ∆S(ω, t) = O( 1ω ), such that ∆S is killed at finite time for large ω argument, regardless of the initial value. An exact solution of the full non-linear equation (109) can also be obtained using the hodograph transform, see Appendix E.

Back to the discrete recursion

Let us give a short discussion of the discrete case, e.g. finite .

a. Corrections to the stationary S-transform. Although it is difficult to solve the recursion relation (102) for finite , we can easily obtain the next orders in of the stationary S transform S ∞ (ω), which show how the stationary distribution deviates from the inverse Wishart ensemble. With the same choice for S ξ (ω) in ( 104) we obtain upon expanding (103)

S ∞ (ω) = -m - σ2 2 ω + m 2 + 7 4 mσ 2 ω + 3σ 4 ω 2 4 + 2 -m 3 - 31 8 m 2 σ2 ω - 53 12 mσ 4 ω 2 - 25 16 σ6 ω 3 + O( 3 ) (115)
One can also expand (102) in powers of ω at fixed for a general noise and one obtains

S ∞ (0) = S ξ (0) -1 , S ∞ (0) = S ξ (0)S ξ (0) S ξ (0) + , S ∞ (0) = S ξ (0) S ξ (0)(S ξ (0) + 1) 2 S ξ (0) + 2S ξ (0) 2 (S ξ (0) + 1) 2 (S ξ (0) 2 + S ξ (0) + 1) (116)
In the case of the noise considered here in (104) this leads to

S ∞ (0) = - m m + 1 , S ∞ (0) = - σ2 (m + 1) 3 (m + 2) , S ∞ (0) = 2σ 4 (m (3m + 10) + 9) (m + 1) 5 (m + 2) 2 (m (m + 3) + 3) (117)
b. Stationary moments. The moments of the spectrum, denoted ϕ k = lim

N →∞ 1 N N i=1
λ k i can be extracted from the S transform. For instance the first moment is obtained from S(0) = 1 ϕ1 . Denoting ϕ k,∞ the moments of Z ∞ and ϕ k,ξ the moments of ξ n the first equation in (116) gives the first moment

ϕ 1,∞ = ϕ 1,ξ 1 -ϕ 1,ξ = - m + 1 m (118)
where the second equation is specified to the model (104). Similarly the second moment reads

ϕ 2,∞ = 2 ϕ 2,ξ -ϕ 4 1,ξ (1 -ϕ 1,ξ ) 3 (1 + ϕ 1,ξ ) = m(m + 1) 2 - σ2 m +2 m 3 (119)
Higher moments can be obtained from the series expansion of the S-transform, see Eq. (B21). One can compare with the case N = 1. In that case the stationary moments can be computed iteratively, by simply averaging the moments of the recursion relation. The result coincides with (118) for the first moment (which can be understood from the fact that E[ξ n ] ∝ I). For the second moment however it reads ϕ 2,∞ =

2 ϕ 2,ξ (1+ϕ 1,ξ ) (1-ϕ 2,ξ )(1-ϕ 1,ξ )
which is markedly different from (119). See more details on the moment computations in the scalar and matrix cases in App. F. c. Edges and stationary density corrections. We turn to the edges of the spectrum, that can be retrieved from the S-transform as showed in App. B by solving the system (B25). At order in the expansion (115), the two edges of the spectrum are found at positions:

λ (1) ± = λ (0) ± + -1 ± σ 2 √ σ2 -2m + O( 2 ) ( 120
)
where the order-zero terms λ (0)

± are the edges of the scaled inverse-Wishart stationary distribution given in Eq. ( 112). We note that the correction terms are negative, such that both edges are pushed to the left under the order-correction.

Finally, the stationary density can be obtained from the S-transform (115) by solving for g and using the Stieltjes inversion formula or Sokhotski-Plemelj formula (B3). At first order in , one obtains:

ρ(λ) = ρ (0) (λ) + 6 + (11m -9σ 2 )λ + 2m(3m -σ2 )λ 2 + m 3 λ 3 2πσ 2 λ 3 |m| (λ -λ (0) -)(λ (0) + -λ) + O( 2 ) ( 121
)
where ρ (0) is the stationary distribution given in Eq. ( 112). One can check that the integrated correction vanishes as required by normalization of ρ. A plot of the density at zeroth and first order in is showed in Fig. 5. 121), as the blue and orange lines respectively. Parameters were fixed as m = -10, σ = 1, = 2.10 -3 in this plot. We see that the correction of the density at first order in induces a leftward translation consistent with the negative correction to both edges in Eq. ( 120). We also note that the expansion fails close to the left edge λ (0) -.

B. Stieltjes transform approach

A complementary way to approach the large-N problem is to study the evolution of the Stieltjes transform. We return here to the continuum limit. Let us denote g the fluctuating Stieltjes transform defined as follows

g(z, t) := 1 N N i=1 1 z -λ i (t) . ( 122
)
Note that we recover the expected Stieltjes transform defined earlier through g(z, t) = E[g(z, t)]. Note also that the Stieltjes transform is self-averaging such that the fluctuating transform g(z, t) converges almost surely to g(z, t) for N → ∞. We emphasize that we will treat the large-N limit by dropping negligible terms, without attempt at mathematical rigor. For a fixed argument z, the stochastic evolution of g is, from the perturbative evolution of eigenvalues [START_REF] Nadal | Nonintersecting Brownian interfaces and Wishart random matrices[END_REF]:

dg(z, t) = 1 N N i=1 1 (z -λ i ) 2 dλ i + 1 2N N i=1 2 (z -λ i ) 3 dλ i .dλ i (123) = - ∂ ∂z 1 N N i=1     1 z -λ i     1 + mλ i + σ 2 1 j N j =i λ i λ j λ i -λ j     + σ 2 λ 2 i (z -λ i ) 2     dt + 1 N N i=1 √ 2σλ i (z -λ i ) 2 dW i (124)
up to O(dt 3/2 ). With the scaling of σ as σ √ N , the noise term is negligible for large N , such that:

dg(z, t) = - ∂ ∂z   g + m(zg -1) + σ2 N 2 i   1 2 j =i λ i λ j λ i -λ j 1 z -λ i - 1 z -λ j + λ 2 i (z -λ i ) 2     + O( 1 √ N ) (125) = - ∂ ∂z g + mzg + σ2 2 (zg -1) 2 + σ2 2N 2 i λ 2 i (z -λ i ) 2 + O( 1 √ N ) ( 126
)
To go from the first to the second line, we have added and subtracted the term i = j in the double sum, and removed the constant term -m in the derivative. For large N replacing g by the averaged Stieltjes transform and neglecting the last term which is of order 1 N , the evolution of g(z, t) is given by [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF]:

∂ ∂t g(z, t) = ∂ ∂z -g + (σ 2 -m)zg - 1 2 σ2 (zg) 2 (127)
An exact solution of this non-linear partial differential equation can be obtained using the hodograph transform, see Appendix E. It is furthermore shown in this Appendix that (127) is equivalent to the evolution (109) of the S-transform, as expected but not directly obvious. We also note that the stationary Stieltjes transform g stat (z) is

g stat (z) = |m| g W (|m|z) (128)
where

g W (z) = 1 z 2 (1 + κ) z -κ -κ (z -z -)(z -z + ) (129)
is the Stieltjes transform of W following the inverse-Wishart law with coefficient κ = |m| σ2 , with z ± = 1 κ (1+κ± √ 2κ + 1). Note that the branch has been chosen so that g W (z) 1/z at large z. This shows once again as in Eq. ( 111) that the stationary distribution for the Kesten evolution is the law of W |m| , with W as above. Note that the evolution equation for the Stieltjes transform in Eq. ( 127) can be studied in the Hamilton-Jacobi framework recently introduced by Grela, Nowak and Tarnowsky in [START_REF] Grela | Eikonal formulation of large dynamical random matrix models[END_REF]. In the notations of [START_REF] Grela | Eikonal formulation of large dynamical random matrix models[END_REF], where the Stieltjes transform g is associated to a momentum variable p and the variable z is thought of as a position variable, the Hamiltonian related to our model is:

H HJ (p, z) = p -σ2 -m zp + σ2 2 z 2 p 2 (130)
Our model and its Hamiltonian are an example of application of the Hamilton-Jacobi framework in large dynamical random matrix models, which we add to the list of examples given by the authors.

C. Expected characteristic polynomial and Stieltjes transform for finite N

Here we study Π(z, t) = N i=1 (z -λ i ) the characteristic polynomial of the evolving matrix at time t. It is an interesting quantity because its average satisfies a closed equation for any finite N . This property has been also used to study the standard Dyson Brownian motion [START_REF] Blaizot | Universal shocks in random matrix theory[END_REF][START_REF] Blaizot | Universal shocks in the Wishart random-matrix ensemble[END_REF][START_REF] Blaizot | Universal shocks in the Wishart random-matrix ensemble. II. Nontrivial initial conditions[END_REF].

Since it is linear in the λ i 's its stochastic evolution is simply dΠ =i Π z-λi dλ i . Using [START_REF] Nadal | Nonintersecting Brownian interfaces and Wishart random matrices[END_REF] and the identities

∂ z Π = i 1 z-λi Π and ∂ 2 z Π = i =j 1
(z-λi)(z-λj ) Π, one finds after some manipulations the evolution of its expectation Π(z, t) := E[Π(z, t)] as:

∂ t Π(z, t) = N m -(N -1) σ 2 2 Π(z, t) + -1 + z(-m + (N -1)σ 2 ) ∂ z Π(z, t) - σ 2 2 z 2 ∂ 2 z Π(z, t) (131)
One can look for the stationary solution of this equation. By inspection we find that it is given by the degree N polynomial

lim t→+∞ Π(z, t) = a N z N L -1-2m σ 2 N ( 2 σ 2 z ), a N = N -2m σ 2 -1 N -1 (132)
It is interesting to note that the dynamics of Π(z, t) can be related to the same Morse potential as studied in section II. Indeed one can perform a change of function and write

Π(z = e x , t) = f (x) ψ(x, t) , f (x) = e (N -1 2 -m σ 2 )x+ 1 σ 2 e -x (133)
and one finds that ψ(x, t) evolves according to minus the Schrödinger Hamiltonian in the Morse potential, i.e.

∂ t ψ = σ 2 2 -∂ 2 x + V (x) -N ψ (134)
where V (x) is defined setting β = 2 in (47) and N = -

(2m+σ 2 ) 2 4σ 4
is the N + 1th energy level of the Morse potential. We recall that the energy level and the eigenfunctions are given in [START_REF] Rider | Matrix Dufresne identities[END_REF]. Because the initial condition for Π(z, t = 0) is a polynomial in z of degree N , the initial wavefunction ψ(x, t = 0) belongs to the subspace spanned by the N + 1 lowest levels of the Morse potential. Since the time evolution is reversed, it converges to the highest energy accessible which recovers (132).

It is interesting to note that the characteristic polynomial is linked to the fluctuating Stieltjes transform as:

g(z, t) = 1 N ∂ z ln Π(z, t) (135)
which upon averaging over the noise leads to

g(z, t) = 1 N ∂ z E[ln Π(z, t)] 1 N ∂ z ln Π(z, t) := g(z, t) (136)
where the approximation is valid in the large N limit where the expectation value and the log can be exchanged. We introduce the notation g for simplicity. The evolution of the expected characteristic polynomial in Eq. ( 131) gives, after some manipulations, an equation valid for any N for the evolution of g defined in (136):

∂ t g = ∂ z -1 + z(-m + (N -1)σ 2 ) g - σ 2 N 2 z 2 g 2 - σ 2 2 z 2 ∂ z g (137)
Taking N large such that g g with σ 2 = σ2 /N and neglecting terms of order 1/N we recover by a different method the evolution for the Stieltjes transform, which was obtained in the previous section in Eq. ( 127):

∂ ∂t g(z, t) = ∂ ∂z -g + (σ 2 -m)zg - 1 2 σ2 (zg) 2 (138)

IV. LINKS WITH RIDER-VALK Ó AND GRABSCH-TEXIER

Our contruction of Kesten matrices bears some similarities, but also some differences, with some matrix stochastic diffusions studied by Rider and Valkó on the one hand, and by Grabsch and Texier on the other. We discuss the precise relation between these models in the present section.

A. Rider-Valkó

Rider and Valkó define in [START_REF] Rider | Matrix Dufresne identities[END_REF] the following geometric N × N drifted matrix Brownian motion M :

dM t = 1 2 -µ M t dt + M t dH t , M 0 = I (Itô) (139)
where the elements of the Brownian matrix H t are N 2 independent real Brownian motions, in contrast with the Hermitian symmetry assumed in this work, with real (β = 1) or complex (β = 2) elements. Since M is not symmetric, the authors study the law of (M M T ) s . To this aim, they define the process

Q t = M -1 -t ∞ -t M s M T s ds M -T -t
for negative time t, and show that it obeys the Itô SDE:

dQ t = ((1 + tr Q t )I + (1 -2µ)Q t ) dt -dH t Q t -Q t dH T t (Itô) (140)
The stationary distribution of this SDE, i.e. the law of

Q 0 = ∞ 0 M s M T s ds, is proved to be the following inverse- Wishart distribution: γ -1 2µ (X) ∝ (det X) -µ-1-N -1 2 e -1 2 tr X -1 (141)
As mentioned in the introduction, their work is an extension to the matrix realm of the Dufresne identity which gives the distribution of the exponential Brownian functional ( 14) as an inverse-gamma law.

B. Grabsch-Texier

In [START_REF] Grabsch | Topological phase transitions in the 1D multichannel Dirac equation with random mass and a random matrix model[END_REF], Grabsch and Texier study topological phase transitions occurring in a multichannel random wire modelled by a random-mass Dirac equation. They introduce the Riccati matrix Z t which follows the SDE:

dZ t = -Z 2 t + k 2 I -µGZ t -µZ t G dt + σ β (Z t dB t + dB t Z t ) (Stratonovich) (142)
where both symmetry cases β = 1, 2 are considered, such that the matrix Z is symmetric for β = 1 and Hermitian for β = 2. B t is a matrix Brownian motion with correlations given by G. In the case where G = I, B t is identical to the one used in this work, see definition below Eq. ( 27), provided one sets σ β = β/2. We note here the time variable as t in order to respect the conventions of our paper, but the evolution parameter is in fact the position x along the longitudinal direction of the disordered wire. Notice that this stochastic equation is given with the Stratonovitch prescription, contrary to the Itô prescription chosen for all SDEs in this work and in [START_REF] Rider | Matrix Dufresne identities[END_REF]. We recall that the Itô and Stratonovitch prescriptions are two different ways to interpret a SDE with multiplicative noise such as Z t dB t . Z t and dB t are taken to be statistically independent under Itô, whereas they are not independent under Stratonovitch, see [START_REF] Van Kampen | Itô versus Stratonovich[END_REF] for a pedagogical introduction.

Using the Fokker-Planck operator acting on the entries of the matrix Z, Grabsch and Texier demonstrate that the matrix diffusion for Z admits the following stationary distribution (for σ β = β/2):

f (Z) ∝ (det Z) -µ-1-β N -1 2 e -1 2 tr{G -1 (Z+k 2 Z -1 )} (143)
Let us take the special case G = I. After the change of variables {Z, k 2 } → { Z, } the SDE becomes in the → 0 limit:

dZ t = (I -2µZ t ) dt + σ β (Z t dB t + dB t Z t ) (Stratonovich) (144)
which admits as a consequence the following stationary distribution, restoring the parameter σ β for future comparisons:

f (Z) ∝ (det Z) -β 2σ 2 β µ-1-β N -1 2 e -β 4σ 2 β tr Z -1 (145)
In a recent work [START_REF] Grabsch | Wigner-Smith matrix, exponential functional of the matrix Brownian motion and matrix Dufresne identity[END_REF], Grabsch and Texier encountered the SDE (144) in a different problem: the evolution of the (symmetrized) Wigner-Smith time delay matrix Q of a multi-channel disordered wire, where the size L of the disordered region plays the role of the time parameter in the SDE. More precisely, the scattering process is defined by a Schrödinger equation with random potential coupling the N channels. Their work generalizes a problem which is known to coincide in the case N = 1 with the Brownian exponential functional that we presented in equation ( 14). In the general N case, they show that the symmetrized Wigner-Smith matrix Q can be written as:

Qt = X t t 0 dsX -1 s X † s -1 X † t (146)
where X t follows the following Stratonovitch flow:

dX t = (-µ + σ β dB t )X t (Stratonovich) (147)
such that Qt satisfies the SDE (144) where Z t is replaced by Qt . They rewrite (146) by defining Λ † s = X t X -1 t-s , which satisfies dλ s = (-µ + σ β d Bs )Λ s in Stratonovitch prescription, where d Bs = dB t-s . Using a reordering s → t -s in (146), similar to the one explained for the scalar case in Footnote 1, they obtain:

Qt law = t 0 Λ † s Λ s ds (148)
As a result they show that both Qt and the r.h.s. of Eq. ( 148) are distributed as (145) in the limit of large time.

As pointed out by Grabsch and Texier, this is very close to the results of Rider and Valkó in the real case (β = 1) as can be seen by comparing the stationary distributions (141) and (145) with σ β = β/2. Fixing β = 1, it is worth noting that the two SDEs (140) and (144) differ by a drift term (Q + I tr Q)dt which can be shown to be exactly the drift obtained by transforming a Stratonovich SDE with noise term ZdB + dBZ into an Itô SDE. The only difference is the absence of structure in the noise matrix dH t used by Rider and Valkó. Grabsch and Texier show in [START_REF] Grabsch | Wigner-Smith matrix, exponential functional of the matrix Brownian motion and matrix Dufresne identity[END_REF] that this is transparent at the level of the distribution. Indeed, the non-symmetric structure can be dealt with by gauging out the non-Hermitian part of dH t in a way that leaves the law of ∞ 0 M s M T s ds unaltered. Since the expression of Q as a matrix Brownian functional and the stationary distribution for SDE (144) were both established also in the complex Hermitian case (β = 2) by Grabsch and Texier, their work is in fact an extension of the matrix Dufresne identity to the complex case.

C. Connections with the present work

The works by Rider-Valkó and Grabsch-Texier are closely interconnected since they study the same matrix diffusion, as we have seen from the SDEs or from the Brownian functional expression. This latter point of view is naturally related to the Kesten recursion as was shown for the scalar case in I B. However, instead of defining the matrix Kesten problem from the Brownian functional, we have chosen in this work to study the recursion [START_REF] Forrester | Log-Gases and Random Matrices[END_REF]. We showed that, in the continuous limit, this recursion yields the SDE [START_REF] Lacroix-A-Chez-Toine | Non-interacting fermions in hard-edge potentials[END_REF]:

dU t = (I + mU t ) dt + σ U t dB t U t (Itô) (149)
This is close to the diffusion of Rider-Valkó and Grabsch-Texier, as the drift terms are essentially the same, but there is a crucial difference in the noise term, which is √ U dB √ U instead of U dB + dBU . Remarkably, although these two matrix processes are different, the eigenvalue processes are actually identical provided one compares Eq. ( 149) in Itô prescription and Eq. ( 144) in Stratonovitch prescription, with some correspondence between the parameters. This is showed in Appendix A where we study the joint evolution of the eigenvalues in the Grabsch-Texier model (which was not studied in [START_REF] Grabsch | Topological phase transitions in the 1D multichannel Dirac equation with random mass and a random matrix model[END_REF][START_REF] Grabsch | Wigner-Smith matrix, exponential functional of the matrix Brownian motion and matrix Dufresne identity[END_REF]). We find that we can identify the eigenvalues of U t and those of Z t if we set

σ β = σ 2 and µ = -m 2 + σ 2 4 (N + 2 -β).
Inserting these values in the stationary distribution (145) reproduces exactly our result [START_REF] Kesten | The limit distribution of Sinai's random walk in random environment[END_REF]

for β = 1, 2.
The evolution processes of the eigenvectors, however, have no obvious reason to be related. Since the stationary matrix measure is isotropic in both cases, the two models converge to the identical inverse-Wishart distribution.

We see that there is a subtle interplay between the form of the noise and the prescription used, such that terms of the Stratonovitch-Itô drift combine with the perturbative eigenvalue drift to recover the same eigenvalue flow.

As the last note on relations to other works, let us finally mention that our Fokker-Planck equation ( 29) appears in a very recent work of Ossipov [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] related to the Wigner-Smith time-delay matrix evolution in a multichannel disordered wire. Indeed, see Eq. ( 8) of this work, where the interaction term τ 2 i τi-τ k is equal to the interaction term of Eq. ( 29) plus a term proportional to τ i . In our notations, the work in Ref. [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] is restricted to the case m = 0 and σ = 1/ √ 2N , with null initial condition. In this special case and for large N , the author obtains the solution of the Stieltjes transform evolution equation, by a mapping to a Burgers equation. We extend this result by solving Eq. ( 127) with arbitrary parameters m, σ and arbitrary initial condition in Appendix E. We stress that the finite N evolution equation solved in [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] lacks a noise term, such that it does not qualify as a solution of the finite N diffusion problem.

D. Stochastic generalization : the matrix Bougerol identity

Bougerol's identity is a probabilistic result related to Dufresne's identity, which gives the law of the same Brownian exponential functional of Eq. ( 14) where one replaces the integral by a stochastic integral with respect to an independent Brownian motion w ν t (with drift ν). It states that the law of

V ∞ law = ∞ 0 e γt+Wt dw (ν) t (150)
is characterized by the density

P (V ∞ ) ∝ e 2ν arctan(V∞) 1 + V ∞ 2 1 2 -γ . (151)
Note that the heavy tail at V ∞ has the same power-law exponent as U ∞ in the Dufresne case, see Eq. ( 18). This result has recently been extended to the matrix case by Assiotis in [START_REF] Assiotis | A matrix Bougerol identity and the Hua-Pickrell measures[END_REF], where it was showed that the matrix analogue of (150) is distributed according to a Hua-Pickrell measure, the natural matrix extension of the generalized Cauchy measure of the scalar case (151). The authors of the present paper believe that a generalized Kesten recursion

Z n+1 = ν I + √ C n + Z n ξ n ν I + √ C n + Z n (152)
where C n is a standard GOE or GUE, would allow to explore this stochastic generalization. We conjecture that the Hua-Pickrell measure is the stationary distribution of the generalized Kesten recursion in the limit → 0 and leave the study for a future work.

V. CONCLUSION

In this paper we have proposed and studied an extension of the celebrated Kesten recursion and random variable 1 + z 1 + z 1 z 2 + . . . to N × N positive defined real (β = 1) and Hermitian (β = 2) matrices, the scalar case being recovered for N = 1. By studying the evolution of the matrix eigenvalues we have shown that the salient feature of the Kesten variable, i.e. the spontaneous generation of a heavy tail distribution, persists in the case of matrices for any N , while at large N the effects of the correlations and of the spectral rigidity, familiar in random matrix theory, also become important.

We have studied in details the continuum limit of the Kesten recursion, which is related to the exponential functional of the Brownian motion. In the matrix case, we have derived the Fokker-Planck equation which describes the flow of the eigenvalues and obtained the stationary distribution limit. It shows that for the continuum matrix Kesten recursion the stationary matrix distribution is given by the inverse Wishart ensemble for any N , which generalizes the classical results of Dufresne and of two of the present authors for N = 1. By relating the Fokker-Planck equation to an imaginary time Schrödinger problem, we have shown that the evolution of the eigenvalues, under a simple transformation, maps onto the quantum evolution of N fermions in a Morse potential. Although both cases are presumably integrable, in the Hermitian case β = 2 the fermions are non-interacting, leading to a simple solution for the dynamics of the original eigenvalues. We showed that the above mentioned stationary state for the Kesten recursion can be retrieved from the ground state of the fermion system.

Since it is an interesting problem in its own right, we have explored in details the physics of the ground state of N non-interacting fermions in the Morse potential. In principle this system can be realized in cold atom experiments, where the fermions are trapped in potentials of tunable shape. Although the Morse potential has only a finite number of bound states, and can accommodate only a finite number of fermions, one can tune its parameters so that this number is large. In the large N limit the Fermi gas develops two edges and we have computed the scaling forms of the density and correlations in the bulk and near the edges. These can be related, via some (a priori non obvious) transformation, to the Marcenko-Pastur distribution and to the correlations in the Wishart-Laguerre ensemble. If we specify the parameters of the Morse potential to be relevant for the matrix Kesten problem, we find two distinct large N scaling limits, depending on whether the variance of the noise is scaled as 1/N or remains O(1). In the first case (weak noise) the two spectral edges are "soft" and the correlations at these edges are described by the Airy kernel. In the second case (strong noise) the right edge is pushed to infinity. The largest eigenvalues (which correspond to the fermions at the far right) are described by the Bessel determinantal point process, i.e. the physics of the "hard edge" in random matrix theory. Remarkably, this provides a precise description of the many body correlations of the heavy tail of the matrix Kesten recursion. Note that these results for non-interacting fermions were obtained from the case β = 2. In the case β = 1, the fermions have in addition a Sutherland type interaction, and we have obtained the explicit form of their ground-state wavefunction.

Finally, we have shown that the matrix Kesten recursion proposed in the present work, although distinct from those studied by Rider and Valko and by Grabsch and Texier, share the same distribution of eigenvalues in its continuum limit (modulo some subtleties with Itô and Stratonovich prescriptions that we elucidated). It is thus rewarding that a universal behavior seems to emerge in that problem.

The present work leaves open several questions and suggests some interesting extensions. First, concerning the discrete matrix recursion we have been able to derive exact results only in the large N limit, using free probability theory. A self-consistent equation for the S-transform of the Kesten matrix was obtained in the stationary limit. We have not been able to solve this equation, except in the case of the continuum limit. It would thus be nice to study this equation further, and also to extend the analysis to finite N . In particular, a problem which remains open is whether the remarkably simple condition E(z ν ) = 1, which determines the heavy tail exponent ν for the scalar case N = 1, has any analog in the matrix case.

Let us finish by recalling our initial motivation for studying the present model, which is to explore possible matrix (non-commuting) generalizations of the famous directed polymer problem (which is related to the KPZ stochastic growth equation). The directed polymer problem amounts to study sums over a set of paths of products of scalar weights along these paths. The idea is to replace these weights by random matrices. Such objects already appeared, e.g. in Chalker-Coddington type of models for localization in d = 2 [START_REF] Cardy | Quantum network models and classical localization problems[END_REF]. In the case of positive definite matrix weights however, it has only been studied very recently [START_REF] O'connell | Interacting diffusions on positive definite matrices[END_REF]. Note also an operator generalization, which also recently appeared in the context of stochastic quantum spin chains [START_REF] Jin | From stochastic spin chains to quantum Kardar-Parisi-Zhang dynamics[END_REF]. We hope that the present study will also stimulate progress in these directions.

Symmetric case β = 1

In the symmetric case, the noise matrix is dB t =

dHt+dH T t √ 2
. Let us truncate the matrix evolution of the model studied in this paper and keep only the noise term:

dU t = σ U t dB t U t = σ U t dH t + dH T t √ 2 U t (A4)
The evolution of eigenvalue λ i , with corresponding eigenvectors {|v i,t }, can be expressed perturbatively [START_REF] Tao | Topics in Random Matrix Theory[END_REF] as:

λ i,t+dt = λ i,t + σ v i,t | U t dH t + dH T t √ 2 U t |v i,t + 1 j N j =i v j,t | σ √ U t ( dHt+dH T t √ 2 ) √ U t |v i,t 2 λ i,t -λ j,t (A5)
λ i,t+dt = λ i,t + σλ i,t v i,t | dH t + dH T t √ 2 |v i,t + σ 2 1 j N j =i λ i,t λ j,t v j,t | dHt+dH T t √ 2 |v i,t 2 λ i,t -λ j,t (A6)
From the independence relations of the A i,j variables, let us define independent η i,j variables following a N (0, 1) law:

v i,t | dH t + dH T t √ 2 |v j,t = √ dt( A i,j + A j,i √ 2 ) = √ dt √ 2 η i,i if i = j η i,j if i = j (A7)
The eigenvalue λ i then evolves according to:

λ i,t+dt = λ i,t + σλ i,t √ 2dt η i,i + σ 2 1 j N j =i λ i,t λ j,t η 2 i,j λ i,t -λ j,t dt (A8) dλ i,t = σλ i,t √ 2dt η i,i + σ 2 1 j N j =i λ i,t λ j,t λ i,t -λ j,t dt + σ 2 1 j N j =i λ i,t λ j,t (η 2 i,j -1) λ i,t -λ j,t dt (A9)
Following [START_REF] Tao | Topics in Random Matrix Theory[END_REF], the last term can be discarded in the stochastic equation as it has zero mean and order dt 3 in the third moment. Finally, defining N independents Brownian motions (B i ) 1 i N , the evolution of the eigenvalues of U is:

dλ i = √ 2σλ i dW i + σ 2 1 j N j =i λ i λ j λ i -λ j dt (Itô) (A10)
The moment generating function, or T -transform, is then:

T (z) = z g(z) -1 = b a u ρ(u) z -u du = ∞ k=1 ϕ k z k (B7)
and has the following scaling property:

T aM (z) = T M ( z a ) (B8)
The B (or Blue) transform is defined as the functional inverse of the Stieltjes transform B(g(z)) = z for all z in C\[a, b], and the R-transform is defined as:

R(ω) = B(ω) - 1 ω (B9)
From the scaling property of g, we obtain scaling properties for B and R:

B aM (ω) = aB M (aω) (B10) R aM (ω) = aR M (aω) (B11)
It can be shown that the R-transform admits a Taylor expansion as ω → 0:

R(ω) = ∞ k=1 κ k ω k-1 (B12)
where the coefficients κ k (M ) are the free cumulants of M , which can be expressed in terms of the moments ϕ k in a different way from cumulants in standard probability theory. The definition of free cumulants and their combinatorial interpretation in terms of non-crossing partitions can be found in [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Mingo | Free Probability and Random Matrices[END_REF][START_REF] Speicher | Free probability theory and non-crossing partitions[END_REF]. An important result from free probability theory is the additivity of free cumulants for the sum of two free matrices A and B, i.e. for all k:

κ k (A + B) = κ k (A) + κ k (B) (B13)
Additivity then extends to the R-transfoms of free matrices A and B:

R A+B (ω) = R A (ω) + R B (ω) (B14)
This is usually referred to as Voiculescu's free addition. A similar result holds for free multiplication, that we now turn to. In this regard, we define the S-transform as:

S(ω) = ω + 1 ωT -1 (ω) (B15)
where T -1 is the functional inverse of the T -transform. From the scaling property of T , we obtain:

S aM (ω) = 1 a S M (ω) (B16)
We state the following relations between S and R transforms:

S(ω) = 1 R(ωS(ω)) (B17) R(ω) = 1 S(ωR(ω)) (B18)
The main free probability result that we use in this paper is the following: for two free self-adjoint matrices A and B, the self-adjoint product M = A 1/2 BA 1/2 verifies the following factorization property on the S transforms [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Voiculescu | Limit laws for random matrices and free products[END_REF]:

S M (ω) = S A (ω) × S B (ω) (B19)
This is the free multiplication, or free multiplicative convolution, theorem, expressed here for the self-adjoint product of free random matrices.

Note that the S-transform can be expressed as a function of the moments ϕ k as

S(ω) = ω + 1 ω ∞ k=1 ϕ k ω k <-1> (B20)
where subscript < -1 > denotes the operation of taking the inverse with respect to composition of formal power series. To invert that relation and obtain the moments from the expansion of S(ω) in powers of ω one can write, denoting z = T -1 (ω):

1 z = ωS(ω) ω + 1 = series in ω (B21)
One then inverts this series, and obtains ω as a series in 1/z whose coefficients are the moments, i.e. ω = k 1

ϕ k z k = T (z).
It is interesting to note that the knowledge of the S-transform allows to find the position of the edges λ e of the spectrum. In generic cases these edges are determined by the roots g e of the equation B (g e ) = 0 , λ e = B(g e ) (B22)

From their definitions, one sees that the relation between the B and S transforms is obtained by eliminating ω in the following system:

B(g) = ω + 1 ωS(ω) (B23) ω = gB(g) -1 (B24)
As a consequence we find that the position of the edges are the solutions of the following system: S(ω e ) + ω e (ω e + 1)S (ω e ) = 0

λ e = ωe+1 ωe S(ωe) (B25)
One can check for instance that for S(ω) = 1 -ω 2κ one recovers as solutions the two edges of the inverse-Wishart ensemble given in (C9). At order in the expansion (115), one recovers the edges given in (120).

(z) = z 2 1 -1 -4 z 2
and its R-transform is: R B √ N

(ω) = ω. By free sum of a scaled identity and a scaled Wigner matrix, the R-transform of ξ is:

R ξ (ω) = 1 + m + σ2 ω (B26)
We deduce from (B17) the S-transform of the noise matrix ξ:

S ξ (ω) = -1 -m + (1 + m ) 2 + 4σ 2 ω 2σ 2 ω (B27)
where the correct square-root branch is chosen here such that S has a regular behaviour as ω → 0. In this region, we have then

S(ω) = 1 1+ m - σ 2 (1+ m) 3 ω + O(ω 2 ).
C -1 if β = 1 (such that X † = X T ), or complex centered gaussian law of variance C -1 if β = 2. By construction, M is a symmetric (resp. Hermitian) positive semidefinite matrix with real (resp. complex) entries if β = 1 (resp. β = 2). The matrix distribution of X is P (X) ∝ e -C β 2 tr(X † X) , and taking into account the Jacobian of the transformation X → M we obtain the Wishart matrix distribution of M :

P w (M ) = Q T,N,C,β det(M ) β 2 (T -N +1)-1 e -C β 2 tr M (C1)
where the normalization constant Q T,N,C,β is:

Q T,N,C,β = Cβ 2 β 2 N T π -β 2 N (N -1) 2 N j=1 Γ β 2 (T -N + j) = (Cβ/2) β 2 N T Γ N,β ( β 2 T ) (C2)
where the generalized multivariate Gamma function is Γ

N,β (a) = π βN (N -1) 4 N j=1 Γ a -β j-1 2
. The joint eigenvalue distribution of M on (R + ) N follows:

P w ( λ) = K T,N,C,β ∆( λ) β N i=1 λ β 2 (T -N +1)-1 i e -C β 2 N i=1 λi (C3)
where the normalization constant K T,N,C,β is:

K T,N,C,β = Cβ 2 β 2 N T N j=1 Γ 1 + β 2 Γ 1 + β 2 j Γ β 2 (T -N + j) = 1 N ! Cβ 2 β 2 N T Γ β 2 N π β 2 N (N -1) Γ N,β β 2 N Γ N,β β 2 T (C4)
In the limit of large N and T with N T = q fixed, and with the usual convention C = 1, the Empirical Spectral Distribution (ESD) of the rescaled matrix 1 T M converges to the Marcenko-Pastur distribution:

ρ 1 T M (λ) = 1 N N i=1 δ(λ -λ i ) → ρ MP (λ) = 4λq -(λ + q -1) 2 2qπλ , λ ∈ (1 - √ q) 2 , (1 + √ q) 2 (C5)
This limit distribution is characterized by the following Stieltjes, R and S transforms:

g MP (z) = (z + q -1) -z -(1 - √ q) 2 z -(1 + √ q) 2 2qz R MP (ω) = 1 1 -qω S MP (ω) = 1 1 + qω (C6)

The inverse-Wishart ensemble

The inverse-Wishart distribution is the distribution of M = W -1 , where W follows the Wishart distribution, with the parameters defined above. Taking into account the Jacobian of the M → M -1 transform, we obtain the inverse-Wishart matrix probability density, with the same normalization constant Q T,N,C,β :

P iw (M ) = Q T,N,C,β det(M ) -β 2 (T +N -1)-1 e -C β 2 tr M -1 (C7)
The corresponding joint eigenvalue distribution is, with the same normalization constant K T,N,C,β :

P iw ( λ) = K T,N,C,β ∆( λ) β N i=1 λ -β 2 (T +N -1)-1 i e -C β 2 N i=1 1 λ i (C8)
In the limit of large N and T with N T = q fixed, and with the usual convention C = 1, the ESD of the rescaled matrix (1 -q)T × M then converges to the inverse Marcenko-Pastur distribution:

ρ (1-q)T ×M (λ) = 1 N N i=1 δ(λ -λ i ) -→ ρ IMP (λ) = κ πλ 2 (λ + -λ) (λ -λ -), λ ∈ [λ -, λ + ] (C9)
Replacing now u = u(ω, S) = σ2 ω 2 + 2ω(S + m) we finally obtain the result announced above in (E2). The function φ is determined from the initial condition ∆S 0 (ω) = ∆S(ω, t = 0)

φ(u) = - 2 σ u + ( m σ ) 2 arctanh σΩ(u) + m σ u + ( m σ ) 2 (E10)
where the function Ω(u) is obtained by inverting (E3). Putting all together we find

t = 2 σ 2ω∆S(ω, t) + ( m σ ) 2 arctanh    σ(ω -Ω(2ω∆S(ω, t))) 2ω∆S(ω, t) + ( m σ ) 2 -(σω + m σ ) σΩ(2ω∆S(ω,t))+ m σ √ 2ω∆S(ω,t)+( m σ ) 2    (E11)
which is finally equivalent to the result given in (E2).

Stieltjes transform

The same method can be applied to the Stieltjes transform. Starting from Equation (127) for g(z, t), we write a partial differential equation on the function t(z, g) using ∂ g t = 1

∂tg and ∂ z t = -∂zg ∂tg . One obtains the linear equation

1 + (m -σ2 )z + σ2 z 2 g ∂ z t -(m -σ2 )g + σ2 zg 2 ∂ g t = 1 (E12)
One homogeneous solution is found as :

u(z, g) = g + (m -σ2 )zg + σ2 2 z 2 g 2 (E13) Inverting u(z, g) gives g(z, u) = -1 -(m -σ2 )z -z (1/z + (m -σ2 )) 2 + 2σ 2 u σ2 z 2 (E14)
where the branch is chosen in order to satisfy g(z) ∼ 1 z for z → ∞, with m < 0. Changing variables as t(z, g) → t(z, u) in (E12) we obtain 1 + (m -σ2 )z + σ2 z 2 g(z, u) ∂ z t = 1 (E15) thus :

∂ z t = -1 z (1/z + (m -σ2 )) 2 + 2σ 2 u (E16)
This integrates as :

t(z, u) = -1 2σ 2 u + (m -σ2 ) 2 arcsinh (2σ 2 u + (m -σ2 ) 2 )z + m -σ2 σ√ 2u + φ(u) (E17)
where φ is an arbitrary function of u that can be determined from the initial condition at time t = 0. Inserting the definition of u = u(z, g) from (E13) we obtain the exact solution for the Stieltjes transform g = g(z, t) in the implicit form by solving t = -1 where we denote A(u) = m -σ2 2 + 2σ 2 u and we recall that Z(u) is defined in (E19) from the initial condition g 0 . One can check that the second equation in (E21) gives z = Z(u) for t = 0, which is equivalent to U (z, 0) = u(z, g 0 (z)) which, plugged in (E20), gives back g(z, t) = g 0 (z). In the large time limit, the r.h.s. in the second equation diverges with corresponds to z going to ∞, hence we recover the stationary value U (z, t → ∞) = m -σ2 2 which in turn yields the stationary Stieltjes transform (128).

2σ 2 g + (σ 2 zg + (m -σ2 )) 2 arcsinh (2σ 2 g + (σ 2 zg + (m -σ2 )) 2 )z + m -σ2 2σ 2 (1 + (m -σ2 )z)g + σ4 z 2 g 2 + φ(u(z, g)) (
This solution reproduces the one obtained by Ossipov in Eq. ( 18) of [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] in the particular case m = 0, σ = 1/ √ 2 and r = t/2 plays the role of the time variable, with his z 0 = z 0 (z, t) related to our u = U (z, t) by 1/z 2 0 = u + 1 4 . Note that in [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] the solution is only given for the initial condition g 0 (z) = 1/z while we have given here the solution for a general initial condition. In the limit t → ∞ (equivalently r 1), a scaling form was found in [START_REF] Ossipov | Scattering Approach to Anderson Localization[END_REF] in the case m = 0. Note however that a well-defined stationary solution does not exist then, since m = 0 is the critical case.

Time evolution equivalence of S and Stieltjes transforms

Let us show that the evolution equations obtained for the S-transform (109) and for the Stieltjes transform (127) are indeed equivalent. For ease of notation, let us drop the dependence on t for all functions and denote f = T -1 the functional inverse of the T -transform of U t such that S(ω) = ω+1 ωf (ω) . Then, we obtain from (109):

- ω + 1 ωf 2 ∂f ∂t = ∂ ∂t S = - ω + 1 ωf (ω) m + σ2 ω + ω + 1 ωf (ω) + ω ωf (ω) -(ω + 1)(f (ω) + ωf (ω)) ω 2 f 2 (ω) (E22)
Simplifying, this leads to:

∂ ∂t f = 1 + f (ω) m + σ2 ω - (ω + 1)f (ω) f (ω) (E23)
Let us deduce the evolution of T = f -1 , between two consecutive instants t 2 = t 1 + dt. We know that : From the definition of T as T (z) = zg(z) -1, we finally obtain the evolution of g as in Eq. (127):

T -1 t2 (
∂ ∂t g = ∂ ∂z -g + σ2 -m zg - 1 2 σ2 z 2 g 2 (E27)
1. Discrete scalar case

The scalar Kesten recursion, with the regularization parameter , is Z n = ξ n ( + Z n-1 ) with the initial condition Z 0 = 0. Denoting ϕ k,n = E[Z k n ] the kth moment of the scalar random variable Z n and ϕ k,ξ the kth moment of ξ, we have the following recursion for the first moment:

ϕ 1,n = ϕ 1,ξ ( + ϕ 1,n-1 ) ⇒ ϕ 1,n = (ϕ 1,ξ + ϕ 2 1,ξ + • • • + ϕ n 1,ξ ) = ϕ 1,ξ ϕ n 1,ξ -1 ϕ 1,ξ -1 ----→ n→∞ ϕ 1,ξ 1 -ϕ 1,ξ (F1)
If we choose, as in Eq. ( 15) of the main text, ξ = 1 + m + √ 2σ √ X with X ∼ N (0, 1), we obtain

ϕ 1,n=+∞ = - m + 1 m (F2)
The limit n = +∞ is finite only for ϕ 1,ξ < 1, i.e. m < 0. The recursion for the second moment is ϕ 2,n = ϕ 2,ξ ( 2 + 2 ϕ 1,n-1 + ϕ 2,n-1 ) and we obtain the second moment for every step n as: The limit n = +∞ is finite only for ϕ 2,ξ < 1, i.e. 2σ 2 + m(2 + m ) < 0. In this work we are interested in random variables ξ n which are restricted to be positive. The Eqs (F1) and (F3) are still valid. The Gaussian choice for ξ n violates this positivity condition, but only with exponentially small probability (as exp(-1/ )) for small . Hence we expect that the explicit formula (F2) and (F4) are useful in that case within an expansion in powers of at small . The same applies below in the matrix case.

ϕ 2,n = 2 ϕ 2,ξ (ϕ 1,ξ -1) (ϕ 1,ξ + ϕ 2,ξ ) ϕ n 2,ξ -ϕ 1,ξ 2 (ϕ 2,ξ -1) ϕ n 1,ξ + ϕ 1,ξ -ϕ 2,ξ + 1 + ϕ 2,ξ (1 -ϕ 1,ξ ) (1 -ϕ 2,ξ ) (ϕ 2,ξ -ϕ 1,ξ ) ----→ n→∞ 2 ϕ 2,ξ (1 + ϕ 1,ξ ) (1 -ϕ 1,ξ )(1 -ϕ 2,ξ ) (F3) Injecting ξ = 1 + m + √ 2σ √ X,

Continuous scalar case

In the continuous limit, the scalar Kesten evolution is the Itô stochastic equation ( 16): dU t = (1 + m U t ) dt + √ 2σ U t dW t , and we choose here for simplicity U 0 = 0. The evolution of the moments ϕ k,t = E[U k t ] can be computed recursively from this stochastic evolution. The first moment is found as:

∂ t ϕ 1,t = (1 + mϕ 1,t ) ⇒ ϕ 1,t = - 1 -e mt m ---→ t→∞ - 1 m (F5)
where the large time coincides as expected with the → 0 limit of the discrete result (F2). The evolution of the squared process U 2 t is obtained from the Itô rules as dU 2 t = 2U t 1 + (m + σ 2 ) U t dt + 2 √ 2σ U 2 t dW t . Taking the expectation yields the evolution equation of ϕ 2,t :

∂ t ϕ 2,t = 2(ϕ 1,t + (m + σ 2 )ϕ 2,t ) ⇒ ϕ 2,t = m e 2t(m+σ 2 ) -2e mt + 1 + 2σ 2 (1 -e mt ) m (m + σ 2 ) (m + 2σ 2 ) ---→ t→∞ 1 m(m + σ 2 ) (F6)
where the large time limit coincides as expected with the → 0 limit of the discrete result (F4). Note that the stationary variable U ∞ exists for m < σ 2 , but that the first and second moments are finite at t = +∞ only for m < 0 and m < -σ 2 respectively. The above time dependent solutions however are valid in all cases, and show how the moments grow with time in the cases where their limits do not exist.

Discrete matrix case

The N × N matrix Kesten recursion Z n = I + Z n-1 ξ n I + Z n-1 makes it possible to obtain the moments ϕ ] and assume that the distribution of ξ is invariant under rotations, i.e. conjugation by an element of the orthogonal (resp. unitary) group for β = 1 (resp. 2). This assumption imposes that E[ξ] = ϕ 1,ξ I and we thus have a similar recursion for the first moment as in the scalar case:

E[Z n ] = ( I + E[Z n-1 )]ϕ 1,ξ ⇒ E[Z n ] = ϕ (N ) 1,n I with ϕ (N ) 1,n = ϕ 1,ξ ϕ n 1,ξ -1 ϕ 1,ξ -1 ----→ n→∞ ϕ 1,ξ 1 -ϕ 1,ξ (F7)
With ξ defined as in Eq. ( 23) we obtain ϕ 2,t=+∞ can be verified to agree with the first and second moments of the matrix inverse-Wishart distribution obtained as the stationary distribution of the continuous model in Eq. [START_REF] Kesten | The limit distribution of Sinai's random walk in random environment[END_REF].

FIG. 1 .

 1 FIG. 1. The Morse potential V (x) in (53), plotted with parameters x0 -2.565 and g = 13, which correspond to σ = 1, m = -3 and N = 10 in the Kesten model. The energy of the highest occupied state by the N = 10 fermions (i.e. the Fermi energy F given in Eq. (63)) is indicated by the dashed line. For these parameters the Morse potential contains N b = 13 bound states, it is thus almost full.

  a. Scaling σ ∼ 1/ √ N . In the case where σ = σ/ √ N with σ kept constant, one has g = N (1 -m σ2

FIG. 2 .

 2 FIG. 2. In the regime where σ = σ/ √ N , we plot the relative positions of the two edges x± -x0 (left) and the derivative of the potential at the edges V (x±)/N 2 (right) as a function of the positive variable -m σ2 . Note that x+ -x0 diverges when -m σ2

FIG. 3 . 2 √ 10 ,

 3210 FIG. 3. Scaling σ ∼ 1/√ N : Plot of the 1-particle density P1(x) given in[START_REF] Gautié | Non-crossing Brownian paths and Dyson Brownian motion under a moving boundary[END_REF] for N = 10 particles with σ = 2, such that σ = 2 √ 10 , and m = -3 (solid blue line). For these parameters, we have x--2.59, x0 -1.95 and x+ 0.39. The bulk density P bulk (x) given in (69) obtained from the LDA is the dashed gray line. The left and right edge densities P edge ± (x) obtained from the Airy kernel in[START_REF] Voiculescu | Symmetries of some reduced free product C*-algebras[END_REF] are the dashed red lines. The inset shows the region near the right edge for N = 50, with σ = 2, m = -3 as above. We see that P edge + (x) is closer to the exact density for this higher number of particles.

FIG. 4 .

 4 FIG.4. Constant σ : Plot of the 1-particle density P1(x) given in[START_REF] Gautié | Non-crossing Brownian paths and Dyson Brownian motion under a moving boundary[END_REF] for N = 10 particles with σ = 2, m = -3 (solid blue line). For these parameters, we have x--4.45, x0 -3.76 and x+ 2.26. The bulk density P bulk (x) given in (69) obtained from the LDA is the dashed gray line. The density at the left edge P left edge (x) (red dashed line) is given by the Airy kernel, see Eq. (79). The density P far right (x) on the far right near x+ ∼ log N (orange dashed line) is obtained from the Bessel kernel, and given in Eq. (84). The inset shows the region near the right edge for N = 50, with σ = 2, m = -3 as above. For these parameters, we have x+ 3.81. We see that the exact density is very closely approximated by P far right (x) in that region.

FIG. 5 .

 5 FIG.5. This figure shows the plot of the stationary density at orders 0 and 1 in , see Eq. (121), as the blue and orange lines respectively. Parameters were fixed as m = -10, σ = 1, = 2.10 -3 in this plot. We see that the correction of the density at first order in induces a leftward translation consistent with the negative correction to both edges in Eq. (120). We also note that the expansion fails close to the left edge λ

2 .

 2 S-transform of ξ We obtain the S-transform of the noise matrix ξ = (1 + m )I + σ √ B √ N considered in the text, where we assume 1 + m > 0. Here B √ N is a symmetric Wigner matrix with semi-circle eigenvalue density with support [-2, 2]. Its Stieltjes transform is g

2 2σ 2 u

 22 E18)It is possible to rewrite it using the initial condition g(z, t = 0) = g 0 (z). Let us define the function Z(u) such thatu(z, g 0 (z)) = g 0 (z) + (m -σ2 )zg 0 (z) + σ2 2 z 2 g 0 (z) 2 = u ⇔ z = Z(u) (E19)Then the solution isg(z, t) = -1 -(m -σ2 )z -z (1/z + (m -σ2 )) 2 + 2σ 2 U (z, t) σ2 z 2 (E20)where U (z, t) is obtained by invertingU (z, t) = u ⇔ u)Z(u) + m -σ2 ) sinh t A(u) + A(u)Z(u) + m -σ2 cosh t A(u) -m + σ2

  ω) -T -1 t1 (ω) = dt 1 + f (ω) m + σ2 ω -(ω + 1)f (ω) f (ω) (E24) Evaluating at ω = T t2 (z) = T t1 (z) + dt ∂T ∂t : z -z + dt ∂T ∂t × T -1 (T (z)) = dt 1 + z m + σ2 T (z) -(T (z) + 1)(T -1 ) (T (z)) z (E25) Using (T -1 ) (T (z)) = 1 T (z), we obtain:∂ ∂t T = T (z) + 1 z -T (z) 1 + z(m + σ2 T (z)) (E26)

  we obtainϕ 2,n=+∞ = (2 + m )(1 + (2σ 2 + m(2 + m ))) m(2σ 2 + m(2 + m )) (F4)

= 1 N 1 N

 11 Tr E[Z k n ] from the distribution of ξ and the initial condition that we fix to Z 0 = 0. We denote ϕ (N ) k,ξ = Tr E[ξ k

1 N 1 N 2 E 1 N 2 E

 11212 coincides with the large-N result(118). Let us turn to the second moment. The isotropy condition on the ξ distribution imposes that there exist four coefficients a, b, c, d such that for any symmetric (resp. Hermitian if β = 2) matrix M [102]:E[ξM ξ] = aM + b N Tr(M )I , E[(Tr M ξ) 2 ] = N c Tr(M 2 ) + d(Tr M ) 2 (F8) Note that ϕ 2,ξ = a + b and E[( Tr ξ) 2 ] = c + d.These two sets of parameters are not independent however, indeed choosing the matrix M to be M ij = δ i1 δ j1 , one obtains the relationEξ 2 11 = a + b N = N c + d (F9)We can then write the expectation of the matrix Z 2 n as:E[Z 2 n ]= E[ I + Z n-1 ξ n ( I + Z n-1 ) ξ n I + Z n-1 ] (F10) = aE ( I + Z n-1 ) 2 + b N E [( I + Z n-1 ) Tr( I + Z n-1 )] (F11)Taking the trace and dividing by N yields:(Tr Z n-1 ) 2 (F12)The correlations of ξ in Eq. (F8) give furthermore:E (Tr Z n ) 2 = E (Tr(ξ( I + Z n-1 )) 2 = N c E Tr(( I + Z n-1 ) 2 ) + d E (Tr(( I + Z n-1 )) 2 (Tr Z n-1 ) 2 ) (F14)Equations (F12) and (F13) form a coupled recursion system which is fully solvable. For simplicity, we simply give the stationary values ϕ(N ) 2,n=+∞ and E (Tr Z ∞ ) 2 : ϕ (N ) 2,n=+∞ = 2 (ϕ 1,ξ + 1)(a + b + bc -ad) (ϕ 1,ξ -1)(a + d + bc -ad -1) (F15) E (Tr Z ∞ ) 2 = N 2 2 (ϕ 1,ξ + 1)(c + d + bc -ad) (ϕ 1,ξ -1)(a + d + bc -ad -1) (F16)We check that fixing N = 1, a = d = ϕ 2,ξ and b = c = 0 recovers the scalar result (F3). In the large-N limit and for an arbitrary matrix M , the random variable 1 N Tr(M ξ) is self-averaging. The vanishing of its variance yieldsd = ϕ 2 1,ξ and c = o( 1 N ). Recalling that ϕ 2,ξ = a + b such that b must remain O(1), relation (F9) then yields that a = ϕ 2 1,ξ in this limit. With these substitutions, we verify that ϕ (N ) 2,n=+∞ of Eq. (F15) matches the general formula for the second moment obtained for large N from free probability tools in Eq. (119).With ξ defined as in Eq. (23) we obtain a= (1 + m ) 2 + σ 2 δ β=1 , b = N σ 2 , c = σ 2 (1 + δ β=1 )/N, d = (1 + m ) 2 such that: ϕ (N ) 2,n=+∞ = 2 + 2 m + (2 + m )(m(2 + m ) -N σ 2 ) m (m(m + 2) -σ 2) (σ 2 (1 + δ β=1 ) + m(m + 2)) (F17) E (Tr Z ∞ ) 2 = N (m + 2) -σ 2 (1 + δ β=1 ) + N -σ 4 (1 + δ β=1 ) + (mσ + σ) 2 δ β=1 + m(m + 1) 2 (m + 2) m (m(m + 2) -σ 2 ) (σ 2 (1 + δ β=1 ) + m(m + 2))

Note that Ut is the continuum analogue of Eq. (2) and can be written as t 0 e A t -As ds which has the same law (at fixed t) as t 0 e As ds in Eq. (12) by a reordering of variables analogous in the continuum to the discrete reordering employed from Eq. (2) to (3). However, note that these processes are not equivalent as can be seen by differing stochastic evolutions. In particular, the evolution of Ut is closed in Ut while it is not the case for the quantity on the r.h.s. of Eq. (12).

ACKNOWLEDGMENTS

We thank A. Krajenbrink, M. Potters, G. Baverez and A. Lafay for useful discussions as well as S. N. Majumdar and C. Texier for pointing out connected references. We are grateful to a referee for bringing the Bougerol identity and its matrix generalization to our attention. P.L.D. thanks the LPTMS Orsay for hospitality while this work was completed and acknowledges support from ANR grant ANR-17-CE30-0027-01 RaMaTraF.

Appendix A: Perturbative evolution for eigenvalues Let dH t a matrix with entries √ dt ξ i,j with ξ i,j ∼ N (0, 1) iid. Letting {|v i } an orthonormal basis of vectors, let us define A i,j as √ dtA i,j = v i | dH t |v j , which is a Gaussian variable as a linear combination of independent Gaussian variables. Furthermore A i,j and A i ,j with a different pair of indices are independent, indeed we have:

Using these notations, let us compute the perturbation theory for eigenvalue evolution in the Matrix Kesten model.

Hermitian case β = 2

The previous computation extends to the Hermitian case by taking dB t = dHt+dH T t +i(d Ht-d HT t ) 2

where dH t and d Ht are independent samples of the same distribution. The perturbative evolution of eigenvalues can be computed as above. Indeed, defining η i,j and ηi,j independent sets of standard gaussian variables:

We conclude that the evolution of the eigenvalues of U is, in the Hermitian case:

We conclude with the full evolution of the main text, by adding the drift that was left aside in this section. The perturbative evolution of eigenvalues of U t evolving according to [START_REF] Lacroix-A-Chez-Toine | Non-interacting fermions in hard-edge potentials[END_REF] is given, in symmetric and Hermitian cases, by:

Eigenvalue evolution of the Grabsch-Texier model

We derive the perturbative evolution of the set of eigenvalues of the matrix Z t evolving as the Grabsch-Texier model:

where the matrix is real symmetric (β = 1) or complex Hermitian (β = 2), and the noise matrix is structured as in Eq. ( 22):

where H t and Ht are, as above, both filled with N 2 independent standard real Brownian motions. This reproduces the noise correlations given in Eqs. [START_REF] Saichev | Theory of Zipf 's law and beyond[END_REF][START_REF] Kesten | Random difference equations and Renewal theory for products of random matrices[END_REF] of [START_REF] Grabsch | Wigner-Smith matrix, exponential functional of the matrix Brownian motion and matrix Dufresne identity[END_REF] by setting σ β to β 2 times a constant. In order to derive the evolution of the eigenvalues, it is now convenient to take the matrix SDE (A14) to the Itô prescription, which gives:

This is obtained as follows. First one recalls that in Stratonovitch prescription, the noise part of (A14) has the form

). Then, expanding (at order dt

2 (Z t dB t + dB t Z t ) and using the noise structure given in (A15) leads to (A16).

Let us truncate the evolution as above and keep the noise term only dZ t = σ β (Z t dB t + dB t Z t ). The perturbative evolution of eigenvalue λ i is:

From the noise structure (A15), we obtain with 2N 2 independent η i,j and ηi,j ∼ N (0, 1) variables:

Adding finally the drift terms, the perturbative evolution of the eigenvalues λ evolving under the Itô SDE (A16) (or equivalently the Stratonovitch SDE (A14)) is:

Reordering terms, we obtain:

We recover the evolution obtained for our model in Eq. (A13) by fixing

Appendix B: Free probability

RMT transforms and main results

Free probability theory is the study of free random variables, where freeness is a generalization of independence to non-commuting variables. We will not detail the definition of freeness further than stating that two large matrices are free if their eigenbasis are randomly rotated with respect to one another. For detailed introductions to the field and its application to random matrices, see [START_REF] Potters | A First Course in Random Matrix Theory for Physicists, Engineers and Data Scientists[END_REF][START_REF] Bun | Cleaning large correlation matrices: Tools from Random Matrix Theory[END_REF][START_REF] Mingo | Free Probability and Random Matrices[END_REF][START_REF] Tulino | Verdú Random Matrix Theory and Wireless Communications[END_REF][START_REF] Novak | Three lectures on free probability[END_REF]. In this Appendix section, we detail the main RMT and free probability tools that we use in this paper. Let us consider a large N × N matrix M , which has a one-cut spectrum with a density of eigenvalues ρ = 1 N λ∈Sp(M ) δ λ such that:

A standard RMT object is the Stieltjes transform, defined on C \ [a, b]:

We note that the density can be retrieved from the Stieltjes transform as:

We also give the following scaling and inverting properties of g (if M is invertible):

The behaviour of g at infinity is clear from the definition, as well as the general series expansion in terms of the moments ϕ k of M , with The Wishart ensemble consists of N × N matrices M = X † X where X is a T × N matrix (T N if T and N are integers, T > N -1 in the general case) with independent entries following a real centered gaussian law of variance where q = 1 2κ+1 and

In addition to the rescaling by T , which is obvious from the previous subsection, we insert an additional factor 1 -q such that the mean of the limit spectrum is 1. In the main text, we denote the rescaled matrix by W = (1 -q)T × M . This limit distribution is characterized by the following Stieltjes, R and S transforms:

From the joint eigenvalue density of the inverse-Wishart distribution, we can obtain the marginal density for the highest eigenvalue λ 1 in the limit where it is very large, with a factor N for the choice of the eigenvalue pushed to large values:

Inserting in this equation the values for C, T given in [START_REF] Derrida | Velocity and diffusion constant of a periodic one-dimensional hopping model[END_REF] one obtains the equation [START_REF] Doussal | First-passage time for random walks in random environments[END_REF] given in the text. For some results on the distribution of the smallest eigenvalue of the Wishart matrix, see [START_REF] Edelman | Eigenvalues and Condition Numbers of Random Matrices[END_REF].

We note that recent mathematical interest in the Wishart and inverse-Wishart ensembles has unveiled connections with the combinatorics of Hurwitz numbers [START_REF] Cunden | Integer moments of complex Wishart matrices and Hurwitz numbers[END_REF][START_REF] Gisonni | Laguerre Ensemble: Correlators, Hurwitz Numbers and Hodge Integrals[END_REF] and an explicit ergodic decomposition [START_REF] Assiotis | Ergodic Decomposition for Inverse Wishart Measures on Infinite Positive-Definite Matrices[END_REF].

Appendix D: Morse potential

We review the single-particle quantum system in a Morse potential [START_REF] Morse | Diatomic molecules according to the wave mechanics. II. Vibrational levels[END_REF][START_REF] Dong | Factorization methods in quantum mechanics[END_REF]. We are looking to solve the following eigenproblem, where g > 0 and x 0 are two constants:

Let us take z = 2ge -(x-x0) and define φ(z) = z g-k-1 2 e -1 2 z f (z). Inserting this in the differential equation gives:

The eigenproblem is solved if f solves the generalized Laguerre equation 0 = zf (z) + (2g -2k -z)f (z) + kf (z) such that regularity conditions impose that f = L (2g-2k-1) k is a generalized Laguerre polynomial. The k-th eigenfunction is then:

where

. This is thus a state in the discrete spectrum, with corresponding eigenvalue:

The maximum number of fermions that can be stacked in the discrete spectrum of the Morse potential is the largest integer N such that the last level exists i.e. N -1 < g -1 2 ⇐⇒ N < g + 1 2 . We denote this maximum occupation number by a strict integer part N = g + 1

2

-.

In addition to the discrete spectrum of bound states at negative energies k , there is a continuum of positive energy eigenstates = p 2 > 0. One can check that the associated eigenfunctions are of the form ψ p (x) ∼ e (x-x0)/2 W g,ip (2ge -(x-x0) ) in terms of the Whittaker W function, such that ψ p is a real function decreasing fast for x < 0 and oscillating for x > 0.

Putting all together we can write the Green's function as follows [START_REF] Zhang | Morse Potential, Contour Integrals, and Asian Options[END_REF]:

We point to [START_REF] Duru | Morse-potential Green's function with path integrals[END_REF] for an expression of this Green's function in path integral form.

Appendix E: Exact solutions for the evolution of S and Stieltjes transforms

S-transform

In this Appendix we study the evolution equation ( 109) for the S transform, first to linear order near the stationary solution and next in an exact manner using an hodograph transformation.

Recalling the stationary solution S stat (ω) = -m -σ2 2 ω and writing S(ω, t) = S stat (ω) + ∆S(ω, t), we have the following equation for ∆S :

Expanding to linear order in ∆S, one obtains a linear differential equation whose solution is given in ( 114). An exact solution of the non-linear equation ( 109) can be obtained using the hodograph transform, we point to [START_REF] Polyanin | HandBook of nonlinear partial differential equations[END_REF] for a general presentation of the technique. The result is given in an implicit form for ∆S(ω, t) as

where the function Ω(u) is obtained by inverting

This result is obtained as follows. Let us rewrite the equation ( 109)

We now perform the hodograph transformation, that is we take S and ω as variables and consider t as a function t(ω, S). The rules for the partial derivatives are such that ∂ S t = 1 ∂tS and ∂ ω t = -∂ωS ∂tS . The equation (E4) leads to

which is a first-order linear differential equation for the function t(ω, S). Let us denote u(ω, S) a solution of the homogeneous equation

A general method in first-order partial differential equations, see [START_REF] Polyanin | HandBook of first-order partial differential equations[END_REF], consists in choosing now ω and u as independent variables, i.e. we look for a solution of (E5) in the form

Using that u is a solution of (E5) setting the r.h.s. to zero, we obtain the following equation for the function t(ω, u)

where S(ω, u) = 1 2ω (u -σ2 ω 2 ) -m is obtained from inverting u(ω, S). The general solution of this equation is

Appendix F: Moments

In this Appendix, we compute the evolution of the first moments in the discrete and continuous versions of the scalar and matrix Kesten models. For clarity, we recall the relatively standard results for the scalar case and in the matrix case we show that there are considerations of invariant theory, which are interesting as they constrain the stationary measure.

We check that the large N limit of ϕ (N ) 2,n=+∞ , while fixing σ = σ/ √ N , coincides for β = 1, 2 with the r.h.s. of (119). The corrections are O(1/N β ) for β = 1, 2. We also find in this limit that

where each term of the l.h.s. are O(N 2 ) but the difference is O(1).

Continuous matrix case

Let us consider the Itô matrix stochastic differential equation ( 27): dU t = (I + mU t ) dt + σ √ U t dB t √ U t , with U 0 = 0, and study the evolution of moments ϕ

We have the following for the process E[U t ]:

where ϕ

(N )

1,t=+∞ coincides with the → 0 limit of the discrete result (F7). For the second moment, the evolution of the squared matrix is, according to Itô rules:

where we used that E[dB t U t dB t ] = E[Tr(U t )I + U t δ β=1 ], as seen from the following correlations:

Taking the expectation value one finds Equations (F23) and (F25) form a closed first-order differential system. For the sake of simplicity, we only give the stationary solutions ϕ We see that ϕ ) coincide as expected with the → 0 limit of the discrete results (F15) and (F16). In the limit of large N with σ = σ/ √ N , we find that

which coincides with the limit → 0 of the discrete result (F19). We recover the scalar continuous result (F6) by injecting N = 1 and β = 1 (or equivalently by also rescaling σ → √ 2σ in the complex case β = 2). Finally, we stress