A Dyson Brownian motion model for weak measurements in chaotic quantum systems

We consider a toy model for the study of monitored dynamics in a many-body quantum systems. We study the stochastic Schrodinger equation resulting from the continuous monitoring with a rate $\Gamma$ of a random hermitian operator chosen at every time from the gaussian unitary ensemble (GUE). Due to invariance by unitary transformations, the dynamics of the eigenvalues $\{\lambda_\alpha\}_{\alpha=1}^n$ of the density matrix can be decoupled from that of the eigenvectors. Thus, stochastic equations are derived that exactly describe the dynamics of $\lambda$'s. We consider two regimes: in the presence of an extra dephasing term, which can be generated by imperfect quantum measurements, the density matrix has a stationary distribution, and we show that in the limit of large sizes the distribution of $\lambda$'s is described by an inverse Marchenko Pastur distribution. In the case of perfect measurements instead, purification eventually occurs and we focus on finite-time dynamics. In this case, remarkably, we find an exact solution for the joint probability distribution of $\lambda$'s at each time $t$ and for each size $n$. Two relevant regimes emerge: at small times $t\Gamma= O(1)$, the spectrum is in a Coulomb gas regime, with a well-defined continuous spectral distribution in the limit of $n\to\infty$. In that case, all moments of the density matrix become self-averaging and it is possible to characterize the entanglement spectrum exactly. In the limit of large times $t \Gamma = O(n)$ one enters instead a regime in which the eigenvalues are exponentially separated $\log(\lambda_\alpha/\lambda_\beta) = O(\Gamma t/n)$, but fluctuations $\sim O(\sqrt{\Gamma t/n})$ play an essential role. We are still able to characterize the asymptotic behaviors of entanglement entropy in this regime.


INTRODUCTION
In recent years, the study of quantum many-body systems in the presence of continuous monitoring of its local degrees of freedom has received much attention.The motivations for this case are of various kinds: first, from a practical point of view, the combinations of measurements and quantum gates are essential ingredients of quantum computation [1,2].In addition to this fact, recent studies have shown that quantum measurements (both projective and weak) can induce peculiar phase transitions, denominated measurement-induced phase transitions (MIPTs) [3][4][5], which arise exclusively by looking at the statistics of trajectories via nonlinear observables [6][7][8][9].Prime examples are transitions in entanglement dynamics [3][4][5] or purification time [10] induced by increasing the strength of the measurements (e.g., the probability with which each site is measured in the case of projective measurements or the rate of measurements in the case of weak measurements).In this case, it was found that going from weak to strong measurements, one can move from a volume-law entangled phase [11][12][13] to an area-law one [14][15][16][17][18].Results along these lines were initially obtained by using random circuits [19], in which quantum evolution occurs through unitary gates chosen from a uniform distribution (Haar) possibly within an appropriate subset of the unitary group.Then, in the limit of a large local Hilbert space dimension, it has been shown that a phase transition is present, and the corresponding critical point lies in the directed percolation universality class.However, it has been clarified that this simplification is not accurate in general [20]: the critical point is described by a peculiar conformal (scale invariant) theory not exactly solvable in general [21][22][23].Following these results, much theoretical [24,25] and numerical [26][27][28][29] work has been devoted to thoroughly characterizing MIPTs by linking them to various other phenomena, such as quantum error correction [30], purification [31][32][33], state preparation, quantum communication [34], and the complexity of classical simulation of quantum systems [35,36].Recent results have attempted to study this phenomenology in simplified toy models [37,38].However, in the case of noninteracting systems (e.g., free fermions), it has been shown that the transition has a radically different nature: the volume-law phase is immediately unstable for arbitrarily weak quantum measurements [39][40][41][42][43][44][45], and yet, in some cases, a subvolume phase can survive before the onset of the area-law phase [46][47][48][49][50], with a universality class more akin to problems of Anderson localization and disordered conductors, studied in the context of nonlinear sigma models [51,52].In the interacting case, one of the main difficulties in studying individual trajectories of monitored systems lies in the fact that the probability of each trajectory depends on the state itself, in accordance with Born's rule.This aspect produces an inevitable nonlinearity when considering statistical averages over measurement outcomes, similar to that faced in disordered systems.In analogy to that approach, it is possible to proceed through the replica trick, thus studying N identical copies of the system and eventually considering the relevant N → 1 limit (in contrast with the more common N → 0 limit relevant in disordered systems [53,54]).The fundamental ingredient for the replica limit is the possibility of performing an analytic continuation in N by means of an exact formula.Some studies avoid this problem by replacing the N → 1 limit with the simpler and more explicit N = 2 case, where in fact two copies are considered.Alternatively, recent results have been obtained using mean-field models [55][56][57][58], where analytical progress can be achieved.In particular, in [59], a discrete fully connected model was introduced in which the dynamics in the presence of measurements can be mapped onto a variant of branching Brownian motion and then studied in terms of the Fisher and Kolmogorov-Petrovsky-Piskunov reaction-diffusion Equation [60][61][62].However, the nature of these results in more general problems is not completely clear.More recently, a model of spin in the presence of noise and weak measurements has been introduced in which it is possible to study the dynamics in the presence of measurements analytically and in particular explicitly consider the N → 1 limit [63]: this showed that the replica limit can be highly nontrivial, with a MIPT present for any integer N > 1, but disappearing in the N → 1 limit.
Beyond the study of the transition, it becomes interesting to characterize the specific behavior of the two phases.In particular, the volume-law phase is hard to simulate from the classical point of view because of the strong presence of entanglement surviving very long times.In this paper, we focus on the characterization of the dynamics in the presence of weak measurements far from the critical point, within the volume-law phase.To do so, we introduce a model of the evolution of the density matrix of the system purely based on random matrices, akin to the one considered in Ref. [64].We choose a dynamics based on the stochastic Schrödinger Equation [65], invariant under unitary rotations in the n-dimensional Hilbert space, such that the dynamics of the eigenvalues of the density matrix decouples from that of the eigenvectors, as is the case in Dyson Brownian motion (DBM).In this way, unitary dynamics becomes inessential since it does not affect the eigenvalues.Once we have derived a stochastic equation for the evolution of the eigenvalues, we study two relevant limits of it: in the first case, we consider dynamics induced by "imperfect" measurements, i.e., a case in which a fraction of the measurement results is unknown.This is effectively equivalent to introducing a dephasing term that on each trajectory tends to bring the density matrix closer to the identity.In this case, we prove the existence of a steady state that describes the distribution of the eigenvalues: at large n, the resulting ensemble is described by an inverse-Wishart matrix.Next, we explore the dynamics for perfect measures: in this case, the steady state at long times reduces trivially to a pure random state uniformly distributed in the steady state; however, this happens with a characteristic transient at small and intermediate times.We present the exact solution for the joint distribution of eigenvalues at all times and analyze its effects on the entanglement entropy.

PRELIMINARIES ON UNRAVELINGS AND TRAJECTORIES
The combined dynamics of a quantum system undergoing both unitary evolution and measurements can always be modeled as a quantum channel [66], whose output contains both the state of the quantum system and the classical encoding of the results of the measurements.According to Choi's theorem [67], a quantum channel Φ(ρ) can always be expressed in terms of a set of Kraus operators {K a } a , first discussed in Ref. [68].Those provide an explicit operator-sum decomposition of the quantum channel as Φ(ρ) = a K a ρK † a , and describe the dynamics by means of linear maps from the system's Hilbert space to itself.Because of trace-norm conservation, Kraus operators satisfy the condition a K † a K a = 1, with 1 the identity operator in the Hilbert space of the quantum system.However, the specific decomposition of Φ in terms of the Kraus operators is not unique as {K ′ a } a with K ′ a = b U a,b K b , for an arbitrary unitary transformation U a,b , define the same quantum channel Φ(ρ).Different decompositions of the same quantum channel are referred to as unravelings.In the context of repeated measurement, the Kraus operators are factorized K a = K a1 K a2 . . .K a T where a = (a 1 , a 2 , . . ., a T ) labels the collection of all measurement outcomes performed at each time step.For MIPTs, one is interested in the single trajectory where the initial density matrix ρ(0) evolves as where the probability of a specific trajectory has been introduced according to Born's rule as p a = tr[ρ a ].Given any functional of the state F [ρ], we define the average over trajectories as If the results of the measurements are not known, one only has access to linear functionals of the state, such as the quantum expectation of any observable, e.g., ⟨tr[ Ôρ]⟩ = tr[ ÔΦ[ρ(0)]], which depends solely on the quantum channel Φ[ρ] and is independent on the specific unraveling.However, this is not true for nonlinear functionals such as the Renyi's entropies, where one considers the functionals and in particular, the Von Neumann entropy defined as the limit k → 1, or, explicitly, where {λ α } n α=1 denotes the eigenvalues of ρ.Quantities like ⟨S 1 [ρ]⟩ can be used as order parameters for MIPTs.By discarding the Born's weight for trajectories, one can define a different average where all measurement outcomes have equal weight and are statistically independent at different times.Unbiased outcomes can be obtained via post-selection by artificially attributing the same probability to each trajectory.In some contexts, ⟨. ..⟩ 0 is used as an approximation for the more physical ⟨. ..⟩.We will see that both cases emerge naturally in our framework.

THE MODEL
We consider a toy model for continuous monitoring in a quantum system with a Hilbert space of dimension n.Let us explain the model by splitting the unitary evolution and the monitoring parts.For the unitary evolution, we have the Hamiltonian increment where hermiticity is ensured requiring that dh αβ = dh * βα .Otherwise, the increments are chosen as a Hermitian Brownian motion with covariance dh αβ dh * γδ = δ αγ δ βδ dt/n.Equivalently, in the limit dt → 0, we can write dh αβ = √ dt h αβ , where h αβ is an n × n Hermitian matrix, drawn from the Gaussian unitary ensemble (GUE) and the measure is over the independent entries of a Hermitian matrix (with a semi-circle spectrum of support [−2, 2] in the large n limit).A crucial property of the GUE distribution is its invariance under unitary conjugation, P (h) = P (uhu † ) for any unitary u ∈ U (n).The dynamics induced by such a unitary evolution is simply described by ρ + dρ = e −idH ρe idH , where the Itô calculus conventions are implied.For measurements, let us first introduce the standard form of the stochastic Schrödinger Equation [65,69,70] for an imperfect monitoring of a Hermitian operator Ô.Within a small time step ∆t, the system is evolved coupling the operator Ô with an extra spin 1/2 (ancilla) whose value along a reference direction is projectively measured with an outcome a = ±1.As explained in the previous section, this evolution can be encoded into the Kraus operators [63] for a summary of this derivation), where Γ is a rate quantifying the measurement strength.In the limit ∆t → 0, the collection of measurements a can be used to define a standard Wiener process Y (t) (dY 2 = dt) and one obtains the stochastic evolution equation with the dephasing superoperator The parameter x ≥ 0 quantifies an extra source of dephasing.It can be seen as due to a fraction of measurements whose outcomes are not known or more generally to a coupling with an external dephasing bath.The limit x = 0 corresponds to perfect measurements.In our model, within each infinitesimal time step dt, we choose Ô(t) = n α,β=1 o αβ (t) |α⟩ ⟨β| to be a random operator, with components o αβ drawn from the GUE distribution (7).We can thus introduce a Hermitian Brownian motion setting do αβ = o αβ (t)dY , satisfying do αβ do * γδ = δ αγ δ βδ dt/n.Using this, according to the rules of stochastic calculus, we can rewrite the dephasing part as do αβ |α⟩ ⟨β|.Note that at finite ∆t, o αβ (t)dY is the product of two Gaussian distributions that is not Gaussian.However, this is irrelevant in the ∆t → 0 limit, where only the covariance is relevant to define a Wiener process.
In the following, we will assume that the density matrix is initially prepared in the infinite-temperature state ρ = 1/n.Thus, because both the Hamiltonian increment and the observable Ô are always chosen from the GUE, the distribution of the density matrix at any time is itself invariant under unitary transformation.In other words, diagonalizing ρ = uΛu † , u will be Haar distributed in the unitary group U (n) and will completely decouple from the dynamics of Λ = diag(λ 1 , . . ., λ n ).In the following, we will focus on the dynamics and the distribution of the eigenvalues {λ α } n α=1 .As they are unaffected by the unitary dynamics, in the following we will ignore the latter.
Before dealing with the general n case, it is worth considering specifically the n = 2 case as an instructive warm-up.We can parameterize the density matrix in terms of a vector ⃗ r = {r 1 , r 2 , r 3 }, within the Bloch sphere (|⃗ r| ≤ 1), as where ⃗ σ = {σ 1 , σ This can be recast into a closed stochastic equation for the modulus r = |⃗ r| (as required by unitary invariance), which takes the form dY being the standard Wiener process satisfying dY 2 = dt defined before Equation (8).In order to recast the previous equation into the standard Langevin form, let us perform the change in variables r = tanh(ω), ω ∈ R + .In terms of ω, we indeed have dω = −ΓV ′ (ω)dt + √ ΓdY , where the potential is given by the following: As a consequence, the evolution of the probability distribution P (ω, t) can be described through the associated Fokker-Planck (FP) equation which admits a stationary probability distribution P stat (ω) = 1 Z e −2V (ω) .Coming back to the original coordinate, one has a result already derived in [63] in the context of mean-field approximation.Since tr(ρ 2 ) = (1 + r 2 )/2, it follows that the purity has a nontrivial stationary distribution as well for any x > 0. As x → 0, the distribution becomes more and more peaked around r = 1, eventually collapsing to P stat (r) = δ(r − 1) for x = 0: this is consistent with the fact that perfect measurements eventually lead to purification.For x = 0, the finite-time solution of Equation ( 13) can be explicitly worked out as follows: normalized in ω ∈ R + .The origin of this form will be clarified for general n in Section .We can use these results to compute the Von Neumann entropy (4).In particular, on each trajectory, it can be related to the variable r with S max = log n = log 2, the entropy of the maximally mixed state.We can compute the behavior of ⟨S 1 ⟩ in different regimes; in the stationary state at x > 0, one has For x = 0, the system purifies, and ⟨S 1 ⟩ → 0. For large times, we have See Appendix for the details of the derivation and for the higher orders.

DYNAMICS OF THE SPECTRUM Stochastic Evolution of the Eigenvalues
In order to derive an evolution equation for the eigenvalues in the general case, we follow the standard approach [71] and make use of second-order perturbation theory, as follows: Now, one has to substitute the infinitesimal variations ρ α as determined in Equation (8).To further simplify the calculation, we take advantage of unitary rotational invariance, so that ρ at time t is assumed to be diagonal in an appropriate basis ρ αβ = λ α δ αβ (no sum over α).Thus, we have whereas for the off-diagonal elements, we simply have dρ α̸ =β = √ Γdo αβ (λ α + λ β ).We can consider the diagonal elements of d Ô and rescale them as dB α = √ n do αα , so that dB α dB β = δ αβ dt are standard Wiener increments.Then, one finally finds where we rescaled γ = Γ/n.Let us notice that, as it should be, the evolution of Equation ( 21) preserves the trace, as α dλ α = d tr(ρ) = 0. Additionally, for x = 0, the configuration λ α = 1, λ β̸ =α = 0 is a fixed point, since as expected, in the absence of any dephasing (x = 0), measurements eventually lead to purification and the density matrix reduces to a rank-1 projector.It is useful to further manipulate Equation ( 22) by rewriting λ α − 1/n = 1 n β̸ =α (λ α − λ β ), where we used β λ β = 1.In such a way, we obtain For x = 0, Equation ( 22) presents some analogy with the eigenvalue flow studied in [72], see Equation (28) therein with (1 + mλ i )dt → 0, σ 2 = 4γ and β = 2.The sum over β ̸ = α is clearly a manifestation of the typical Coulombic repulsion between eigenvalues.On top of that, the form of the noise is non-diagonal, as a combined effect of exactly preserving the trace and attributing Born's probability to trajectories.Then, in order to simplify the dynamics, we introduce a new set of unconstrained variables.

Mapping to Unconstrained Variables
As the dynamics of the eigenvalues exactly preserves the traces, it is useful to write them as Of course, the mapping between the λ's and the y's is not one-to-one, as a global rescaling of all the y's does not affect the mapping.This freedom can be used to simplify the evolution of the y's, in particular, by writing we can look for a force term F α (⃗ y) which is homogeneous of degree 1.From standard Itô calculus, one can verify that where g(⃗ y) is an arbitrary homogeneous function (g(c⃗ y) = g(⃗ y)), which correctly reproduces Equation ( 22) for the evolution of the λ's.Equation ( 23) has the advantage of involving only a diagonal noise term, which is nonetheless multiplicative.We can further simplify its form setting w α = log y α , which leads to where Fα ( ⃗ w) = e −wα F α (e w ).Setting g(⃗ y) = (n − 1)/2, one can express where whereas the non-conservative force is given by f α = β e w β −wα .We will now explain how the dynamics induced by Equation ( 25) can be solved in different regimes.
We now consider the dynamics induced by imperfect measurements at finite x = O(1) in the large n limit.When the dephasing term dominates the right-hand side of Equation ( 22), the eigenvalue distribution is expected to be peaked around λ = 1/n.Therefore, it is convenient to rescale the eigenvalues via with the trace condition becoming α λα = 2n/x.With this scaling, we can analyze the noise term in Equation ( 22), which is proportional to λα (dB α −x/(2n) β λβ dB β ) ∼ λα dB α .Indeed, by summing the last term in quadrature, we see it scales like n −1/2 and can thus be neglected for n ≫ 1.This can be understood as the sum α λα being extensive and in the n → ∞ limit its fluctuations are subleading with respect to its average.In this scaling, Equation ( 22) becomes This associated stochastic dynamics is exactly solvable even at finite n.Indeed, by introducing again logarithmic variables wα = log λα , we can reduce this to a Langevin equation, analogous to Equation ( 25), as follows: We stress that this equation involves a completely conservative force where the potential V x ( ⃗ w) now reads as follows: It is helpful to comment between the relation of the variables wα introduced here and w α in Section .We can see that they match, up to an inessential shift of the center of mass.The initial value of the sum α e wα can be chosen arbitrarily and in the limit n → ∞ it becomes a conserved quantity.Since we implicitly chose α e wα = α λα = 2n/x, we can rewrite λ α = e wα / α e wα coherently with Section .Thanks to the Langevin form of Equation ( 30), we can easily deduce the joint distribution which is the stationary solution to the associated FP equation.In terms of the rescaled eigenvalues λα which is, upon rescaling λα → λα /2n = λ α /x, the joint distribution function for the eigenvalues of a matrix in the inverse-Wishart ensemble, with m = n(1 + x/2).We refer to Appendix for a brief discussion of such an ensemble.
Let us now introduce the one-particle density function For matrices in the inverse-Wishart ensemble, the eigenvalue density takes the inverse-Marchenko-Pastur (IMP) expression reported in Equation (123) of Appendix .It is thus immediate to verify that the distribution for the rescaled eigenvalues defined above reads and it is shown in Figure 1.The endpoints λ± are given by The distribution in Equation ( 35) is valid in the n → ∞ limit which is also needed for Equation (29) to become a valid approximation of the original Equation (22).The behavior of f IMP ( λ, x) as x is varied and defines a sharply peaked distribution around 2/x for x ≫ 1, as the interval λ+ − λ− = 16 x 2 1 + x 2 collapses to a point.Conversely, in the limit of nearly perfect measurements x ≪ 1, the rescaled eigenvalues spread on the positive semi-axis between λ− → 1/4 and λ+ → ∞, whereas f ( λ) remains peaked around λ = 1/3.As this corresponds to λ → 0, this signals the purification dynamics induced by perfect measurements.For any finite n, it is important to notice, however, that even in this regime, one has to assume x ≫ n −1 for our approximation to be valid, so that the limits x → 0, n → ∞ do not commute.The Von Neumann entropy, defined in Equation ( 4), can be computed in the stationary state, in terms of the inverse-Marchenko-Pastur eigenvalue distribution above, as follows: with S max = log n the entropy of the maximally mixed state.The last integral on the right-hand side can be evaluated by means of contour integration so that one has the remarkably simple exact formula, as follows: In particular, this allows us to evaluate its asymptotic behavior in the opposite regimes x ≫ 1 and n −1 ≪ x ≪ 1, as follows:

Finite Time Dynamics
The finite time evolution of the set of unconstrained { wα } n α=1 used in this section has been studied in detail by one of the authors in Ref. [72].Our Langevin equation (30) indeed reduces to Equation (42) of Ref. [72] upon the identification t → Γxt, m = −1, σ = 4/xn, β = 2. Consistent with our result in Equation ( 35), a stationary inverse-Wishart distribution was also obtained there.Here, we briefly summarize the procedure and the results with our notation and parameters.The associated FP equation can be mapped to an imaginary time quantum problem via the transformation with P stat ( ⃗ w) given in Equation ( 32) (we better discuss the quantum mapping from FP to Schrödinger in Section ).The ensuing Schrödinger equation reads with the Hermitian Hamiltonian and the finite energy shift E 0 , which is the n-particle ground state of H.As the trajectories of the n particles are non-crossing, the problem is mapped onto fermions.The one-particle potential U ( wα ) is of Morse type: The single-particle spectrum of the Morse potential features a finite number of bound states ψ k ( w) at eigenenergy so that the lowest eigenvalue is ε 0 (n) = 0.The single-particle finite time evolution from an initial w 0 to a final w coordinate is determined by the Euclidean propagator G M using the spectral decomposition over the bound and the scattering states Then, the n-fermion propagator is obtained in terms of the Karlin-McGregor determinant formula [73] for the probability evolution of n non-intersecting particles from ⃗ w0 to ⃗ w.From the quantum Euclidean propagator, one can write the solution to the original FP problem using Equation (41) as given the initial condition P ( ⃗ w, t = 0) = α δ( wα − w0 α ).

THE PERFECT MEASUREMENT DYNAMICS
In this section, we present the full solution of the joint eigenvalue distribution for the measurement problem at arbitrary n and time t, setting x = 0.

Exact Solution at Finite Time
The stochastic Equation (25) takes the Langevin form with the potential in Equation (27).The joint distribution P ( ⃗ w, t) for the variables w's satisfies the FP equation Formally, one can obtain a stationary solution to the FP equation, taking P stat = e −2V .Of course, this is not normalizable, consistent with the fact that at t → ∞ the density matrix simply purifies.Nonetheless, it can once again be used to convert the FP into a quantum problem, similarly to what was conducted in Section with Equation ( 41).One has which does satisfy a Schrödinger evolution with Hamiltonian operator Remarkably, after using some algebra (see Appendix and also [74,75]), one can show that implying that the Schrödinger evolution in Equation ( 50) amounts to free diffusion with the additional non-crossing constraint when two w's collide, which maps on free fermions, as the external Morse potential Equation (44) present in [72], see Equation (47), is absent here.The propagator for such dynamics can be obtained using the Karlin-McGregor determinant formula [73], also used in Equation (47).For a given initial condition P ( ⃗ w, t = 0) = α δ(w α − w 0 α ), we can obtain the solution which leads to In this expression, we introduced the free diffusion propagator, namely G 0 (w, w ′ , τ ]/ √ 4πτ .In the following, we focus on the case of the infinite-temperature initial state, where all the initial conditions coincide ⃗ w 0 → 0 (the specific value is inessential).In that limit, Equation (54) further simplifies as P ( ⃗ w, t) → 1 Zt P 0 ( ⃗ w, t) ( α e wα ), where and the time dependent constants Z 0 t , Z t are both determined by enforcing the normalization of P 0 ( ⃗ w, t) and P ( ⃗ w, t), respectively.

Relation between the Two Averages
The distribution P 0 ( ⃗ w, t) has appeared in [75,76] in studying the dynamics of the so-called isotropic Brownian motion in the context of May-Wigner instability.In contrast, the expression for P ( ⃗ w, t) has the extra factor α e wα , whose origin can be traced back to Born's rule, expressing (up to an irrelevant normalization) the probability in Equation (1) as p a ∝ α e wα .It is useful to identify the where ρ is the unnormalized density matrix introduced in Equation ( 1).In analogy with other problems of a multiplicative process involving random matrices [77][78][79][80], we will refer to the variables w α as the Lyapunov exponents.In the following, consistent with the definition in Equation ( 5), we will add a subscript 0 to the averages computed with P 0 , with the general relation for an arbitrary function F ( ⃗ w).For instance, for the Renyi's entropies we can write the following: Then, using Equation ( 56), we can rewrite the Von Neumann entropy as In particular, introducing the moments of the eigenvalues and of the trace as we can express Equation ( 60) as Thus, below we will study these moments with the measure P 0 .We observe that in Equation ( 55), one can recognize two Vandermonde determinants, since where we set δ k = (n + 1)/2 − k.This implies that P 0 describes a determinantal point process, which allows us to obtain several exact results, as discussed below.

Average of Schur's Polynomials
The calculation of several quantities, including Renyi's entropies (3), requires the expressions of correlation functions of the Lyapunov exponents w's.Here, we explain how they can be expressed systematically.First, because of the symmetry under exchange in the w's, we can restrict to symmetric functions.Then, because of the determinantal structure in Equation (), following [81] a complete and useful basis of symmetric functions is given by Schur's polynomials and each correlation function can be expressed once the expectation of Schur's polynomials is known.To a partition κ = (κ 1 , . . ., κ n ) of the integer m = j κ j , with κ 1 ≥ κ 2 ≥ . . .≥ κ n ≥ 0, one associates the corresponding Schur polynomial via [82] s κ (y) = det(y Setting y α = e wα and denoting h j = κ j + n − j, we can now express the average where C is the normalization and the constant A t raised to the power m accounts for the shift of the center of mass and will be fixed below.We can use Andreief identity [83] to express it in terms of a single determinant where we defined We can thus express the coefficients I k,h in terms of the Hermite polynomials H n (x) = (−1) n e x 2 ∂ n x [e −x 2 ] as where again we absorbed some extra constants by redefining C. The fact that at large x, H ℓ (x) = 2 ℓ x ℓ + O(x ℓ−1 ) can be used to express the determinant This last determinant is once again a Vandermonde one which can be expressed via (64).We can now plug this back into Equation ( 67) and fix the constant C using s κ=0 (y) = 1, where h j → n − j.We finally obtain where we recognized the equality which expresses the number of semistandard Young diagram of shape κ and n entries [82].We can now use n j=1 h 2 j − (j − 1) 2 = (2n − 1)m + 2ν(κ) [82], with where κ ′ = (κ ′ 1 , κ ′ 2 , . ..) denotes the partition dual to κ, with κ ′ i = #{κ j |κ j ≥ i}.The above average can then be expressed as follows: Power-Law Symmetric Polynomials For later convenience, we also introduce the power-law symmetric polynomials and additionally, for any integer partition µ = (µ 1 , . . ., µ n ), we can set which form a complete linear basis for symmetric polynomials.For instance, we can express Since Schur's polynomials are also a complete basis, one can find a change in basis between the two.It is given as (see Equation (3.10) in [81,84]) where χ κ µ represents a character of the symmetric group S m , with m = α κ α , on the irreducible representation and the conjugacy class labeled, respectively, by the integer partitions κ and µ.

Calculation of the Moments
Thanks to the relation between power-law symmetric polynomials and Schur's ones, given in Equation ( 79), we can use the previous result to express the moments of eigenvalues, c.f. Equation (61) with y α = e wα : which corresponds to Equation (79) in the case µ = (m), i.e., a single cycle of length m.In this case, the only nonvanishing characters χ κ (m) are those related to L-shaped partitions κ = (m − r, 1 r ) = (κ 1 = m − r, κ 2 = 1, . . ., κ r+1 = 1), with r = 0, . . ., m − 1.In particular, one has χ = (−1) r .Therefore, from Equation ( 74), we obtain 2ν(κ) = m(m − 2r − 1), and The sum in Equation ( 80) thus becomes which can be expressed for integer n in terms of the hypergeometric function 2 F 1 (a, b; c; z) [85], obtaining the following: In both expressions, we have fixed the value of A t = e −2γt(n−1) using ∂ m M (0) = ⟨ α w α ⟩ 0 = 0 in our conventions.Note that M (0) = n as required, and that the lowest moments have simple expressions, e.g., one has M (1) = ⟨ α e wα ⟩ 0 = ne 2nγt .Moreover, it is also possible to compute the moments of the trace defined in Equation ( 62), with y α = e wα , by making use of the same relation.We can express where the sum is over the partition κ of the integer m.

Equivalent Formulations
Equation ( 55) can appear in different contexts that provide interesting interpretations.Following [86][87][88], and using ( 64) and (65), at a fixed time t, one can identify the Lyapunov exponents w's with the spectrum of the matrix where H is drawn from the GUE distribution (7) and D = diag( n−1 2 , n−3 2 , . . ., − n−1 2 ) = diag(δ 1 , . . ., δ n ), see also [91].The scaled distribution of the largest eigenvalue of W was computed at large n in [89], see Section III there, and found to be GUE Tracy-Widom, since no localization transition takes place when the elements of D are equispaced; see also [90].Besides, equation ( 85) is particularly effective for numerical sampling from the distribution Equation (55).
Equivalently, the w's can be seen as performing a DBM with β = 2 in inverse time ∼ 1/t.For time inversion of the DBM we refer to, e.g., Appendix B in [89] with equally spaced initial condition.Indeed, one can rewrite (85) as with s = 1/(4γtn), and the eigenvalues ⃗ x(s) of X have the same joint law as a DBM at time s with initial condition ⃗ x(0) = ⃗ δ/n.This correspondence is given in Appendix .Note that the ⃗ x(s) form a determinantal point process [88], whose kernel is given in Appendix .The interpretation of Equation ( 85) already indicates a qualitative behavior for the spectrum of the w's: at initial times, the first term dominates and the distribution will resemble the one of a GUE in with eigenvalues in the support [−4 √ tΓ, 4 √ tΓ].At larger times, the second term becomes more important and the distribution becomes uniformly spread in [−2tΓ, 2tΓ].It is also useful to write the distribution P 0 ∼ e −n 2 Et(f ) as a functional of the single particle density with the functional One can interpret this as the energy of gas particles in the presence of an external harmonic trap that repel each other.When Γt = O(1), the different terms are of the same order: at large n, the particles will have a separation ∝ tγ/n and a continuous description emerges.We will analyze this regime in the next section.We will then discuss the long-time regime which emerges when t ∼ n/Γ = 1/γ.

Coulomb Gas Regime Γt ∼ O(1)
To quantitatively analyze the crossover between the semicircle and the uniform distribution predicted by the form Equation (85), we rescale the time in the following way and consider the limit n → ∞ at fixed τ .The above correspondence in (86) shows that s = 1/τ , and W = τ X.The interpretation as a DBM in inverse time allows us to compute the limiting resolvent g(z, τ ) := lim using the complex Burgers equation (see details in Appendix ).One finds that g = g(z, τ ) satisfies the parametric equation We can set z = w − i0 + and g = g r + iπf , where f is the spectral density.Eliminating the real part g r = Re(g) one finds (see Appendix ) that the density f = f τ (w) is the solution of Its support is an interval with edges at w = ±w e = ± arccosh(1 As anticipated from the qualitative analysis of Equation ( 85), the density interpolates (see Figure 2) between a semi-circle at small τ , with 12 , and a square distribution at large τ with w e = τ 2 + log(eτ ) + 1 τ + O(1/τ 2 ).At all times τ > 0 the density vanishes at the edges as a square root f (w) ∼ B τ |w e − w| 1/2 , with B τ = √ 2/[πτ 3/4 (4 + τ ) 1/4 ], and Wigner-Dyson gap statistics near the edges.Although the density is known in parametric form, one can also use Equation ( 82) to extract the exponential moments.After some manipulation on the hypergeometric function, we find m (x) denotes the generalized Laguerre polynomial, with lim m→0 µ(m) = 1.A numerical check confirms that the density implicitly defined in Equation ( 92) satisfies The fact that the variables become dense and described by a continuous distribution ensures that Ω(1)/n → µ(1) in Equation ( 62) is self-averaging and thus Universal Regime Γt = O(n)

Scaling of the Edge
Now, we consider the situation where time is scaled as t = O(n/Γ).Using the rescaled rate γ, we can equivalently say that γt = O(1) whereas n → ∞.In this regime the repulsion due to the log sinh term in becomes more relevant.As a result, the Lyapunov exponents w are well separated one from the other.More precisely, assuming the ordering w 1 < w 2 < . . .< w n , one expects that the separation between two consecutive variables w α+1 − w α = O(γt), so that the position of the edge we = ⟨w n ⟩ 0 = O(tγn).However, each variable w α will have fluctuations of O( √ γt) so that the interactions are relevant and complicate the analysis.A similar lattice structure was already demonstrated in a related model involving a log sinh-potential [92].Here, the Coulombic repulsion log |w − w ′ | is an additional ingredient that modifies that short-distance behavior.We can obtain a preliminary understanding of this regime by taking the n → ∞ limit from the general formula of moments (83).In this limit, the moments are dominated by the largest w α ∼ we near the edge.We can determine their value expanding Γ(n + m − r)/Γ(n − r) ≃ n m in Equations ( 81) and ( 82) and performing the sum, obtaining where we also fixed we = 2nγt + log n.This formula suggests the existence of non-trivial statistics governing the fluctuations of the Lyapunov exponents w's at the edge.However, we note that it does not define a normalized distribution since M(m → 0) = ∞.This is due to a crossover at small m between a regime dominated by the bulk with M (m → 0) = n and one dominated by the edge M (m) n≫1 = O( wm e ).An identical formula can be obtained in the analogous long-time limit considering the product of Wishart random matrices (see Equation (3.41) in [93]).This is a manifestation of universality expected at large times in monitored systems and more generally in products of random matrices [94].
The calculation of the moments of the trace in this regime is much more involved, but it can be expressed as in Equation ( 84).Indeed, in the limit of large n, we can replace s κ (1) → n m χ κ 1 m /m! [82].One can understand this using s κ (y 1 = 1, . . ., y n = 1) counts the number of semistandard Young tableau of shape κ and involving n entries, whereas χ κ 1 m counts the number of standard Young tableau with entries 1, . . ., m.At large n, the difference between semistandard and standard Young tableau becomes irrelevant and s κ (1)/χ κ 1 m ∼ n m ∼ n m /m!.We thus obtain which was obtained in [94] in a different model as an additional manifestation of the universality.

Asymptotics at Large γt
When γt ≫ 1, the separations are |w α − w β | ≫ 1 ∀α, β.Additionally, their fluctuations are much smaller than their separation.In this regime, we can approximate in the potential which is the 1-dimensional Coulombic repulsion.This kind of potential was recently studied in [95] in the context of ranked diffusion.This potential induces a force f α = 1 2 β̸ =α sign(w α − w β ), which simply counts the number of slower particles in the back of w α , minus the faster ones in its front.Assuming a given ordering w 1 ≪ w 2 ≪ . . .≪ w n , we obtain the simple dynamics where the drift velocity v α := 1 2 (2α − n − 1).This equation is an example of ranked diffusion (RD), namely, a diffusion process in one dimension where the n particles undergo a drift term proportional to their respective rank [95,96].In the RD problem, however, the particles can cross freely, whereas here they cannot cross, which leads to different types of fluctuations at finite t.Nevertheless, in both models in the regime γt ≫ 1, the equations decouple and can be solved separately, namely where we recall that the B α (t)'s are n independent standard Brownian motions, each of variance t at time t.Equation (101) gives some characterization of the joint probability distribution of the Lyapunov exponents w in this large time γt ≫ 1 regime.It is not complete however, as there are O(1) contributions of the joint cumulants of the w α which persist at infinite time.These were computed in Section III-D of Ref. [95] by a saddle point method, and that calculation is easily extended to the present model in Appendix .As a result, Equation ( 101) must be treated with care when computing exponential moments.

ENTANGLEMENT ENTROPIES FOR CONTINUOUS MONITORING
Now, we make use of the results obtained in the previous section for the unbiased ensemble ⟨. ..⟩ 0 to characterize the behavior of the entanglement entropies.

Short Time Regime
First of all, let us consider the Coulomb gas regime.In this case, we can use the fact that the moments are self-averaging, i.e., n −1 tr ρk in law → µ(k).Thus, from Equation ( 58), we obtain where the Von Neumann entropy can be recovered in the limit k → 1.Using L (1) ) for τ ≫ 1, we also obtain the asymptotic expansions from which we can extract the k → 1 limit where γ E is the Euler-Mascheroni constant (not to be confused with the rescaled rate γ).Note that in this regime, the Born rule (2) and unbiased outcomes (5) give the same result, since the weight of each trajectory factorizes from the rest.

Universal Regime
As we discussed in Section , this regime is much harder to address as the w's are strongly correlated but the moments M (m) and Ω(m) are not independent: they are dominated by the behavior around the edge (see Equations ( 97) and ( 98)).However, we can still access the Von Neumann entropy by making use of Equation (63).Equation ( 97) admits a simple analytic continuation for m → 1 and we obtain Unfortunately, the dependence on m in Equation ( 98) is much less transparent.A way to perform the analytic continuation was developed in [94].Here, we analyze the asymptotic behavior.At small γt's, because of the symmetry ν(κ ′ ) = ν(−κ), where κ ′ is the integer partition dual to κ, it is clear that Ω(m) is an even function of time for every m .Thus, at small γt ≪ 1, we can conclude that Ω(m) = O(γt) 2 .Up to this order, the dynamics of Von Neumann entanglement is fully captured by Equation (105).Thus, we have which connects nicely with the behavior for τ ≫ 1 obtained in Equation ( 104), recalling that τ = 4γnt.Although the regime of intermediate times is hard to access, we can still estimate the asymptotic expansion γt ≫ 1.As discussed in Section , the Lyapunov exponents separate linearly in time because of the drift Equation (101).Returning to the eigenvalues λ α = e wα / β e w β , it is evident that if w n ≫ w β̸ =n , then λ α ∼ e wα−wn ∼ δ α,n , and the first n − 1 eigenvalues are suppressed by purification.In particular, using Equation ( 101), we also have that where v α are the drift velocities in Equation (100) and ⟨. . .⟩ B denotes the average over the independent Brownians.The sum of the first n − 1 averages gives α<n ⟨λ α (t)⟩ ∼ e −4γt , yielding Moreover, the ratio between any two eigenvalues ⟨λ α̸ =n (t)⟩, ⟨λ β̸ =n (t)⟩ which is independent of number n of diffusing particles, and indicates that the first n−1 eigenvalues are logarithmically equispaced, as we show in Figure 3.Because of the exponential separation of the eigenvalues, from Equation (108), we obtain the estimate for the Von Neumann entropy which is in agreement with the n = 2 result of Equation ( 18) if one replaces Γ = 2γ.A more detailed discussion of this approximation is given in Appendix .
Although this suggests that for any n only the two largest eigenvalues matter for the calculations, the prefactor of the asymptotic expression of ⟨S 1 ⟩ does not reduce to the one from the n = 2 model.Indeed, a more careful calculation (see Appendix ) shows that the saddle-point evaluation of the numerator of ⟨S 1 ⟩ is shifted with respect to the one of the normalization, so that the two prefactors do not cancel out.In particular, we have which once again agrees with the result of the n = 2 case, whereas the prefactor (n − 1)/n tends to 1 in the n → ∞ limit.

CONCLUSIONS
In this work, we studied the purification dynamics induced by weak measurement, where the measurement operators are random matrices drawn from the GUE.The ensuing is written in terms of a stochastic Schrödinger equation, invariant under unitary transformation.Because of rotational invariance, we were able to determine the evolution of the density matrix eigenvalues {λ α } n α=1 , both for imperfect and perfect monitoring.In both cases, no MIPT takes place, as the system is far away from the critical point and is characterized by volume-law scaling of the entanglement entropy.
In the former case of imperfect measurements x > 0, we were able to thoroughly characterize the stationary distribution of the density matrix in the limit n → ∞.There, the stationary joint probability distribution of the eigenvalues of the density matrix takes the typical expression of matrices from the inverse-Wishart ensemble, and the spectral density attains the corresponding inverse-Marchenko-Pastur form.Anyway, we showed that a solution for the finite-time stochastic dynamics also exists via a mapping to a quantum problem in imaginary time.In this context, the eigenvalue evolution at x > 0 and for finite n is still unknown.The relevant n = 2 case seems to suggest that the stationary state does not belong to the inverse-Wishart ensemble, and that a more sophisticated stratagem must be sought in order to treat the non-conservative forces due to additional dephasing.
In the case of perfect measurements x = 0, we were able to fully determine the joint distribution function for the eigenvalues, at both finite time and finite n.Again, this can be conducted by mapping the stochastic dynamics onto fermions in imaginary time.Surprisingly, those undergo free diffusion as the external potential due to weak measurement disappears.Moreover, this allows us to factorize the joint probability distribution with respect to the tr ρ = 1 constraint, thus making the computation of average quantities accessible.We identify two different regimes.The first is a short-time one, where a Coulomb gas distribution appears with a time-dependent continuous density characterizing the spectrum of Lyapunov exponents.The density exhibits a characteristic crossover between the GUE semicircle and the uniform distribution.At later times, the Lyapunov exponents separate further and a continuous description is not possible anymore.Our calculations are perfectly consistent with the prediction of [94] and provide a new example supporting the universality of this regime.Subsequently, the Lyapunov exponents separate linearly over time, and the interactions between them become less important.In this regime, the dynamics resemble that of a ranked diffusion.We use this information to characterize the dynamics of entanglement in the various time regimes.Some comments are important.It would certainly be intriguing to characterize the crossover at x → 0 small and large times.Furthermore, from a technical point of view, it would be interesting to tie the regime emerging at long times with the spectral distribution of Lyapunov exponents near the edge.
Finally, interesting questions endure, characterizing the differences induced by studying the problem within other random matrix models, such as the Gaussian orthogonal ensemble.
Note added.Recently, Ref. [97] appeared, where a similar approach is independently studied in the case of unbiased outcomes.Equation ( 55) is also obtained and the Coulomb gas regime is studied, also identifying the crossover between a GUE and a uniform distribution.

Identities
To perform the calculation in the text we use the following two identities: and The first one is well known from Calogero's papers (see references and generalizations in, e.g., Appendix A of [74]) and the second is specific to the present case.We will now derive an expression for ⟨S 1 ⟩ (t) assuming n = 2.In this case, the finite-time probability density is given by Equation ( 15), whereas we have Replacing the explicit form of P (ω, t) we obtain (setting Γ = 1 for simplicity) Expanding for large times, and considering erf(x) = e −x 2 / √ 2πx 1 − x −2 /2 + O(x −4 ) , we have hence making use of the identity we find the result (18) present in the main text.

Inverse-Wishart Ensemble
The unitary (β = 2) Wishart ensemble is defined by the distribution for any squared n × n matrix A, with integer m ≥ n [98,99].A Wishart-distributed matrix A can be written as A = HH † , with H a rectangular n × m matrix with complex Gaussian entries.Equation (118) leads to the following joint distribution for the eigenvalues {a i } n i=1 , as follows: The marginal distribution of the eigenvalues, i.e., the spectral density, can be obtained from the previous equation.
In the limit of infinitely large matrices n → ∞, a standard result is the Marchenko-Pastur distribution with the rescaled eigenvalues ã = a/2n.The endpoints of f MP are defined by ã± = (1 ± m/n) 2 .It is then possible to derive the inverse-Wishart ensemble of matrices B = A −1 .In terms of the eigenvalues b i = a −1 i , one simply has to find the distribution as follows: or, in matrix terms So it is easy to see that for n = 2, we have m = 1 which does not fulfill the condition m ≥ n.The inverse-Marchenko-Pastur distribution for the eigenvalues of an inverse-Wishart matrix simply reads, following from Equation (120): with the inverse rescaled eigenvalues b = ã−1 = 2nb.The new endpoints are b± = 1/ã ∓ .Clearly, the spectral density in Equation ( 123) is not well defined in the limit m = n.

Equivalent Dyson Brownian Motion
Consider the DBM ⃗ x(s) for β = 2 where the db i (s) are independent standard Brownian motions, with db i (s)db j (s) = δ ij ds, and we choose a fixed ordered initial condition x 0 1 > x 0 2 > • • • > x 0 N .At fixed time s, ⃗ x(s) has the same law as the spectrum of the random matrix X = diag(x 0 i ) + √ sH.Its propagator takes the form Until now, the initial condition is arbitrary.We now choose x 0 j = δ j /n, regularly spaced.In this case, det G 0 (x i , x 0 j , s/n) can be explicitly evaluated and one finds for x 1 > x 2 > • • • > x N as follows: To make the connection with the distribution P 0 ( ⃗ w, t) of the main text, we set s = 1/(4γtn) and x j = w j /(4γtn).Then, one recovers P 0 ; namely, one has This is in agreement with (86) in the main text.The above correspondence is exact and valid for any n.We now consider the large n limit.In the text, we scaled time as t = τ 4γn and considered the limit n → ∞ at fixed τ .As mentioned in the text around (86), we can focus on the DBM ⃗ x(s) in inverse time s = 1/τ .Let us denote µ s (x) as its density, and h(z, s) = 1 N i 1 z−xi(s) its resolvent.We will now compute both quantities, and from them we will deduce which are displayed in the text.The DBM density µ s (x) interpolates between being uniform in [−1/2, 1/2] at small s, and a semi-circle of support √ s[−2, 2] at large s.To compute it at all times, we recall that its resolvent h(z, s) satisfies in the large n limit the complex Burgers equation The general solution is Here, we have a uniform initial density Hence, one finds the parametric solution Upon rescaling (128) one obtains (91) in the text.To compute the DBM eigenvalue density µ s (x), one sets z = x−i0 + and h = h r + iπµ.Taking the imaginary part of (132) gives cosh h r = a s (µ) := cos πµ + sin πµ 2sπµ .
Inserting h r within the real part of (132), one finds where the two branches correspond to h r > 0 and h r < 0. Equation (134) determines µ s (x) for a given s.Upon the rescaling (128), one obtains (92) in the text.Note that an analogous calculation was performed in [76] in a different context.

Kernel
From the relation to the DBM described in the previous Appendix, and from Ref. [88], one knows that the ⃗ w form a determinantal point process.This means that both their joint PDF and their correlation functions (obtained by integrating over some of the w α ) are equal to determinants involving a kernel.Here, we obtain a convenient form for this kernel using bi-orthogonal polynomials.We follow the method of Ref. [93], c.f. Equations (3.11)-(3.13)therein.To simplify the expressions, in this section we fix γ = 1/4.Equivalently, one can recover the general form by replacing t → 4γt.One defines the two sets of polynomials n is the Stirling number of the first kind.These polynomials are not monic, but they are normalized to have unit integral.Indeed, they satisfy the orthogonality relation From these polynomials, one defines a first kernel as To obtain a centered process, however, one defines the shifted kernel Both kernels are self-reproducing, i.e., one has K 2 n = K n .The density of the ⃗ w for any n and any time t is then given by and it is normalized to unity.One can check that one has indeed We can also use a similarity transformation to re-write the Kernel as with with y k = 1 − k.This last form is analogous to Proposition 2.3 and Equation (2.18) in [88].We can further expand the last term, using One can check that the exact expression for the moments M (m) in (83) are recovered from M (m) = e −m(n−1)t/2 dwJ n (w, w)e wm (since the quadratic term in w cancels with the one in w ′ , the integration over w leads to a delta function), where the prefactor is related to the shift with respect to the edge position.
This form can be used to derive the limit of the kernel when n → ∞, focusing on the neighborhood of the edge.We set which is defined so that the largest Lyapunov is moved at ω = 0. From the explicit limit, we obtain the expression Exponentiating the denominator with an integral, we can also express with This is the kernel describing the Lyapunov exponents in the universal regime Γt = O(n).A similar kernel was obtained in formula (1.15) in [100] and in formula VI.22 in [101] in the context of the universal edge statistics of products of Ginibre matrices.The connection to the Dyson Brownian motion was also discussed in [101].

Large Time Moments from a Saddle Point
We perform here a calculation analogous to the one in Section III-D of Ref. [95].It is valid for any n.Let us denote O the ordered sector w 1 < • • • < w n .For ⃗ w ∈ Ω, one can rewrite exactly (not keeping track of time-dependent normalizing constants) P 0 ( ⃗ w, t) ∼ e −S , with S = S 0 + S int as follows: For γt ≫ 1 the term e −4γt S has a saddle point for The saddle point remains in the ordered sector Ω, providing m α − m α−1 + 1 > 0 for all α.That case corresponds to the particle crossing being irrelevant.The interaction term takes the form (up to a time-dependent constant) For γt ≫ 1, the last term is irrelevant compared to the first (provided the particle crossings are irrelevant).The saddle-point method then leads to Taking derivatives of log G[ ⃗ m] gives all the O(1) joint cumulants of the variables w α in the large time limit.This estimate is valid as long as the saddle point is in the O sector.
Specializing m α = δ α,n , i.e., the largest of the w's, gives In these results, we see that the prefactor includes the interactions with all the particles (not just the neighbor).Furthermore, we see from the condition that the saddle point remains in the Ω sector that (i) Equation ( 154) is valid for all m > −1 (in agreement with ( 83)) (ii) Equation (155) for α < n is valid only for −1 < m < 1.Indeed, for m > 0, the rightmost particle is pulled out of the gas, and hence particle crossings are irrelevant.An "internal" particle, however, cannot be pulled too strongly without crossing its neighbors.Finally, note that evaluating Ω(m) = ⟨( α e wα ) m ⟩ 0 by the same method leads to an additional term −mz n in S, plus another term − log(1 + α<n e −4γt(zα−zn) ), irrelevant at large times.Hence, we find Ω(m) ≃ M (m) by this saddle-point method in the γt ≫ 1 limit.

Long-Time Entanglement Entropy
Here, we derive the long-time asymptotic expression for the entanglement entropy, for arbitrary value of n.We start from the expression of the entanglement entropy as computed in Equation ( 63), namely where M (m) and Ω(m) are defined in Equations ( 61) and ( 62), respectively, and the average ⟨•⟩ 0 is taken on the probability distribution in Equation (63).Let us set 4γt = t for the rest of the section, and let us define w α = u α t, so that we have Let us notice that It is now convenient to introduce new variables u α = v α + ξ α / √ t.The integration domain is now up to O(t −1 ) terms, where we denoted with the symbol * α>β the sum over couples such that α > β and (α, β) ̸ = (n, n − 1), and with ∆(u) the Vandermonde determinant restricted to all pairs but the (n, n − 1) one.Let us notice that, by parity, the O(t −1/2 ) terms vanish once integrated on the Gaussian measure.Finally, we have that, up to O(t −1 ) terms, the numerator becomes t −n/2 e −t(I(v)−vn−1) ∆(v) ( Introducing the rotated variables ξ ± = (ξ n ± ξ n−1 )/ √ 2 the last integral becomes trivial, so we have the following for the numerator: (2π) n/2 t −n/2 e −t(I(v)−vn−1) ∆(v)(1 + (πt) −1/2 ) . (170) Denominator: imposing again the condition ∂ uα (I(u) − u n )| v * α = 0, one finds the solutions In terms of the new variables u α = v * α + ξ α / √ t the integration domain is now so that in the limit t → ∞ the integral will run on the whole R n up to exponentially small correction.In this case, the (n, n − 1) pair in the Vandermonde determinant should not be singled out.Analogously to the numerator case, moreover, the O(t −1/2 ) corrections coming from the Vandermonde determinant vanish once integrated.Finally, for the denominator we have the following: (2π) n/2 t −n/2 e −t(I(v up to O(t −1 ) corrections.
Entanglement Entropy: replacing the results (170) and (173) in the expression (165), we find Taking into account the explicit form of v α and v * α , we have whereas the prefactor reads as follows: Restoring t = 4γt, we reach the result of Equation (111) of the main text.

Figure 1 .
Figure 1.Inverse-Marchenko-Pastur distribution.The spectral density f ( λ) for the rescaled eigenvalues λ defined in Equation (28) is shown.In the large n limit, it takes the inverse-Marchenko-Pastur form given in Equation (35).In both plots, the orange line displays the theoretical curve, whereas the blue histogram bars are computed after a numerical simulation (with n = 50) of the weak measurement protocol.(a) The spectral density at x = 0.2 in the range [0, 4/x].At small x, most eigenvalues are located in vicinity of λ− → 1/4, as an effect of purification.(b) The spectral density at x = 5.0 in its domain [ λ−, λ+].For larger values of x, the rescaled eigenvalues take finite values around their average ⟨ λ⟩ = 2/x.

Figure 2 .
Figure 2. Short time behavior.The density fτ (w) in the short time t = τ /4Γ ∼ 1/Γ regime, with τ finite.fτ (w) features a crossover between a semi-circle at small τ and a square distribution at larger τ .(a) The small-τ semi-circle distribution is displayed for τ = 5 and increasing n from n = 10 to n = 100.The red solid line shows the theoretical n → ∞ curve.(b) The large-τ square distribution is shown for τ = 50 and increasing n, with the red solid line showing the n → ∞ curve.(c) The crossover from semi-circle towards square distribution is shown for increasing τ from τ = 5 to τ = 50.All solid lines represent the theoretical density wew • fτ (w) on the rescaled w/we axis, where we is the edge coordinate for fτ (w) in the n → ∞ limit.

Finite-
Time Von Neumann Entropy for the n = 2 Case
so that in the limit t → ∞, we only have the constraint ξ n−1 < ξ n , up to exponentially small corrections.Let us now consider the Vandermonde determinant.The pair (n − 1, n) yields a factor t −1/2 (ξ n − ξ n−1 ).Instead, the other terms give the following: β)̸ =(n,n−1)
2 , σ 3 } denotes the set of Pauli matrices.The modulus r = |⃗ r| coincides with the difference |λ 1 − λ 2 | between the two eigenvalues of ρ, whereas the normalization of the trace, forcing λ 1 + λ 2 = 1, is implicit in the parameterization.It is also convenient to similarly parameterize the measurement operator as O = o 0 1 + ⃗ o • ⃗ σ.We can assume Ô to be traceless so that o 0 = (o 11 + o 22 )/2 = 0, as this contribution would be inessential in either case.The three remaining one-indexed random variables o 1 = Re(o 12 ), o 2 = −Im(o 12 ), o 3 = (o 11 − o 22 )/2 are real-valued and satisfy o α o β = δ αβ /4.In terms of these variables, the infinitesimal variation dr α = tr(dρσ α ) reads and a continuum of scattering states ϕ p ( w), such that E p = p 2 , with p ∈ R. Let us notice that the number of bound states is always greater than the number of particles, i.
e., ⌊n(1 + x/4) + 1/2⌋ > n ∀x > 0, and the scattering states are empty.The n-fermion ground state is thus E 0 = k E k .Consequently, the n-fermion energy states of the original FP operator are simply given by