We study the two-dimensional t-J model at finite temperature directly in the thermodynamic limit using purification represented by an infinite projected entangled-pair state (iPEPS). We reach temperatures down to T/t = 0.1 and hole concentrations up to 1 - n 0.25, and provide benchmark thermodynamic-limit results for the specific heat, uniform susceptibility, and charge compressibility. We identify a susceptibility maximum T* that tracks the buildup of short-range antiferromagnetism and a shallow compressibility enhancement upon cooling in the same doping window. To expose the underlying microscopic mechanism, we introduce dopant-conditioned multipoint correlators that quantify how holes reorganize nearby exchange: Single holes weaken adjacent antiferromagnetic bonds, while nearest-neighbor hole pairs produce a cooperative response that reinforces antiferromagnetism on the parallel plaquette edge. Over the same parameter window, d-wave pairing correlations remain short-ranged. These results provide experiment-compatible thermodynamic-limit benchmarks and establish dopant-conditioned correlators as incisive probes of finite-temperature spin-texture reorganization in doped Mott insulators.

Quantum phase transitions provide a controlled route for generating many-body excitations, but the dynamics after the critical point can be as important as the initial defect creation. Recent progress in quantum annealing has made it possible to access coherent nonequilibrium dynamics in programmable Ising systems with thousands of superconducting qubits. Here we use this capability to study a longitudinally biased quantum Ising chain, where Kibble–Zurek defect creation is followed by nonintegrable post-critical dynamics. The longitudinal bias confines kink–antikink excitations into mesonic bound states, so that the final spin configurations encode both the production of defects near the critical point and the subsequent evolution of the confined excitations. Using energy-scale rescaling and zero-noise extrapolation, we find that the defect density follows the expected biased Kibble–Zurek/Landau–Zener crossover and agrees with matrix-product-state simulations with uniform bias. In contrast, magnetization, spatial profiles, and minority-domain statistics reveal that mesonic evolution is interrupted by localization of the post-critical domain pattern. Matrix-product-state simulations with disorder reproduce this separation between robust defect creation and localized post-critical dynamics. Our results show that large-scale quantum annealers can probe the fate of critical excitations beyond defect counting.

The emergent practical applicability of the Quantum Approximate Optimization Algorithm (QAOA) for approximate combinatorial optimization is a subject of considerable interest. One of the primary limitations of QAOA is the task of finding a set of good parameters, which is usually done using a variational optimization loop. Parameter transfer, or parameter concentration, is a phenomenon where QAOA angles trained on problem instances that are self-similar tend to perform well for other problem instances from that similar class. This suggests a potentially highly efficient and scalable nonvariational learning method for QAOA angle finding. In this work, we systematically study QAOA parameter transferability from small problem sizes (16 and 27 decision variables) onto large problem instances (up to 156 qubits) for heavy-hex graph Ising models with geometrically local higher-order terms using the Julia based QAOA simulation tool JuliQAOA to perform classical angle finding for up to 49 QAOA layers (p). Parameter transfer of the fixed angles is validated using a combination of full state vector, projected entangled pair states, matrix product state, and LOWESA numerical simulations. We find that the QAOA parameter transfer from single instances applied to other (unseen) problem instances does not in general provide monotonically improving performance as a function of p-there are many cases where the performance temporarily decreases as a function of p-but despite this the transferred angles have a general trend of improved expectation value as the QAOA depth increases, in many cases converging close to the true ground-state energy of the 100+ qubit instances. We also sample the hardware-compatible Ising models using the ensemble of transfer-learned QAOA parameters on several superconducting qubit IBM quantum processors with 127, 133, and 156 qubits. We find continuous solution quality improvement of the hardware-compatible QAOA circuits run on the IBM noisy intermediate-scale quantum processors up to p = 5 on ibm_fez, up to p = 9 on ibm_torino, and up top = 10 on ibm_pittsburgh.

This work introduces SpinGlassPEPS.jl, a software package implemented in Julia, designed to find low-energy configurations of generalized Potts models, including Ising and QUBO problems, utilizing heuristic tensor network contraction algorithms on quasi-2D geometries. In particular, the package employs the Projected Entangled-Pairs States to approximate the Boltzmann distribution corresponding to the model's cost function. This enables an efficient branch-and-bound search (within the probability space) that exploits the locality of the underlying problem's topology. As a result, our software enables the discovery of low-energy configurations for problems on quasi-2D graphs, particularly those relevant to modern quantum annealing devices. The modular architecture of SpinGlassPEPS.jl supports various contraction schemes and hardware acceleration.

The quantum Fisher information (QFI) is a geometric measure of state deformation calculated along the trajectory parametrizing an ensemble of quantum states. It serves as a key concept in quantum metrology, where it is linked to the fundamental limit on the precision of the parameter that we estimate. However, the QFI is notoriously difficult to calculate due to its nonlinear mathematical form. For mixed states, standard numerical procedures based on eigendecomposition quickly become impractical with increasing system size. To overcome this limitation, we introduce a numerical approach based on Lyapunov integrals that combines the concept of symmetric logarithmic derivative and tensor networks. Importantly, this approach requires only the elementary matrix product states algorithm for time evolution, opening a perspective for broad usage and application to many-body systems. We discuss the advantages and limitations of our methodology through an illustrative example in quantum metrology, where the thermal state of the transverse-field Ising model is used to measure magnetic field amplitude.

Optimization problems pose challenges across various fields. In recent years, quantum annealers have emerged as a promising platform for tackling such challenges. Here, to provide an additional perspective, we develop a heuristic tensor-network-based algorithm to reveal the low-energy spectra of Ising spin-glass systems with interaction graphs relevant to present-day quantum annealers. Our deterministic approach combines a branch-and-bound search strategy with an approximate calculation of marginals via tensor-network contractions. Its application to quasi-two-dimensional lattices with large unit cells of up to 24 spins, realized in current quantum annealing processors, requires a dedicated approach that uses sparse structures in the tensor-network representation and GPU hardware acceleration. We benchmark our approach on random problems defined on Pegasus and Zephyr graphs with up to a few thousand spins, comparing it against the D-Wave Advantage quantum annealer and the simulated bifurcation algorithm, with the latter representing an emerging class of classical Ising solvers. In addition to examining the quality of the best solutions, we compare the diversity of low-energy states sampled by all the solvers. For the largest considered independent and identically distributed problems with over 5000 spins, the state-of-theart tensor-network approaches lead to solutions that are 0.1% to 1% worse than the best solutions obtained by Ising machines while being 2 orders of magnitude slower. We attribute these results to approximate contraction failures. For embedded tile planting instances, our approach reaches approximately 0.1% from the planted ground state, a factor of 3 better than the Ising solvers. While all three methods can output diverse low-energy solutions-e.g., differing by at least a quarter of spins with energy error below 1%-our deterministic branch-and-bound approach finds sets of a few such states at most. In contrast, both Ising machines prove capable of sampling sets of thousands of such solutions.

Recent demonstrations on specialized benchmarks have reignited excitement for quantum computers, yet their advantage for real-world problems remains an open question. Here, we show that probabilistic computers, co-designed with hardware to implement Monte Carlo algorithms, provide a scalable classical pathway for solving hard optimization problems. We focus on two algorithms applied to three-dimensional spin glasses: discrete-time simulated quantum annealing and adaptive parallel tempering. We benchmark these methods against a leading quantum annealer. For simulated quantum annealing, increasing replicas improves residual energy scaling, consistent with extreme value theory. Adaptive parallel tempering, supported by non-local isoenergetic cluster moves, scales more favorably and outperforms simulated quantum annealing. Field Programmable Gate Arrays or specialized chips can implement these algorithms in modern hardware, leveraging massive parallelism to accelerate them while improving energy efficiency. Our results establish a rigorous classical baseline for assessing practical quantum advantage and present probabilistic computers as a scalable platform for real-world optimization challenges.

Quantum computers hold the promise of solving certain problems that lie beyond the reach of conventional computers. However, establishing this capability, especially for impactful and meaningful problems, remains a central challenge. Here, we show that superconducting quantum annealing processors can rapidly generate samples in close agreement with solutions of the Schrödinger equation. We demonstrate area-law scaling of entanglement in the model quench dynamics of two-, three-, and infinite-dimensional spin glasses, supporting the observed stretched-exponential scaling of effort for matrix-product-state approaches. We show that several leading approximate methods based on tensor networks and neural networks cannot achieve the same accuracy as the quantum annealer within a reasonable time frame. Thus, quantum annealers can answer questions of practical importance that may remain out of reach for classical computation.

We identify persistent oscillations in a nonintegrable quantum Ising chain. In the integrable chain with nearest-neighbor interactions, the nature, origin, and decay of post-transition oscillations are tied to the Kibble-Zurek mechanism. Remarkably, when coupling to the next-nearest neighbor is added, the resulting nonintegrable ”zigzag” chain (still in the quantum Ising universality class) supports persistent oscillation: Topological defects (kinks) appear as a result of the quantum phase transition. However, in a ”zigzag” Ising chain defects can form Cooper pairs. The oscillation of the Cooper-pair condensate has a frequency that depends on the binding energy gap between the paired and the unpaired defects, so it can be excited by resonant driving. While one might have expected that the integrability-breaking ”zigzag” coupling causes relaxation, the oscillations we identify are persistent: Their longevity in our simulations is likely limited only by numerical accuracy. This oscillation of the Cooper-pair condensate of kinks is a manifestation of quantum coherence and should be experimentally accessible.

Doped Mott insulators host intertwined spin-charge phenomena that evolve with temperature and can culminate in stripe order or superconductivity at low temperatures. The two-dimensional t-J model captures this interplay yet finite-temperature, infinite-size calculations remain difficult. Using purification represented by a tensor network - an infinite projected entangled-pair state (iPEPS) ansatz - we simulate the t-J model at finite temperature directly in the thermodynamic limit, reaching temperatures down to one tenth of the hopping rate and hole concentrations up to one quarter of the lattice sites. Beyond specific heat, uniform susceptibility, and compressibility, we introduce dopant-conditioned multi-point correlators that map how holes reshape local exchange. Nearest-neighbor hole pairs produce a strong cooperative response that reinforces antiferromagnetism on the adjacent parallel bonds, and single holes weaken nearby antiferromagnetic bonds; d-wave pairing correlations remain short-ranged over the same window. These results provide experiment-compatible thermodynamic-limit benchmarks and establish dopant-conditioned correlators as incisive probes of short-range spin-texture reorganization at finite temperature.

We analyze the correlation between the energy, momentum, and spatial entanglement produced by two luminal jets in the massive Schwinger model. Using tensor network methods, we show that for m/g>1/π, in the vicinity of the strong- to weak-coupling transition, a nearly perfect and chargeless effective fluid behavior appears around the midrapidity region with a universal energy-pressure relationship. The evolution of energy and pressure is strongly correlated with the rise of the spatial entanglement entropy, indicating a key role of quantum dynamics. Some of these observations may be used to analyze high multiplicity jet fragmentation events, energy-energy and energy-charge correlators at current collider energies.

In a recent preprint [1] (arXiv:2503.05693), Tindall et al. presented impressive classical simulations of quantum dynamics using tensor networks. Their methods represent a significant improvement in the classical state of the art, and in some cases show lower errors than recent simulations of quantum dynamics using a quantum annealer [2] (King et al., Science, eado6285, 2025). However, of the simulations in Ref. [2], Ref. [1] did not attempt the most complex lattice geometry, nor reproduce the largest simulations in 3D lattices, nor simulate the longest simulation times, nor simulate the low-precision ensembles in which correlations grow the fastest, nor produce the full-state and fourth-order observables produced by Ref. [2]. Thus this work should not be misinterpreted as having overturned the claim of Ref. [2]: the demonstration of quantum simulations beyond the reach of classical methods. Rather, these classical advances narrow the parameter space in which beyond-classical computation has been demonstrated. In the near future these classical methods can be combined with quantum simulations to help sharpen the boundary between classical and quantum simulability.

We investigate the relation between static correlation functions in the ground state of local quantum many-body Hamiltonians and the dispersion relations of the corresponding low energy excitations using the formalism of tensor network states. In particular, we show that the Matrix Product State Transfer Matrix (MPS-TM) - a central object in the computation of static correlation functions - provides important information about the location and magnitude of the minima of the low energy dispersion relation(s) and present supporting numerical data for one-dimensional lattice and continuum models as well as two-dimensional lattice models on a cylinder. We elaborate on the peculiar structure of the MPS-TM's eigenspectrum and give several arguments for the close relation between the structure of the low energy spectrum of the system and the form of static correlation functions. Finally, we discuss how the MPS-TM connects to the exact Quantum Transfer Matrix (QTM) of the model at zero temperature. We present a renormalization group argument for obtaining finite bond dimension approximations of MPS, which allows to reinterpret variational MPS techniques (such as the Density Matrix Renormalization Group) as an application of Wilson's Numerical Renormalization Group along the virtual (imaginary time) dimension of the system.

We examine the stationary-state equations for lattices with generalized Markovian dephasing and relaxation. When the Hamiltonian is quadratic, the single-particle correlation matrix has a closed system of equations even in the presence of these two processes. The resulting equations have a vectorized form related to, but distinct from, Lyapunov's equation. We present an efficient solution that helps to achieve the scaling limit, e.g., of the current decay with lattice length. As an example, we study the super-diffusive-to-diffusive transition in a lattice with long-range hopping and dephasing. The approach enables calculations with up to 104 sites, representing an increase of 10 to 40 times over prior studies. This enables a more precise extraction of the diffusion exponent, enhances agreement with theoretical results, and supports the presence of a phase transition. There is a wide range of problems that have Markovian relaxation, noise, and driving. They include quantum networks for machine-learning-based classification and extended reservoir approaches for transport. The results here will be useful for these classes of problems.

We present an open-source tensor network Python library for quantum many-body simulations. At its core is an abelian-symmetric tensor, implemented as a sparse block structure managed by logical layer on top of dense multi-dimensional array backend. This serves as the basis for higher-level tensor networks algorithms, operating on matrix product states and projected entangled pair states, implemented here. Using appropriate backend, such as PyTorch, gives direct access to automatic differentiation (AD) for cost-function gradient calculations and execution on GPUs or other supported accelerators. We show the library performance in simulations with infinite projected entangled-pair states, such as finding the ground states with AD, or simulating thermal states of the Hubbard model via imaginary time evolution. We quantify sources of performance gains in those challenging examples allowed by utilizing symmetries.

Resonant transport occurs when there is a matching of frequencies across some spatial medium, increasing the efficiency of shuttling particles from one reservoir to another. We demonstrate that in a periodically driven, many--body titled lattice, there are sets of spatially fractured resonances. These ``emanate'' from two essential resonances due to scattering off internal surfaces created when the driving frequency and many--body interaction strength vary, a scattering reminiscent of lens flare. The confluence of these fractured resonances dramatically enhances transport. At one confluence, the interaction strength is finite and the essential resonance arises due to the interplay of interaction with the counter--rotating terms of the periodic drive. We discuss the origin and structure of the fractured resonances, as well as the scaling of the conductance with system parameters. These results furnish a new example of the richness of open, driven, many--body systems.

Adiabatic preparation of a critical ground state is hampered by the closing of its energy gap as the system size increases. However, this gap is directly relevant only for a uniform ramp, where a control parameter in the Hamiltonian is tuned uniformly in space towards the quantum critical point. Here, we consider inhomogeneous ramps in two dimensions: initially, the parameter is made critical at the center of a lattice, from where the critical region expands at a fixed velocity. In the 1D and 2D quantum Ising models, which have a well-defined speed of sound at the critical point, the ramp becomes adiabatic with a subsonic velocity. This subsonic ramp can prepare the critical state faster than a uniform one. Moreover, in both a model of $p$-wave paired 2D fermions and the Kitaev model, the critical dispersion is anisotropic -- linear with a nonzero velocity in one direction and quadratic in the other -- but the gap is still inversely proportional to the linear size of the critical region, with a coefficient proportional to the nonzero velocity. This suffices to make the inhomogeneous ramp adiabatic below a finite crossover velocity and superior to the homogeneous one.

The Minimally Entangled Typical Thermal States (METTS) are an ensemble of pure states, equivalent to the Gibbs thermal state, that can be efficiently represented by tensor networks. In this article, we use the Projected Entangled Pair States (PEPS) ansatz as to represent METTS on a two-dimensional (2D) lattice. While Matrix Product States (MPS) are less efficient for 2D systems due to their complexity growing exponentially with the lattice size, PEPS provide a more tractable approach. To substantiate the prowess of PEPS in modeling METTS (dubbed as PEPS-METTS), we benchmark it against the purification method for the 2D quantum Ising model at its critical temperature. Our analysis reveals that PEPS-METTS achieves accurate long-range correlations with significantly lower bond dimensions. We further corroborate this finding in the 2D Fermi Hubbard model at half-filling. At a technical level, we introduce an efficient \textit{zipper} method to obtain PEPS boundary matrix product states needed to compute expectation values. The imaginary time evolution is performed with the neighbourhood tensor update.

The simulation of open many-body quantum systems is challenging, requiring methods to both handle exponentially large Hilbert spaces and represent the influence of (infinite) particle and energy reservoirs. These two requirements are at odds with each other: larger collections of modes can increase the fidelity of the reservoir representation but come at a substantial computational cost when included in numerical many-body techniques. An increasingly utilized and natural approach to control the growth of the reservoir is to cast a finite set of reservoir modes themselves as an open quantum system. There are, though, many routes to do so. Here, we introduce an accumulative reservoir construction-an ARC-that employs a series of partial refreshes of the extended reservoirs. Through this series, the representation accumulates the character of an infinite reservoir. This provides a unified framework for both continuous (Lindblad) relaxation and a recently introduced periodically refresh approach (i.e., discrete resets of the reservoir modes to equilibrium). In the context of quantum transport, we show that the phase space for physical behavior separates into discrete and continuous relaxation regimes with the boundary between them set by natural, physical timescales. Both of these regimes "turnover" into regions of over-and underdamped coherence in a way reminiscent of Kramers' crossover. We examine how the range of behavior impacts errors and the computational cost, including within tensor networks. These results provide the first comparison of distinct extended reservoir approaches, showing that they have different scaling of error versus cost (with a bridging ARC regime decaying fastest). Exploiting the enhanced scaling, though, will be challenging, as it comes with a substantial increase in (operator space) entanglement entropy.

Sampling a diverse set of high-quality solutions for hard optimization problems is of great practical relevance in many scientific disciplines and applications, such as artificial intelligence and operations research. One of the main open problems is the lack of ergodicity, or mode collapse, for typical stochastic solvers based on Monte Carlo techniques leading to poor generalization or lack of robustness to uncertainties. Currently, there is no universal metric to quantify such performance deficiencies across various solvers. Here, we introduce a new diversity measure for quantifying the number of independent approximate solutions for NP-hard optimization problems. Among others, it allows benchmarking solver performance by a required time-to-diversity (TTD), a generalization of often used time-to-solution (TTS). We illustrate this metric by comparing the sampling power of various quantum annealing strategies. In particular, we show that the inhomogeneous quantum annealing schedules can redistribute and suppress the emergence of topological defects by controlling space-time separated critical fronts, leading to an advantage over standard quantum annealing schedules with respect to both TTS and TTD for finding rare solutions. Using path-integral Monte Carlo simulations for up to 1600 qubits, we demonstrate that nonequilibrium driving of quantum fluctuations, guided by efficient approximate tensor network contractions, can significantly reduce the fraction of hard instances for random frustrated 2D spin-glasses with local fields. Specifically, we observe that by creating a class of algorithmic quantum phase transitions, the diversity of solutions can be enhanced by up to 40% with the fraction of hard-to-sample instances reducing by more than 25%.