Efficient circular Dyson Brownian motion algorithm

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 $\mathcal{Z}$ is a normalization constant that has the interpretation of a partition function, and $\beta>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 non-trivial dynamics. The equilibrium as well as the non-equilibrium 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 long-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 $\beta=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.

One distinguishes between orthogonal ($\beta=1$), unitary ($\beta=2$), and symplectic ($\beta=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 time-reversal symmetry. This section considers the unitary class, which is arguably the most convenient one.

Let $H(t)$ be an $N\times 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\times N$ matrix $M$ sampled from the Gaussian unitary ensemble can be constructed as

where $A$ is an $N\times 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 $\tilde{M}=U^{\dagger}(t)MU(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 $\tilde{M}$ can be replaced by a new sample from the Gaussian unitary ensemble. The increments $d\lambda_{n}(t)=\lambda_{n}(t+dt)-\lambda_{n}(t)$ of the eigenvalues $\lambda_{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\to\infty$, $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\times 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\times 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 $\mathbbm{1}+i\sqrt{dt}M$, 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 orthonormality-preserving 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 $\mathbbm{1}+i\sqrt{dt}M$ is invariant under unitary transformations.

Let $\tilde{M}=U^{\dagger}(t)\,M\,U(t)$ with the time-dependent unitary matrix $U(t)$ chosen such that it diagonalizes $Q(t)$. As before, $\tilde{M}$ can be replaced by a new sample from the Gaussian unitary ensemble. Circular Dyson Brownian motion of the eigenvalues $e^{i\theta_{1}(t)},e^{i\theta_{2}(t)},\dots,e^{i\theta_{N}(t)}$ entails that the increments $d\theta_{n}(t)=\theta_{n}(t+dt)-\theta_{n}(t)$ of the eigenphases $\theta_{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+\mathcal{O}(x^{2})$. For $t\to\infty$, $Q(t)$ converges to a sample from the circular unitary ensemble irrespective of the initial condition $Q(0)$.

The Gaussian random matrix ensembles have the property that the sum of $n$ independent samples is a sample again, although with a prefactor $\sqrt{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}\neq e^{A+B}$ when $A$ and $B$ do not commute. Time-evolution 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\times 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 upper-triangular 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 upper-triangular matrix, for example, by requiring them to be non-negative. Let

Then, $Q\to Q\Lambda$ and $R\to\Lambda 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\times N$ unitary matrix with time-dependent 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\to\infty$) 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\tau$ is some yet undetermined function of $dt$, which for small enough $d\tau$ will be found to be equal to $dt$. The aim is to show that $U(d\tau)$ corresponds to $U(dt)$ of Eq. (9). In Eq. (10), the columns $u_{n}$ of $U(d\tau)$ 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 $Nd\tau\ll 1$, these substractions do not significantly alter the directions of the columns, which are then rotated randomly since $\mathbbm{1}+i\sqrt{d\tau}A$ is invariant under unitary transformations. As the columns of $U(d\tau)$ are rotated randomly, $U(d\tau)$ corresponds to $U(dt)$ of Eq. (9) for a proper choice of $d\tau$, provided that $Nd\tau\ll 1$.

The limitation on the maximum value of $d\tau$ 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 $\Sigma$ is a diagonal matrix with real-valued nonnegative entries. Löwdin symmetric orthonormalization gives the unitary matrix $U=U_{1}U_{2}^{\dagger}$, 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\to U$ by Löwdin symmetric orthonormalization, then $V^{\dagger}MV\to V^{\dagger}UV$ for unitary matrices $V$. The Ginibre unitary ensemble is invariant under unitary transformations. These two facts combined guarantee the rotations generated by $U(d\tau)$ to be random. Thus, $U(d\tau)$ corresponds to $U(dt)$ of Eq. (9) for a proper choice of $d\tau$, without $d\tau$ required to be small.

The relation between $dt$ and $d\tau$ can be established by requiring $U(dt)$ [Eq. (9)] and $U(d\tau)$ [Eq. (14)] to be identically distributed. For $U(d\tau)$ [Eq. (14)], let $u_{1}(d\tau)$ 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(\infty)=0$. Figure 1 shows that the ensemble average of $F$ is almost perfectly described at all times already for $N=10$ by $F=(1-Nd\tau/2)^{2}$ before and $F=\big{(}(1+Nd\tau)^{-2}+2\sqrt{2}Nd\tau\big{)}^{-1}$ after the intersection at $Nd\tau\approx 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^{-Ndt}$ since $e^{-Ndt_{1}}e^{-Ndt_{2}}=e^{-N(dt_{1}+dt_{2})}$. Equating $F(dt)$ to the piecewise expression $F(d\tau)$ introduced above gives

which can be inverted numerically to find $d\tau$ as a function of $dt$. Up to first order, the approximation $dt=d\tau$ can be made.

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)=\operatorname{diag}(1,1,\dots,1)$ is taken such that $Q(dt)=U(dt)$. The values of $d\tau$ corresponding to these values of $dt$ are given in the caption. One observes that the matrices on the left and right show identical characteristics.

A sample from the unitary equivalent of the Rosenzweig-Porter model, considered next, can be obtained as $Q(dt=N^{-\gamma})$ by taking $Q(0)=\operatorname{diag}(e^{i\theta_{1}^{(0)}},e^{i\theta_{2}^{(0)}},\dots e^{i\theta_{N}^{(0)}})$ with the phases $\theta_{n}^{(0)}$ sampled independently from the uniform distribution ranging over $[0,2\pi)$ [21]. See Refs. [42, 43] for an introduction to the Rosenzweig-Porter model and its relation to Dyson Brownian motion. Level statistics are here quantified by the average ratio $\langle r\rangle$ of consecutive level spacings [44, 45]. For unitary matrices with ordered eigenphases $\theta_{n}$, the $n$-th ratio is defined as

The average is taken over all $n$ and a large number of realizations. Wigner-Dyson level statistics are characterized by $\langle r\rangle\approx 0.600$, while Poissonian level statistics obey $\langle r\rangle\approx 0.386$. The Rosenzweig-Porter model shows a transition (at finite dimension, a crossing) from Wigner-Dyson to Poissonian level statistics at $\gamma=2$. When plotted as a function of $(\gamma-2)\ln(N)$, the average ratio is numerically found to be independent of $N$ (finite-size collapse) [46, 21]. 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\tau$.

Circular Dyson Brownian motion describes the Brownian dynamics of particles interacting through a long-range 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 $\beta=2$, physically corresponding to broken time-reversal symmetry). For short times $Ndt\ll 1$, Eq. (14) can be replaced by the computationally cheaper Eq. (10), and the first-order approximation $dt=d\tau$ 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 symplectic 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.