Random-matrix models of monitored quantum circuits

We study the competition between Haar-random unitary dynamics and measurements for unstructured systems of qubits. For projective measurements, we derive various properties of the statistical ensemble of Kraus operators analytically, including the purification time and the distribution of Born probabilities. The latter generalizes the Porter-Thomas distribution for random unitary circuits to the monitored setting and is log-normal at long times. We also consider weak measurements that interpolate between identity quantum channels and projective measurements. In this setting, we derive an exactly solvable Fokker-Planck equation for the joint distribution of singular values of Kraus operators, analogous to the Dorokhov-Mello-Pereyra-Kumar (DMPK) equation modelling disordered quantum wires. We expect that the statistical properties of Kraus operators we have established for these simple systems will serve as a model for the entangling phase of monitored quantum systems more generally.

of the conventional MIPT, with a "trivial" critical point at p c = 1.Nevertheless, we show below that these models yield a plethora of new and exact results, some of which should capture universal features of the entangling phase in spatially local systems, though we will not explicitly consider spatially local systems in this paper.
We achieve this by applying a combination of ideas from random-matrix theory and disorder physics to analyze this problem, which differ from the field-theoretic replica methods that have largely been the analytical tools of choice for analyzing monitored dynamical phases in previous work [12].This approach allows us to determine analytically the behaviour of various physical quantities that were previously mostly accessible by numerical simulations or by heuristic arguments, including the purification time, the dynamics of Rényi entropies and the distribution of Born probabilities.We also introduce an unstructured model with near-identity weak measurements that can be seen as a monitored analogue of Dyson Brownian motion (see Refs. [9] and [11] for related constructions).The Fokker-Planck equation describing this model in the continuous time limit turns out to coincide with a known solvable model of Calogero-Sutherland type [13].
This exact solution grants us a thorough understanding of the dynamics of the full set of singular values of Kraus operators along quantum trajectories, for all times.
Our perspective in this paper differs from most previous work on monitored quantum circuits in another important respect.Monitored dynamical phases and the transitions between them are usually diagnosed by averaging physical quantities over ensembles of quantum trajectories weighted by their Born probabilities, and the main objects of study are the numerical values of these averages [12].The quantities being averaged are furthermore usually nonlinear functions of the density matrix on each quantum trajectory, leading to a "post-selection barrier" that hinders direct comparison with experiments [14][15][16].Post-selection and averaging tend to be advocated on the grounds that they are indispensable for probing "typical" quantum trajectories, which determine the monitored dynamical phase but are invisible at the level of the conventional density matrix.
In this paper, we instead study time evolution along typical quantum trajectories directly (see also Ref. [8] and [17]).Considering time evolution along single trajectories in lieu of Born-rule averaging can be motivated by analogy with the ergodic hypothesis in classical statistical physics, which similarly permits replacing ensemble averaging in chaotic systems with time-averaging along individual phase-space trajectories.However, in order to apply such reasoning to monitored quantum systems, we must first understand the extent to which distinct quantum trajectories are statistically alike.We must also understand the shape of the probability distributions that capture this statistical similarity.The latter will depend on both time and on the specific system under consid-eration, raising the question of how far such probability distributions reflect universal features of the underlying dynamical phase.To our knowledge, these probability distributions have not been studied in detail.Here, we derive them for monitored Haar-random quantum dots and return to questions of universality in the conclusion.
The paper is structured as follows.We first consider projective measurements and introduce the statistical ensembles of Kraus operators that will be the central objects of study.We present exact expressions for the Lyapunov spectrum of these Kraus operators, based on a mapping to so-called "truncated unitary ensembles" in random matrix theory [18,19], from which we determine the purification time analytically.We then derive the exact distribution of Born probabilities, and explain how this generalizes the so-called Porter-Thomas distribution [20][21][22] for random unitary circuits to a two-parameter family of distributions determined by both the measurement density p and the circuit depth t.Finally, we consider time evolution starting from maximally mixed initial states, which clarifies earlier results [7,10,23] on the time evolution of Rényi entropies as reflecting a crossover in time from a narrow Wigner-semicircle-like to a broad (and to a first approximation log-normal) distribution of singular values of Kraus operators.
We next turn to the weakly measured case, for which we derive a DMPK-like equation [24,25] that captures the time evolution of the joint distribution function of singular values of Kraus operators.The resulting Fokker-Planck equation is exactly solvable [13] in a manner analogous to the time-reversal-symmetry-breaking case of the DMPK equation [26].Using this exact solution, we demonstrate explicitly that the joint distribution function of singular values exhibits a "semicircleto-square" crossover between log-GUE statistics and log-normal statistics as a function of time.
This represents a remarkably complete understanding of this model's dynamics that corroborates our results on projective measurements.We conjecture that this behaviour is universal for entangling phases of monitored quantum systems, in the same manner that random matrix theory captures the spectral properties of generic closed quantum systems.throughout this paper.Each measurement layer comprises independent, projective, single-qubit Ẑ measurements acting on a mean number of qubits pL per layer, with 0 ≤ p ≤ 1.The specific spatial probability distribution of these measurements will depend on the specific model of interest, to be fixed below.The difference between the models that we consider here and more standard [12] spatially local Haar random quantum circuits is depicted schematically in Fig. 1.We denote the string of measurement outcomes in the jth measurement layer by m j , and denote the measurement history of the whole circuit along a given quantum trajectory by the tuple m = (m 1 , m 2 , . . ., m t ).
We denote the projection operator corresponding to projection onto the measurement outcome m j by Pm j .In general, the rank of Pm j equals 2 L−|m j | , where |m j | denotes the number of measurements made at time j.
For a pure initial state |ψ(0)⟩, the time-evolved state along a given quantum trajectory m can be written as [12] where the Kraus operator and the Born probability of the string of measurement outcomes m is given by This can be generalized to arbitrary initial density matrices ρ(0), whose time evolution under this dynamics is given by Where applicable, we make the choice [7] to unravel this evolution as an ensemble of singletrajectory density matrices drawn with probabilities p(m) = Tr[ Km (t)ρ(0) K † m (t)].We refer to the review article Ref. [12] for a more detailed discussion of Kraus operators for monitored quantum circuits.When considering mixed initial states, for simplicity we will always restrict our attention to the maximally mixed initial state for which along each quantum trajectory.
The basic quantity of interest in this work is the statistical ensemble of Kraus operators { Km }, defined in general by first sampling over the Haar measure, then sampling over measurement locations and finally sampling over measurement outcomes.Since we will be focusing on Haar random unstructured systems, the latter average will mostly be redundant, and expectation values E will always denote expectation values along a fixed quantum trajectory with respect to the Haar measure on each unitary layer unless specified otherwise.Thus we do not explicitly consider Bornrule averaged quantities in this work, although we address the distribution of Born probabilities with respect to the Haar measure in Section II E.
The main tool at our disposal for studying the ensemble of Kraus operators Km will be their singular value decomposition, where we write Dt = diag(σ 1 (t), . . ., σ N (t)), with the convention that σ 1 (t) ≥ . . .≥ σ N (t) ≥ 0. We emphasize that the matrices of singular vectors Vt , Ŵt and the singular values σ n (t) are random variables that depend sensitively on the circuit realization, and suppress the measurement history in Eq. ( 8) for notational convenience.

B. Rank collapse and dynamical purification
Because the Kraus operators above generically include projective measurements, their rank ] is a random variable that does not increase in time and σ n (t) = 0 for n > r(t).
The decay of the rank defines a "rank-collapse time", given by the mean stopping time which is infinite if the rank does not decay.
In general, the singular values σ n (t) will decay in an average sense as t → ∞; in particular, the Oseledets ergodic theorem [27] guarantees the existence of Lyapunov exponents with probability one, provided that n ≤ r(t) with probability one for all time.To see the physical meaning of the first few singular values and their Lyapunov exponents, suppose that λ 2 is defined and consider time evolution from a maximally mixed state, as in Eq. (7).Then where v n denotes the nth column of Vt in Eq. (8).Thus the Born probability decays at a rate τ −1 = 2λ 1 [8], while the trajectory density matrix It is clear that the decay of the ratio determines the rate of convergence of ρm (t) to a pure state.We thus define the "purification time" For more general initial states, such as pure states, τ −1 P is properly thought of as the expected rate of convergence of the trajectory density matrix to a rank-one projection along the random singular vector v 1 .The specific distribution of v 1 will depend on the specific model and monitored dynamical phase under consideration.Thus a more general perspective on the purification time τ P is that it defines the timescale on which a monitored quantum system "forgets" its initial state [9] and begins to reveal universal properties of its monitored dynamical phase.We return to this point in the conclusion.
The above discussion reveals that there are two important timescales that determine the fate of the monitored quantum dot (and indeed arbitrary projectively measured quantum circuits) at asymptotically long times, namely the rank-collapse time and the purification time.However, these two timescales are in tension, because standard results [27,28] on the existence of Lyapunov exponents λ n for n > 1 require that σ n (t) is generically non-zero for all time.Thus in order to define the purification time or higher Lyapunov exponents rigorously, we require that τ R.C. = ∞.This is not the case for the standard formulation of monitored random circuits, according to which measurements are performed randomly and independently at every site [2] leading to rank collapse in finite time.Previous work on models with rank collapse implicitly assumes a parametric separation of scales in L, and then estimates τ P numerically by simulating the system for times much shorter than τ R.C.
(see [8]).The downside of this approach is that there is an inherent "fuzziness" in the definition of τ P on the order of inevitable statistical fluctuations 1/ √ τ R.C. induced by rank collapse.While this is mostly a mathematical subtlety when it comes to estimating the purification time numerically (since for the canonical spatially local model τ R.C. ∼ 1/(2p) L ≫ 1 by a mapping to bond percolation [4]), it becomes a serious obstacle even for numerical estimates of higher Lyapunov exponents λ n with n > 2, which succumb to rank reduction at much earlier times than τ R.C. .
We resolve this subtlety for monitored quantum dots by noting that the timescales τ R.C. and τ P naturally pertain to two distinct microscopic models, with and without rank collapse respectively, for which the separation of scales Eq. ( 16) can be proved analytically as L → ∞.This then justifies attempting to estimate τ P numerically in the model with rank collapse.These models, which we refer to as "Model I" and "Model II" respectively, differ only in the spatial distribution of measurements in each measurement layer and are defined as follows.
For Model I, the measurements in each measurement layer are performed randomly and independently at each qubit with probability p, as for the standard formulation of the MIPT [2].Then the number of measurements in each layer fluctuates and rank collapse occurs when all the qubits in a given layer are measured, which occurs with probability p L .For Model II, we perform a fixed number of measurements pL ∈ {0, 1, . . ., L} per layer.Note that for the monitored quantum dots under consideration in this paper, the specific location of these measurements does not matter, by Haar randomness of the unitary layers.
For both models, we can determine the rank-collapse time in finite systems analytically, in contrast to the situation for spatially local models, for which closed forms are not easily obtained (even if the asymptotic behaviour is understood [4]).First consider Model I. From Eq. ( 9), the rank-collapse time is the expected value of the stopping time T such that r(T ) = 1 and r(t) > 1 for t < T .Note that P( Next consider Model II.In this case r(t) = 2 (1−p)L > 1 almost surely for p < 1 and it follows that Below, we will restrict our attention to Model II with 0 < p < 1 and write M = 2 (1−p)L for the rank of its Kraus operators.

C. Mapping to truncated unitary and Ginibre ensembles
The "truncated unitary ensembles" of random matrices [18] consist of square submatrices of Haar random matrices drawn from U (N ).We can relate the Kraus operators Eq. ( 2) for Model II to products of M × M truncated unitary matrices drawn from U (N ) as follows.We first note that by Haar randomness and transitivity of U (N ) on measurement outcomes, the statistics of singular values of { Km (t)} is unchanged if we fix the measurement outcomes Pm j = P in each layer to be the same.Writing ∼ for equality of singular value statistics, we have at each time step t.To proceed further, we use the fact that P 2 = P to write P Ût . . .P Û1 = ( P Ût P )( P Ût−1 P ) . . .( P Û2 P ) P Û1 , implying that Km (t) where in the computational basis, Rj = P Ûj P is an M -by-M truncated unitary matrix padded by N − M rows and columns of zeros, with the same configuration of nonzero entries for each j, and Ŝ1 = P Û1 has M unit singular values and N − M zero singular values, since It follows from this analysis that the singular values of Km (t) are distributed as the singular values of products of t − 1 truncated unitary matrices.While various properties of products of truncated unitary matrices have been derived analytically in the literature [19,[29][30][31], including explicit expressions for their Lyapunov spectrum that we will discuss further below, the simplest results are obtained in the limit M/N → 0, in which K is distributed as a product of Ginibre random matrices [29].In the physical context of monitored random circuits introduced above, this regime is realized in the thermodynamic limit L → ∞ in which the number of qubits tends to infinity at any fixed p > 0.
Intuitively, the emergence of Ginibre ensembles corresponds to the fact that as M/N → 0, the non-zero elements of Rj look "random", in the sense that they are asymptotically distributed as the elements of a Ginibre matrix, up to rescaling [32,33].To be precise, let Âj denote an M -by-M complex Ginibre matrix whose elements are i.i.d. with probability distribution function This choice of normalization ensures that in the absence of truncation, M = N , the expected squared norm of each row and each column of Bj coincides with the expected squared norm of each row and each column of an N -by-N unitary matrix.It can then be shown [29] that as has the same probability distribution as the non-zero elements of Rt Rt−1 . . .R2 , and therefore that the non-zero singular values of Km (t) and Ĝ(t) are identically distributed in the thermodynamic limit.We note that the spectral density of the product Ĝ(t) has been studied in some detail [34,35], and will make use of some of these results below.

D. Computation of Lyapunov exponents and the purification time
We now compute the entire Lyapunov spectrum for Model II.This follows immediately from the mapping Eq. ( 21) from Model II to the products of truncated unitary matrices, combined with an expression for the Lyapunov exponents of such products derived recently [19,30,31].Together, these results imply that the Lyapunov exponents for Model II are given by where ψ(x) = Γ ′ (x)/Γ(x) denotes the digamma function.
It follows from Eq. ( 25) that the inverse purification time for Model II by the standard digamma function identity [36] for p < 1, implying that for L ≫ 1 and p < 1, which is one way to make Eq.( 16) precise.The functions of p appearing on either side of Eq. ( 27) are plotted in Fig. 2; it is clear that they coincide only in the approach to the "transition" as p → 1 − .In particular, we have proved that the purification time grows exponentially in the system size, and having established Eq. ( 16) we henceforth suppress the subscript specifying Model II.

E. The distribution of Born probabilities
Now consider applying Model II dynamics to an arbitrary initial pure state |ψ⟩, with a view to determining the distribution of Born probabilities as a function of both the density of measurements  This implies that where Note that so far we have assumed nothing about the distribution of unitary layers Ûj .If we now assume that Ûj are Haar-random matrices on the full Hilbert space, a drastic simplification occurs and the measurement outcomes at distinct time steps become pairwise uncorrelated.Explicitly, we have for all j, where the v j,k are i.i.d.complex Gaussian random variables with probability density Thus let X K,j denote a set of independently distributed gamma random variables, each with p.d.f.f K (x).Then follows a beta distribution with parameters M and N − M and p.d.f The usual Porter-Thomas distribution for random unitary circuits [21] corresponds to measuring every qubit after each unitary layer, and is recovered from Eqs. ( 34) and ( 35) in the limit that p = 1 and M = 1, for which is asymptotically exponentially distributed with mean 1/N for large N .
Let us now return to the original problem of determining the probability distribution function of p(m) for Model II.We showed above that for i.i.d.beta random variables Y j ∼ Beta(M, N − M ).This implies in particular that log p(m) from which log-normality of p(m) at long times is immediate.It follows by the central limit theorem that log p(m) in distribution, where the mean the variance and Z t ∼ N (0, 1) is a unit normal random variable.Thus we have shown that the Born probabilities p(m) for Model II are asymptotically log-normally distributed, with a shape that is determined by the parameters N , M and t.
Exact expressions for the distributions of Born probabilities can be obtained as follows.Denoting the characteristic function of a given random variable W by φ W (θ) = E[e iθW ], it follows by Eq.
(37) that The characteristic function of the beta distribution Eq. ( 35) is given by implying that the PDF for the random variable log p(m) after t time steps is given by Note that the constraint x < 0 closes the contour in the upper half-plane.This contour encloses all the poles of the integrand, since by the recurrence relation for Gamma functions the integrand has N − M distinct poles of order t, which are evenly spaced along the positive imaginary axis, θ j = i(M +j) for j = 0, 1, . . ., N −M −1. Thus with g(θ) given by Eq. (45).
In principle, this calculation yields an exact expression for the p.d.f.f log p(m) (x) of log Born probabilities for all times.In practice, these expressions rapidly become cumbersome.At t = 1 we recover Eq. ( 35) up to the necessary change of variables, while for t = 2 we find that For all times, Eq. (37) implies that in the large system limit, and similarly that Thus we have proved that for large systems, the "typical" Born probability, i.e. the mean of log p(m), coincides with the uniform distribution over all 2 pLt possible measurement outcomes up to time t, while the distribution of p(m) remains narrow on a log scale until times comparable to the purification time.

F. Dynamics of Rényi entropies
We now illustrate how the above results provide an analytical handle on time evolution along quantum trajectories, starting from a maximally mixed initial state as in Eq. ( 7).We will focus on the behaviour of Rényi entropies with α > 0 (including the von Neumann entropy defined as the limit α → 1) along a given quantum trajectory, which can be expressed in terms of the corresponding Kraus operator's singular values as [10] S (α) We discuss the limits of short and long times separately.Our analysis is similar in spirit to that of Ref. [10], with a greater degree of analytical control owing to the results of previous sections.

Short times
We assume throughout this section that L ≫ 1, so that the approximation of Kraus operators by a product Ĝ(t) of t−1 Ginibre matrices as in Eq. ( 24) holds, and let ρ(σ, t) = M n=1 δ(σ −σ n (t)) denote the density of singular values of Ĝ(t).Then the Rényi entropies are determined by the ratios Let us write where ρ(σ, t) = E[ρ(σ, t)], for the mean and fluctuations of ρ(σ, t) about its mean respectively.
Similarly, we write is dominated by the moments of the mean spectral density, provided the leading fluctuation cor- m (1)   is small.Let us assume this at times that are short compared to the purification time (for an analytical argument that such short-time fluctuations are small in the context of weak measurements, see Section III D).Then the Rényi entropies are dominated by their non-fluctuating part In the regime of times 1 ≪ t ≪ M , we have [35] ρ(σ, t) yielding which predicts that for all α > 0, recovering a result that had been previously derived in various specific limits using various distinct physical arguments [7,10,23,37] and is consistent with our precise definition of the purification time in Eqs. ( 15) and (29).
For α = 2, 3, . . ., M (but not the von Neumann entropy) a more refined estimate is available, using the exact result [34] for such integer moments.This yields for the non-fluctuating contribution to the Rényi entropies.The simplest non-trivial case of this formula is the second Rényi entropy for which which is easily verified to recover Eq. ( 57) in the limit t ≪ M .

Long times
We next consider the dynamics of Rényi entropies at times that are long compared to the purification time, for which the leading singular values of the Kraus operators are well separated from one another with high probability, so that the semicircle law Eq.( 55) no longer provides a good approximation to their distribution.First note that as t → ∞ the squared singular values are distributed as [19,30] where the Y n are independent normal variables with mean λ n and variance In particular, the Rényi entropies and von Neumann entropies are dominated by the largest two singular values with high probability as t → ∞, implying that in terms of the ratio ν(t) defined in Eq. ( 14) (similar expressions were obtained in Appendix E of Ref. [10]).By our definition of the purification time Eq. ( 15), it follows that as t → ∞.Thus we have constructed a solvable model that confirms earlier proposals for the late-time dynamics of entropy in entangling phases of spatially local monitored systems [7,10,23], together with the idea that at long times t ≫ τ P , the purification time τ P captures the typical [10] behaviour of a broad (and to a first approximation log-normal [30,38]) distribution of S α m (t).Unfortunately, a more detailed characterization of the late-time distributions of ν(t) and S α m (t), including their mean (rather than typical) values, is beyond the scope of the long-time prediction Eq. ( 61) and related expressions [38], because these results do not contain information about the joint distribution of σ 2 1 (t) and σ 2 2 (t) that would be necessary to accurately model ν(t).We therefore turn to a model with weak measurements, for which the joint distribution function of singular values can be obtained analytically.

III. WEAK MEASUREMENTS A. Model and Fokker-Planck equation
We now consider a generalization of Model II defined in Section II B above, whose only distinction from Model II is that its measurement layers consist of independent weak measurements on pL qubits per layer, which act on a single qubit as where we fix 0 ≤ ϵ ≤ 1.We again write the Kraus operators along a single quantum trajectory at time t as Km (t) = Pmt Ût . . .Pm 1 Û1 , where Pm j now denotes the tensor product of pL weak measurement operators as in Eq. ( 65), and the identity on the remaining (1−p)L qubits.Regardless of the specific set of measurement outcomes, which we mostly suppress for economy of notation, we can write where Λj has eigenvalues with respective multiplicities It follows from these expressions that We would like to understand the evolution of the singular values σ 1 (t) ≥ . . .≥ σ N (t) ≥ 0 of Km (t).
Thus consider the singular value decomposition Km (t) = Vt Dt Ŵ † t as in Eq. ( 8).Then Letting Û = Ût+1 V † t , which inherits Haar randomness of Ût+1 at time t + 1, we can write Introducing the operators B = Û † Λt+1 Û and Xt = 2 2pLt Km (t) Km (t) † , we have Xt+1 where as above ∼ denotes a common distribution of singular values (which are eigenvalues in this case).Let x 1 (t) ≥ . . .≥ x N (t) ≥ 0 denote the eigenvalues of Xt , which are given by x n (t) = 2 2pLt σ 2 n (t) in terms of the singular values σ n (t) of the Kraus operators.Then, if pLϵ 2 ≪ 1 so that typical eigenvalues of Λ are small assuming large pL, Eq. ( 72) can be expanded perturbatively in ϵ to yield Collecting terms yields a perturbative model for the time evolution of the singular values x n (t), namely This equation splits naturally into multiplicative "noise" and "drift" contributions, given by Thus To eliminate multiplicative noise, we write x n (t) = e yn(t) and, again assuming that ϵ is small, find that y n (t + 1) − y n (t) = log (η n (t + 1) + ∆ n (t + 1)) = η n (t + 1) + D n (t + 1) + O(ϵ 3 ), where the new drift term To proceed further, we note that the matrix elements and are amenable to Haar averaging.In particular, the two-point function and the four-point functions by standard results [39], where angle brackets ⟨...⟩ denote the Haar average over Û ∈ SU (N ), implying that and where the noise strength , implying that the average value of the drift term (conditioned on the circuit realization at time t) is given by m̸ =n e ym(t) e yn(t) − e ym(t) (87) Thus, assuming that fluctuations of D n (t + 1) about its Haar-averaged value are negligible, we obtain the Langevin equation m̸ =n e ym(t) e yn(t) − e ym(t) . ( This can be written in a more symmetric form by making the change of variables 2z n (t) = y n (t) + ΓN 4 t, which yields Finally, letting Γ → 0 yields the continuous-time Fokker-Planck equation where P (⃗ z, s) denotes the joint PDF of ⃗ z(t) = (z 1 (t), z 2 (t), . . ., z N (t)) at time and In particular, the timescale Γ −1 determines the dynamics of singular values, and should be thought of as the purification time for these models, an interpretation that will be confirmed by our results for the short and long time dynamics of Rényi entropies below.
The continuous-time Fokker-Planck equation Eq. ( 90) was first derived [13] in the context of a stochastic matrix model known as "isotropic Brownian motion" [40].In that work, it was also pointed out that Eq. ( 90) is exactly solvable via a connection to Calogero-Sutherland models.In order to be self-contained, we now derive this exact solution.

B. Exact solution of the Fokker-Planck equation
To solve Eq. ( 90), we first note that D n can be written in gradient form where the "prepotential" arises naturally in the theory of the classical hyperbolic Calogero-Sutherland model [41].This implies that Eq. ( 90) can be written as [42] ∂ Letting ψ = e Φ/2 P then yields an imaginary-time Schrödinger equation for ψ, namely [42,43] with an effective potential To proceed further write z nm = z n − z m and note that 1 4 By an identity apparently first published by Calogero and Perelomov [44], we have yielding a complete cancellation of interactions in the effective potential which is reminiscent of (but simpler than) the cancellation of interactions that occurs for the β = 2 case of the DMPK equation [26].Making the change of variables ψ = e − N (N 2 −1) s ψ finally reduces the Fokker-Planck evolution Eq. ( 90) to a free diffusion equation (albeit with non-trivial boundary conditions to be discussed shortly).
It remains to solve Eq. (102) in the "ordered sector" z 1 ≥ z 2 ≥ . . .≥ z N , subject to the initial condition where the regulators ε 1 > ε 2 > . . .> ε N > 0 will be taken to zero at the end of the calculation and are needed to avoid the singularities in Eq. ( 90) at collision planes z j = z k .We further impose boundary conditions for all s ≥ 0 that are sufficient for the vanishing of probability flux at collision planes, and therefore guarantee conservation of probability within the ordered sector.
Our method of solution follows that of Refs.[26,45] for the β = 2 DMPK equation but differs in its details.It is first useful to note that for every solution to Eq. ( 102), there exists a solution to (90), given by In terms of the "wavefunction" ψ in Eq. ( 105), we find that the initial condition and the vanishing of the wavefunction at collision planes, for s ≥ 0, are sufficient for the initial and boundary conditions Eqs. (103) and Eq. ( 104) on P to hold.To satisfy the boundary conditions in the ordered sector, we extend ψ to the whole space and consider a "fermionic" initial condition where sgn : S N → ±1 denotes the sign of the permutation κ ∈ S N .We can write this more suggestively as a Slater determinant of single-particle factors where det[A jk ] denotes the determinant of the matrix A with elements A jk .Then the Slater determinant of heat kernels (110) satisfies both the diffusion equation Eq. ( 102) for s ≥ 0 and the initial and boundary conditions Eqs. ( 103) and (104) in the ordered sector.The final nontrivial step is to set the regulators ε j → 0.
To this end, we first note that the leading behaviour of the denominator of Eq. (105 is a Vandermonde determinant of order N − 1 in each of the ε j .We next expand the heat kernels in terms of Hermite polynomials [36] e to the same order in ε, where zj = z j / √ 4s, εk = ε k / √ 4s.However, at this order, the sum in Eq. ( 112) is nothing but a product of square matrices and so where in the second line we applied Gaussian elimination to the matrix of Hermite polynomials to obtain another Vandermonde determinant.Combining the above expressions, we deduce that solves the Fokker-Planck equation for all s ≥ 0, satisfies the initial condition P (⃗ z, 0) = N j=1 δ(z j ), and conserves probability in the ordered sector.This recovers the solution to Eq. ( 90) presented in Ref. [13].We note for future reference that the random variables z n (t) modelled by Eq. ( 114) are related to the singular values σ n (t) of Kraus operators at time t by the change of variables and recall that s is related to t by Eq. (91).

C. Very short times: emergence of log-GUE
Let us first consider the probability distribution Eq.(114) as s → 0 (note that s ≪ 1 corresponds to t ≪ Γ −1 , i.e. times that are short compared to the purification time).We have Thus the leading contribution to P (⃗ z, s) is from the distribution function P v.s.t.(⃗ z, s), given by This is strongly reminiscient of the distribution function of the Gaussian Unitary Ensemble (GUE) on the full space [46], and indeed the two distributions are related by a simple change of variables where the prefactor of N !arises because we restricted the domain of P (⃗ z, s) to the ordered sector.
We will call the regime of validity of the approximation Eq. ( 117) the "very-short-time regime", to be determined below.In the very-short-time regime, the singular values of Kraus operators follow a "log-GUE" distribution, which is related to the conventional GUE in the same sense that the log-normal distribution is related to the normal distribution.Defining the random variable , it follows by the Wigner semicircle law for the conventional GUE [46] that the non-fluctuating part of the level density for N ≫ 1 and very short times.The boundary of the very-short-time regime is set by the requirement in Eq. ( 119) that |z j − z k | ≪ 1 for all j, k: from Eq. (120), this is the case for We can use this observation to extend the analysis of Rényi entropies in Section II F 1 to very short times as follows.First note that in terms of ρ(z, t), the Rényi entropies along a given quantum trajectory can be written as In particular, we have where I 1 (w) denotes a modified Bessel function of the first kind [36], so that the non-fluctuating part of each Rényi entropy is given by The small-argument asymptotic behaviour implies perturbative time dependence of the Rényi entropies at very short times.Meanwhile, log-GUE statistics and the semicircle law Eq.( 120) break down at times 1 N ≪ Γt ≪ 1 that are short but not very short, as we now discuss.
D. Short times: semicircle-to-square crossover Let us now consider how the semicircle law corresponding to the exact distribution function P (⃗ z, s) differs from Eq. (120) beyond the regime of very short times.We first note that P (⃗ z, s) can be written as a Boltzmann weight for a harmonically trapped gas of particles on a line at positions where the effective potential interpolates between Dyson's "log-gas" [46] in the limit of small interparticle separations [13] |z j − z k | → 0 that was discussed in the previous section, and a one-dimensional Coulomb gas at large interparticle separations [47] |z j − z k | → ∞, as follows from the asymptotic behaviour log (z sinh z) ∼ log (|z|) as |z| → ∞.The latter model is also known in the literature as the "onedimensional jellium model", various properties of which can be derived analytically [48][49][50], such as its uniform density of states [48].
It will be instructive to understand this crossover to uniformity in terms of the density of states ρ(z).Suppressing explicit time dependence, we have Note also that ρ(z) satisfies the constraint Thus for large N ≫ 1, the non-fluctuating part of the density of states ρ(z) minimizes the functional Finally differentiating with respect to z yields the singular integral equation [51,52] which is distinguished from the usual semicircle law by the presence of the coth (z − w) term in the kernel, that breaks the spatial rescaling symmetry of the former.This breaking of scaling symmetry allows ρ to flow from a semicircle law to a uniform distribution as s increases.We note that an exact but implicit solution to the integral equation Eq. ( 130) was proposed recently in the literature [52].
We instead proceed numerically and solve for solutions to Eq. ( 130) supported on a finite interval [−a, a] with a > 0 for s > 0 (note that the microscopic distribution function Eq. ( 114) is inversion symmetric in z).Numerical results obtained by discretizing the principal value integral in Eq. ( 130) and inverting the resulting matrix to obtain ρ(z, s) are shown in Fig. 3 and confirm that the semicircle law discussed in the previous section quickly converges to a uniform distribution beyond the regime of very short times.The resulting late-time profile is sometimes called a "square law" [13] and in this sense Eq. (130) captures a "semicircle-to-square" crossover with increasing s.
Expanding the right-hand side of Eq. ( 132) perturbatively in z yields a linear estimate N z/a.We will show that Eq. ( 132) is close to this linear estimate over a large "bulk region" |z ± a| > ε, where ε > 0 is an order one constant.To this end, we define the function which quantifies the deviation of the right-hand side of Eq. ( 132) from linearity.The function η(z) is odd and strictly increasing in z on the interval (−a, a).It is also convex on the interval [0, a − ϵ].
Together these observations imply an upper bound valid over the bulk region.Finally, we note that to leading order in a and for ε of order one,  131) yields an accurate bulk solution provided we make the identification corresponding to the regime of times Γt ≫ 1 N , which strictly excludes the very-short-time regime identified in the previous section.Thus we have shown that from the short-time regime onwards, the uniform ansatz Eq. (131) yields a good approximation to the solution of Eq. ( 130), excluding a boundary region of width 2ε that has negligible measure compared to the bulk region.The latter point in particular means that the uniform ansatz Eq. ( 131) can reliably be used to estimate integrals over singular values when the separation between neighbouring singular values is small.
For example, at short times the uniform ansatz Eqs. ( 131) and (137) predicts that whose leading time dependence perfectly matches Eq. ( 57) for projective measurements, provided that Γ −1 is interpreted as the purification time.That Γ −1 indeed defines the purification time in the sense of Eq. ( 15) will be checked carefully below, but can be seen here by noting that the timescale for the characteristic spacing between adjacent z n to become comparable to their longtime standard deviation, invalidating the treatment of ρ(z, t) as a uniform distribution, is defined by the condition 2a(s)/N = O(s 1/2 ), i.e.Γt = O(1).
Finally, we note that the energy functional Eq. ( 128) appearing in the Dyson gas formulation of our problem can be used to estimate the strength of fluctuations of the Rényi entropies, thereby validating our estimation of short-time Rényi entropies by their non-fluctuating values as in Eq.
(138).We first expand the functional W perturbatively about the mean field ρ(w) to yield a quadratic effective action whose kernel It follows that the expected value of the correlation function δρ(z)δρ(w) with respect to the Boltzmann weight ∝ e −S is given by where the Green's function satisfies We refer to Refs.[53,54] for a careful justification of analogous arguments for the DMPK equation.
To proceed further, we note that the bulk short-time behaviour described above is consistent with the linear approximation to the kernel K, for which the Green's function defines an inverse in the bulk region.We note that Eq. ( 144) implies the results of Ref. [50] on the variance of linear statistics for the one-dimensional Coulomb gas, which is a nontrivial test of its validity.Let us define the mean moments n (α) = a −a dz ρ(z) e 2αz and their fluctuations δn (α) = a −a dz δρ(z) e 2αz .Combining the short-time approximations Eq. (131) and Eq. ( 144) and neglecting boundary effects yields the estimates for N Γt ≫ 1.From these expressions, we deduce that the leading contribution to the covariance of Rényi entropies grows quadratically in time for N Γt ≫ 1.In particular, since (Γt) 2 ≪ | log Γt| for Γt ≪ 1, this guarantees that fluctuations in the Rényi entropies will be small compared to the mean Rényi entropies S(α) until the end of the short-time regime Γt ≈ 1, at which the S(α) are order one by Eq. ( 138) and thus have magnitude comparable to their fluctuations Eq. ( 146).

E. Long times: log-normality, Lyapunov exponents and level repulsion
Anticipating that the z n are well-separated for s ≫ 1, we define the asymptotic drift velocities Noting that N n=1 c 2 n = 4N (N 2 − 1)/3 and j<k (z j − z k ) = 1 2 N n=1 c n z n , we can rewrite Eq. (114) as Thus as s → ∞, where the dominant contribution at asymptotically long times is given by since the Vandermonde determinant and The latter manifestly exhibits level repulsion as ζ → 0 + .With respect to the improved long-time measure Eq. ( 157), typical values of ν(t) behave as implying that which recovers Eqs. ( 158) and (159) at times t ≫ Γ −1 , while the mean value of ν(t) also decays exponentially, (albeit at a slightly different rate from the typical value), implying for example that the mean decay exponentially for t ≫ Γ −1 .

IV. CONCLUSION
We have presented various solvable models of monitored quantum circuits ("monitored Haarrandom quantum dots") that realize the entangling phase of monitored quantum dynamics.While these models have antecedents in the literature [9][10][11], the full extent of their analytical tractability does not appear to have been exploited until now.This analytical tractability has allowed us to derive the first exact expressions that we are aware of for quantities such as the purification time, the Lyapunov spectrum and the distribution of Born probabilities in the entangling phase of a monitored quantum system.By constructing explicit mappings from monitored Haar-random quantum dots to well-studied models in random matrix theory, such as products of truncated unitary matrices [18,29] and isotropic Brownian motion [13,40], we have further provided a template for realizing random-matrix universality in monitored quantum systems.
Our proposed notion of universality for such systems is analogous to the use of random matrix theory (RMT) as a description of the spectra of closed quantum systems; in the same way that RMT captures spectral correlations of generic chaotic quantum systems, we conjecture that the models discussed in this paper capture certain features of the entangling phase of generic monitored quantum systems, including the spatially local quantum circuits in which the entangling phase was first identified [2,4].
Before describing our proposal in detail, let us explain how universality for monitored quantum systems must differ from random-matrix universality and the related eigenstate thermalization hypothesis (ETH) [55] in closed quantum systems.Most obviously, the qualitative behaviour of the unstructured models considered in this paper depends on three parameters: the Hilbert space dimension N , the purification time τ P and the time t under consideration.The latter two parameters have no natural counterpart in closed quantum systems, whose dynamics is Hamiltonian and therefore time-translation invariant.In particular, we expect that monitored quantum systems will exhibit different facets of universal behaviour depending on the dimensionless parameter t/τ P , with qualitative differences between the short-time (t/τ P ≪ 1) and long-time (t/τ P ≫ 1) regimes.
Second, instead of the matrix elements of local observables and energy-level-spacing statistics that are usually of primary interest in formulating and testing ETH [55], the monitored setting suggests new arenas for universality, such as the statistics of single-trajectory Born probabilities and Rényi entropies discussed in this paper, both of which reflect the singular-value statistics of the underlying Kraus operators.
First consider the short-time regime, t/τ P ≪ 1.Our results on the distribution of Born probabilities in Section II E and on the sample-to-sample fluctuations of Rényi entropies in Section III D suggest that for the models studied in this paper, both the Born probabilities and the Rényi entropies will be narrowly distributed about their means (possibly on a log scale) in the short-time regime.We expect this conclusion to hold more generally for the entangling phase in arbitrary spatially local models.This means that on timescales that are short compared to the purification time (excluding initial transients as in Section III C), the behaviour of the system is to a first approximation decoupled from its measurement history, and to this extent, all quantum trajectories are statistically alike.
By contrast, our results on both Born probabilities and Rényi entropies at long times (Sections II E, II F 2 and III E) indicate that these quantities exhibit a broad and normal (or normalderivative-like) distribution on a log scale at long times t/τ P ≫ 1.Thus the observed behaviour in the long-time regime will depend strongly on the measurement history, and distinct quantum trajectories will exhibit very different properties.We expect similar behaviour for the entangling phase of spatially local systems at long times, even at the level of the shapes of these probability distributions, which should reflect universal properties of large products of identically distributed random matrices [35].
Finally, we note that the state-vector of the system at long times should exhibit some degree of universality by Eq. ( 13), which predicts that the long-time density matrix along any given trajectory is a random pure state.This is reminiscent of Berry's conjecture [56], which posits that the eigenvectors of a chaotic Hamiltonian are distributed as Gaussian random vectors.For the projectively-measured Haar-random quantum dot, the long-time state-vector v 1 in Eq. ( 15) is distributed as a Gaussian random vector in the M -dimensional image of the most recent measurement layer.More generally, for spatially local systems the long-time state-vector v 1 is no longer expected to be perfectly Gaussian for entangling phases, as can be seen from the presence of logarithmic in L corrections to the bipartite entanglement entropy of the long-time state in spatially local systems [57,58], which would otherwise exhibit purely volume-law dependence on L as predicted by the Page curve [59].
Note that in this paper, we do not consider the Born-rule weighted averages of quantities such as Rényi entropies.We emphasise that the possibility of Born-rule weighting does not arise for the Lyapunov exponents, which by definition characterize a product of independently and identically distributed random matrices.On the other hand, for Rényi entropies we expect that introducing Born-rule averaging will leave our results essentially unchanged at times t ≪ τ P but will have an effect at times t ≫ τ P .Accounting for Born-rule averaging within the approach in Section III would lead to a distinct Fokker-Planck equation from the one studied above, as discussed in Ref. [9].
Important goals for future work include directly testing the above predictions of universality in spatially local realizations of the entangling dynamical phase (for which there is currently a dearth of analytical understanding) and developing a similarly universal characterization of disentangling phases.Another interesting question is whether the kinds of quantities that we have computed analytically in this paper, such as Lyapunov spectra and distributions of Born probabilities, can shed further light on monitored phases of "non-interacting" or Gaussian quantum circuits [9,10,[60][61][62][63], for which it is known that the purification time scales quadratically with the system size [9,10].
In the Haar-random setting, it seems worth understanding how far the generalizations of the Porter-Thomas distribution that we identify in Section II E imply hardness-of-sampling results analogous to what is known for random unitary circuits [22].Such results would both complement the prediction of a computational-complexity transition in monitored random circuits [64] and lend theoretical support to proposals [65,66] for diagnosing monitored dynamical phases that avoid post-selection by computing cross-entropies instead.
Shortly after this work was completed, related results appeared in two papers [67,68].
collapse and dynamical purification C. Mapping to truncated unitary and Ginibre ensembles D. Computation of Lyapunov exponents and the purification time E. The distribution of Born probabilities F. Dynamics of Rényi entropies 1.Short times 2. Long times arXiv:2312.09216v2[quant-ph] 5 May 2024 FIG. 1.A schematic illustration of how the Haar-random monitored quantum circuits that we consider in this work (right) differ from more standard spatially local examples (left).The legs at the bottom of each picture correspond to individual qubits, blue rectangles depict Haar random unitary gates acting on the Hilbert space of the incoming legs, and red dots indicate single-qubit projective measurements.

p
and the circuit depth t.The Born probability of measuring a string of measurement outcomes m after time t can be written explicitly as p(m) = ∥ Pmt Ût . . .Pm 1 Û1 |ψ⟩∥ 2 .Introducing rotated projection operators P ′ m j = Û † 1 . . .Û † j Pm j Ûj . . .Û1 , we can write this as

FIG. 2 .
FIG. 2. The rank-collapse time for Model I (dashed red line) versus the purification time for Model II (solid blue line).It is clear from this plot that these timescales are exponentially well separated as L → ∞, except in the immediate vicinity of p = 1.

1 .FIG. 3 .
FIG.3.Numerical solutions to the singular integral equation Eq. (130) on intervals [−a(s), a(s)] with the endpoint position a(s) determined implicitly by number conservation.We set N = 50 and discretize the integral over one thousand points.The accuracy of our scheme is confirmed by its recovery of the semicircle law discussed in the previous section as s → 0, and we observe a clear semicircle-to-square crossover from a semicircle law to uniform behaviour between the very-short-time (s ≪ 0.01) and short-time (0.01 ≪ s ≪ 1) regimes defined by Eq. (137).
Let us now attempt to understand this crossover to uniformity analytically.Thus consider the uniform ansatz ρs.t.(w) = z) sinh (a + z) (a − z) sinh (a − z) .
bulk region.Comparison with the original integral equation (130) reveals that the uniform ansatz Eq. (