Efficient circular Dyson Brownian motion algorithm

Circular Dyson Brownian motion describes the Brownian dynamics of particles on a circle (periodic boundary conditions), interacting through a logarithmic, long-range two-body potential. Within the log-gas picture of random matrix theory, it describes the level dynamics of unitary ("circular") matrices. A common scenario is that one wants to know about an initial configuration evolved over a certain interval of time, without being interested in the intermediate dynamics. Numerical evaluation of this is computationally expensive as the time-evolution algorithm is accurate only on short time intervals because of an underlying perturbative approximation. This work proposes an efficient and easy-to-implement improved circular Dyson Brownian motion algorithm for the unitary class (Dyson index $\beta = 2$, physically corresponding to broken time-reversal symmetry). The algorithm allows one to study time evolution over arbitrarily large intervals of time at a fixed computational cost, with no approximations being involved.


I. INTRODUCTION
Brownian motion describes the stochastic dynamics of microscopic particles in a thermal environment [1,2].It connects a broad variety of topics, including thermal physics, hydrodynamics, reaction kinetics, fluctuation phenomena, statistical thermodynamics, osmosis, and colloid science [3].Brownian motion is intimately related to random matrix theory, which plays a key role in the understanding of quantum statistical mechanics and quantum chaos [4][5][6].Random matrices have eigenvalue statistics that typically can be studied using the so-called log-gas picture [7,8].For matrices with real eigenvalues, the joint probability distribution P of the eigenvalues is then written as a Boltzmann factor where Z is a normalization constant that has the interpretation of a partition function, and β > 0 is a parameter known as the Dyson index that has the interpretation of an inverse temperature.The Hamiltonian H describes a collection of classical massless particles on a line (the eigenvalues) repelling each other over long ranges through a logarithmic two-body potential, held together by a confining background potential.The log-gas picture describes long-range interacting particles.It has been found, for example, to accurately describe the level statistics across the many-body localization transition [9].As the Hamiltonian in Eq. (1) does not contain a kinetic term, the particles obey nontrivial dynamics.The equilibrium as well as the nonequilibrium dynamics of the particles ("level dynamics") are described by a phenomenon referred to as Dyson Brownian motion [10,11].Dyson Brownian motion turns out to provide a good description rather generically when * buijsman@pks.mpg.delong-range interactions are involved.As such, these dynamics (as well as the corresponding stochastic evolution of the eigenstates [12,13]) have found applications in studies on, for example, disordered systems [14][15][16][17], random matrix models [18][19][20][21], many-body localization [22,23], quantum information dynamics [24][25][26], and cosmological inflation [27][28][29][30].
Circular Dyson Brownian motion for unitary ("circular") matrices describes Dyson Brownian motion on a circle (periodic boundary conditions) and without background potential.A common scenario is that one wants to know about an initial configuration evolved over a certain interval of time, without being interested in the intermediate dynamics.Dyson Brownian motion can be evolved over a time interval of arbitrary length at a fixed computational cost, with no approximations being involved (see below for a more detailed explanation).Circular Dyson Brownian motion, however, requires extensively many evaluations over small intermediate intervals because of a perturbative approximation underlying the time-evolution algorithm.Circular Dyson Brownian motion is thus a process that is computationally expensive to simulate, which moreover is subject to a loss of accuracy with progressing time.Despite significant recent [31,32] and less recent [33][34][35][36] analytical progress on circular Dyson Brownian motion out of equilibrium, improved numerical capabilities are thus desired.
This work proposes an improved, easy-to-implement circular Dyson Brownian algorithm for the unitary class (Dyson index β = 2, corresponding to systems with broken time-reversal symmetry).The algorithm does not require intermediate evaluations, and thus operates at dramatically lower computational cost compared to the currently used algorithm.Moreover, it does not involve approximations, and is thus not subject to a loss of accuracy with progressing time.In short, it constructs the desired unitary matrices by orthonormalizing the columns of certain non-Hermitian matrices for which the elements perform Brownian motion.Similar to Dyson Brownian motion for Hermitian matrices, this Brownian motion process can be time evolved at a computational cost independent of the length of the time interval.

II. DYSON BROWNIAN MOTION FOR HERMITIAN AND UNITARY MATRICES
One distinguishes between orthogonal (β = 1), unitary (β = 2), and symplectic (β = 4) random matrix ensembles [7].These names reflect the type of transformations under which the ensembles remain invariant.Physically, the type of invariance determines the behavior of a system under time reversal.For example, the orthogonal class correspond to time-reversal systems, whereas the unitary class correspond to systems with broken timereversal symmetry.This section considers the unitary class, which is arguably the most convenient one.
Let H(t) be an N ×N Hermitian matrix with elements depending on time t [10].The initial condition H(0) can be either random or deterministic.Dyson Brownian motion for Hermitian matrices of the unitary class is a stochastic process described by where the time step dt, in order for the eigenvalue dynamics to obey Dyson Brownian motion, is supposed to be small enough such that the eigenvalues of H(t + dt) can be obtained accurately by second-order perturbation theory.Here, M is a sample from the Gaussian unitary ensemble that is resampled at each evaluation.An N ×N matrix M sampled from the Gaussian unitary ensemble can be constructed as where A is an N × N matrix with complex-valued elements A nm = u nm + iv nm with u nm and v nm sampled independently from the normal distribution with mean zero and variance 1/2.Let M = U † (t)M U (t), where the time-dependent unitary matrix U (t) is chosen such that it diagonalizes H(t).The Gaussian unitary ensemble is invariant under unitary transformations, meaning that M can be replaced by a new sample from the Gaussian unitary ensemble.The increments dλ n (t) = λ n (t + dt) − λ n (t) of the eigenvalues λ n (t) when evolving from time t to t + dt obey where terms of order three and higher have been ignored.It can be shown that this time-evolution indeed describes a Brownian motion process, for example by writing down the corresponding Fokker-Planck equation.For t → ∞, H(t) converges to a (scaled) sample from the Gaussian unitary ensemble irrespective of the initial condition H(0).
Dyson Brownian motion can also be studied for unitary matrices [10,33].Let Q(t) be an N × N unitary matrix with time-dependent elements.Similar to the above, the initial condition Q(0) can be either random or deterministic.Circular Dyson Brownian motion for the unitary class is generated by where again M is an N × N sample from the Gaussian unitary ensemble that is re-sampled at each evaluation.
For small enough dt, the matrix exponent can be approximated by the first-order expansion 1 + i √ dtM , which is invariant under unitary transformations (the second and higher-order terms are not).The matrix Q(t + dt) is thus obtained by applying infinitesimal orthonormalitypreserving random rotations on the columns of Q(t)."Random" here means that rotations in each direction are equally likely, which agrees with the observation that As before, M can be replaced by a new sample from the Gaussian unitary ensemble.Circular Dyson Brownian motion of the eigenvalues e iθ1(t) , e iθ2(t) , . . ., e iθ N (t) entails that the increments dθ n (t) = θ n (t + dt) − θ n (t) of the eigenphases θ n (t) when evolving from time t to t + dt are given by where terms of order three and higher have been ignored.Equations ( 4) and ( 6) describe similar dynamics on a microscopic scale since 2 tan(x/2) = x + O(x 2 ).For t → ∞, Q(t) converges to a sample from the circular unitary ensemble irrespective of the initial condition Q(0).

III. THE ALGORITHM
The Gaussian random matrix ensembles have the property that the sum of n independent samples is a sample again, although with a prefactor √ n.Equation ( 2) and its equivalents for the orthogonal and symplectic classes thus do not require the time step dt to be small.This implies that numerically obtaining H(T ) from H(0) can be done in a single instance, at a computational cost independent of T .Equation ( 5) for the evolution of unitary matrices does not allow for a similar argument since e A e B ̸ = e A+B when A and B do not commute.Timeevolution for unitary matrices can thus naively only be accomplished by subsequently evolving over infinitesimal time intervals.Equation ( 5) moreover is subject to a loss of accuracy with progressing time as it describes the desired dynamics only up to first order.The starting point in establishing an improved algorithm is the observation that a random unitary matrix (circular unitary ensemble) can be obtained by orthonormalizing a set of random vectors [37,38].Let A be an N × N matrix with elements A nm = u nm + iv nm with u nm and v nm sampled independently from the normal distribution with mean zero and unit variance.Such a matrix is known as a sample from the Ginibre unitary ensemble [8,39].The QR decomposition decomposes A in a unitary matrix Q and an uppertriangular matrix R with real-valued diagonal elements.This decomposition is not unique.It can be made unique by fixing the signs of the diagonal elements of the uppertriangular matrix, for example, by requiring them to be non-negative.Let Then, Q → QΛ and R → ΛR is the QR decomposition with the upper-triangular matrix having non-negative diagonal entries.One can prove that the resulting unitary matrices Q obey the distribution of the circular unitary ensemble.Algorithmically, such unitary matrices are obtained by performing Gram-Schmidt orthonormalization (discussed below) on the columns of A. A sample from the circular unitary ensemble can thus be obtained by orthonormalizing a set of random vectors.Let U (dt) be an N × N unitary matrix with timedependent elements.The goal is to express Q(t + dt) of Eq. (5) as where dt is not necessarily small.Equation ( 5) indicates that the dynamics of Q(t) are generated by orthonormality-preserving random rotations of the columns.Thus, U (dt) interpolates between an identity matrix (dt = 0) and a sample from the circular unitary ensemble (dt → ∞) in a way such that U (dt) is invariant under unitary transformations.In other words, it generates a finite orthonormality-preserving random rotation of the columns of Q(t).Generalizing the above algorithm generating random unitary matrices, consider the QR decomposition Here, A is again a sample from the Ginibre unitary ensemble.This ensemble is invariant under unitary transformations.The parameter dτ is some yet undetermined function of dt, which for small enough dτ will be found to be equal to dt.The aim is to show that U (dτ ) corresponds to U (dt) of Eq. ( 9).In Eq. ( 10), the columns u n of U (dτ ) result from Gram-Schmidt orthonormalization of the columns m n of the left-hand side, In words, the n-th column is obtained by substracting the projections on the first n − 1 columns, followed by normalization.Columns with a higher index undergo more substractions than columns with lower indices.For N dτ ≪ 1, these substractions do not significantly alter the directions of the columns, which are then rotated randomly since 1 + i √ dτ A is invariant under unitary transformations.As the columns of U (dτ ) are rotated randomly, U (dτ ) corresponds to U (dt) of Eq. ( 9) for a proper choice of dτ , provided that N dτ ≪ 1.
The limitation on the maximum value of dτ can easily be overcome by adapting a different, appropriate, orthonormalization procedure.Löwdin symmetric orthonormalization is a procedure for which the columns are treated symmetrically, that is, the outcome is independent of the ordering [40].For M denoting some matrix, consider the singular value decomposition Here, U 1,2 are unitary matrices and Σ is a diagonal matrix with real-valued nonnegative entries.Löwdin symmetric orthonormalization gives the unitary matrix U = U 1 U † 2 , which can be shown to be optimal in the sense that the distance between the columns m n of M and u n of U acquires the minimal possible value [41].This invites to consider the SVD-decomposition with A denoting a sample from the Ginibre unitary ensemble.If M → U by Löwdin symmetric orthonormalization, then V † M V → V † U V for unitary matrices V .The Ginibre unitary ensemble is invariant under unitary transformations.These two facts combined guarantee the rotations generated by U (dτ ) to be random.Thus, U (dτ ) corresponds to U (dt) of Eq. ( 9) for a proper choice of dτ , without dτ required to be small.The relation between dt and dτ can be established by requiring U (dt) [Eq.( 9)] and U (dτ ) [Eq. ( 14)] to be identically distributed.For U (dτ ) [Eq. ( 14)], let u 1 (dτ ) denote the first column (the choice for the first column is arbitrary), and consider the overlap The overlap is shifted and scaled such that F (0) = 1 and F (∞) = 0. Figure 1 shows that the ensemble average of F is almost perfectly described at all times already for N = 10 by F = (1 − N dτ /2) 2 before and after the intersection at N dτ ≈ 0.66211710937.These expressions have been found empirically.Next consider U (dt) [Eq.( 9)].Equation (5) dictates, as can be verified numerically, that the product of two independent samples U (dt 1 ) and U (dt 2 ) is from the same distribution as U (dt 1 dt 2 ).For F defined similar as above, this means that F (dt) = e −N dt since e −N dt1 e −N dt2 = e −N (dt1+dt2) .Equating F (dt) to the piecewise expression F (dτ ) introduced above gives which can be inverted numerically to find dτ as a function of dt.Up to first order, the approximation dt = dτ can be made.

IV. NUMERICAL VERIFICATION
This Section provides a numerical verification of the algorithm proposed above.First, the focus is on the structure of the resulting matrices.Figure 2 shows density plots of |Q(dt)| 2 for matrices of dimension N = 50 at short (dt = 0.02) and longer (dt = 0.05) times obtained through Eq. ( 5) [left, "naive"] and Eqs. ( 9), (14), and ( 16) [right, "efficient"].The initial condition Q(0) = diag(1, 1, . . ., 1) is taken such that Q(dt) = U (dt).The values of dτ corresponding to these values of dt are given in the caption.One observes that the matrices on the left and right show identical characteristics.
The average is taken over all n and a large number of realizations.Wigner-Dyson level statistics are characterized by ⟨r⟩ ≈ 0.600, while Poissonian level statistics obey ⟨r⟩ ≈ 0.386.The Rosenzweig-Porter model shows a transition (at finite dimension, a crossing) from Wigner-Dyson to Poissonian level statistics at γ = 2.When plotted as a function of (γ − 2) ln(N ), the average ratio is numerically found to be independent of N (finite-size collapse) [21,46].Figure 3 shows that the algorithm proposed in this work leads to the same results, and illustrates the capability of the algorithm proposed in this work to operate at large matrix dimensions (here, up to N = 10 000).Reference [47] (Fig. 1) shows a visually indistinguishable plot obtained using Eqs.( 9), (10) with the first-order approximation dt = dτ .

V. CONCLUSIONS AND OUTLOOK
Circular Dyson Brownian motion describes the Brownian dynamics of particles interacting through a longrange two-body potential in a one-dimensional environment with periodic boundary conditions.This work proposed an easy-to-implement algorithm [Eqs.( 9), ( 14), and ( 16)] to simulate circular Dyson Brownian motion for the unitary class (Dyson index β = 2, physically corresponding to broken time-reversal symmetry).For short times N dt ≪ 1, Eq. ( 14) can be replaced by the computationally cheaper Eq. ( 10), and the firstorder approximation dt = dτ can be used instead of the more complicated relation (16).The latter approach is a generalization of a commonly used algorithm generating samples from the circular unitary ensemble, proposed in Refs.[37,38].In contrast to the currently used circular Dyson Brownian motion algorithm [Eq.( 5)], here the time step dt does not have to be small, and no approximations have been involved.This allows one to study time-evolution over arbitrarily large time intervals at a computational cost independent of the length of the time interval, without loss of accuracy.In typical settings, this algorithm dramatically reduces the computational costs, thereby for example opening the possibility to perform detailed studies without the need for high-performance computing facilities.
An arguably interesting follow-up question would be how to modify the algorithm for the orthogonal and sym-plectic classes.From a sample Q of the circular unitary ensemble, a sample S from the circular orthogonal ensemble can be obtained as S = Q T Q [48].It is thus tempting to hypothesize that circular Dyson Brownian motion for the orthogonal class can be simulated by the algorithm proposed in this work, and by taking the product of the transpose of the resulting unitary matrix and the resulting unitary matrix itself as the output.
Circular Dyson Brownian motion can be used to numerically generate non-ergodic unitary matrices ("unitaries") with fractal eigenstates and a tunable degree of complexity [21,47].Next to what is mentioned above, this work can thus be expected to be relevant for future studies on the emergence and breakdown of statistical mechanics in the context of unitary (periodically driven) systems.It also relates to recent developments on algorithms generating random rotations [49].Dyson Brownian motion recently attracted a spurge of interest in the context of the Brownian SYK model [50][51][52][53][54][55][56].Unitary Brownian quantum systems are of current interest in the context of Brownian quantum circuits [57][58][59][60][61][62].This work finally can be expected to provide new opportunities in the context of the non-trivial dynamics of Brownian quantum systems.