Unitary long-time evolution with quantum renormalization groups and artificial neural networks

In this work we combine quantum renormalization group approaches with deep artificial neural networks for the description of the real-time evolution in strongly disordered quantum matter. We find that this allows us to accurately compute the long-time coherent dynamics of large, many-body localized systems in non-perturbative regimes including the effects of many-body resonances. Concretely, we use this approach to describe the spatiotemporal buildup of many-body localized spin glass order in random Ising chains. We observe a fundamental difference to a non-interacting Anderson insulating Ising chain, where the order only develops over a finite spatial range. We further apply the approach to strongly disordered two-dimensional Ising models highlighting that our method can be used also for the description of the real-time dynamics of nonergodic quantum matter in a general context.

Introduction. The understanding of emergent behavior in quantum many-body systems is largely based on the discovery of effective descriptions of analytically unsolvable models [1]. An essential toolkit to find the former constitute renormalization group (RG) methods. They are traditionally applied on systems in thermal equilibrium, thereby explaining many collective phenomena including structured phases, phase transitions, critical scaling and universality. In the past decade, real-space RGs have been developed that aim to explain analogues of these well-known phenomena also in systems where a thermodynamic treatment breaks down due to strong quenched disorder [2][3][4][5][6][7][8][9][10].
Whereas real-space RGs successfully operate in the stationary setting at the level of individual eigenstates [11][12][13][14][15][16], reaching a quantitative description of the dynamical properties of quantum many-body systems appears even more challenging. So far, coherent dynamics of quantum matter far from equilibrium has been mostly simulated using tensor networks methods [17][18][19] or exact diagonalization [20,21] with recent developments targeting dynamical descriptions in terms of machine learning methods by utilizing Restricted Boltzmann Machines (RBM) [22] or, more general, Artificial Neural Networks (ANN) [22][23][24]. Still, accessing quantitatively the long-time dynamics for large quantum many-body systems, especially in spatial dimensions beyond one, represents a major challenge [25][26][27].
In this work, we show how ANNs can be utilized in a different way for numerically exactly time-integrating effective descriptions of generically interacting systems, generated by RG methods. As a concrete example, we explore the temporal build-up of MBL spin-glass order out of a simple polarized state for a large, disordered spin chain, see Fig. 1b, among other long-time dynamics in 1D-and 2D-lattices. We begin by formulating a prototypical strong-disorder RG (SDRG) for spin-1/2 systems of arbitrary spatial dimension and map its transformations into the time-domain. As a result, we obtain a quantum In the course of the RG, long-distance and higher-order couplings emerge. Adding time-dependence leads to a further broadening with increasing time. (b) Spatiotemporal build-up of MBL spinglass order in a random quantum Ising chain with 64 lattice sites after quenching a paramagnetic initial condition into the symmetry-broken phase at J = 5h, J (x) = h/5. Dashed lines indicate emergence of light-cones except for the non-interacting case, J (x) = 0, where the order stops developing at a finite distance. The numerical data is obtained from an average over 25 disorder realizations.
circuit, see Fig. 1a, as an effective description of the timeevolution operator. Hereafter, we show that this circuit can be encoded efficiently into deep ANNs associated with typical initial conditions for quantum real-time dynamics. This allows us to quantitatively represent time-evolved many-body quantum states not only at short but also long times. We note that our method avoids a discretization of time but relies on a renormalized Hamiltonian that is assumed to effectively describe the relevant physics up to some finite but nevertheless long time scale.
The scheme we introduce in the following can be applied for any generic, but strongly disordered, spin-1/2 system. Concretely, we will apply it to a paradigmatic interacting disordered quantum Ising model [28] of the form with next-neighbor couplings J ij ∈ [−J, J] and local magnetic fields h i ∈ [−h, h] drawn randomly from uniform distributions. We use periodic boundary conditions. For the 1D case we also add random transverse couplings J ] to obtain a generic and interacting model.
Solving the time evolution. Before describing the utilized renormalization procedure and the training of the ANN in detail, let us start by outlining the general scheme for solving quantum real-time evolution utilizing strongdisorder RGs. Such an RG generates a sequence of local unitary transformations U k in order to iteratively obtain a simplified effective description of the considered quantum many-body system. In the time domain, we will show that this leads to the following representation of time-evolved quantum many-body states: where a time-dependence is added to the RG-transformations U k through a generalized interaction picture, see the derivation below. The above equation maps quantum dynamics onto a quantum circuit generated by the local unitaries U k (t). As we assume that the effective description in terms of the final Hamiltonian H (n) 0 after the end of the RG procedure can be solved exactly, the complexity of the quantum circuit emerges solely unitaries U k (t). We find that such quantum circuits can become a non-perturbative object, as the spatial support of the U k (t) typically grows over time developing longdistance and higher-order couplings with large overlaps, see Fig 1a. A central contribution of this work is to outline a numerically exact scheme to encode |ψ(t) QC and therefore the RG transformation itself into an ANN using machine-learning techniques. The numerical learning effort in obtaining |ψ(t) QC , as well as its memory requirement, scales at most quadratically with system size while being independent on the targeted time t or the spatial dimension.
Dynamical strong-disorder Renormalization Group. In principle, quantum circuits such as in Eq. (2) can be generated using a variety of standard SDRGs. In the following we introduce a variant of a SDRG, which as we find improves the quantitative accuracy of the resulting scheme.
As other SDRGs, the dynamical variant we introduce is based on a local separation of energy scales. Consequently, at the beginning of each iteration k we pick the strongest coupling, also called "fast mode", whose corresponding term in the Hamiltonian we call H 0 . For the first iteration this could be either a spin interaction J ij , J (x) ij , or transverse field h i , see the Hamiltonian in Eq. (1). Those terms in the Hamiltonian which are not commuting with H 0 we denote by V . These can be eliminated perturbatively using a Schrieffer-Wolff transformation (SWT) [29] by applying a unitary transformation W k = e S k on the Hamiltonian with a generator S k satisfying [H 0 , S k ] = V and S † k = −S k [12], at the expense of the renormalization In general, this modifies existing couplings and leads to the generation of new terms in the Hamiltonian. After the SWT the fast mode is decoupled from the remainder and can then be faithfully removed from the system as a second-order local integral of motion (LIOM) [10,30,31]. After n such iterations, an unperturbed Hamiltonian H (n) 0 is obtained, formed by the set of LIOMs.
The newly generated couplings after each iteration are, of course, not known a-priori, especially if the SDRG is designed regardless of details of the model like range of interaction, dimensionality etc. We approach this problem by represent at each stage of the RG the Hamiltonian as a sum of arbitrary Pauli-strings σ α1 l1 . . . σ α M l M with a real coefficient λ l1,...,l M each. Certainly, this approach can entail a costly handling of numerous generated higherorder couplings, see below, but it opens the possibility to take into account many-body resonances, which are neglected using earlier SDRGs [11,32] and related, socalled flow equation approaches [27,33,34].
In addition the accuracy of the RG can be further increased by splitting the SWT into infinitesimal unitary transformations, closely resembling in spirit the flow equation framework. This turns out to be particularly helpful in the vicinity of a critical point, h ≈ J here for the 1D model, where the SWT is least controlled and in order to capture many-body resonances to an arbitrary degree. For a detailed presentation of the technical details, see the appendix. To control the exponential number of couplings {V i } generated during the RG, we first neglect those terms where |V i | t * −1 which are much smaller than the inverse of the targeted time scale t * say, as they do not influence physics up to t * . Secondly, we perform the continuous renormalization only w.r.t. those V i whose relative magnitude lies above a fixed threshold, |V i |/|H 0 | > 1. Therefore we have a tradeoff, that is controlled by , between exactness and total number of couplings within the RG-generators S k and the renormalized Hamiltonian H (n) 0 . In our computations, typically ranges from 10 −4 . . . 10 −2 , depending on the closeness to the critical point, h ≈ J, or the ergodic transition, J (x) ≈ J. Later we will present a quantitative analysis of our RG w.r.t. the dynamics of local observables.
To derive the timedependence of U k (t) we express the time-evolution operator in the renormalized basis, which yields We achieve a much more robust learning of the ANN upon successively commuting each factor e S k to the left until its counterpart e S † k (t) is reached. Identifying U k (t) = eS † k (t) eS k gives then the desired form as in (2). Here, S k denotes the total application of all rotations from e S l , l < k on S k , see the appendix for details.
Training the artificial neural network. Utilizing ANNs as a variational ansatz for many-body wavefunctions has seen an active development recently [22,35] becoming competitive with or partially even superior to other state-of-the-art methods. [23,36,37]. In contrast to the commonly used time-dependent variational principle (TDVP), we introduce another way of training an ANN.As such, the scattering operators U 1 (t), . . . , U n (t) are consecutively used during n iterations to train the network. As the U k (t) are still local operators with a finite support in real space, we perform for each iteration k a supervised learning procedure to find the set of complex network parametersW (k) that minimize the Fubini-Study metrics, given by Notice that we assume always properly normalized wave functions. The network H ANN (W, s) can be considered as a deep extension of a complex-valued RBM with up to three hidden layers, see the appendix for details. After convergence, the "learned" solutioñ W (k) is passed to the next iteration as W (k+1) . To complete the learning procedure, we write L[W (k) ] at the k-th iteration, while omitting the index, as with ψ( s) = s|ψ W and U loc ( s, t) being the equivalent of the local energy known from TDVP methods. We access the above sum with a Markov chain Monte-Carlo (MCMC) algorithm which, as we can confirm, is sign-problem free in all our computations. In the same way, we calculate the gradient ∂L/∂W i with the backpropagation algorithm and pass the result to a stochastic gradient-descent optimizer referred to as PADAM [38].
Benchmarking. In order to quantify the overall accuracy of our approach we first benchmark the RG-component and the machine learning part individually. For the former task, we calculate Eq. (2) for small system sizes exactly using a matrix representation of the quantum circuit.  2 shows a comparison of the local magnetization with the result obtained from exact diagonalization for a system of L = 12 spins. The plot reveals that the accuracy of the dynamics depends crucially on the inclusion of many-body resonances, which is tuned by the only free RG-parameter , see above. For practical purposes, we set indirectly by imposing a maximum total number n of couplings within all RG-generators S k . Here, n = 3(10)L corresponds to the label of excluded (included) many-body resonances and matches (exceeds) the number of original couplings. Already for n = 10L we observe a very good agreement even for the longest times. Importantly, the result can be systematically improved by increasing n.
Next, let us benchmark the training of the ANN. For this purpose we ideally would like to check the overlap F = | ψ W (n) (t)|ψ QC (t) | of the final ANN-state to the one obtained from an exact application of the quantum circuit, which is impossible for large system sizes. Nevertheless, we can offer a lower bound | denotes the partial overlaps measured at the end of each iteration k, which are a by-product of the training procedure. We plot F * in Fig. 2 as a function of time for 1D-and 2D-lattices. It shows a high, macroscopic overlap even for large sys- tem sizes and a systematic improvement on adding more units and hidden layers to the ANN. From this finding we conclude that the quantum circuit can be applied essentially numerically exactly on the ANN. For comparison, we also plot the result of a perturbative treatment |ψ QC (t) ≈ e −iH (n) 0 t s k exp ( s| S k − S k (t) |ψ 0 ) | s , i.e. a cumulant expansion of the quantum circuit, that neglects (higher-order) commutators between different S k . It shows a rapid decay and thus confirms the circuit's non-perturbative nature. In the appendix we show further benchmarks of the whole framework for a large integrable system.
Numerics. As an application of our framework we now explore non-equilibrium dynamics involving global quenches that has been difficult to access so far in the large system size and long-time limit. It is known from previous RG-studies that a symmetry-broken state will keep a non-zero Edwards-Anderson order parameter in the long-time limit starting from symmetry-broken states if the system is in the MBL-spin glass (MBL-SG) phase [32]. Here, we aim to address the build-up of spatiotemporal order starting from a Z 2 -symmetric state upon quenching into the MBL-SG phase. We detect the spatiotemporal dynamics of the MBL-SG order via [39], where ij denotes the reduced density matrix of two lattice sites i, j, while ν enumerates its four eigenvectors | ij and eigenvalues (probabilities) p (ν) ij . Fixing a distance |i − j| we average χ ij (t) across all associated pairs and disorder realizations. This quantity can be interpreted as a local version of the Edwards-Anderson order parameter, which is otherwise mostly used to detect MBL-SG order in a static context, but which doesn't exhibit a natural extension to the dynamical regime considered here. Figure 1b shows χ d (t) both for an interacting MBL (J (x) = h/5) and a non-interacting Anderson localized (J (x) = 0) case for a 1D chain of 64 spins. At short times tJ≤J/J (x) = 25, an almost identical light-cone for the buildup of MBL-SG correlations is visible, which appears consistent with a logarithmic growth. On longer time scales we observe a fundamental difference between the Anderson and MBL cases. For the non-interacting Anderson-localized limit the growth of MBL-SG order stops, while for J (x) > 0 a second light-cone arises at a timescale that we estimate as ∼ 1/J (x) . Interestingly, we find that all light-cones do not become more open as we quench deeper into the MBL-SG phase but rather the more close we quench to the critical point. This behavior is reminiscent of the l-bit picture, where LI-OMs become more extended on approaching criticality. We will draw a connection to this picture below. As expected, a quench within the MBL-PM phase does not show any SG-order. Right at criticality, J = h, even without interaction, we find that the order becomes genuinely long-range as it decays algebraically with distance within the light-cone. For the interacting case, inside the SG-phase, we observe an exponential decay with distance, but having an essential difference to the non-interacting case: the order at any fixed distance does not saturate, but increases strictly monotonically for all observed times within the light-cones. This is a drastic non-perturbative effect of the interacting model. It is particularly obvious for next-neighboring spins, see Fig.1b). The important question whether this growing will eventually lead to a finite plateau for |i − j| → ∞ requires access to even much later times, which we currently cannot access. When initializing the system in a symmetry-broken state, as studied in previous works, the stability of MBL-SG order originates from the large overlap with the LI-OMs. The mechanism for the build-up of long-range order from symmetric states as targeted in this work is of fundamentally different origin, as the initial state is oriented orthogonal to the LIOMs. Here, it is essential to generate long-distance quantum correlations between LIOMs. This is not possible in the Anderson localized limit because the LIOMs are independent, as we also see from our results in Fig. 1. Only in the interacting MBL limit the MBL-SG order can develop. Quantum correlations between two lattice sites i and j can emerge on a time scale [J (x) ] −1 e |i−j|/ξ where ξ denotes a typical localization length. Consequently, at a given time t MBL-SG order can be generated over distances d ∼ ξ log[J (x) t] explaining the appearance of the logarithmic light-cone in Fig. 1.
As a closing point, we now turn briefly to quantum many-body dynamics in two dimensions. Whether a nonergodic phase due to strong disorder exists there has remained an outstanding challenge [40]. Its difficulty originates from the percolation of many-body resonances [41,42]. We find that at least for sufficiently small or large external fields, the latter can be effectively captured using our framework up to an unprecedented long timescale. Fig.  3b shows the temporal evolution of the local magnetization in a quadratic, rectangular lattice, using essentially the same quench protocol as above. In contrast to the glassy dynamics of a chain, see Fig. 3a, the lattice exhibits a rapid decay of magnetization at h J, consistent with thermalization. On the other hand, for h J a stable non-thermal plateau is reached. Our result therefore numerically confirms a presumed quasi-localization [41,42] in the disordered 2D transverse-field Ising model at infinite temperature.

Conclusion.
We have demonstrated how many-body quantum dynamics can be simulated for generic spin-1/2 systems up to exponentially long times given that sufficiently strong disorder breaks ergodicity at least up to the targeted timescale. Importantly, this includes an unbiased treatment of many-body resonances, which allowed us to obtain quantitative results in general and to go beyond one-dimensional systems. We could show that our proposed framework does not fundamentally rely on any specific details of the model and scales up to systems sizes far beyond of what is possible with exact diagonalization. This opens up for broad investigations e.g. of non-thermal behavior and quantum aging dynamics in higher dimensions [43,44], long-range interacting systems [45][46][47] or localization in lattice gauge theories [48]. Since this work has shown that deep ANNs are able to apply the proposed quantum circuit numerically exact, the ansatz could also be well suited for random unitary circuit models e.g. to study operator spreading [49] [50] [51] or measurement induced localization transitions [52] [53].

Structure of the artificial neural network
We use a complex-valued feed forward network taking a spin-configuration s as input layer and returning the activation of a single output unit as H ANN (W, s) = log[ψ W ( s)]. In between those, there are one or more hidden layers, each passing the previous, weighted and biased activations W (ν) v (ν−1) + b (ν) through a nonlinear activation function f (z) to the next layer, v which justifies our choice, although we note that a formal way of deriving an optimal f (z) does not exist in machine learning. In the above definition, W = ( a, b, W ) summarizes all network parameters. Unfortunately, however, by using f (z) = log cosh(z) we frequently observe the occurrence of numerical instabilities during training, caused by two poles located at ±iπ/2. These instabilities are triggered whenever z comes close to those poles. This manifests itself in sudden jumps of L[W (k) ], which can ultimately make convergence impossible. To fix this problem, we use an approximationf (z) ≈ log cosh(z) that "smooths" the poles while