{'fraction_non_alphanumeric': 0.03207814923599881, 'fraction_numerical': 0.01062671198970331, 'mean_word_length': 4.754272692745917, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analyzing non-equilibrium quantum states through snapshots with artificial neural networks

Current quantum simulation experiments are starting to explore non-equilibrium many-body dynamics in previously inaccessible regimes in terms of system sizes and time scales. Therefore, the question emerges which observables are best suited to study the dynamics in such quantum many-body systems. Using machine learning techniques, we investigate the dynamics and in particular the thermalization behavior of an interacting quantum system which undergoes a dynamical phase transition from an ergodic to a many-body localized phase. A neural network is trained to distinguish non-equilibrium from thermal equilibrium data, and the network performance serves as a probe for the thermalization behavior of the system. We test our methods with experimental snapshots of ultracold atoms taken with a quantum gas microscope. Our results provide a path to analyze highly-entangled large-scale quantum states for system sizes where numerical calculations of conventional observables become challenging.

-After a global quench in a thermalizing system, local observables approach a value which corresponds to their expectation value in a typical microcanonical many-body eigenstate of the system [1][2][3]. Depending on the properties of the system and the initial state, the path to thermal equilibrium can vary. For example, conserved quantities can slow down the equilibration process [4][5][6] or a quasi-stationary prethermal state can form, which exhibits properties different from the true thermal equilibrium state [7]. Quantum simulation experiments can enable the observation of the time-evolution of a quantum many-body system starting from a non-equilibrium state with almost perfect isolation from the environment. In the past decade, a variety of non-equilibrium phenomena has been observed with examples ranging from exotic phases realized through Floquet driving [8][9][10] to many-body localization [11] and prethermalization [12]. In many cases, theory can provide a clear prediction which observables should be studied, such as a given order parameter for a well-known phase transition. For some problems, however, it is not as clear which observable to look at, and by making a choice for one specific quantity, valuable information might be discarded. In many platforms with microscopic readout, Fock space snapshots of the quantum many-body state are the measured data set. Fock space snapshots provide a wealth of information about the quantum many-body state by providing access to both local observables and non-local, high-order correlations.
In order to address the challenge of finding suitable observables, artificial neural networks have recently emerged as a valuable tool in quantum many-body physics [13][14][15][16][17], and in nonequilibrium statistical me-chanics [18]. Previous machine learning approaches to study non-equilibrium systems have focused on quantities such as the entanglement spectrum [19][20][21] or full eigenstates [22], which are, however, experimentally inaccessible. In this work we study the dynamics of an interacting quantum many-body system in terms of experimental Fock space snapshots with the help of neural networks, Fig. 1a). We find this analysis to have two main advantages: (i) these snapshots are directly measured in many quantum simulation platforms, and large numbers of snapshots can be routinely obtained. (ii) Raw data is used, where no analysis for specific quantities has taken place and all available information can be used without any bias. We consider the one-dimensional Bose-Hubbard model (1) Here,â ( †) i annihilates (creates) a boson on site i and n i =â † iâ i is the particle number operator. The first term corresponds to hopping between neighboring sites, the second term is the interaction, here fixed at U/J = 2.9, and the last term is the quasi-periodic potential mimicking on-site disorder with amplitude W , which can be created in a cold atom setup with an incommensurate lattice as h i = cos(2πβi + φ). In this work, we consider 1/β = 1.618. This system exhibits a many-body localized (MBL) phase, where thermalization breaks down as the disorder strength is increased beyond a critical value. The transition from an ergodic to a many-body localized phase is fundamentally different from the well-studied  The Bose-Hubbard model with a quasi-periodic disorder potential exhibits a many-body localized (MBL) phase, where thermalization breaks down, as the disorder strength is increased beyond a critical value. a) We study the dynamics of the system after a quench for different disorder strengths by evaluating snapshots from a quantum gas microscope with neural networks. b) A neural network is trained to distinguish exact diagonalization snapshots at W/J = 0.3 and W/J = 11 for U/J = 2.9 and a system with 8 and 12 sites at time tJ = 100 after a global quench. After the training process is finished, snapshots at intermediate values of the disorder strength are used as input. The plot shows the resulting classification for numerical data (shaded band) as well as experimental snapshots (symbols). As the system size is increased, the fraction of snapshots classified as MBL begins to increase at larger values of W , indicating the transition in the finite size system. The accuracies are averaged over two independent runs and the errors denote one s.e.m. case of equilibrium phase transitions, as it describes a non-equilibrium setting [23][24][25][26][27][28][29][30][31]. Finding the transition point is numerically challenging, because it is usually obtained from entanglement properties or the level statistics, which can only be obtained for small system sizes where full diagonalization of the Hamiltonian is possible. Here, we focus on Fock space snapshots of the manybody quantum state as input data, which are the direct output of quantum gas microscopy experiments and thus experimentally readily accessible for the systems of interest. This approach has the advantage that significantly bigger system sizes can be reached experimentally. We consider the dynamics of two one-dimensional sys-tems of 8 and 12 sites, which are initialized in a Mottinsulating state with exactly one particle per site. In Fig. 1, we first train the network to distinguish snapshots of the many-body quantum state, obtained from exact diagonalization calculations, for low (W/J = 0.3) and high (W/J = 11.0) disorder strength for an interaction strength of U/J = 2.9 in the comparatively long-time limit at time tJ = 100. We average over ten different disorder realizations, obtained by varying the phase φ in the potential. After the network has learned to label the extremal cases correctly with sufficiently high accuracy (> 90%), we input snapshots for intermediate values of the disorder strength. After training the neural network on numerically simulated snapshots, we use experimental data as input, where each snapshot stems from a different disorder realization. As output, for each disorder strength we obtain the fraction of snapshots labeled as many-body localized and thermalizing, see Fig. 1. Based on these results, we conclude that the many-body localization transition is located within the range of W/J ≈ 4 − 8 with strong finite-size drifts. This result is in agreement with previous experiments [32,33], which considered conventional observables such as the local entropy. Notably, the local entropy exhibits volume law scaling both in the thermal and the MBL phase and is thus by itself not sufficient to locate the transition without exact numerics [32]. Our results, in contrast, are able to distinguish the two phases without any theoretical input, which suggests that the network learned a more suitable observable to distinguish the two phases. In [34], we show the level statistics for system sizes L = 6, 7, 8 for comparison. Similar to the machine learning analysis of a disordered spin chain based on the entanglement spectrum in [19], the transition found by the neural network is as sharp as the level statistics, but exhibits a small shift to larger disorder strengths.
While we have only compared two extremal disorder strengths in the long-time limit, the full dynamics of the system contain much more information. We proceed by analyzing the time-and disorder-strength dependence of the system after the global quench.
Learning thermalization.-We now investigate the system's approach to thermal equilibrium by comparing each time step to a thermal state of the same Hamiltonian. The performance of the network in distinguishing dynamics from equilibrium can then be used as a probe of thermalization. In order to compare the time evolved state to thermal equilibrium, all conserved quantities of the model should be considered [3]. In our experiment, both the energy density and the particle number are conserved during the many-body evolution. The energy density of the initial state is matched by choosing the temperature of the thermal state accordingly. We take the conservation of the total particle number into account by calculating the thermal state within a fixed particle number sector. We numerically generate snapshots from such a state in thermal equilibrium as well as from the time-evolved state for each time step under consideration. For each time step, we train the network to label the snapshots from the thermal equilibrium distribution as equilibrium, and the snapshots from the numerically time-evolved initial state as dynamics. The neural network parameters optimized for each time step seperately. We then test the network's performance by inputting experimental data with different evolution times. In Fig. 2a) the resulting classification into the categories dynamics versus equilibrium is shown as a function of time. Here, we average over 12 different disorder realizations and take snapshots at the corresponding effective temperatures. For small W/J, the system thermalizes comparably fast: for times tJ > 10, the network reaches an accuracy of 50%, equivalent to guessing between the two classes. This means the network fails to distinguish snapshots from the time-evolved state from the corresponding thermal state. For high values of W/J the system fails to thermalize on the time-scales accessed here, and the network is able to distinguish the current timestep from the thermal equilibrium state with a high accuracy.
Using an interpretable network architecture [45], we find that for intermediate disorder strengths, higher order correlations play a role in the classification task, see [34]. We study the long time limit at tJ = 100 for a range of values of the disorder strength. As shown in Fig. 2c), the fraction of snapshots classified as dynamics rises strongly between W/J ≈ 4 and W/J ≈ 8 and reaches values close to 1, indicating that the system has not reached thermal equilibrium. We benchmark our experimental results by testing the network with theoretical snapshots not used during training and find good agreement throughout the range of the covered parameters. This procedure has the advantage that the features used to make the classification can vary for different time steps and the network specifically searches for differences between the current time and thermal equilibrium. It is therefore in principle capable of identifying specific observables that have not yet reached their thermal equilibrium value and thus find, for example, (almost-) conserved quantities. Indeed, with this method we find deviations from thermal equilibrium already in the range of W/J ≈ 2 − 5, in contrast to the results from the classification scheme in Fig. 1b). This indicates an improved sensitivity of our method. Here we consider a system which exhibits a transition from thermalizing behavior to many-body localization, which constitutes a canonical example in the study of non-equilibrium phenomena. Note, however, that our scheme is not limited to the system considered here and can be applied to a variety of models. This method also allows to detect, for example, prethermal behavior and the existence of conserved quantities that keep their value during the dynamics and therefore never reach a generic thermal equilibrium value. Another canonical model to study equilibration behavior is the transverse field Ising model, which has an extensive number of conserved quantities. In [34], we show that a neural network performs significantly worse in distinguishing the time-evolved state from an approximative generalized Gibbs ensemble, where a few conserved quantities are taken into account, than the simple thermal state discussed above, where only the energy density is considered. This highlights the capability of our approach to identify conserved quantities, which can drastically alter the thermalization process. Our method comes at the expense that one needs snapshots from the thermal density matrix for training, whichespecially in the case of a non-thermalizing phase such as MBL -may need to be generated numerically. In the following, we overcome this limitation by analyzing the transition in the dynamics with an unsupervised scheme that, in principle, does not rely on theory data.
Confusion learning.-Several unsupervised learning schemes that use the network performance to probe whether and where a phase transition or more general, a qualitative change in the data, exists have been proposed [37][38][39]. Here, we adapt a scheme termed "confusion learning" introduced in Ref. [37]. In brief, the scheme works as follows: We have a dataset of snapshots for values of the disorder strength 0.3 ≤ W/J ≤ 11.0. The goal is to test whether a value W * exists at which the data changes qualitatively. We start with a guess for W * and label all snapshots for W ≤ W * as phase A and correspondingly all snapshots with W > W * as phase B. Assuming the snapshots are qualitatively different for W ≤ W * as compared to W > W * , the network should achieve a high accuracy in assigning the correct labels. However, if there is no qualitative change at the W * under consideration, there will be confusion about the correct labels and the accuracy will thus be lower. Therefore, if there is a qualitative change in the data, the accuracy as a function of W * will be maximal if W * corresponds to the transition point. Trivially, the test accuracy is expected to approach unity when the guessed W * corresponds to the minimum or maximum value of W , because all data are labelled equally and no confusion occurs. In total, the presence of a critical point is therefore signalled by a characteristic W -shape of the test accuracy as a function of the control parameter. We train the neural network with numerical snapshots in the long-time limit (tJ = 100) in order to test for the presence of a phase transition. Subsequently, we use experimental data as input to the network, Fig. 3a). The data shows the onset of a maximum around W * /J = 7, indicating the presence of a critical point in agreement with Fig. 1b). The contrast in the W -shape achieved here is comparable to the signal seen for a spin model in [37], where instead of snapshots the entanglement spectrum is used as input to the neural network. In order to isolate the signal of the phase transition from the trivial part of the W -shape, we subtract the accuracy obtained when training on randomly labeled data. The resulting difference, shown in the inset of Fig. 3a), exhibits a clear peak at W * /J = 7, that indicates the transition between the different dynamical phases. We also check with theoretical snapshots not used during training and find qualitatively similar behaviour. We attribute the slight deviation in the maximum to the coarse resolution in the disorder strength for the experimental data. Since differences in the thermalization behavior only present themselves in the course of the dynamics, we expect the phase transition to remain hidden at short evolution times. In order to reveal this effect, we perform the same method with theoretical snapshots at different evolution times. In Fig. 3b), the resulting accuracy achieved by the network is shown as a function of W * . These results have several advantages compared to the previous methods: as opposed to Fig. 1b), we do not a priori assume that there is a transition. Moreover, we specifically train the network to find differences between the snapshots at all available values of the disorder strength, thus avoiding bias from the choice of training data.
Summary and Outlook.-In this work, we used machine learning techniques to study the non-equilibrium dynamics after a global quench in the one-dimensional Bose-Hubbard model with a quasi-periodic disorder potential. We used supervised as well as unsupervised machine learning methods to probe for a qualitative change in experimental snapshots as the disorder strength is tuned. Comparing the results for systems with 8 and 12 sites, we find that the critical value of the disorder strength increases with the system size, proving the need for methods applicable in large -experimentally accessible -systems. In contrast to standard tools to locate the MBL transition, the methods used here can be directly applied to experimental data taken with a quantum gas microscope and are not limited to small system sizes. We furthermore studied the approach to thermal equilibrium -or lack thereof -by training a neural network to distinguish snapshots from the current time step from snapshots from a thermal ensemble at the same energy and particle density. The accuracy achieved by the network indicates how non-thermal the time-dependent quantum many-body state is. An exciting future research direction consists of applying the same scheme to identify conserved or almostconserved quantities in experimentally accessible data, for example by using a generalized Gibbs ensemble for comparison. Apart from the concrete system studied here, it would be interesting to consider other models and phenomena, for example quantum scars [40,41] and Hilbert space fragmentation [42][43][44]. In order to gain additional physical insights, interpretability is an extremely important direction for future work and it would be interesting to study which observables the network uses to make the classifications considered here [45], and how those observables change during the time evolution of the many-body system.

SUPPLEMENTARY INFORMATION TRANSITION
Upon increasing the disorder strength, the level statistics of the Hamiltonian evolves from the Gaussianorthogonal ensemble to Poisson statistics as the system enters the MBL phase [25,26]. We consider the level spacings where E (n) φ is the n-th eigenenergy of Hamiltonian (1) in the main text with disorder realization given by φ. The ratio of adjacent gaps is then given as (3) For a given system size L, we fix the particle density to one particle per site and for each disorder strength W/J obtain the average value of this ratio over 30 disorder realizations, i.e. different values of φ. In Fig. 4, the resulting average value of the ratio of adjacent energy gaps r is shown as a function of the disorder strength for system sizes L = 6, 7, 8. The shift of the drop in r to larger values of W/J agrees well with the results presented in Fig.1b) of the main text. In order to relate to previous work, in particular Ref. [33], we directly evaluate observables from the snapshots and calculate the transport distance ∆x, defined as with FIG. 5. Observables. The transport distance ∆x and onsite fluctuations F, see text, are evaluated from the same snapshots. Shaded bands correspond to exact diagonalization snapshots, symbols are based on experimental data. We simultaneously evaluate snapshots from ten different disorder realizations in the numerical data. In the experimental data, each snapshot is from a different disorder realization. The errors denote one s.e.m. and the on-site fluctuations F, defined as in Fig. 5. Comparing the output of the neural network in Fig.1b) of the main text with the transport distance ∆x shows a similar behavior, from which one might conjecture that the network uses a similar observable to make the distinction. Note that with this approach, we are able to make a quantitative prediction on the basis of single or few snapshots, for which the observables shown in Fig. 5 are not clearly converged to their average value.
Confusion learning: experimental data In Fig. 6, we show the same analysis as presented in Fig.3 of the main text, but using a network that has been solely trained on experimental data. A clear "W"shape does not emerge. The experimental result agrees qualitatively with the result using exact numerical data, where we only used data at the same values of disorder strength as available from the experiment. We hence attribute the lack of a clear "W"-shape to the rough grid of values of the disorder strength -consisting of eight values between W/J = 2 and W/J = 10.

Unsupervised learning of the transition
We use the unsupervised scheme introduced in [38,39] to locate the transition to the many-body localized phase as a function of the disorder strength W . In Fig. 7  Snapshots of the many-body quantum state of a system with 6, 8, and 12 sites, U/J = 2.9, and various disorder strengths W/J are analyzed using the unsupervised learning scheme introduced in [38,39]. A neural network is trained to label a given snapshot with the corresponding value of W/J. If a qualitative change in the data occurs, the derivative δW pred /δW label will exhibit a maximum. The plot shows W pred as a function of W label . Shaded bands are exact diagonalization and symbols experimental data. The results are averaged over 25 independent runs and the error bars corresponds to one s.e.m. . text.

LEARNING THERMALIZATION
In order to compare the time evolved state to thermal equilibrium, all conserved quantities of the model should be taken into account. The energy density of the initial state can be matched by choosing the temperature of the thermal state accordingly. In particular, the energy density of the initial state |ψ 0 is given by E i = ψ 0 |Ĥ|ψ 0 . The effective temperature T eff is then determined such that the density matrix of the system,ρ β = 1 Z exp(−β effĤ ), with the inverse temperature β eff = 1/T eff and Z = tr(exp(−β effĤ )) fufills The energy density E(β) = tr(Ĥρ β ) is calculated for a range of values β until the effective temperature is determined such that Eq. 7 is fulfilled. Due to the disorder potential, this effective temperature varies for different values of φ, where φ determines the disorder realization. In Fig. 8, the effective inverse temperature is shown as a function of φ for a system with L = 6 sites at unity filling for interaction strength U/J = 2.9 and disorder strength W/J = 4.5.
Interpretability: higher-order correlation functions Following the same approach as in main text Fig. 2, but using the correlator convolutional neural network (CCNN) architecture, we can gain insights into the information used to solve the classification task. In particular, the order of correlation considered enters as a  9. Accuracy as a function of order of correlations. Accuracy obtained for comparing snapshots from time step tJ = 7.3 to a thermal ensemble using the correlator convolutional neural network (CCNN) architecture [45]. The order of correlation functions used by the neural network serves as a hyper parameter in the CCNN architecture. The top plot shows the obtained accuracy for different disorder strengths for correlation orders 2,3,4, and 5. The bottom plot shows the difference to the accuracy obtained when taking second order correlation functions into account to allow for a better comparison of the results between the different disorder strengths.
hyper parameter of the network architecture. The network can thus be trained to distinguish dynamics from thermal equilibrium taking into account for example only correlations up to second order. Note that given the Fock space snapshots we consider here, all correlations considered are density-density correlations. In Fig. 9, we compare the accuracies obtained by the network when taking into account correlations up to second, third, fourth, and fifth order for different disorder strength, for comparing the time step tJ = 7.3 to thermal equilibrium. Since the overall scale increases significantly with increasing disorder strength, as the state becomes less and less thermal, we show in Fig. 9 bottom the accuracies obtained for correlations of order 3 − 5 with the accuracy for order 2 subtracted. For low and high disorder, taking into account correlation functions of order higher than two does not significantly increase the accuracy obtained by the CCNN. However, for intermediate disorder strengths, higher order correlations play a role in the classification task and the accuracy improves by up to 5% when considering higher order correlators. This result is in accordance with Ref. [33], which showed sizable higher order correlations in the critical regime at intermediate disorder strength.

Transverse field Ising model and generalized Gibbs ensemble
The transverse field Ising model, has an extensive number of local conservation laws, which are known to be [35,36] I (n,+) = −J j (S xx j,j+n +S yy j,j+n+2 )+h(S xx j,j+n−1 +S yy j,j+n−1 )  10. Thermalization and conserved quantities. We consider the dynamics after initializing the system in a product state in the transverse field Ising model (TFIM). In the TFIM, there is an extensive number of conserved quantities, which have to be considered to correctly describe the longtime limit. Top: If only the conserved energy is taken into account (red), the long-time limit is still comparably easy to distinguish from the thermal density matrix. The performance of the network drops significantly if the first (i.e., most local) two conserved quantities are additionally taken into account in an approximation of the generalized Gibbs ensemble (green). Bottom: Difference between the accuracies shown in the top plot.
with I (1,+) = H and We numerically simulate the dynamics for J = 1, h = 2 starting from the initial product state with θ = 3.0 using exact diagonalization for a system of size L = 12. We then train a neural network to distinguish snapshots sampled from the time-evolved state  11. Learning thermalization. The system is initialized in a uniform state (|111111111111 ) and the ensuing time evolution is investigated. In each time step, the neural network is trained to distinguish snapshots from the current time step from snapshots from the long-time limit. The plots show the resulting accuracy for W/J = 6.4, U/J = 2.9.
from a thermal ensemble. For the thermal ensemble, we consider two different choices: • a thermal ensemble ρ = 1/Z exp(−βH), where the inverse temperature β is determined by the energy ψ 0 |H|ψ 0 as discussed above • an approximation to the generalized Gibbs ensemble, ρ GGE = 1/Z GGE exp(−βH − λ 2 I 2,+ − λ 3 I 3,+ ) with Lagrange multipliers λ 2,3 determined in the same way for the conserved quantities defined above, with Z GGE the corresponding partition sum for the approximation of the GGE, such that trρ GGE = 1. We determine the parameters to match the expectation values in the initial state, yielding β = 0.304, λ 2 = 0.336, and λ 3 = 0.1.
As shown in Fig. 10, the performance of the neural network in distinguishing dynamics from thermal/generalized Gibbs ensemble is significantly worse in the latter case, indicating that the long-time dynamics is better described by the GGE. This emphasizes the importance in taking conserved quantities into account in order to correctly describe the equilibration behavior, and in particular the ability of a neural network to capture the differences between the time evolved state and the standard thermal ensemble.

Distinguish from long-time limit
In order to study the dynamics of the quantum manybody system, we here compare snapshots from the current time step to the long-time limit. In a thermalizing system, the long-time limit corresponds to a thermal equilibrium state and the scheme is thus basically the same as the thermalization learning scheme introduced in the main text. This is, however, not the case for the MBL phase. In Fig. 11, the accuracy achieved on a test set not used during training is shown as a function of time. In each time step, the neural network parameters are optimized to enable the classification of snapshots into the categories current timestep versus long-time limit. This procedure has the advantage that the features used to make the classification can vary for different time steps and in particular, the network specifically searches for differences between the current time and the long-time limit. It is therefore in principle capable of identifying specific observables that have not yet reached their longtime value. In Fig. 12, the accuracy achieved on a test set not used during training is shown as a function of time when starting from the product state for W/J = 0, 4.8, 8.9. We compare the resulting accuracy to the imbalance, defined as where N s is the number of snapshots, the first sum runs over all snapshots, n s i is the occupation of site i in snapshot s, and n ref i = 0(2) for i even (odd). In all three cases, the first tunneling events cause a sharp decay in the imbalance on a time scale of one hopping time. The accuracy with which the neural network can distinguish the current time step from the long-time limit is always larger than the difference of the imbalance to its long-time limit. In the many-body localized case, Fig. 12c), the accuracy shows qualitatively the same behavior as the imbalance. However, in the critical phase, Fig. 12b), there is no fast initial decay in the accuracy and instead, it is still higher than 50%, which corresponds to its lower limit, after several hundred hopping times. Without disorder, Fig. 12a), the imbalance has almost reached its long-time value after one hopping time. The accuracy with which the network can distinguish the current time step from the long-time limit decays on a slower time-scale of about ten hopping times.