{'fraction_non_alphanumeric': 0.04060838747784997, 'fraction_numerical': 0.011311281748375665, 'mean_word_length': 4.3807405053233754, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Typicality of nonequilibrium (quasi-)steady currents

The understanding of the emergence of equilibrium statistical mechanics has progressed significantly thanks to developments from typicality, canonical and dynamical, and from the eigenstate thermalization hypothesis. Here we focus on a nonequilibrium scenario in which two nonintegrable systems prepared in different states are locally and non-extensively coupled to each other. Using both perturbative analysis and numerical exact simulations of up to 28 spin systems, we demonstrate the typical emergence of nonequilibrium (quasi-)steady current for weak coupling between the subsystems. We also identify that these currents originate from a prethermalization mechanism, which is the weak and local breaking of the conservation of the energy for each subsystem.

Introduction.-Daily experience teaches us that when two large systems prepared in different equilibrium states are put into contact, a long-lasting current, which we will refer to as (quasi-)steady, emerges between them. Furthermore, after this possibly long intermediate time, which depends on the size of the two objects, the overall system relaxes to an equilibrium state. To understand and characterize the emergence of the (quasi-)steady current, generally one considers each system to be described by an equilibrium, canonical or microcanonical, ensemble. However, here we ask ourselves if the systems driving the (quasi-)steady state current can be described by single pure states. Advances in "pure state quantum statistical mechanics" have shown that equilibration can emerge in isolated quantum systems [1,2]. The seminal papers by Deutsch [3] and Srednicki [4] brought to the eigenstate thermalization hypothesis (ETH), which postulates that, for a nonintegrable system, eigenstates that are close in energy have similar local properties. ETH has been tested in a variety of systems [5] and reviews can be found in Refs. [1,[6][7][8][9][10]. Meanwhile, the attention on the foundation of statistical mechanics has led to the development of the notion of typicality [2,[11][12][13], within canonical formalism [14,15] and for unitary dynamical processes [12,16] with more specific scenarios such as prethermalization [17] and perturbation [18,19]. For example, dynamical typicality states that pure states with the same initial expectation value for some observables will likely have similar expectation values at any later time [12,16]. This was used to show that a "weak" version of ETH [20][21][22][23][24] is both necessary and sufficient for the vast majority of states to thermalize [17].
Within the thermalization dynamics, the system may experience prethermalization [25][26][27][28][29][30][31]. A key signature of this phenomenon is a separation of timescales during the relaxation, for instance, a fast initial dynamics in which the system relaxes to an intermediate (often called prethermal) state, and a slower relaxation toward the true thermal state. This separation of time scales can be seen in systems with a weakly broken conserved quantity, and it has been explored in ultracold atoms experiments [32]. In Ref. [33] it was shown that prethermalization is typical in the presence of weak coupling.
Much less is known regarding the emergence of nonequilibrium (quasi-)steady current when coupling two systems in pure states. Here we show the emergence of such typical quasi-steady current based on the notion of dynamical typicality. We highlight that there have been remarkable studies on the emergence of typical dynamics for nonequilibrium systems [34][35][36], and also insightful works to extend ETH to open quantum systems [37,38]. However, in our work, we do not assume a priori that a steady current can be reached, and we consider a unified framework for both the emergence of (quasi-)steady currents and thermalization.
For the two "baths" we take two nonintegrable spin chains coupled at one of their edges, and we initialize the baths either in single eigenstates or random pure states taken within an energy shell. For weak coupling between the systems, we show that the resulting (quasi-)steady current is typical in the sense that it converges, when the system size increases, to what would be obtained from initializing the baths in microcanonical states. We also verify that the value of the current converges towards the prediction from the ETH ansatz. Furthermore, we are able to show that the dynamics that leads to the formation of a long-lasting current, and the eventual thermalization of the two coupled chains, can be understood in the framework of prethermalization. In fact the dynamics of each bath, when decoupled, conserves their own energy, while the coupling between the baths breaks this conservation law. For weak coupling between the two chains one thus expects a slow "prethermalizationlike" dynamics, where the prethermal state actually ap-proaches a nonequilibrium steady state. Finally, we numerically show that the relaxation dynamics is proportional to the square of the coupling between the chains, and inversely proportional to their length, which can be derived within the prethermalization paradigm [30,31]. We thus expect, in the thermodynamic limit, that the current will exist indefinitely.
Model .-Two finite-size quantum systems H L and H R are treated as the left and right baths, which are coupled via an interaction term V , i.e., the total Hamiltonian of the system is H = H L + H R + V . Each bath is chosen to be a nonintegrable spin chain with bond and site Hamiltonian (h b n and h s n ) given by while for H R the site labelling ranges from N L +1 to N . Here N = N L +N R is the total number of spins, while N L and N R are the length of the left and right bath respectively. The interaction term V is given by where γ is the coupling strength. We consider J yz = J zz , h x = −1.05J zz , h z = 0.5J zz , for which values each bath is nonintegrable with Gaussian unitary ensemble (GUE) level statistics [39]. The Hamiltonians H L and H R , which have no conserved quantities apart from energy, have also been used to study out-of-time-ordered correlators [40]. In the following we will work in units for which J zz = = 1. Initial conditions.-We are interested in the currents generation when the baths are prepared with different local equilibrium properties, e.g., the baths are prepared at different initial energies per spin. More specifically, for each bath, we consider an energy shell Ξ L/R = [E L/R − ∆ L/R /2, E L/R + ∆ L/R /2] where E L/R is the average energy and ∆ L/R the width of the shell.
To this end we consider three different scenarios: i) eigenstate pairs -each bath is prepared at an arbitrary eigenstate corresponding to the shell Ξ L/R , i.e., where |i L (|j R ) is the i(j)-th eigenstate of H L (H R ); ii) typical-state pairs -each bath is a typical superposition of states within the shell, with random complex numbers c L i and c R j drawn from the Haar measure, i.e., iii) microcanonical ensemble pairs -each bath is prepared in a microcanonical state of the energy shell Ξ L/R , with ρ mic Ξ L/R and d L/R is the number of states in the shell Ξ L/R .
Energy currents.-The energy current operator with respect to the left environment is defined as The expectation value I exact = Tr I L ρ(t) gives the exact current from evolving the initial condition via the full Hamiltonian H, and it takes three forms I eig exact , I typ exact and I mic exact depending on whether the initial condition ρ(0) is |ψ eig ψ eig |, |ψ typ ψ typ | or ρ mic [41]. A true bath, by definition, should be infinitely large with a continuous spectrum. In our simulation, we would like to consider systems as large as possible to avoid finite-size effects. This poses a great numerical challenge due to the exponential growth of the space and time complexity. For this purpose, various strategies have been proposed [42,43]. We developed a highly optimized time evolution algorithm that allows to study system sizes up to N = 28 on a personal computer. Two main numerical techniques are used: i) a Suzuki-Trotter based time evolution algorithm which only requires the storage of a single quantum state; ii) a highly parallelized and cache-friendly implementation of the gate operations (See Ref. [44] for details).
Since such exact simulations are numerically demanding, we complement these calculations with a perturbative approach. On top of significantly reducing the time and memory demands of the calculations, the perturbative approach also allows us to gain analytical insights on the typicality of the currents as well as the importance of weak coupling between the baths. It also shares the same spirit of the weak coupling master equations formalism [45][46][47][48][49]. In the weak coupling limit γ → 0, the perturbative current expression with respect to different initial conditions can be written as [44,50] whereḂ and we have defined the two-time correlation func- We have also defined B L/R (t) = e iH L/R t B L/R e −iH L/R t and the same applies toḂ L (t). We call I typ perturb , I mic perturb , I eig perturb as typical, microcanonical, and eigenstate currents respectively depending on the initial conditions used. Note that when equilibrium states or single eigenstates are used, the expression for the current Eq. (5) can be significantly simplified and the first term vanishes [44].
Typical (quasi-)steady currents.-We first consider weak coupling γ = 0.01 where we expect the perturbative results to be consistent with the exact ones. In Fig. 1 we consider baths each of size N L = N R = 12 and with energy windows Ξ L/R given by [4.95, 5.05] and [−5.05, −4.95], respectively. In Fig. 1(a), we depict I both for exact and perturbative calculations (respectively for empty symbols and lines), and for the three different initial conditions. We observe perfect agreement between exact currents and perturbative currents for all the initial conditions considered. For the typical initial conditions, the results present larger oscillations, due to the dynamics of initial coherence for typical states. These results validate the use of perturbative methods for computing currents as all the currents computed are close to each other. In addition, the microcanonical current becomes constant after a time of order 1 indicating the emergence of a (quasi-)steady current.
Another question to be answered is whether also the currents I typ exact or I typ perturb can be considered as typical, i.e., the vast majority of the currents using different typical initial conditions show similar dynamics. This means that if we consider an ensemble of typical currents, their average E [I typ ] ≈ I mic with a bounded variance. Perturbatively, one can show this for the (quasi-)steady current in a similar fashion to the study of typicality in isolated quantum systems [11,34]. In fact, the average typical correlations functions E [C typ ] is equivalent to the microcanonical ones by considering, for example, where  delta. Since the left and right correlation functions are independent, it is thus straightforward to conclude from Eq. (5) that the average typical currents is equal to the microcanonical one, i.e., E I typ perturb = I mic perturb . This is an indicator of the typicality of (quasi-)steady currents. However, one should also consider the variance of the typical currents. To this end, we perform numerical computations of the exact currents I typ exact and I eig exact at a time t = 10, and provide a histogram of these values in Fig. 1(b). Already for bath sizes N L/R = 12, we observe a clear peak near the prediction from microcanonical initial states (continuous black line) which corresponds to the average of the typical states too (solid blue square).
We then analyze the dependence of the variance with the size of the baths. Relying on the extensive character of H L and H R , we linearly shift the energy shell for the initial conditions as the baths' size increases, i.e., Ξ L/R are such that E L/R = N L/R ε L/R and ∆ L/R = N L/R δ L/R , where ε L/R is the energy per spin in each bath and δ L/R is the shell width per spin chosen to be 0.1/12. In Fig.  2 (a), we show the standard deviation for different realizations of typical currents and eigenstate currents std[I] (rescaled over the value of the current E [I]) versus system size N . Such decaying behavior serves as numerical evidence of the typicality of nonequilibrium currents [51]. Fig. 2(b), which depicts the value of I mic perturb (to which the typical current converges) for different energy per spins, also shows a clear convergence of the current with the system size. In fact, larger systems sizes, each denoted by different markers, result in smaller oscillations and the current versus energy curve converges towards a smooth line. To have a deeper understanding of the emergence of this smooth curve, we compare the numerical perturbative results with predictions from ETH. The ETH ansatz, which each bath follows [44], states that for a local observable O, its matrix elements in the energy basis are [4,52] where is a mildly varying function of E and ω and R ij is a matrix with normal distributed random elements. We thus estimate the smooth function f (E, ω) for each bath from a system of size N = 24 and evaluate the correlations stemming from Eq. (5), 44,53], from which we can determine the current expected from a large system. In Fig. 2(b) we thus show with a continuous black line the value from this ETH-based approach for a system with N = 2000. Indeed, the microcanonical currents converge towards ETH predictions showing that the emergence of a smooth dependence of the current versus energy can be deduced from ETH.
Prethermalization and thermalization.-To address how the system goes from the presence of a (quasi-)steady current to the absence of current, we consider larger magnitudes of the coupling strength γ and we perform exact simulations. We stress here that perturbative calculations do not allow one to observe thermalization because, implicitly, they assume that one takes the thermodynamic limit first, and then the infinite-time limit. In Fig. 3(a) we depict the exact typical current versus time, for a single random typical initial condition, divided by γ 2 (because of the dependence on γ of the perturbative current) for different values of the coupling strength with N = 24. In particular we have γ = 0.05, 0.1, 0.3, 0.5 and 1.0 for darker shades of blue. As a reference, we add the microcanonical perturbative current I mic perturb as a black dot-dashed line. The typical current of a sin-gle random realization and for small γ follows closely the perturbative current before deviating from it. As γ increases, the typical current deviates earlier, and more markedly from the perturbative value. In Fig. 3(b) we consider different coupling strengths γ, and system sizes N . We show a convergence at long times of the current when plotted against γ 2 t/N . This is a consequence of the perturbative character of the local coupling for larger systems; in fact the energy of the baths is extensive while the two baths are only coupled at a single site.
When analyzing the current evolution over these time scales we can clearly recognize three regimes: a first transient regime, a regime in which one can observe (quasi-)steady currents for weak coupling (which we have shown to be typical), and a thermalization regime in which the current goes to zero. The emergence of the intermediate regime with a (quasi-)steady current is analogous to the phenomenon of prethermalization in isolated systems in which a thermalizing system with a weakly broken conserved law would relax slowly due to the underlying presence of conservation laws (the prethermalized regime), and then slowly relax toward the final states. In our case the conserved quantity is the energy of each bath, which is broken by the coupling between them. This dynamical correspondence with prethermalization is strengthened by the γ 2 relaxation rate [30,31].
Conclusions.-We studied the emergence of longlasting currents between differently prepared nonintegrable systems in analogy to the emergence of prethermalization in isolated systems. Each bath has its own thermalizing dynamics which conserves its energy. This conservation law is weakly broken due to the coupling between the baths, and this results in a typical fast relaxation towards a nonequilibrium scenario with, in the weak-coupling regime, a (quasi-)steady current. This is then followed by a slow relaxation towards thermalization. We have shown that such dynamics is typical in the sense that it is quantitatively the same for any typical state chosen from two different energy shells in the respective baths, and specifically we have shown that this current is consistent with the prediction from ETH, and with preparing the baths in microcanonical states. The convergence towards ETH also results in a smooth dependence of current with energy difference. To obtain these results, we have used both perturbative and exact numerical methods. For the latter we used a highly parallelized algorithm which allowed us to simulate spin chains of the size of 28 spins (14 spins per bath) and more. We have thus studied a model which makes it possible for the observation of nonequilibrium steady current as well as the long time thermalization. This connects the advances in equilibrium statistical mechanics with those in nonequilibrium open quantum systems.
We would like to stress here that the study of transport and thermalization properties of systems by preparing two subsystems in different states and coupling them has already been undertaken in the past years [54][55][56][57][58][59][60][61]. With this work we lay a stronger basis on the application of these approaches, highlighting that the evolution of a relatively small number of different pure-state initial conditions can lead to accurate estimates of the currents expected from the microcanonical initial conditions. In the future we will consider adding an intermediate system between the baths so as to study conditions for the emergence of a (quasi-)steady current within this intermediate system. Furthermore, it would be important to consider whether the conditions for the emergence of the (quasi-)steady current can be loosened, for example considering also integrable baths [55,[62][63][64][65][66].

I. PERTURBATIVE AND EXACT CURRENTS
In this section we derive the perturbative current expression (Eq. (5) of the main paper) similar to Ref. [1] which focused on transport problems through a weakly coupled interface. In this work, we consider two finite quantum spin chains with a total length N , where B L and B R are operators acting respectively on the left and right baths and γ is the coupling strength. We work in units for which = 1. Since the total energy is a conserved quantity for each region, the energy current operators can be defined as the rate of change of the Hamiltonian from the left region given by where we definedḂ L ≡ −iγ B L , H L . To find the current I L = Tr ρI L , we consider the Dyson series expansion of the total density operator in the interaction pictureρ(t) = e i(HL+HR)t ρe −i(HL+HR)t of the composite system H with respect to the coupling strength γ,ρ The corresponding current operator in interaction picture is thus I L (t) = e i(HL+HR)t I L e −i(HL+HR)t and I L = Tr ρ(t)I L (t) . By considering a separable initial condition where the initial density operator is the tensor product of two operators acting only on the left of right bath, ρ(0) = ρ L (0) ρ R (0), it follows from Eqs. (S3, S4) that a perturbative current expression can be found as Here C L and C R are the correlation functions defined as Up to this point, we have not yet specified the density operator ρ L and ρ R . When the density operators are such that [H L , ρ L ] = [H R , ρ R ] = 0, e.g., when the density operators are given by some equilibrium states, the correlation functions can be simplified in a time translational invariant form (in the sense that it only depends on the difference of t and τ ) One can further simplify the current expression by considering a change of variable t − τ → τ and the current can then be expressed as because − → C L (τ ) = − ← − C L (τ ) andḂ L (t)B R (t) = 0 for equilibrium states. The first remark is that the above correlation functions Eqs. (S9, S10, S11) and the current expression Eq. (S12) only apply to equilibrium states such as eigenstates, microcanonical distributions, and canonical distributions and they do not apply to typical initial states. Typical initial states in general do not commute with the Hamiltonian H L or H R . As a consequence, the time translational invariance condition does not hold and the free evolution terṁ B L (t)B R (t) does not vanish. Thus, in Fig. 1 in the main paper, the perturbative results for typical initial conditions are obtained via Eq. (S5) and those for equilibrium states are obtained via Eq. (S12). The second remark is that these perturbative current expressions assume that the environments are not affected by the interaction, in a similar manner commonly used in perturbative quantum master equations such as Born-Markov master equation [2].
For numerical exact calculations to measure the current, we follow the current operator definition given by Eq. (S3) and the full expression of the current operator I L , with respect to the left environment Hamiltonian in Eq. (1) of the main paper, and here in Eq. (S28), is given by Correspondingly, the current operator for right bath can be defined as the energy entering the right part, that is (S14)

II. DYNAMICAL TYPICALITY AND TYPICAL CORRELATION FUNCTIONS
The perturbative current expression given by Eq. (S12) allows us to extend the notion of typicality to currents via the correlation functions. We first show the typicality of the correlation functions in Eqs. (S6, S7, S8) in a similar manner for isolated systems with respect to some typical pure states ρ typ drawn from a given ρ mic . If follows that one can show that the typical behavior will be inherited by the perturbative expression for the currents given in Eq. (S5).
More in detail, if dynamical typicality holds, then one can write that for vast majority of the typical initial states, where the superscripts refer to typical-state initial condition and microcanonical initial condition. One can examine the expectation value of these typical correlation functions and it can be shown that they are equivalent to the microcanonical ones, i.e., E [C typ ] = C mic . For example, considering − → C typ L (t, τ ) given by Eq. (S6) and some typical initial states, ρ typ L ij = c * i c j defined in Eq. (3) of the main paper, and C typ R . One can thus infer that the typicality for currents holds, at least in the weak coupling limit for which perturbative expressions are accurate, as the correlation functions for left and right environments are independent to each other. Hence, one can expect that Furthermore, as shown in Fig. 2 of the main paper, and as can be seen in Fig. S5, we have checked numerically that the variance of the currents is bounded and decreases with the system size.

A. Correlation functions via eigenstate thermalization hypothesis
Since each of the environment are nonintegrable and thus follows the eigenstate thermalization hypothesis, one can evaluate the correlation functions − → C L (τ ), ← − C L (τ ), and C R (τ ) in a manner similar to Refs. [3][4][5] and thus a perturbative current expression. Considering the ETH ansatz for the operator B, where B ij are the elements of B in the energy basis, E i an energy eigenvalue, E = (E i + E j )/2, ω = E j − E i and R ij is a matrix whose elements follow a normal distribution with zero mean and unit variance. We consider the correlation function for a single eigenstates with eigenvalue E i − → C eig L (τ ) which, using Eq. (S9), can be written as It follows that the corresponding correlation expression via ETH is By taking the continuum limit where To obtain the last expression in Eq. (S24) we have defined and we have also set |R ij | 2 = 1 because this is the value of the variance of the elements of the matrix R ij . Finally, we can safely drop the ω/2 shift in f L , and write f L (E i , ω), as f L varies mildly with respect to E i and quickly decays with respect to ω. This results in Following the same steps, we can also write a simpler to evaluate expression for C ETH R (τ ), which is

B. Bath Hamiltonian and local operators properties: diagonal and off-diagonal
In this work we have chosen a nonintegrable local Hamiltonian which has no conserved quantities except for the energy. In the specific, although the results do not depend on this particular choice of Hamiltonian, we have considered the following one for the left bath and the same for the right bath's Hamiltonian H L , except that the site number n goes from N L + 1 to N L + N R . Dynamical typicality requires common initial expectation values for the observables of interest, a property that readily occurs in systems following ETH. In Fig. S1(a) we show i|σ x NL |i L versus the energy per spin for the eigenstate |i L which, as the system size increases, follows a smooth curve as expected from ETH.
For the off-diagonal elements, in Fig. S1(b) we show the off-diagonal elements i|σ x NL |j L multiplied by the square root of the Hilbert space size 2 NL/2 , versus the energy difference ω. To make this plot more relevant to the main paper, here we have considered a window of energies for the left bath (we can transfer these results to the right bath since it follows the same Hamiltonian) Ξ L = N L [ε L − δ L /2, ε L + δ L /2] where ε L is an average energy per particle, while δ L indicates the linear scaling factor of the width of the shell. Furthermore, the off-diagonal elements i|σ x NL |j L are selected in the energy window such that |(E i + E j )/(2N L ε L ) − 1| < N L δ L [6,7]. Also in Fig. S1(b) we see the expected result from a system which follows ETH, i.e. that the (rescaled) off-diagonal elements tend to a smooth curve f L (E i , ω).

C. Scaling analysis
In order to evaluate Eqs. (S26, S27) we need to further analyze the scaling of the different quantities included. The energy windows Ξ L/R considered shift both in average value and in width linearly with the system size. This is due to the extensive nature of the Hamiltonians H L and H R used. One natural numerical check that one can do is to verify the scaling of the maximum and minimum value of the energy, respectively E max and E min . These quantities are shown in Fig. S2 (a, b), which clearly depicts a linear scaling. In Fig. S2 (c), after fitting the many-body density of states i δ(E − E i ) with a Gaussian function, we show that the standard deviation of such a Gaussian-like density of states scales as a function √ N L , which confirms that for this local Hamiltonian, the density of states indeed scales as the square root of the system size. This allows us to find the ratio G(E i , ω) = exp [S(E i + ω)]/ exp [S(E i + ω/2)] defined in Eq. (S25), which is depicted in Fig. S2(d) for different sizes. Note that there exists a fixed point for ω = 0 and the lines converge towards a single line near it.
In Fig. S3, we show the function f (E, ω) for finite size systems with fixed energy per spin. The function f (E, ω) is extracted by considering the moving average of the off-diagonal elements. Fig. S3(a) shows the off-diagonal elements which are the same as Fig. S1(b) with the full range of ω. Fig. S3(b) is the moving averages of the off-diagonal elements of panel (a) and we interpolate them as the envelope function f (E, ω) for a fixed energy per spin. The moving average is computed by considering averaged absolute values of the off-diagonal elements in a window δω/E L [7]. Different system sizes gives a consistent f (E, ω) function profile upon a scaling with 2 NL/2 given by the prefactor in the ETH ansatz in Eq. (S21).
An important consequence from the function f profile, as demonstrated in Fig. S3 (a, b), is that the contribution near ω = 0 is the most significant one. This means that the fixed point in Fig. S2 contributes most significantly to the many-body density of states ratio and thus the correlation integral is also invariant under a scaling with respect to system size N .

IV. EQUIVALENCE OF CURRENTS FOR LEFT AND RIGHT BATHS
In the main paper, we show the emergence of long-lasting steady currents, however, the data we show refers solely to the energy current leaving the left bath. For completeness, here we show that the current entering the right region, see Eq. (S14), coincides with the one leaving the right bath (up to small deviations due to energy being temporarily accumulated at the interface). This is shown in Fig. S4 where the left current is represented by continuous lines and the right current by dashed ones. In particular, we show that the correspondence of the two currents is verified both in weak (orange curves) and strong (blue curves) coupling regimes. This is due to the boundness of the local interface operator V in the overall Hamiltonian Eq. (S2).

V. EXAMPLES OF CURRENTS FROM TYPICAL INITIAL CONDITIONS
To give further insight into the dynamics of different typical initial conditions |ψ typ , in Fig. S5 we plot 500 currents versus time, each corresponding to a different typical initial condition (see Eq. (3) of the main paper). These currents are computed numerically in an exact manner, without perturbative approximations. Of these 500 initial conditions, we single out a few random ones in different colors to allow the reader to appreciate their dynamics. The black line shows the average of all initial conditions considered. Each panel corresponds to different system sizes, specifically N = 22, 24, 26, 28 from (a) to (d). We have randomly selected five typical initial states to demonstrate the current dynamics as indicated by the markers with different colors and the background grey shades includes all the realizations.

VI. DETAILS OF THE OPTIMIZED TIME EVOLUTION ALGORITHM
Since the time evolution algorithm used in this work is vital for our large-scale simulations, the details of the numerical implementation are presented here as we believe that they could be useful for future works. The evolution for a short time δt, i.e. U (δt) = exp(−iHδt), can be evaluated using two-site evolution operators U σ n σ m σnσm after Suzuki-Trotter decomposition [8,9]. Here, with σ n we indicate the spin configuration at site n, hence a state |ψ can be represented using the basis state vectors |σ 1 , σ 2 , . . . , σ N with the coefficients c σ1, σ2,...,σ N as |ψ = σ1, σ2,...,σN c σ1, σ2,...,σ N |σ 1 , σ 2 , . . . , σ N . (S29) Applying the two-site evolution operator U  S6. Schemes to aggregate a fixed size matrix for the inner for loop occurred in the two-body gate operation U for case (a) both the two tensor indices are within the first 5 sites, (b) both tensor indices been contracted are beyond the 5-th site. The tensor legs inside the blue dashed box mean that they should be placed in the inner loop, while the rest should be placed in the outer loop.
One can thus think of c σ1σ2···σ i ···σ j ···σ N as a multi-leg tensor, with each leg taking the possible values σ i , as depicted in Fig. S6 . At this point one can make two considerations that are the key to significantly speed up the computations: 1) the memory access complexity and computational complexity are essentially the same (both of the order O(2 N )), which is unfriendly to modern computers and 2) the memory access to the quantum state is in general random, which will result in a further slowing down. We remark that the memory complexity is here manageable, even on an ordinary laptop, for system sizes of N = 28 or more because the Suzuki-Trotter algorithm requires the storage of one single vector only. Hence the memory required for N = 28 is 8 gigabytes for complex double precision. More importantly, we use a fused tensor permutation and contraction technique, which is demonstrated in Fig. S6. Eq. (S30) then means to contract the two tensor indices of the gate operation U σ i σ j σiσj with the two corresponding indices of the quantum state, which are marked in orange in Fig. S6. This could be straightforwardly implemented as an outer for loop over the black (uncontracted) legs followed by an inner for loop over the orange (contracted) legs where a small matrix-vector multiplication should be done. However, a naive implementation of this approach would in general be very inefficient due to the low computational density and random memory access. The central idea of our technique is to aggregate a fixed number of indices, which is determined by the cache size of the computer. In this case, we choose 5 indices as a compromise in using the space in the cache and considering the time needed to load and unload the cache. By doing so, the inner loop has an increased computational density because there we implement a multiplication between a 4 × 4 matrix, which is U σ i σ j σiσj , and 4 × 8 matrix which represents the values of the coefficients c σ1σ2···σ i ···σ j ···σ N when varying only 5 indices. One further important point, which we depict in Fig. S6, is that the performance is improved if one considers always sites 1 to 3 to be in the inner loop, while the other two sites will depend on where the two-site unitary operator U σ i σ j σiσj is acting upon. Two examples of the five sites used are enclosed by a dashed blue line in Fig. S6. On top of this, the outer for loop needed to evaluate Eq. (S30) can be perfectly parallelized since the inner loop only accesses non-overlapping data of the quantum state.