Wave function network description and Kolmogorov complexity of quantum many-body systems

Programmable quantum devices are now able to probe wave functions at unprecedented levels. This is based on the ability to project the many-body state of atom and qubit arrays onto a measurement basis which produces snapshots of the system wave function. Extracting and processing information from such observations remains, however, an open quest. One often resorts to analyzing low-order correlation functions - i.e., discarding most of the available information content. Here, we introduce wave function networks - a mathematical framework to describe wave function snapshots based on network theory. For many-body systems, these networks can become scale free - a mathematical structure that has found tremendous success in a broad set of fields, ranging from biology to epidemics to internet science. We demonstrate the potential of applying these techniques to quantum science by introducing protocols to extract the Kolmogorov complexity corresponding to the output of a quantum simulator, and implementing tools for fully scalable cross-platform certification based on similarity tests between networks. We demonstrate the emergence of scale-free networks analyzing data from Rydberg quantum simulators manipulating up to 100 atoms. We illustrate how, upon crossing a phase transition, the system complexity decreases while correlation length increases - a direct signature of build up of universal behavior in data space. Comparing experiments with numerical simulations, we achieve cross-certification at the wave-function level up to timescales of 4 $\mu$ s with a confidence level of 90%, and determine experimental calibration intervals with unprecedented accuracy. Our framework is generically applicable to the output of quantum computers and simulators with in situ access to the system wave function, and requires probing accuracy and repetition rates accessible to most currently available platforms.

Programmable quantum devices are now able to probe wave functions at unprecedented levels.This is based on the ability to project the many-body state of atom and qubit arrays onto a measurement basis which produces snapshots of the system wave function.Extracting and processing information from such observations remains, however, an open quest.One often resorts to analyzing low-order correlation functions -that is, discarding most of the available information content.Here, we introduce wave function networks -a mathematical framework to describe wave function snapshots based on network theory.For many-body systems, these networks can become scale free -a mathematical structure that has found tremendous success and applications in a broad set of fields, ranging from biology to epidemics to internet science.We demonstrate the potential of applying these techniques to quantum science by introducing protocols to extract the Kolmogorov complexity corresponding to the output of a quantum simulator, and implementing tools for fully scalable cross-platform certification based on similarity tests between networks.We demonstrate the emergence of scale-free networks analyzing experimental data obtained with a Rydberg quantum simulator manipulating up to 100 atoms.Our approach illustrates how, upon crossing a phase transition, the simulator complexity decreases while correlation length increases -a direct signature of build up of universal behavior in data space.Comparing experiments with numerical simulations, we achieve cross-certification at the wave-function level up to timescales of 4µs with a confidence level of 90%, and determine experimental calibration intervals with unprecedented accuracy.Our framework is generically applicable to the output of quantum computers and simulators with in situ access to the system wave function, and requires probing accuracy and repetition rates accessible to most currently available platforms.

I. INTRODUCTION
Harnessing and probing many-body systems at the single particle/qubit level are hallmark features of presentday quantum simulators and computers [1][2][3][4].One of the most drastic demonstrations of these tools is the possibility of taking a large number of 'photos' of a manybody system, obtained via projective measurements of the full many-body wave function.While this flood of available observations could be seen as a blessing, it immediately encounters practical as well as conceptual challenges: how can this large amount of information be processed, without a priori discarding some (in the data science language, before performing a dimensional reduction)?What can one learn, that is, e.g., not available by utilizing low-order correlation functions?Answering these questions requires a structured, mathematical understanding of the experimental wave function snapshots, that addresses the information limbo between traditional many-body theory based on few-points correlation functions [5], and full-fledged -but experimentally limited to few-particle systems -tomographic methods [6].
Here, we develop a theoretical framework to characterize and classify experimentally accessible collections of wave-function snapshots utilizing network theory, that is scalable and allows one to retain all available information.The backbone of our method is a mapping between collections of wave-function snapshots, and a 'wave-function' network, schematically depicted in Fig. 1, that is applicable to spin-, bosonic, and fermionic systems.Utilizing well-established tools in network theory we unravel several key characteristics of the underlying quantum wave function, that are inaccessible by conventional means.
The pivotal finding is that the resulting quantum wave function networks can become scale free -a mathematical structure that has found wide-spread application in several fields, ranging from power distribution and internet networks to epidemics [7][8][9].We demonstrate this property using experimental snapshots obtained on a Rydberg quantum simulator operating with more than 100 atoms [2,10] and with large-scale numerical simulations using neural quantum states [11,12].We then argue about its generic applicability to state preparation protocols, and discuss how other types of networks -Erdos-Renyi [13] -can instead emerge if the resulting dynamics describes uncorrelated states.In terms of observables, required resources and applicability regimes, our approach is complementary to other methods aimed at fully characterizing quantum states via snapshots such as those based on classical shadows [14], randomized measurements [15][16][17], and chaotic dynamics [18,19].Its main distinctive features, that we elaborate upon below, are direct interpretability and straightforward scalability for strongly correlated, low-temperature states.
The correspondence between quantum simulator outputs and conventional network theory immediately enables a transfer of methods and concepts from previously disconnected fields.We leverage this connection to address two challenges in the field of quantum simulation.Firstly, we show that we are able to characterize the complexity of the quantum simulator output by determining its Kolmogorov complexity -the accepted absolute measure of information content of finite objects [20,21], that quantifies the (in-)compressibility of the quantum wave function information as contained in the snapshots.This allows us to demonstrate the emergence of critical behavior at the level of information complexity, directly probing at the wave function level the emergent simplicity dictated by renormalization group theory.
Secondly, we introduce a method to perform crossplatform verification of quantum simulators [22,23].The method is based on the full network information without the need to perform an exponentially increasing number of measurements for increasing system size, which is the case for generic cross-verification based on the density matrix [22,23].By means of the Epps-Singleton test [24] we identify, with statistical significance, a time scale beyond which cross-verification falters due to experimental imperfections not covered by our theoretical description.In addition, we provide statistically rigorous bounds for previously observed time-delay effects, that demonstrate the capability of our methods to identify systematic effects that are invisible to low-order correlation functions.Beyond these two demonstrative tools, the quantum wave function networks introduced in this work provide a new generically applicable framework to probe and characterize the quantum many-body wave function accessible in a variety of atomic and solid state quantum hardware, solely requiring in situ imaging of the many-body wave function.
In this section, we describe how data sets generated by a collection of wave function snapshots can be represented by a network structure with nodes and links.For the sake of simplicity, we consider a many-body system composed of spin-1/2 degrees of freedom defined on a two-dimensional lattice: the approach can be straightforwardly generalized to continuum theories, as well as to different types of local Hilbert spaces.
Snapshot data set.-Each wave function snapshot, labeled by an index j, takes the form: where s j m is the measured value of the spin at position m.N is the total number of sites in the system, while w are the external parameters related to the snapshot -in our case the Hamiltonian couplings.Each of these configurations corresponds to a single data point embedded in a data space whose embedding dimension is N .This is depicted in Fig. 1(a) with the three examples of green, orange and blue dots.
The data set we are interested in is formed by the collection of all available snapshots: where N r is the number of available snapshots, that is, the number of realizations.The data set might, in principle, include repetitions -e.g., X l = X f for some l = f -in particular, at very small volumes.It is possible to take care of them, as we detail in [25].However, to simplify the remainder of the discussion here, we assume no repetitions are present.
From data sets to wave function networks.-We now discuss how to translate the wave function snapshot data sets into a network structure.There are two key choices that have to be made: (i) the selection of a proper metric in the embedding space, that allows one to compute distances between data points; and (ii) a criterion to activate links between data points, based solely on their distances.
The choice of a proper metric is an important aspect of the approach.Taking inspiration from recent results in the context of classical and quantum statistical mechanics models, we use the Hamming distance [26,27].Given two configurations X i , X j , such distance counts the number of spins that are aligned differently and reads The statistics of Hamming distances are related to arbitrary rank correlation functions between local degrees of freedom (i.e., s k ) [26].Hence, they are sensitive to shortrange and long-range correlations alike, which justifies their use as a similarity measure to define links between As an example, taking snapshots of a classical antiferromagnet below its critical temperature will feature the antiferromagnetic state as a hub (top cartoon), while doing so well above its critical temperature will lead to a graph with no hubs and random connections (bottom cartoon).
nodes.Specifically, we define a (geometric) network from our data sets by adopting the following procedure: 1.Each point X i in the data set represents a node.
2. If two nodes are at distance d < R, we draw a link.
3. The distance R is chosen in a way that is dependent on the number of samples taken and reflects the typical value of distances for a given set of external parameters w.In particular, we define R as where r c (i) is the distance between point X i and its c-nearest neighbor.
An essential aspect of our approach is the choice of a cutoff R that avoids overcounting isolated nodes while keeping a non-trivial network structure (i.e., in which all the nodes are not simply fully connected).In particular, we show that our main conclusions are independent of the choice of R for a certain range of values of c in Eq. (4).We discuss this issue in detail in Appendix VII A.

A. Network representation and correlations
At a naive level, one could expect that such WFNs simply reflect the intrinsic randomness of the wave function sampling -that is, they are ultimately generated by a Poissonian process.It turns out that this intuition is fundamentally incorrect.
In order to underpin the relation between network representation and correlations, we start by schematically illustrating the above procedure in Fig. 1

(a) (i) -(iv).
A graphical example of a network with spin-1/2 systems and cutoff radius equal to 1 (that is, only configurations differing by a single spin flip are connected) is depicted in Fig. 1(b): there, the black circle represents the Neel state that is connected to several other states by a single spin-flip -and, thus, minimal distance -while it is not connected to any other states.This example allows us to intuitively connect physical properties to network ones: a wave function network that carries correlations will feature 'hubs', that is, few states with many connections, and a lot of states with few connections.Conversely, a random "infinite-temperature" state will likely feature a majority of states with an intermediate number of neighbors, and won't feature either hubs or states with very few links.
The simple picture defined above is, per se, not particularly informative: however, it crucially sheds light on the classes of networks we can expect depending on how correlated the system is.This description of correlated states is reminiscent of what happens in several classes of scale-free networks: typically described by a scaling relation of the probability distribution P k associated to the number of connections k of each node (or, as is more commonly called, the degree distribution), that follows a power law distribution Such a function monotonically decreases with k, and allows us to distinguish between the majority of nodes that have few links, and the minority of those that have many links (see Fig. 1(b)).While the prominence of hubs seems to be mostly relevant to ordered states, it is in fact a property that is even more robust in the presence of very strong correlations -such as, e.g., those emerging at quantum critical points.Conversely, networks representing random states will not be scale free, and can be construed as Erdos-Renyi (ER) networks -where the probability P k of a node having k neighbors is approximately given by a Poisson distribution [13].
We emphasize that in the many-body regime, the number of snapshots, N r , available from an experiment, is typically insufficient to tomographically reconstruct the wave function, i.e., N r 2 N .The WFN construction aims at a characterization of a state that focuses solely on the most important (yet unknown) degrees of freedom in the system, and not its entire data structure.So, our method is conceptually different from tomographic methods, including those based on specific ansatze.

B. An illustrative example: quantum Ising model at equilibrium
Before discussing the experimental relevance of WFNs, we illustrate the emergence of scale-free networks in many-body systems by utilizing an example borrowed from equilibrium statistical mechanics.In Fig. 2, we show the degree distribution P k obtained via sampling the partition-function snapshots of the 2D quantum Ising model on a square lattice, with samples in the z-basis.The Hamiltonian reads: It features a quantum phase transition at g c ≈ 3.04 separating a non-correlated disordered phase (for g > g c ) from a ferromagnetically ordered state.The corresponding P k is obtained by taking snapshots of the partition function, calculated via stochastic series expansion Monte Carlo simulations [28,29] for a system of N = L × L = 8 × 8 sites, at inverse temperature β = 2L, which in our calculations was high enough to observe convergence within statistical uncertainty of energy and squared magnetization, i.e., to reach the ground state regime.Hence, the generated data sets correspond to the ground-state WF snapshots described above.Fig. 2 (a) displays the results from N r = 10 5 realizations.Deep in the paramagnetic phase, g = 5.0, there are only weak correlations: the corresponding network is very well described by a Poisson distribution with k = 1.In the correlated regime, which is also the most entangled one, close to the phase transition, the network is described by a scale-free structure.We note that such a scale-free structure is unrelated to the absence of scale at criticality (scale-free networks can still be compatible with the presence of real-space finite scales [7]).
Once the degree k of neighbors becomes a sizeable fraction of the total size of the network, we observe deviations from a scale-free profile, as expected.Above this size, the network properties are influenced by limited sampling [30].To further inspect this, we plot in the lower panel P k against k for various N r .We observe that indeed, the origin of the bending is due to the finite number of samples, and that the curves for various N r are with k = 1.As expected, in the vicinity of the critical point, the WFN becomes scale free, with α 2.4 (dashed line).For comparison we compute P k using both linear (triangles) and logarithmic (circles) histograms.Panel (b) shows P k for different values of Nr for g ≈ gc using logarithmic histograms, again with the scale free ditribution shown as a dashed line.In all cases we build the network using a cutoff R = r1 (where r1 is defined in Eq. ( 4)).
all compatible with a single power law, in this case, with exponent α 2.4.Changing the cutoff distance R used to build the WFN does not affect the power-law scaling behavior of P k for k above a certain threshold k c ; see Appendix VII A and VII B for more details.
At equilibrium, we expect the same dichotomic structure to appear generically in models that feature both weak-and strong-correlated regions.The key question we address below is, whether such structures are purely theoretical constructions, or if they can indeed be representative of the intricate dynamics taking place in quantum simulators, that are (i) off-equilibrium, open, anda key difference from simulations -(ii) inherently probed with very high but not 100% fidelity.In the paramagnetic (PM) regime, one expects a network description compatible with a Renyi-Erdos network, while in the vicinity to the antiferromagnetic (AF) region, that contains the Kibble-Zurek regime, a scale-free network structure is expected with power-law degree distribution, P k , as illustrated by the network structures.The panel (b) presents the P k vs k of the experimentally observed wave-function networks for a square lattice with L = 8.At short times, i.e., before crossing the phase transition, the distribution decreases exponentially (similar to Erdos-Renyi degree distribution with k 1, represented by the dashed lines in the graphs).At later times (t > 3 µs), we observe a power law decay over two orders of magnitude, limited only by a bending that is due to a finite value of Nr.Graph (c) shows NQS simulations of this quasi-adiabatic protocol for the same square lattice.The scale-free behavior of P k is again observed until one is sensitive to the effects of finite sampling.We note that the value of the decay exponent is α < 2, signifying very stable wave function network properties, that will be discussed later in the presence of defects.For all cases we considered WFNs with Nr = 2500 nodes.

III. EXPERIMENTAL OBSERVATION OF ERDOS-RENYI AND SCALE-FREE WAVE FUNCTION NETWORKS
A. Experimental data and analysis of the network We now discuss the network structure of quantum simulation experiments.We analyze a recent experiment, that focuses on the quasi-adiabatic state preparation of a large antiferromagnetic state using a Rydberg quantum simulator [10].This protocol plays a fundamental importance in quantum simulation and computing, and is very widely employed in atomic physics platforms.In addition, it typically features both regimes of no-correlations (short times), and of strong correlations, enabling us to test predictions based on both Erdos-Renyi and scale-free networks.Below, we summarize the main features of the experiment, that have been reported in [10].
The experiment consists of arrays of laser-cooled Rb atoms, individually trapped into optical tweezers separated by a distance a.Each atom can be considered as a pseudo-spin, the ground state being |↓ and a Rydberg state state being |↑ .Initially all the atoms are prepared in |↓ .The atoms are then laser-excited to Rydberg states via a two-photon transition, so that the effective time-dependent Hamiltonian describing the dynamics reads: with n i = (σ z i + 1)/2, and σ α i Pauli matrices at the site i.Here, we have that J ij = C 6 /r 6 ij as the atoms interact via the van der Waals interaction.This quantum spin model exhibits both paramagnetic and antiferromagnetic phases in its ground state, for a schematic phase diagram see Fig. 3(a).
In the experiment a dynamical process has been implemented which, upon varying slowly Ω(t) and δ(t) over time, transforms an initial paramagnetic state into an antiferromagnetic one, as depicted in Fig. 3(a).The adiabatic theorem guarantees that such a transformation is possible for ground states of systems with a nonzero gap whenever the parameter variations are sufficiently slow.Close to a continuous quantum phase transition, however, the gap closes for a thermodynamically large system and excitations are generated unavoidably.Importantly, the celebrated quantum Kibble-Zurek mechanism (QKZ) predicts that this defect generation, and on a more gen-eral level the dynamical properties itself of crossing such a transition, displays universal behavior controlled by the underlying quantum phase transition [31][32][33].In the context of two-dimensional systems, this has been recently described at the theoretical level [34], and signatures have been observed in Rydberg experiments [35].For a finitesize system, such as the ones we deal here, the gap remains always finite.Because of this, a crossover from a QKZ regime towards an adiabatic regime emerges upon lowering the velocity of the ramp [34].In the experiment an antiferromagnetic ordering pattern has been achieved with a correlation length of the order of the system diameter, so that it is to be expected that the system resides in the crossover regime between QKZ and adiabaticity.
In what follows we will support the experimental data with numerically exact theory calculations, which will be key in a later stage in the cross-certification of the quantum simulator output.For that purpose we will use neural quantum states (NQSs), which have been recently introduced as novel class of variational wave functions for the quantum many-body problem [11].Most importantly for the purpose of this work, recent paramount advances have pointed out a route to numerically calculate quantum many-body dynamics in interacting two-dimensional quantum matter beyond what is achievable with other state-of-the-art methods [12,34].For details on the utilized numerical method we refer to Refs.[12,34] and to Appendix C.
Contrarily to the work in [10] we consider for our network analysis here two types of data sets: in the first one we use post-selected data without any defects in the array, i.e., each trap contains exactly one atom.In the second one we instead consider data sets including a mean number of defects of ∼ 3%, coming from an imperfect assembly of the atomic array [36].The purpose of this second choice is that it will allow us to make quantitative statements on the resilience of scale-free structures, and most importantly, on their significance in terms of information -and, thus, complexity -content.
Scale-free and Erdos-Renyi networks.-In Fig. 3(b), we plot the distribution P k for defect-free experimental data for square lattices of size 8 × 8 and N r = 2500 at different times.We identify two regimes: (A) At short times t = 1.52µs,P k decays exponentially with k, and its distribution resembles the one of a random ER network with k 1.This indicates that only limited correlations in the z-basis measurements are present in the system.
(B) Upon approaching the quantum phase transition (t ∼ 2.6µs) and at later times, the distribution changes drastically.In particular, we observe the emergence of a stable power-law profile with α < 2 over almost two orders of magnitude, until at large k finite sampling with N r < ∞ introduces an inevitable cutoff in the form of an exponential decay.This phenomenology is characteristic of scale-free networks.
In Fig. 3(c) we include as a comparison numerically exact theoretical results for P k by means of NQS sim-ulations.We utilize the same system parameters and number of samples as for the experimental data.The simulations capture the exact same qualitative pattern described by the experiment, already indicating that, for the depicted timescale of the experiments, the effect of dissipation on the full many-body wave functions are likely to be negligible, and validating the microscopic modelling at a quantitative level.
As depicted in Fig. 3, at large k, deviations from a power-law scaling become appreciable.Such eventual deviations appear to be a sole effect of working with a finite number of samples N r < ∞, which in turn implies that the range of the power-law behavior in P k can be extended by increasing N r .We show in Fig. 4 the distribution P k for three reference cases of times t by means of data obtained using NQS.Both qualitatively and quantitatively, P k exhibits the same features in all regimes: for ER graphs (Fig. 4(a)), increasing the number of nodes N r yields essentially the same structure of the network (keeping k 1).For the case of scale-free networks, see Figs. 4(b,c), increasing the number of samples has the effect to enlarge the regime in k of power-law behavior shifting the eventual bending, i.e., the deviation from the scale-free structure at large k, to larger and larger k.
Robustness of quantum simulator outputs.-Weobserve that at late times t > 3µs, the exponent α of the power-law tail in P k satisfies the condition α < 2 (see Fig. 4 (c)).As is known from network theory, scale-free networks with such an exponent α exhibit very robust information content with respect to perturbations.We identify such a robustness also in the experimental data.Specifically, as can be seen in Fig. 5, the experimental data sets with defects in the atomic array capture the same scaling behavior as without defects.In analogy to network theory, this analysis provides an interesting tool to characterize the robustness of quantum simulators based solely on their outputs, whenever they are described by scale-free or ER networks.An important comment is in order: such small values of the power law exponents are typically characteristic of finite networks: this is compatible with our theory, since we know that, in the infinite sampling limit N r → ∞ our network becomes infinitely large and it will be unavoidable to generate repetitions of the same snapshot.Such repetitions, however, we have excluded from the beginning and would require an adaption of our approach by means of certain weighted networks.

B. Theory of wave function networks evolution over quasi-adiabatic state preparation
The scale-free and ER WFN phenomenologies we observe in both experiment and numerical simulations are not tied to the specific problem we explore here, but, as we argue in this section, are generic features of quasiadiabatic state preparation protocols.Starting from an uncorrelated product state, it is natural to expect that at short times one may typically find random networks of wave function snapshots, i.e., networks with ER type structures.An example of such an instance is the case covered in this work, where we start from a product state with spins aligned in the z-direction.At short times the unitary dynamics will generate some weak but noticeable superposition of other configurations with a few flipped spins, which we expect to appear similar to the presence of local fluctuations such as those caused by dissipation or thermal fluctuations.These are inherently random, and should therefore yield an ER network, with Poissonlike degree distribution.For N r 2 N , such a process is expected to generate a very sparse network with k 1, due to the fact that the average distance between configurations is roughly a constant.
Upon approaching the quantum phase transition, we observe the emergence of a scale-free network structure.The basic mechanism behind this can be understood upon the inspection of the introduced metric in Eq. ( 3), which is used to impose the fundamental underlying structure on our data sets.The network structure, which we probe through P k , is generated by correlations in distances between different snapshots.Only in the case where such distances are correlated is it possible to find a power-law distribution P k of nodes having a connectivity k.As we discuss in the following, these correlations in distances between nodes in the network might be linked to the real-space correlation length in the system.Upon entering the quantum phase transition regime the system develops a large correlation length of the order of the system diameter due to the almost adiabatic dynamics generated through the experimental protocol.
From previous works on data analysis of snapshot measurements, it is expected that such large real-space correlation lengths yield Pareto distributions, i.e., powerlaw distributions, of distance measures in the data set [26,27].In this light our observations of a scalefree network structure in P k appears natural in particular because P k quantifies correlations between distances of network nodes.We note that the scale-free property of a scale-free network solely concerns P k -indeed, other network properties may carry information that reflects the presence of a finite correlation length.
Once the quantum phase transition regime is reached, the system is effectively described by a large real-space correlation length.Let us note that the dynamical behavior of the considered quantum spin model is expected to be potentially much richer as compared to the mostly studied case of dimension D = 1.Upon entering the broken-symmetry phase, the system will eventually thermalize implying an infinite correlation length.In analogy to classical systems, the temporal process of generating a long-range ordered state is typically associated by coarsening and phase-ordering kinetics [37], which also comes along with universal power-law behavior.In turn this means that upon crossing the quantum phase transition it is expected that the correlation length grows further in time also deep in the broken-symmetry phase.When  linking large real-space correlation lengths with scale-free network structures, this would imply that the scale-free network structure might survive also beyond the quantum phase transition region.This underlines the universal character of the data structure dynamics observed in the experiments.The reasoning above applies also to first order phase transitions, as long as the correlation length at the critical point is larger than the system diameter (so that, in fact, the correlation functions in the system are not able to discern differences with respect to a continuous transition).

IV. APPLICATION 1: KOLMOGOROV COMPLEXITY OF WAVE FUNCTION SNAPSHOTS
The output of a quantum simulator obtained via wavefunction snapshots is, per se, a classical object.How complex must a classical computer program be, in order to reproduce such output?This is quantified by the socalled Kolmogorov Complexity (KC) [20,21].
For generic strings, computing the KC is a NP-hard problem.The same holds true for generic graphs, where the KC is quantified by the Haussdorf dimension [38].This implies that computing the Kolmogorov Complexity of wave function snapshots is an extremely challenging task that cannot be undertaken in general.
However, as noted in the previous sections, quantum simulators often generate scale-free networks: for these, there exist known non-parametric learning algorithms that allow us to estimate the intrinsic dimension of the data points, and thus, the KC, in a manner that does not depend on scale.In particular, we utilize the 2-NN algorithm [39,40], that has already been applied in the determination of critical properties of both classical and quantum statistical mechanics partition functions [26,27].
The starting point is to consider, for each point X j in our data set, the distances to its first and second nearestneighbor, r 1 (X j ) and r 2 (X j ), respectively.Under the condition that the data set is locally uniform in the range of second nearest-neighbors, it has been shown in Ref. [39] that the cumulative distribution function F emp of µ = r 2 (x)/r 1 (x) obeys: where I d is the intrinsic dimension of the data set.The intrinsic dimension quantifies the number of degrees of freedom required to capture the information content of the data set.While this is in principle a length-scale dependent property, our estimator directly focuses on the physically relevant distance that is determined by the sampling of the many-body wave functions.In Fig. 6(a), we depict the relation between F emp and µ obtained from (a1) experiments and (a2) NQSsimulations.In both cases, and for all times considered, the distribution is compatible with Pareto (additional oscillations appear at short times, likely due to the very simple structure of the network).These results guarantee the applicability of the 2-NN approach [39].In Fig. 6(b), we show the time-dependence of the KC as measured by the intrinsic dimension across the ramp.Both experimental and simulation data clearly display two regimes: (i) up to 2 µs, the complexity increases.This effect is trivial: the initial state is very close to a product state along the z-direction, so that at short times there is just one single dominant snapshot as the measurement outcome.The unitary evolution will necessarily generate additional correlations afterwards, thus increasing complexity.(ii) From 2 µs onward, the complexity becomes a monotonously decreasing function of time.This second regime is a manifestation of the emergence of universal behavior whilst crossing a phase transition.Following quasi-adiabatic dynamics, the correlation length monotonously increases as a function of time: this implies that, in order to describe network properties, fewer variables are actually required -at equilibrium, these would just be the critical exponents and the amplitudes of correlation functions.The observations above are thus a direct manifestation of the emergent simplicity associated to universality at critical points [27], and represent, to the best of our knowledge, the first experimental demonstration of the link between complexity and quantum critical behavior.
We note that, after some time, NQS simulations predict a faster decrease of complexity with respect to the experimental data.We attribute this to the fact that the simulations can only partly keep track of the time evolution: the neural network structure utilizes a smaller number of effective variables, compatible with a decrease of KC, with respect to those that are describing the time evolution realized in experiment.

V. APPLICATION 2: CROSS-CERTIFICATION BASED ON NETWORK PROPERTIES
One of the key challenges for quantum computers and simulators is to verify their correct functioning or to certify the validity of their outcome.One basic idea in the field is cross-certification, which consists of directly comparing the output of one quantum machine with another -either quantum or classical.Recent protocols based on random unitary circuitry, aiming to compare the full ground-state wave functions, have been experimentally demonstrated to be superior to tomographic methods [16,41].However, resources still scale exponentially with system size, making the present methods inapplicable to large devices.
Here, we take a complementary angle and focus on a comparison based on wavefunction snapshots that take into account the maximum amount of extractable information with currently available resources.At the formal level, our goal is to compare two distributions in the limit where N r 2 N , i.e. which is relevant to experiments exploring many-body problems (oppositely, for N 10, it is possible to reach by brute force the regime N r 2 N [41]).Clearly, a configuration-by-configuration comparison sampled by two distributions is meaningless unless the states are very close to a product state; for generic states, the probability of sampling the same set of configurations will scale exponentially to zero with system size.
The network representation we use allows us to bypass this limitation.Specifically, we wish to compare two WFNs obtained either by two experiments or by an experiment and a simulation, see Fig 7 (a).For the concrete case considered here, let us point out that the numerical simulation by itself is a formidable challenge, which we target again by means of the NQS approach [11,12,34], see also Appendix C for details.Finding and quantifying similarities between two networks is a problem largely explored in different applications of network theory and is particularly useful for data sets that cannot be distinguished by direct inspection or low-order correlations [42].In our case, such comparisons between networks are directly related to the choice of metric used to define the WFN.For scale-free WFN, this is particularly suitable, as we are guaranteed to have chosen a metric distance capturing correlations in the system.
As a simple and efficient way to compare experimental and simulated WFNs, we check the hypothesis that the corresponding degree distributions are equal by employing a nonparametric test, known as Epps-Singleton (ES) test [43].The latter allows us to identify, with statistical significance, when two WFNs are different.In the following, since we will employ this test to establish the identity of two WFNs, we will take as statistically significant proof of our claim, cases in which the p-value of the test takes values p value > 0.1, i.e., in which both experimental and simulation data are compatible with a common probability distribution.
The results are summarized in Fig. 7.In panels (a1-a6), we present the corresponding p value of the ES tests obtained by comparing experimental data at a given time t exp with the simulation data over a given time window (i.e., 1 < t sim < 4.2 µs).This allows us to identify time windows where the quantum and classical simulators can be cross-certified with statistical significance in terms of the maximum amount of information available from their wave function networks.Interestingly, we note that the cross-certification agreement occurs at time windows that are shifted from the actual experimental times presented in the legends of Figs. 7 (a1-a6).
The fact that the cross-certification agreement of quantum systems occurs at a time t exp that is different from times considered in the simulations t sim can be attributed to miscalibrations of the parameters of the Hamiltonian (e.g., Ω(t) and δ(t)).Similar observations are found when comparing experiments with simulations of physical observables based on matrix product states [10].Although such miscalibrations do not affect the actual physics, quantifying the corresponding time shift is essential for the cross-validation of the quantum simulators.
In general, we find that the ES test can provide, for a given t exp , multiple candidate simulation times t sim for which p value > 0.1, see Fig. 7.In order to finally select across these multiple candidates, we perform a second test by computing for each of the candidates an independent second quantity.Here, we consider the staggered magnetization m stg = ix,iy (−1) ix+iy σ z ix;iy , where we included the results for all the candidate simulation times (see Fig. 7(b)).Finally, we choose the candidate t * sim , in which the simulated results for the order parameter are closest to the experimental data.As one can see from Fig. 7, up to a time of t exp = 3.0 µs, we find that this procedure is capable of cross-certifying the experimental and theoretical data within the achievable accuracy, which is limited, for instance, by the finite time grid of the theoretical data.For intermediate times 3.0 t exp 4.0 µs small deviations start to emerge, whereas for times t exp 4.0 µs the cross-certification fails, which could be caused by dissipation effects in the experiment that are not included in the theory calculation, or by a decreasing accuracy of our variational computation (similarly to what is observed in the complexity scaling).
This scheme further defines an optimal time-shift ∆t * = t * sim − t exp for the experimental data.Fig. 7 (c) shows the estimated values of ∆t * .Importantly, in the time interval that the quantum simulator can be crosscertified, we identify a small time dependence of ∆t * , which has not been addressed previously.We note that the procedure does not work well for t < 1.5µs, as expected: there, the network is not scale-free yet, so a direct comparison can only provide some rough qualitative guidance.

VI. CONCLUSIONS AND OUTLOOK
We have introduced a network theory framework to interpret the maximum amount of information extractable from quantum simulators -wave function snapshots.Remarkably, such networks can become scale-free for strongly correlated states of matter, and are of direct ex-perimental relevance, as we demonstrate with data from a large-scale Rydberg atom array experiment.We have illustrated the power of network description with two applications: demonstrating the scaling of complexity across a quantum phase transition during Kibble-Zurek scaling, and cross-certifying the wave function of a quantum and classical simulator up to system sizes that have never been attained previously.
Our work opens up a series of research directions based on a transfer of methods and concepts between network and quantum science.At the big picture level, it would be important to determine to which extent Erdos-Renyi and scale-free networks are able to characterize quantum simulators and computers.While our framework provides strong evidence that it works at and close to equilibrium, the structure of wave function networks in genuinely out-of-equilibrium situations is presently completely unknown.Understanding network properties corresponding to such dynamics might provide qualitative insights into how equilibrium is established at the wave function level, complementing current efforts focusing on observables, and providing direct links between dynamics and Kolmogorov complexity.Going beyond the case of unitary dynamics, understanding the role of dissipation might help characterize the stability of quantum dynamics to noise, which will ultimately always kick in andvery likely -imprint an Erdos-Renyi structure onto the system wave function.
In addition to conceptual insights, our framework is ideally suited to developing scalable quantum information tools.Examples range from improving crosscertification methods, and, most relevantly, applying them to data sets from large scale experiments, for which computing direct wave function overlap is hopeless.On a broader level, we believe that the parallelism between two very active, but so far disconnected, fields could be an ideal playground for developing new insights into how information is associated to many-body phenomena.
the Gauss Centre for Supercomputing e.V. (www.gausscentre.eu)for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) [44] As described in Sec.II, the structure of the wave function network (WFN) is defined by choosing a cutoff distance, R, in an embedded space defined by the Hamming distances, which allow us to define links between nodes.We now discuss in more detail the influence of R on the observation that WFNs can exhibit a scale-free structure.
Let us consider the list of all pair of distances d(X i , X j ) between nodes X i and X j .One crucial aspect to consider is that the choice of R is bounded by the minimum and maximum distances on such a list (let us call d min and d max , respectively): R < d min would generate a network with all the nodes isolated, while R > d max would generate a featureless, fully connected network.The choice of R introduced on Eq. ( 4) naturally takes into account the typical scale distance in the embedded space, which depends on the N r or the Hamiltonian parameters.
Another important aspect is that we deal with distances in a "high-dimensional" embedded space (the embedded dimension is equal to the number of spins, N ), where the so-called curse of dimensionality is expected to play a fundamental role.For instance, we could expect that the difference between the minimum and the maximum distance (i.e., d min − d max ) would become indiscernible compared to any reasonable choice of R [45] given that the volume of a high-dimensional space increases so fast that the available data becomes sparse when N r 2 N .If this were the case, we would have observed just a featureless, fully connected network.In the correlated regime, however, we observe non-trivial network structures, which can be attributed to the fact that, in reality, the intrinsic dimension of the WFNs is much lower than the dimension of the embedded space [26,27].
Let us discuss how changes on R influence the scalefree WFN. Figure 8 shows the degree distribution P k associated with the WFN generated at the quantum critical point of the quantum Ising model.By increasing R, we observe two main effects.First, the P k is shifted by larger values of k.Second, the threshold k c above which P k starts to behave as power law increases, see Fig. 8 (a).Specifically, we observe that our data for different values of R collapse in a same curve when we rescale the x-axis; see Fig. 8 (b).This result indicates that for the scale-free WFN, the main effect of increasing R is to enlarge the cutoff k c below which the network is not scale-free.In addition, we show the fraction of isolated nodes, f N0 by increasing R.For R = r 10 less then 1% of the nodes are isolated, however we still observe a power law behavior for almost one decade in k.Overall, we notice that for g ≈ g c we can always observe a scale-free WFN for a wide range of choices of R.   4)).We also show the fraction of isolated nodes, fN 0 .In panel (b) we consider a phenomenological collapse for the different P k .

B. Power law fitting of the degree distribution
One of our key results in this work is that wave function networks (WFNs) emerging in the vicinity of the quantum critical points are scale-free networks.As we discuss in this section, our conclusion is based on the empirical observation that it is very likely that a power law function describes the corresponding degree distribution, i.e., First, given a WFN, we count the number of links k of each point in the network.The corresponding list of values is used to generate the histogram P k .Second, we fit the P k by employing the approach proposed in Ref. [46].Specifically, in such an approach, one (i) defines an optimal scaling range, i.e., k > k min , and the value for the power-law exponent α by selecting an optimal fit.In particular, the one that minimizes the Kolmogorov-Smirnov (KS) distance between the empirical data and a set of trial fits.(ii) The power-law distribution (i.e., with the selected k min and α) is used to generate many synthetic data sets, which are compared with the powerlaw form by using the KS distance.The fraction of the synthetic KS distances that are larger than the empirical KS distance is then defined as the p-value of the test statistics.If the resulting p-value is greater than 0.1, the power law is a plausible hypothesis for P k ; otherwise, it is rejected [46].To implement the described approach, we use the Python power law package of Ref. [47].
As shown in Fig. 9, when considering a fitting of P k within an interval between k min = 3 and k max = 120, the results of the ES test satisfy the criterion p value > 0.1.We note that for obtaining p value > 0.1 we have to impose a maximum cutoff k max for P k .The departure from the scale-free behavior for k > k max can be attributed to the finite size of the network, N r .Based on the analysis pre- sented on Fig. 2 (b), we indeed observe the range in k that the power-law behavior is observed increases with N r .For future works, it will be interesting also to investigate the role of system size (i.e., L = √ N ) in the scale-free behavior of P k .In particular, if it is possible to establish a finite-size scaling expression addressing the role of both the network-size N r and system-size L.

C. Simulations with Neural Quantum States
Neural quantum states (NQS) have emerged recently as a new versatile class of variational wave functions [11].The goal is to find an efficient representation of a manybody wave function |ψ in the form of a parameterized function ψ θ (s) that maps a computational basis configuration s = (s 1 , . . ., s N ) to a complex number, such that Here, |s = |s 1 ⊗ . . .⊗ |s N denotes the computational basis states of a system with N degrees of freedom, and for our purposes s i ∈ {↑, ↓}.There is a number of appealing reasons to choose ψ θ (s) in the form of an artificial neural network (ANN) to render the ansatz a NQS.Most importantly, rigorous representation theorems guarantee that any possible wave function can be approximated by an ANN in the limit of large network sizes [48][49][50][51].This means that the approach is numerically exact in the sense that the accuracy of results can be certified selfconsistently by convergence checks.While the general function approximation theorems do not tell us whether the representation in the form of an ANN is efficient, it has been shown that NQS cover some volume law entangled states and correlated states of systems in two spatial dimensions, which are notoriously difficult to capture with established methods [52][53][54][55][56][57].Finally, the complexity of the algorithms involved scales gently with sys-tem size and number of parameters, and large parts are amenable to large-scale parallelization to take advantage of distributed GPU clusters [58].
While the variational ansatz with a limited number of parameters solves the problem of efficient representation, the efficient extraction of information from the wave func-Born probability distribution |ψ θ (s)| 2 ψ θ |ψ θ [59], and the same holds for all quantities of interest appearing in NQS algorithms.Notice that the only way to access the wave function in quantum simulation experiments are projective measurements, which are likewise a sampling of the Born distribution; this is a very useful parallel when attempting a direct comparison of the obtained data, because obtaining samples from the wave function could turn out to be very costly with alternative numerical approaches [10].
An optimal approximate solution of the Schrödinger equation i d dt |ψ θ = Ĥ |ψ θ within the manifold of wave functions |ψ θ is obtained via a time-dependent variational principle (TDVP) [11,12,60].This leads to an ordinary differential equation prescribing the time evolution of the variational parameters, with the quantum metric tensor notice that the imaginary part appears on both sides of the equation as we are considering real parameters [58,60].Hence, the time-evolved wave function starting from a given initial state can be obtained by integrating Eq. (11).In previous works it was found that careful regularization is crucial to achieve state-of-the-art results in this way [12,34].For the present work we developed a new way of phrasing and solving the variational problem, which we call the conditional TDVP.The details of this approach will be described in a separate manuscript [61].All results presented here were obtained in this way.
The network architecture used in our simulation is a variant of the recurrent neural network (RNN) for two-dimensional systems introduced in Ref. [62].The structure of this architecture is depicted schematically in Fig. 10.The starting point is a one-hot encoding σ i,j of the local spin configurations s i,j , i.e., σ i,j = (1, 0) if s i,j =↑ or σ i,j = (0, 1) if s i,j =↓.The neural network is + FIG. 10.Schematic depiction of the used neural network architecture.For evaluation the lattice is traversed along the path indicated by the blue arrow.A hidden state h ij is computed at each site using the one-hot encoded local basis configurations and the hidden states of previously visited neighboring sites as indicated by the pink arrows, which correspond to dense layers.From the hidden state a correlated contribution to the conditional qubit state, qij , is computed and an additional uncorrelated contribution qij added to it to obtain the logarithmic conditional amplitudes χ ij after normalization as noted in Eq. (14).
then evaluated by traversing the two-dimensional lattice in a snake-like manner.Let us denote the k-th lattice site index along the snake path as (i k , j k ) and assume that the linear dimension of the lattice is L. At each lattice site, a conditional single qubit state ψ(s i k ,j k |s 1,1 , . . ., s i k−1 ,j k−1 ) is generated in the way detailed below.From these conditional states, the coefficient of the many-body wave function is obtained as For the conditional states at every lattice site, a local hidden state h (i,j) is computed based on the spin configuration and hidden state of two neighboring sites as Here, f denotes the non-linear activation function and lm denote the weights of the dense layers; double indices indicate summation.At the boundaries, where required neighboring sites do not exist, the corresponding input is replaced by zeros.Next, the hidden state is processed by a dense layer with two-dimensional output (q ij R , qij I ), corresponding to the real and imaginary parts of a complex number qij .This number constitutes the correlated contribution to the logarithmic ↑-coefficient of the conditional local qubit state, up to normalization and a global phase.In addition, we introduced one complex-valued variational parameter qij = qij R + iq ij I for each lattice site, which corresponds to a contribution to the conditional qubit state that is uncorrelated.With q ij = qij + qij we finally produce the logarithmic conditional wave function amplitudes such that ψ(s i k ,j k |s 1,1 , . . ., The uncorrelated contribution qij extends the standard RNN architecture.We introduced it, because we found it difficult with the plain RNN to capture the initial part of the control protocol, where only the orientation of the uncorrelated qubit states is rotated and hardly any correlations are produced.In our architecture qij can fully capture the product state, such that the job of the RNN is just to account for correlations on top of it.Including qij does not affect the autoregressive property of the ansatz introduced by the decomposition into a product of conditionals (12).This means that the architecture allows for direct sampling of uncorrelated configurations at the cost of a single network evaluation per sample [62,63].
For the simulations we incorporate further experimental details, extending the elementary Rydberg atom Hamiltonian given in Eq. ( 7) of the main text.We include spatial laser intensity profiles that were extracted from the experimental setup such that the considered model Hamiltonian reads Here, we introduced the notation (kl) ≡ kL + l to map between double and single indices of the lattice sites; accordingly, the lasers shining in along one of the lattice dimensions exhibit an intensity profile perpendicular to that direction.The spatial and temporal form of the control fields during the considered protocol are shown in Fig. 11a,b.The coupling is U ij = U/∆r 6 ij with nearest neighbor interaction energy U/h = 1.947MHz,where h is Planck's constant, and ∆r (kl)(mn) = (k − m) 2 + (l − n) 2 the Euclidian distance between lattice sites.
At the beginning of the protocol all atoms are prepared in their ground state, meaning that the initial state in the spin language is a polarized state |ψ(t = 0) = |↓, . . ., ↓ .The initial part of the protocol mostly consists of a nearly adiabatic rotation of the polarization.This situation is difficult to address with NQS when using a fixed computational basis, because polarized states that align with the computational basis are hard to encode with NQS.Therefore, we implemented our simulation in a timedependent frame W (t) = exp(−iα(t) i σ y i ) with α(t) as shown in Fig. 11c) such that polarizations that align with the computational basis are avoided throughout the time evolution.and δ k (t).(c) Time-dependence of α(t), which parameterizes the time-dependent choice of the computational basis as described in the text.Initially, the quantization axis aligns with σ x , before it is rotated on the time interval between t0 = 0.8µs and t1 = 1µs with α(t) = − π 2 cos 2 π 2 (t−t0)/(t1 −t0) to align with the σ z quantization axis.

FIG. 1 .
FIG.1.Network description of many-body wave function snapshots.Panel a): construction of the network.First, samples of a wave function are collected (i) and individually mapped onto the target data space (ii).All data are then merged into a single data structure (iii), that defines a set of points in the configuration data space.This data structure is then mapped onto the corresponding wave function network (iv) by drawing links in the network according to a cutoff distance R that is determined by the data structure and the choice of metric (see text).Panel b) physical interpretation of the network structure.Within the network, the number of neighbors of given points follows a specific distribution.Points with a large number of links k (i.e., large number of points within R) are hubs, and are indicated in darker colors.As an example, taking snapshots of a classical antiferromagnet below its critical temperature will feature the antiferromagnetic state as a hub (top cartoon), while doing so well above its critical temperature will lead to a graph with no hubs and random connections (bottom cartoon).

gFIG. 2 .
FIG.2.Degree distribution, P k , for the WFN of the groundstate quantum Ising model.Panel (a) shows P k of the WFN with Nr = 10 5 nodes for g = 5.0 and g = 3.04 ≈ gc.In the paramagnetic region, the resulting network is compatible with a Poisson distribution (solid line, i.e., Erdos-Renyi network) with k = 1.As expected, in the vicinity of the critical point, the WFN becomes scale free, with α 2.4 (dashed line).For comparison we compute P k using both linear (triangles) and logarithmic (circles) histograms.Panel (b) shows P k for different values of Nr for g ≈ gc using logarithmic histograms, again with the scale free ditribution shown as a dashed line.In all cases we build the network using a cutoff R = r1 (where r1 is defined in Eq. (4)).

FIG. 3 .
FIG.3.Observation of scale-free wave-function networks in Rydberg quantum simulators.Panel (a) shows a schematic ground state phase diagram and the quasi-adiabatic state preparation scheme: the inset shows the sweep shape and the corresponding trajectory is represented by the dashed lines in the phase diagram.In the paramagnetic (PM) regime, one expects a network description compatible with a Renyi-Erdos network, while in the vicinity to the antiferromagnetic (AF) region, that contains the Kibble-Zurek regime, a scale-free network structure is expected with power-law degree distribution, P k , as illustrated by the network structures.The panel (b) presents the P k vs k of the experimentally observed wave-function networks for a square lattice with L = 8.At short times, i.e., before crossing the phase transition, the distribution decreases exponentially (similar to Erdos-Renyi degree distribution with k 1, represented by the dashed lines in the graphs).At later times (t > 3 µs), we observe a power law decay over two orders of magnitude, limited only by a bending that is due to a finite value of Nr.Graph (c) shows NQS simulations of this quasi-adiabatic protocol for the same square lattice.The scale-free behavior of P k is again observed until one is sensitive to the effects of finite sampling.We note that the value of the decay exponent is α < 2, signifying very stable wave function network properties, that will be discussed later in the presence of defects.For all cases we considered WFNs with Nr = 2500 nodes.

FIG. 4 .
FIG.4.Dependence of the degree distribution, P k , with the total number of nodes in the WFNs, Nr, at the different values of t.In the scale-free regime the maximum size of the WFN at kmax exhibit a strong dependence with Nr.The WFNs are obtained with data sets generated by NQS simulations of Rydberg experiments.

10 FIG. 6 .
FIG. 6. Complexity scaling in quantum simulators.Panels (a1-a2): cumulative distributions against µ = r2/r1 for selected times.Panels (a1) and (a2) show results for experimental and NQS-simulations data sets, respectively.The quality of a description in terms of Pareto distribution (lines ) increases as a function of time, for both simulation and experiment.Panels (b1-b2): time dependence of the intrinsic dimension, I d , along the quasi-adiabatic time evolution for both experimental and simulated data sets.For t > 2µs, the complexity of the WFN is a monotonously decreasing function of time in both experiments and simulations , capturing the emergent simplicity (decrease of degrees of freedom) that is expected from the emergence critical behavior.

FIG. 7 .
FIG. 7. Comparing the experimental and simulated WFNs, as illustrated in panel (a).In particular, we consider the Epps-Singleton two-sample test to check the hypothesis that the experimental and simulated degree distribution, P k , are equal.For each experimental time texp, we consider ES tests with the different simulated results at times tsim.Both WFNs have Nr = 2500 nodes, and we choose a cutoff distance R = r1 to generate them.Panel (a1-a6) shows the corresponding p value as a function of tsim: results with p value > 0.1 (marked by the dashed lines) are interpreted as statistically significant.To cross-check our analysis, we also consider in panel (b) the order parameter, mstg, as a function of texp.Each simulated result corresponds to the different times tsim for which p value > 0.1.Such an analysis allows us to identify t * sim , the times where the best agreement exists between simulation and experimental data; the results corresponding to t * sim are marked as the black star points in the panels (a1-a6).Finally, panel (c) shows the corresponding time-shift, ∆t * = t * sim − texp, between experiments and NQS simulations (see text).

FIG. 8 .
FIG. 8. Degree distribution, P k , for the WFNs of the groundstate quantum Ising model at the critical point, g = gc.Panel (a) shows P k of the WFNs with different values of R = rc (see Eq. (4)).We also show the fraction of isolated nodes, fN 0 .In panel (b) we consider a phenomenological collapse for the different P k .

FIG. 9 .
FIG. 9. Power law fitting of the degree distribution presented on Fig. 2. We use Komolgorov-Smirnoff statistics to perform the fitting and define the p-value (see text).

FIG. 11 .
FIG. 11. (a,b) Control protocols of the external fields Ω k (t)and δ k (t).(c) Time-dependence of α(t), which parameterizes the time-dependent choice of the computational basis as described in the text.Initially, the quantization axis aligns with σ x , before it is rotated on the time interval between t0 = 0.8µs and t1 = 1µs with α(t) = − π 2 cos 2 π 2 (t−t0)/(t1 −t0) to align with the σ z quantization axis.
. This work is supported by the European Union's Horizon 2020 research and innovation program under grant agreement No. 817482 (PASQuanS), the Agence Nationale de la Recherche (ANR, project RYBOTIN), and the European Research Council (Advanced grant No. 101018511-ATARAXIA).