A Monte Carlo Time-Dependent Variational Principle

Abstract: We generalize the Time-Dependent Variational Principle (TDVP) to dissipative systems using Monte Carlo methods, allowing the application of existing variational classes for pure states, such as Matrix Product States (MPS), to the simulation of Lindblad master equation dynamics. The key step is to use sampling to approximately solve the Fokker-Planck equation derived from the Lindblad generators. An important computational advantage of this method, compared to other variational approaches to mixed state dynamics, is that it is "embarrassingly parallel".

• The paper proposes a Monte Carlo-based stochastic extension of the Time-Dependent Variational Principle to simulate Lindblad master equation dynamics in open quantum systems.
• It demonstrates that ensemble averaging over Matrix Product States recovers physical observables with 1/√N variance scaling, confirmed against analytic benchmarks.
• The method offers excellent parallel scalability and extends variational ansatzes to non-integrable and dissipative models, paving the way for large-scale simulations.

A Monte Carlo Time-Dependent Variational Principle for Open Quantum Systems

Introduction

This paper presents a framework that generalizes the Time-Dependent Variational Principle (TDVP) to dissipative quantum many-body systems, leveraging Monte Carlo sampling to simulate Lindblad master equation dynamics in variational ansatzes for pure states, in particular Matrix Product States (MPS) (1411.5546). The approach constructs an effective variational dynamics for mixed states through the stochastic evolution of variational parameters, thereby extending the applicability of highly successful wavefunction-based methods to general open-system scenarios. A salient computational aspect of the proposed algorithm is its embarrassingly parallel structure, which is pertinent for scalability.

Methodological Foundation

The cornerstone of the approach is the representation of a mixed state $\rho_t$ as an ensemble average over pure states parameterized by a variational ansatz,

$\rho_t = \int_{\mathcal{M}} p_t(\bar{a}, a) \ketbra{\Psi(a)}{\Psi(a)} da d\bar{a},$

where $p_t(\bar{a}, a)$ evolves according to the Lindblad master equation mapped into a Fokker-Planck (FP) equation on the variational parameter manifold. The dynamics of $p_t$ are unraveled into a system of Ito SDEs for the parameters $a$:

$$da^j(t) = (b_Q^j + \langle \bar{L}_\alpha\rangle b_\alpha^j) dt + b_\alpha^j dw_\alpha,$$

where $b_Q^j$ and $b_\alpha^j$ are projections of the Lindblad operators onto the tangent space of the variational manifold.

The algorithm thus amounts to integrating these SDEs for an ensemble of trajectories, each corresponding to a pure-state evolution, then averaging physical observables over the ensemble.

Performance Benchmarks

Benchmarking was performed using MPS as the variational class, with both small and large system sizes and various dissipative regimes. Results for two-qubit systems confirmed convergence of ensemble-averaged observables to exact master equation results with the expected $1/\sqrt{N}$ variance scaling, where $N$ is the number of Monte Carlo samples.

For larger systems, e.g., a 16-site XXZ Heisenberg chain driven at the edges, the method accurately reproduced analytic steady-state magnetization profiles derived in the strong dissipation regime, for sufficiently large bond dimension and sample sizes. Deviations in the bulk for intermediate bond dimensions were attributed to limited entanglement capacity and slow propagation of boundary-induced dissipation.

Figure 1: Quantitative agreement for the XXZ chain with edge driving ($n_\text{sites} = 16$) and low sample numbers, compared to the analytical solution.

Additionally, the method was applied to non-integrable models, including the XXZ chain with bihomogenous dissipation (up-flip and down-flip Lindblad operators on the left and right chain halves, respectively).

Figure 2: Magnetization profiles in the XXZ chain with bihomogenous dissipation, demonstrating the capacity to resolve interference effects in the domain wall region at weak interactions.

One noteworthy outcome is the capacity of relatively moderate bond dimensions (far less than the full Hilbert space dimension) to capture interference and nontrivial correlation effects in dissipative scenarios.

Dissipative Correlation Dynamics

Time-resolved two-point correlation functions, e.g., $\langle \sigma^x_8 \sigma^x_n \rangle$, were computed to probe the dynamical emergence of antiferromagnetic domains and correlation transport in the presence of engineered dissipation.

Figure 3: Left, time-resolved two-point correlation functions in the $K_{XZ}$ model with uniform dissipation; right, the effect of bihomogenous dissipation on the spatial structure of correlations.

The results confirm that the sampling-based TDVP captures both symmetry-breaking effects due to boundaries and the smoothing of correlations across domain walls when competing dissipative processes are present.

Algorithmic Characteristics and Parallel Scaling

The main computational advantage is the parallel independence of sample trajectories. This allows perfect strong scaling in $N$ and enables the simulation of larger systems and longer temporal evolution windows than statevector-based unravelings or variational mixed-state methods constrained by information distance ambiguities.

The convergence rate in observables is governed by $1/\sqrt{N}$, and increasing sample size does not incur cross-sample overhead. This advantage is particularly relevant when exploring system regimes with demanding accuracy requirements or critical dynamics, as in the edge-driven XXZ chain.

Practical and Theoretical Implications

Practically, the method enables physically-motivated variational ansatzes developed for closed systems (such as MPS, PEPS, and other tensor networks) to be extended to arbitrary Lindblad-type open-system dynamics, bypassing the necessity for ad hoc mixed-state parameterizations or complex purifications.

Theoretically, translating the Lindblad equation into a FP equation for variational parameters solidifies the conceptual bridge between open-system stochastic quantum trajectories and variational simulation. The scheme encompasses quantum state diffusion as a special case when the variational class spans the full Hilbert space, and reduces to standard TDVP for unitary evolution in the absence of Lindbladian terms.

Limitations and Future Prospects

While the method demonstrates good accuracy under moderate entanglement and local dissipation, its efficiency is ultimately constrained by the expressiveness of the chosen variational class (e.g., the bond dimension for MPS) and the presence of criticality-induced slow dynamics. However, its amenability to high-performance and grid computing architectures opens up the possibility of addressing far larger and more complex open-system scenarios, such as two-dimensional driven-dissipative lattices, engineered quantum information devices, and quantum thermalization processes.

The method's generic structure allows for the incorporation of more expressive ansatzes (e.g., tree tensor networks, variational neural states), exploration of unbiased (non-product state) initializations, and systematic investigation of dissipative phase transitions.

Conclusion

The paper establishes a practical synthesis of Monte Carlo sampling and variational quantum dynamics in the form of a stochastic time-dependent variational principle for open systems. By mapping Lindblad dynamics onto stochastic flows in parameter space, it significantly extends the applicability of variational tensor network methods to dissipative many-body phenomena, with parallelization properties advantageous for contemporary computational resources. This framework is poised to enable otherwise intractable large-scale simulation and is likely to fuel further advances in the study of quantum non-equilibrium physics, dissipation-engineered order, and variational algorithms for quantum statistical mechanics.