Exact quantum dynamics of Fermi--Hubbard systems using the Gaussian phase-space representation with diffusion gauges

TL;DR

The paper develops and benchmarks a diffusion-gauged Gaussian Phase-Space Representation (GPSR) method for real-time fermionic Fermi–Hubbard dynamics. A randomized SVD-based diffusion gauge constructs noise matrices that extend the practical simulation time, enabling larger 1D, 2D, and 3D Hubbard systems to be simulated beyond the reach of exact diagonalization and mean-field approaches. Compared with analytic gauges, the numerical diffusion gauge significantly delays spiking trajectories and yields accurate two-point correlations (g(2)) relative to ED, while incurring computational costs dominated by diffusion-matrix decompositions that scale as approximately O(n_s^4). The approach demonstrates substantial scalability, including a 36-site 3D system, and provides a path toward simulating more complex, entangled states (e.g., Bell states) and guiding future phase-space-method developments. Overall, the work expands the practical reach of GPSR for real-time quantum dynamics in strongly correlated fermionic lattices and outlines routes for further efficiency and applicability improvements.

Abstract

We use the Gaussian Phase-Space Representation to solve the real-time dynamic of interacting fermions in 1D, 2D, and 3D systems. The method is exact up to a spiking point, which represents a limit on the practical simulation time. The spiking can be delayed, and the practical simulation time extended, by adjusting the gauges of the representation, resulting in different equivalent stochastic differential equations. Here, we work on the so-called diffusion gauge and propose an algorithm to find efficiently new implementations of the noise terms. Compared with our initial results [F. Rousse \textit{et al.} 2024, J. Phys. A: Math. Theor. \textbf{57}, 015303], the new method achieves a significantly longer practical simulation time and can be applied to significantly larger systems.

Figures (8)

• Figure 1: Site occupations of 8-site systems in different dimensions. The initial conditions are global spin waves, Eq. (\ref{['eq:IC8']}), and we are presenting the spin-up occupation, $n_{ii \uparrow}$. The number of curves are different because symmetries cause some occupancies to be equal in the 2D and 3D systems. The black dashed curves represent the ED solutions, the red dash-dotted ones the HF solutions, the solid blue ones the GPSR solutions computed with the analytical noise matrix, and the solid green ones the GPSR solutions computed with the noise matrix from the randomized SVD decomposition. For all systems, switching from an analytical noise matrix (blue curves) to the noise matrix produced by the randomized SVD (green curves) at least doubles the practical simulation time of the GPSR solution.
• Figure 2: Pair-correlations between spin-up particles site occupations. For the 1D 8-site system, correlation between sites 1 and 2, then between sites 1 and 3, the next-to-nearest neighbor.
• Figure 3: Pair-correlations between spin-up particle sites occupation. For the 1D 8-sites system, correlation between sites 1 to 4 and their respective right neighbors (2 to 5). The curves starting at 0 are the correlation between sites 4 and 5 which start without spin-up particles, hence the initial difficulty of the fTWA to compute an accurate correlation. We have here used $10^4$ trajectories, but fTWA can still be improved visibly (closer to the ED solution) if we add up more trajectories.
• Figure 4: On-site correlations on the 4 first sites. In the GPSR plots, the spikes announcing the failure of trajectories are clearer and appear sooner than in the occupation curves (Fig \ref{['fig:FH8_ni']}). fTWA is still struggling to model correlation with near-to-zero occupation. We here have used $10^4$ trajectories, fTWA can still be visibly improved (closer to the ED solution) by adding up more trajectories. We did not plot the HF solution which is a constant 1 here.
• Figure 5: Site occupations of a 3D 36-sites system ($4\times3\times3$). The initial condition is a global spin wave, and we are following the spin-up occupation of two sites. Dynamic simulation of systems of this size is unreachable with the ED method. On this system, the GPSR method with a numerical noise matrix (green line) has a practical simulation time three times larger than the one computed with the analytical noise matrix (blue line). We also observe that fTWA (magenta dashed line) produces an accurate estimation of the occupation until at least the moment GPSR fails.