Adaptive variational low-rank dynamics for open quantum systems

We introduce a novel, model-independent method for the efficient simulation of low-entropy systems, whose dynamics can be accurately described with a limited number of states. Our method leverages the time-dependent variational principle to efficiently integrate the Lindblad master equation, dynamically identifying and modifying the low-rank basis over which we decompose the system's evolution. By dynamically adapting the dimension of this basis, and thus the rank of the density matrix, our method maintains optimal representation of the system state, offering a substantial computational advantage over existing adaptive low-rank schemes in terms of both computational time and memory requirements. We demonstrate the efficacy of our method through extensive benchmarks on a variety of model systems, with a particular emphasis on multi-qubit bosonic codes, a promising candidate for fault-tolerant quantum hardware. Our results highlight the method's versatility and efficiency, making it applicable to a wide range of systems characterized by arbitrary degrees of entanglement and moderate entropy throughout their dynamics. We provide an implementation of the method as a Julia package, making it readily available to use.

Understanding the impact of both incoherent and coherent relaxation processes on quantum systems is therefore crucial for the progress of quantum technologies.Consequently, the development of fast and scalable methods for their classical simulation is equally important.The exact numerical simulation of many-body quantum systems, however, suffers from the curse of dimensionality, namely the exponential growth of computational resources, scaling as d N , required to simulate the dynamics of a system with N identical modes, each with d possible states.As a result, simulating even moderately large quantum systems becomes rapidly an intractable problem [29].
When a system is not isolated, but couples to the degrees of freedom of an external environment, the two become entangled.The system's state is then represented by an Hermitian positive semi-definite operator, the density matrix, resulting from tracing out the environmental degrees of freedom.Under the assumption of a memoryless environment, the time evolution of the density matrix is governed by the Lindblad master equation [30][31][32] which has found application in a very broad range of physical settings [33].In this context, the curse of dimensionality becomes even more daunting, with computational requirements scaling as d 2N .Different approaches for mitigating the complexity of the task exist [34], each one compromising on the degree of accuracy with which different features of the statistical ensemble can be captured.Among these, Monte Carlo quantum trajectory methods [35,36], phase-space methods [37,38], semiclassical approximations [39,40], tensor network techniques [34,[41][42][43][44][45][46][47], and variational approaches [48][49][50][51][52][53][54][55] possibly coupled with Monte Carlo sampling and neural network techniques [56][57][58][59][60][61][62], stand out.
In recent years, another class of methods know as ensemble truncation methods has emerged as a powerful tool for simulating isolated [63][64][65][66][67] and dissipative quantum systems [68][69][70][71].This approach is based on the realization that many quantum systems, particularly those with low entropy, can be effectively represented by a density matrix of significantly lower rank than what the whole Hilbert space would require.This reduction is achieved by focusing on a subset of states that capture the essential structure of the statistical ensemble characterizing the mixed quantum state, thereby reducing computational complexity while maintaining accuracy.The low-rank (LR) hypothesis applies to many systems of physical interest.Among these are systems weakly coupled to the environment, or systems initialized in a pure, i.e. zero-entropy state, in the early stage of their dynamics.Modern quantum-computing platforms, designed to minimize dissipation and noise, often fall in either of the above instances.

arXiv:2312.13676v1 [quant-ph] 21 Dec 2023
A considerable challenge remains, however, determining how the LR states should evolve to optimally represent the density matrix along the time evolution, and in particular how the dimension of the LR subspace should vary to accommodate changes of the entropy over time.Building upon this premise, this paper introduces a novel method encapsulating the most relevant features of the ensemble truncation schemes within the robust framework of the time-dependent variational principle (TDVP) developed for open quantum systems in Refs.[48,[72][73][74].Indeed, our method, which we deem LR-TDVP method, effectively integrates the benefits of both approaches.On one side it leverages the dynamical truncation methods' ability to efficiently represent quantum states with minimal information loss and the dynamically adjustment of the basis' dimension to adapt to changes in the system's entropy.On the other it takes advantage of the variational principle to dynamically modify the LR basis states thus ensuring the optimal fidelity of the evolution.This combination is particularly suited for quantum systems characterized by an arbitrary degree of entanglement but moderate entropy as is often the case in modern quantum hardware.By addressing the challenges faced in the NISQ era, our method offers a significant advancement in the simulation of complex quantum systems, paving the way for new discoveries and applications in quantum computing and simulation.
Our method is designed to be universally applicable, independent of specific system characteristics such as spatial symmetries, particle statistics, geometry of the space, or topology of the interactions.To facilitate its use and integration into various research workflows, we have incorporated it as a method of the QuPhys library [75], a Julia-based [76] quantum physics toolkit.The library, along with our method, is readily accessible and can be found at the repository listed in [77].

II. THE LOW-RANK TDVP
Let us consider an open quantum system whose dynamics is described by the Lindblad master equation where ρ ≡ ρ(t) (for brevity) is the system density matrix at time t, L is the Liouvillian superoperator, and we take ℏ = 1.While the coherent part of the dynamics is described by the Hamiltonian Ĥ acting on a Hilbert space H of dimension N H , the incoherent evolution is described by the dissipator which formalizes the action of the jump operator Γ on the system.
The exact numerical simulation of this dynamics is generally computationally challenging due to the exponential scaling of the required resources, growing as the square of the Hilbert space dimension.The focus of our work is however on low-entropy systems, whose dynamics is known to be well captured, at all times, by a limited number of states {|φ k (t)⟩ ; k = 1, . . ., M } with M ≪ N H .These states span the low-rank subspace H M ⊆ H which is assumed to encapsulate the essential structure of the statistical ensemble throughout the dynamics.In recent years, various heuristic algorithms have been developed that effectively truncate the density matrix ρ to a rank-M matrix, significantly smaller than what the whole Hilbert space would require.Such methods are collectively known as ensemble truncation methods [68][69][70][71].Fundamental to these methods is their ability to dynamically adjust the truncation rank M = M (t) at each time step, adapting to an entropy landscape that varies over time.However, a rigorous and systematic approach for selecting and adapting the low-rank basis states remains a challenge.
Concurrently, variational methods have emerged as a versatile tool for circumventing the exponential space problem by considering trial states from a physically motivated, small subset of the exponentially large Hilbert space [48].In these methods, the density matrix ρ(t) = ρ(θ(t)) is expressed in terms of a set of variational parameters θ(t) which evolve in time to guarantee, within the expressiveness of the ansatz, the optimal approximation of the evolution generated by the action of L on ρ(t).
In this work, we introduce the low-rank TDVP algorithm where the benefits of both classes of methods are united.By employing the McLachlan TDVP we rigorously estimate the optimal subspace H M (t) over which to decompose the evolution at each time step (Sec.II A).Additionally, we develop schemes to dynamically adjust the dimension of this subspace in response to changes in the system's entropy (Sec.II B).Deferring a detailed analysis of the ensemble truncation methods to App.C, we lastly compare our variational method to the heuristic algorithms established in [68][69][70] (Sec.II C).The main features of this method are summarized in Fig. 1.

A. Variational equations of motion
The most general variational parametrization of an arbitrary statistical ensemble can be expressed as:  Schematic representation of the dynamic adjustment of the variational manifold's dimension.Both the manifold and the rank of the truncated density matrix are dynamically modified to maintain a consistent level of accuracy, adapting to changes in the system's entropy.This adjustment reflects the principle that highly entropic systems necessitate a larger number of states for accurate representation, whereas low-entropy ensembles can be effectively described with fewer states.
Throughout the paper we will refer to |φ k (t)⟩ as the variational or low-rank states and to M (t) as the rank of ρ.Importantly, the variational method ensures a dynamical adjustment of the variational states, guaranteeing the optimal set of states is selected at all times to best approximate the evolution under L. To enable this dynamic adjustment, the variational states themselves must be parametrized.In our approach, we adopt a full linear parametrization, decomposing |φ k ⟩ on the computational basis {|e α ⟩} α=1,...,N H as: For clarity, and where it does not lead to confusion, we will omit the time argument in subsequent discussions.While this work focuses on full linear parametrization, alternative parametrizations relying on tensor networks or neural network architectures are indeed possible and will be subject of future investigation.
The dynamical problem outlined in Eq. ( 1) can now be reformulated into a set of equations of motion (EOM) for the variational parameters, comprising the population matrix B = [B ij ] and the coefficient matrix z = [z αk ].To derive these EOM, we apply the McLachlan variational principle, with the additional constraint of trace preservation: where λ is introduced as a Lagrange multiplier.While Eq. ( 5) is general in scope, substituting Eq. (3) for ρ yields the EOM specific for the chosen ansatz (the details of the derivation are provided in App.A): where L = L(ρ)z, L = z † L, S = z † z, while S −1 and B −1 denote the inverse matrices of S and B respectively.Since both S and B may be singular, their inverse is often ill-defined.To overcome this problem, throughout the paper we adopt a smooth regularization criterion based on the singular value decomposition of the two matrices, as proposed in Ref. [78].Details are deferred to App.B. The dynamics of the system can now be obtained by numerical integration of Eq. ( 6) which we perform using high-order adaptive-timestep solvers.Each time step requires the calculation of B −1 and of L. The first we perform using singular value decomposition, while the second only requires matrix multiplications, the most extensive being between (extremely sparse) N H × N H and (dense) N H × M matrices.Notably, since for a linear parametrization Ṡ = 0 (see App. A), the computation of S −1 incurs no additional computational cost after its initial calculation at the beginning of the process.
As previously anticipated, within each integration step, if M ≪ N H (LR hypothesis), both the speed and storage requirements for the calculation are dramatically improved.For additional information on the computational costs and procedures, we refer readers to App.B.

B. Dynamical rank update
A critical aspect of any LR approach is its capability to dynamically adjust the dimension of the LR basis throughout the system's evolution [M = M (t)].This adaptability is essential for accommodating changes in the system's entropy over time.In the ensemble truncation methods [68][69][70], this update is achieved by truncating the full-rank density matrix ρ(t+dt) at each time step with the goal of keeping a control quantity, namely the truncation error ϵ M = 1 − M j=1 p j , below a predefined threshold ϵ max .In this section we detail schemes for efficiently implementing a dynamic update of the dimension of the variational basis in our approach.

Basis inflation
Let's first consider the basis inflation problem, where the system's entropy increases over time.An example of a system presenting this behaviour is developed in Sec.III A. To start, we posit the existence of a control quantity χ(t) that accurately reflects the solution's accuracy and that can be efficiently computed at each time step.Various options for χ will be discussed below.This quantity is monitored through a callback mechanism, continuously evaluated against a predefined upper bound of ϵ max .Initially zero at t = 0, χ grows with the system's entropy, and hence, over time.Let t * denote the time when χ(t * ) = ϵ max .At this time, the rank is incrementally increased from M to M + 1 through the addition of a new state |φ M +1 ⟩ to the LR basis.The population and coherences associated with this new state are initially vanishing [B j,M +1 = (B M +1,j ) * = 0 ], to ensure continuity in the solution.The variational principle ensures that the solution can be made independent of the specific choice of |φ M +1 ⟩.This independence, within the specified tolerance, is guaranteed if the state is added at a time prior to t * , ensuring that |φ M +1 ⟩ (t * ) satisfies the TDVP at that time.To achieve this, we implement periodic checkpoints at intervals of ∆t.Each time a threshold crossing is registered, the integration is restarted with an enlarged basis from the nearest checkpoint before the crossing.This approach also ensures that the bound on 2. Time evolution of the approximate variational error, χ(t), as defined in Eq. (10).The upper bound for χ is set at ϵmax = 10 −4 , which strictly applies for ∆t > 0. The inset displays the overlap between the true solution and the truncated one as obtained from Eq. ( 7).The physical model underlying this example is the XYZ Heisenberg model (detailed in Sec.III A) with N = 9 spins and Jy = 1.The values of all remaining physical parameters [see Eqs. ( 11) and ( 12)] are the same as those used in Fig. 4.
χ remains strict, as the response in χ may not be immediate following the basis inflation.However, we observe that in most situations of interest, an instantaneous basis inflation (∆t = 0) does not significantly compromise the solution's accuracy, suggesting that the variational adjustment of |φ M +1 ⟩ is almost instantaneous.This is illustrated in Fig. 2, where we display χ(t) for both ∆t = 0 and ∆t = 0.2.As shown in the inset, the overlap between the exact and approximated solutions, evaluated using the expression remains virtually unchanged when choosing ∆t = 0.

Basis deflation
For the basis deflation problem, we adopt a similar approach to that used for basis inflation.Consider a case where the system's entropy, after peaking, gradually diminishes as the system approaches the stationary state.An example of such dynamics is discussed in Sec.III C. In this context, to enhance efficiency, it is desirable to reduce the dimension of ĤM(t) when χ falls below a lower threshold ϵ min .Upon crossing this threshold, the rank is incrementally decreased from M to M − 1 by removing a state from the LR basis.The selection of which state to remove from the basis is a critical aspect for the success of the deflation protocol.Intuitively, we should eliminate the state contributing the least to the dynamics, poten-tially the one with the lowest occupancy.However, a low occupancy does not necessarily imply negligible correlations or coherences.Therefore, we first transition to the diagonal basis by diagonalizing ρ(t) as with p 1 ≥ p 2 ≥ . . .≥ p M .Subsequently, we remove the M th state from the basis but maintain the equivalent diagonal representation of the LR basis, resulting in z(t As discussed in Refs.[68][69][70], this matrix shares the same non-vanishing eigenvalues p j as the small M × M matrix σ = Ĉ † Ĉ.The eigenvectors of these two matrices are related through Ĉ.It is important to note that, unlike the basis inflation case, failing to reduce the rank in the deflation scenario does not lead to a less accurate solution.Quite the contrary, the solution will be more accurate, albeit at a higher computational cost.

Control quantities
We now provide a set of possible choices for χ.Given the variational parametrization of the LR states, the natural control quantity for our method is the distance between the true evolution and the variational evolution of the trial state, calculated as [48] The only linear algebra operations involved in the computation of this quantity are matrix multiplications, the largest of which between extremely sparse N H × N H matrices and dense N H ×M matrices.The number of matrix multiplications to be computed scales as D 2 , the square of the number of jump operators.A more efficient choice, which comes at no additional computational cost, is We adopted this choice for all simulations in the paper.Although Tr{L(ρ)} = 0 for any physical density matrix ρ and Liouvillian L, Tr{P L(ρ)} can be non-zero as a result of projecting over an incomplete basis [74,79].Indeed, Tr{P L(ρ)} is only vanishing when ρ is entirely contained in the LR manifold.Assuming the latter to be constructed around the initial state of the simulation, Tr{P L(ρ)} vanishes at t = 0, and in an ideal evolution would remain so throughout the dynamics.Its departure from zero quantifies leakage outside of the variational manifold, and as such, it is a good indicator of the solution's accuracy in real time.This choice for χ is closely related to the truncation error ϵ M of Ref. [68], as can be seen from the perturbative expansion carried out in App. C.
The last criterion we propose takes χ = p M /p 1 where the probabilities p i follow from the diagonalization of ρ(t) according to Eq. ( 8).
All three choices of χ are well suited control quantities as evidenced by Fig. 3 where we display the overlap between the real and variational solution as a function of the threshold ϵ max .As expected, 1 − O [calculated from Eq. ( 7)] and M increase the lower the value of ϵ max .

C. Comparison with ensemble truncation methods
In this section, we outline the significant computational advantages that the LR-TDVP method offers over the ensemble truncation schemes discussed in Refs.[68][69][70]: 1. Ensemble truncation methods execute linear algebra operations on N H ×M (D +1) matrices, leading to larger computational and storage demands than those required by our variational LR density matrix approach, which only requires storing matrices of size at most N H × M .
2. Ensemble truncation methods rely on a Kraus representation of the Liouvillian map.The operators decomposing this map are proportional to √ dt, thereby restricting these methods to the explicit Euler first-order integration scheme.In Ref. [68], this difficulty is partially lifted by replacing the Kraus operator K0 with the full time-evolution operator K0 = exp(−iH eff dt) and using a highorder integration scheme for its application.The contribution of the other Kraus operators to the dynamics is however still limited to an explicit first-order in time Euler scheme.Conversely, our method allows integrating the EOM with a highorder adaptive-timestep solver.This approach not only significantly enhances the efficiency of our method but also offers greater accuracy and flexibility compared to the aforementioned schemes.
3. Checks for basis inflation and deflation can be carried out at no additional costs to the simulation of our dynamics so that the additional diagonalization of an M × M matrix is only necessary when the deflation criterion is triggered.This is not the case in ensemble truncation methods which, for this purpose, require the diagonalization of an M (D + 1) × M (D + 1) matrix at every time step.
4. Unlike ensemble truncation methods, where the selection of states to retain in the dynamics is heuristic in nature, our method's selection process is grounded in the time-dependent variational principle ensuring the adoption of the optimal LR basis at each time step.

III. NUMERICAL SIMULATIONS
In this section, we outline the significant computational advantages that our variational method offers over the ensemble truncation schemes discussed in Refs.[68][69][70], and demonstrate its applicability in a broad range of physical models.
Although a key strength of our method lies in its ability to capture the entire Liouvillian dynamics, we begin our analysis focusing on systems for which the relevant physics occurs in the steady state.Here, we leverage our method only to reach stationarity, disregarding the remainder of the dynamics.In doing so, we showcase the method's effectiveness in analyzing steady-state phenomena in low-entropy systems.In Sec.III A, we investigate the dissipative phase transition in the anisotropic XYZ Heisenberg model.In Sec.III B, we apply our method to bosonic models, specifically to the study of the steadystate of the bosonic simulator of the triangular antiferromagnetic Ising model.This exploration underscores the versatility and general-purposefulness of our approach.
We proceed to showcase the algorithm's capability to accurately capture the full Liouvillian dynamics, particularly in scenarios with non-monotonous entropy profiles.In Sec.III C we provide insight into the mechanisms of basis inflation and deflation by studying the dynamics of the transverse field Ising model in the presence of weak magnetic fields for which the system is known to be lowrank.In Sec.III D, we address the timely problem of simulating the full dynamics of bias-preserving gates in dissipative cat qubit architectures displaying once again the advantage of a LR approximation.

A. Dissipative Heisenberg XYZ model
We consider a two-dimensional lattice consisting of N spin-1/2 particles governed by the Heisenberg XYZ Hamiltonian: where σα i (α = x, y, z) are the Pauli operators acting on the ith spin of the lattice.The incoherent relaxation of each spin is described by the dissipator of σ− j = (σ x j − iσ y j )/2, leading to the master equation: This equation exhibits a Z 2 symmetry, as evidenced by its invariance under the transformation σx,y j → −σ x,y j ∀ j.In the thermodynamic limit, the steady state of the system bearing this symmetry is the fully aligned spindown state, |↓, ↓, . . ., ↓⟩.This state is characterized by having zero magnetization in the xy plane, a defining feature of the so called paramagnetic phase.
In accordance with the conventions established in the relevant literature, all simulations of this model are performed fixing J x = 0.9, J z = γ = 1, and h z = 0.In this section we vary J y between 0.9 and 1.1 to investigate the transition from the paramagnetic to the ferromagnetic phase.Observing the transition in a fully quantum model that retains all correlations in the system is challenging.Indeed, the steady-state magnetization along the y direction, used as an order parameter in cluster mean-field studies [39], is always zero in the full model.Similarly, the homogeneous steady-state structure factor which, in a Gutzwiller approximation is zero in the paramagnetic phase and acquires a positive value only in the ferromagnetic one [95], is not a good order parameter when correlations between different sites are included.Alternatives have been proposed.Specifically, the angle-averaged magnetic susceptibility [94] and the trace-distance susceptibility [97].The latter relies on the distance between steady state matrices evaluated at infinitesimally close values of J y and is thus too sensitive to any form of noise to be of use in approximated solutions.The former, however, would be well suited for our approximation.Its computation is however quite cumbersome as it requires assessing the linear response of the system to small perturbations.Here, we thus prefer to present as an indicator of the transition the variation in J y of that is, the square root of the steady state expectation value of the variance of My .This is a good indicator of the long-range correlations developing in the array across the transition, and its derivative ∂∆M y /∂J y in J y is a clear signature of critical behaviour.We display this quantity in Fig. 4 for increasingly larger system sizes.
A finite-size analysis in the linear lattice length L = √ N of the position of the maximum variance derivative yields an estimate of J c = 1.06 ± 0.03, consistent with the results obtained in Ref. [94].The maxima of the curves themselves follow a power-law growth in L with critical exponent κ = 0.92 ± 0.08.
To demonstrate the increasingly mixed character of the steady-state density matrix ρss , we display in Fig. 5 the steady-state von Neumann entropy and its derivative ∂S/∂J y as a function of J y .As expected, we find that in proximity of the critical point, the entropy sharply rises with a slope that increases with the system size.Concurrently, ∂S/∂J y shows a peaked  4).The inset displays the maximum value of the derivative as a function of the linear lattice length L = √ N .The solid line is a powerlaw fit of the finite-size scaling.We use the same physical parameters as in Fig. 4.
structure which becomes increasingly more pronounced as we move towards the thermodynamic limit.By fitting the maximum entropy derivative with a power law in L [i.e max(∂S/∂J y ) ∝ L λ ], we get an estimate for the critical exponent of λ = 1.6 ± 0.1, again consistent with the state-of-the-art results found in Ref. [94] via the corner space method.
The system's steady state is obtained by numerically integrating Eq. ( 6) for a time γt = 15 long enough for the system to reach stationarity.In all simulations, the rank is dynamically increased starting from an initial choice of M (0) = N + 1.As the initial state of the dynamics we take the pure state |ψ(0)⟩ = |↓, ↓, . . ., ↓⟩ setting to zero the population of the remaining N states making up the initial variational basis.These states are chosen as those with minimal Hamming distance from |ψ(0)⟩.
As S increases, additional states are required to accurately describe ρss , thereby progressively reducing the computational advantage introduced by a LR ansatz.Given that S increases monotonically in J y , the advantage introduced by our method is largest in proximity of the paramagnetic phase, where the entropy remains small.This is corroborated by Fig. 6, where we display M (γt = 15) as a function of J y .Indeed, for highly entropic configurations [c.f.Fig. 5(a)] M → 2 N , i.e the LR simulation becomes equivalent to the full solution of the master equation.
As a figure of merit for the method's computational 0.90 0.95 1.00 1.05 1.10 Rank M of the steady state of the system as a function of the coupling strength Jy.For small values of Jy, i.e small values of S, the rank of the system remains fixed at the initial preset value of N + 1.The rank increases with the entropy of the system according to the discussion in Sec.II B. We use the same physical parameters as in Fig. 4.
performance we take the memory footprint of steadystate simulations similar to those performed above.The free parameters in this analysis are: • The number of spins N , which uniquely determines the Hilbert space dimension.
• The coupling strength J y , which uniquely determines the system's configuration.Recall that all other physical parameters have been fixed.
• The dimension M of the LR space.This can be fixed at the start, as we will do in what follows, or varied dynamically, in which case the free parameter becomes the threshold ϵ max .
For a synthetic two-dimensional representation of the results -memory footprint against Hilbert space dimension -as depicted in Fig. 7, both J y and M need to be fixed.We fix J y = 0.98 and select M , for each value of N , so as to ensure that the relative error in the estimation of M z remains within a 0.1% margin from the actual solution.The distance between the LR estimate and the true value of M z for different lattice sizes is shown in the inset to Fig. 7.While for N < 16 the true value of M z is obtained via numerical integration of the full master equation, for N = 16 we resort to the Monte Carlo trajectory method, averaging over 10 4 trajectories.The results, as illustrated in Fig. 7, highlight the significant computational advantages of our approach.The LR-TDVP method not only accurately captures all the critical features of the XYZ model but also demonstrates improved efficiency in terms of memory usage and computation time compared to both the full solution and the dynamical truncation schemes referenced in [68][69][70].
In concluding this section, it is important to stress that our method is inherently dynamic in nature.As such, it is not expected to always outperform methods specifically tailored for solving steady-state problems.For in- Memory allocations as a function of number of simulated spins N for the full simulation and the LR-TDVP one.We set the value of Jy = 0.98.For each value of N , we choose the rank M ensuring the relative error on the magnetization Mz to be smaller than 0.1%, namely err The inset shows the behaviour of err(Mz) as a function of M .The color scheme in the inset is the same as in Fig. 4, as are all physical parameters aside from Jy.
stance, when compared to the corner space renormalization method, which serves as our reference, the LR-TDVP does not surpass it in terms of the system sizes that can be explored.The corner space method, based on spatial block decimation, is indeed very effective for systems with spatial translational invariance.
One of the key advantages of the LR-TDVP method over the corner space approach, however, is its generality and applicability to a broader range of low-entropy systems, not limited to those characterized by translational invariance.Additionally, this method offers more precise control over error estimation and convergence in the rank.In contrast, the corner space method's truncation at any decimation level introduces errors that, while locally controllable, can accumulate and propagate through different decimation steps, affecting the overall accuracy.This often necessitates running multiple simulations with varying truncation combinations to ensure convergencea requirement not present in our algorithm.
Furthermore, the merging process in the corner space method is not unique and can lead to artificially slow timescales, potentially affecting performance.The fixed but non-unique choice of basis in the corner space method depends on the initial lattice for decimation, meaning that missing a crucial state in the system can lead to erroneous or artificially slow simulations.Conversely, our algorithm allows for a more precise final result, even if part of the dynamics might be initially missed.Error control is much more straightforward, and can be dynamically adjusted by expanding the LR dimension.This flexibility offers a distinct advantage in terms of precision and control over the simulation process.
In summary, while the LR-TDVP method may not always surpass specialized steady-state solvers in terms of system size exploration, it compensates with its generality, error control, and adaptability.This makes it a valuable tool for studying a range of low-entropy quantum systems, including those with non-uniform interactions where methods like the corner space renormalization may not be applicable.

B. Frustrated antiferromagnetism in quadratically driven QED cavities
We consider an array of N = 2, 3 coupled dissipative quadratically-driven QED nonlinear cavities: a bosonic model known to be a valid quantum simulator of the the triangular antiferromagnetic Ising model [93,[98][99][100].The Hamiltonian of the model reads where U is the Kerr nonlinearity, G the two-photon driving field amplitude, ∆ the resonators' detuning frequency, and J the hopping coefficient.We include single-and two-body losses with rates γ and η respectively, so that the system evolves according to the Lindblad master equation In the limit of a strong driving field, the photons in each cavity form a coherent state with phase α or −α.While a positive hopping coefficient J > 0 favours the formation of a statistical mixture of two equiprobable separable coherent states [101] |Ψ ± ⟩ = j |±α⟩ j , a negative one J < 0 favours the emergence of states |Φ ± ⟩ = |±α, ∓α, ±α, ∓α, . ..⟩ with an antisymmetric alignment of phases.The two limits map directly to the ferromagnetic and antiferromagnetic configuration of the effective spin model upon the mapping |α⟩ → |↑⟩ and |−α⟩ → |↓⟩.
Following Ref. [99], we investigate the effective antiferromagnetic spin model, by studying the steady-state properties of the photonic system for varying values of G. Specifically, we focus on the steady-state values of the first-order coherence correlation function and of the von Neumann entropy S defined in Eq. ( 16).
Our results, displayed in Fig. 8, are consistent with the findings in Ref. [99].Indeed, as expected for an antiferromagnetic coupling, the correlation g 1,2 for N = 2 is negative and, for increasing values of G, converges to the asymptotic value g increases from S(G = 0) = 0 to S(G ≫ γ) = ln(2), consistently with the double degeneracy of the ground state manifold of the equivalent spin model.
For an odd number of sites, the presence of geometrical frustration in the system hinders the emergence of |Φ ± ⟩ states in the open system delaying to larger values of G the onset of the antiferromagnetic signatures.The effects of frustrations are visible for N = 3 in the non-monotonous behavior of S as a function of G, and in its asymptotic value of S = ln (6).Indeed, the latter is consistent with the existence of six possible spin configurations, emerging as a result of frustration, and minimizing the energy of the equivalent spin model.
We obtained all results by numerically integrating Eq. ( 6) over a time interval sufficiently long to allow the system to reach the steady state.The Hilbert space of the N cavities is obtained as the N -fold tensor product of truncated single-boson Fock spaces.For numerical convenience, we adjusted the truncation dimension according to the prescription N cut (G) = max{10 , 5G/W }, where W = U 2 + η 2 .For N = 2, a LR subspace of dimension M = (M pm ) N with M pm = 2 (number of states per mode) accurately captures the relevant physics, as evidenced by the indistinguishable results obtained for M pm = 3.For N = 3, a larger LR space is necessary.We present results for M pm = 2, 3, 4, 5. Convergence is observed for M pm ≥ 3.
Direct numerical integration of Eq. ( 18) for determining the steady state of the system quickly becomes impractical, as highlighted by the extensive efforts in Ref. [99].In that approach, the steady state is computed by numerically solving ∂ t ρ = 0, a linear system of up to 10 8 equations.In contrast, our method demonstrates significant potential by efficiently computing the steady state by integration of Eq. ( 18) but retaining only a small number of optimal states.This approach not only reduces computational complexity but also maintains accuracy, showcasing the advantages of our method in handling large-scale quantum systems.In this section, we demonstrate our algorithm's ability to accurately capture the full dynamics of a system while efficiently adapting to an entropy profile that does not monotonically increase over time.Specifically, we examine a scenario where the entropy initially rises, reaches a peak at t max with a maximum value S max , and then decreases towards a steady-state value S ss < S max .
This behavior is exemplified in the transverse field Ising model [J x = J y = 0, J z = γ = 1, and h x > 0 in Eqs. ( 11) and (12)] which is characterized by a dynamics that, under specific initial conditions, remains LR at all times.Indeed, for sufficiently small values of h x in (11), the steady state of the model closely resembles that at h x = 0, which reads |ψ 0 ⟩ = |↓⟩ ⊗N .Consequently, an evolution under small h x starting from the latter state exhibits the desired entropy profile while remaining amenable to a LR representation.In contrast, an evolution starting from the state |↑⟩ ⊗N would saturate the rank, rendering the LR ansatz an inefficient choice.
In Fig. 9(a), we display the time evolution of the structure factor ∆M y as defined in Eq. (15).We remark the perfect agreement between the LR and full simulations.
Figure 9(b) depicts the desired evolution of the system's entropy and highlights the dynamic adaptation of the rank M (t) to its changes throughout the simulation.Specifically, M (t) initially increases, in parallel with S, to maintain the accuracy threshold set by ϵ max .However, for t > t max , as S decreases, M is automatically reduced, thereby enhancing computational efficiency.
In Fig. 9(c), we display the variational error χ, as defined in Eq. (10), which serves as a measure of the error performed at time t.As expected, the automatic adaptation of M ensures that ϵ min < χ(t) < ϵ max .In contrast, when the rank is fixed to either M max = max t [M (t)] or M min = min t [M (t)], there exist regions where χ(t) falls below ϵ min or exceeds ϵ max , respectively.
Generating and stabilizing this manifold is far from trivial, as it requires engineering parity-preserving processes involving exclusively the pairwise exchange of photons between the system and its environment [105,108,110,111].In dissipative cat qubits, for instance, this manifold coincides with the four-degenerate steady-state manifold of the dissipator , modelling engineered two-photon drive and dissipation processes [16,19,20,[103][104][105]122].Indeed, because of the underlying strong Z 2 symmetry of the Liouvillian, the rightkernel of L 0 is [16,107,123,124] As(H) = lim so that, for any choice of ρ0 ρss = lim The coefficients c µ,ν encode the quantum information and are defined as where the operators Ĵµ,ν , spanning the left kernel of L 0 , are defined as [105,124] Here, where I q (•) is the modified Bessel function of the first kind, and n!! = n•(n−2)!! is the double factorial, applied element-wise to all diagonal operators.
The two-photon drive and dissipation processes modelled by L 0 have been realized, for example, on superconducting circuit platforms.In these platforms, the dominant source of errors comes from single photon loss processes, modeled by the dissipator κ 1 D[â] with κ 1 /κ 2 ≪ 1.These processes hinder the code's ability to encode quantum information as well as the performance of logical gates [4,118,122].

Z, ZZ, and ZZZ gates
Let us consider the application of Z rotations of an angle π over N = 1, 2, 3 cat qubits.For a single qubit, these rotations are equivalent to the application of a logical Pauli σz operator, which introduces a phase of (−1) on the logical |1⟩ state [125].In the cat basis, this amounts to exchanging |C + α ⟩ and |C − α ⟩.Similarly, in multi-qubit systems, these rotations correspond to a unitary transformation changing the tensor product Logical Z operations can be approximately implemented by evolving the system for a time T under the Liouvillian L = L 0 +L 1 given by [105,106,118,119,122]: Here, L 1 encompasses both the unwanted single body losses that compromise the accuracy of gate operations, and the ideal Hamiltonian evolution responsible for the rotation.Specifically, the Hamiltonian ĤZ is defined as [117,126]: with If we initialize each qubit in the |+⟩ state, after a time T , ideally, all qubits should be in the |−⟩ state, for which ⟨−| Ĵ++ |−⟩ = 0. Therefore, we can use Ĵ++ t=T (30) as a measure of the phase-flip error probability during the Z rotation.Fig. 10(a) displays P Z as a function of α 2 .Simulations were conducted in a truncated Fock space with dimension N cut (α) = max{20 , ⌈4.5α 2 ⌉}, ensuring the convergence of each data point.Notably, the L simulations, executed with M pm = 3 states per mode, overlay  (30)] as a function of the photon number α 2 .We display as a dashed gray line the analytical approximation PZ ≈ κ1T α 2 + ϵ 2 Z T /α 2 (for N = 1) provided in Ref. [117].(b) Bit-flip error probability PX [Eq.(31)] as a function of the photon number α 2 .To improve visualization, as the curves would largely overlap, we multiply each curve by a different constant, effectively shifting them vertically on the plot.The constants are R = 10 2 , 1, 10 −2 for N = 1, 2, 3 respectively.Full lines are used to display the exponential fit of the LR data.Because of computational constraints, the exact simulations for N = 3 were only possible up to α 2 = 5.The inset shows the von Neumann entropy S at time T as a function of α 2 for N = 1.In both panels, the results obtained from the L simulations (squares) are in agreement with the full solution (yellow crosses).Parameters: κ2 = 1, κ1 = 1/1000, ϵZ = ϵZZ = ϵZZZ = 1/20.perfectly with the results obtained from the full evolution.Simulations for N = 3 were restricted to α 2 < 5 because of prohibitively large memory requirements.
Conversely, since any rotation leaves the states lying on the rotation axis unchanged, if we initialize each qubit in |0⟩ ≈ |α⟩, we expect to find them in that exact same state after evolving under L for a time T .As a quantifier of the bit-flip error rate affecting the system we can therefore take where x = â+â † , and sgn(x) ≈ Ĵ+− .In Fig. 10(b) we display P X as a function of α 2 .Noticeably, our model is able to accurately capture P X even up to the very small error probabilities associated with its exponential suppression in α 2 [19,105,107,108].The suppression coefficients ζ extracted from an exponential fit of the exact simulations are ζ = 2.13 ± 0.01, 2.13 ± 0.01, and 2.16 ± 0.03, respectively for N = 1, 2, 3.The same fit on the LR data yields ζ = 2.13 ± 0.01, 2.14 ± 0.01, and 2.14 ± 0.02.These estimates are in agreement with each other and with the theoretical expectation [117].As we show in the inset of Fig. 10(b), the entropy S of the target state at time T decreases with the number of photons in the system.This entails that progressively less states are necessary to capture the dynamics.To avoid over-representation in our LR ansatz we thus reduce M pm (number of states per mode included in the simulation) with the photon number.Specifically, we obtain good agreement with the full solution by taking M pm = 1 for α 2 > 6.

IV. CONCLUSIONS
In this paper, we introduce a novel, model-independent method designed for the efficient simulation of the dynamics of low-entropy systems.Recognizing that such evolution can be accurately captured by a limited number of states, our method builds upon and advances the previously established ensemble truncation schemes [68][69][70][71], integrating their key features within the framework of the time-dependent variational principle recently developed for open quantum systems [48,[72][73][74].Our approach enhances previous ensemble truncation methods by offering a rigorous and systematic protocol for defining and dynamically modifying the low-rank basis.Furthermore, it extends prior variational descriptions of dissipative systems by introducing a computationally efficient protocol for adapting the size of the variational manifold, thereby dynamically adjusting the rank of the density matrix throughout the simulation.This dynamic adaptation ensures that the low-rank subspace of the Hilbert space optimally represents the system's state at all times.
To ensure easy integration into diverse research workflows, we have implemented our method in Julia [76] and incorporated it into the QuPhys library [75,77], a comprehensive toolkit for quantum system simulation.
Looking ahead, potential extensions of our work include the incorporation of efficient variational representations of the low-rank space, such as tensor-network [47,127] or neural quantum states [128].These advancements could enable the simulation of larger-scale systems.Additionally, a digital version of our method could approximate the simulation of noisy quantum computing hardware, potentially enhancing error-mitigation protocols that rely on estimating hardware-specific noisy quantum channels [129,130].
where the last equality follows once again from Eq. (A8).Upon multiplying by B −1 on the right, and substituting Eq. (A9) for Ḃ, we find that where P = zS −1 z † is the projector on the LR manifold.
Since by definition P ż = 0, the latter reduces to A few comments on the full parametrization are now in order.First, it is evident from Eq. (A13) that τ = z † ż = 0, which consequently reduces Eq. (A9) to the form presented in Eq. ( 6).Second, although Ref. [74] showed that using a Lagrange multiplier method to enforce both energy and trace preservation does not, in general, result in the unitary evolution expected for an isolated system; any linear parametrization as the one in Eq. ( 4) is exempt from this problem.

Appendix B: Computational details
The integration of the evolution problem presented in Eq. ( 6) relies exclusively on linear algebra operations involving matrices of dimensions N H ×M (z and L) and/or M × M (S, L, and B), where M ≪ N .This disparity in dimensions allows for efficient execution of these operations, given the relatively smaller size of at least one dimension of each matrix.The bulk of the computational effort is dedicated to calculating the matrix L, defined as: (B1) By following the order of operations as indicated by the parentheses in the equation, the computation of L involves matrix-matrix multiplications between M × M matrices, N H × M and M × M matrices, as well as between extremely sparse N H × N H and dense N H × M matrices.This approach ensures that no dense N H × N H matrix is ever stored in memory or involved in any multiplication.The process is optimized for efficiency, leveraging the reduced dimensionality to minimize the computational load, and free from the diagonalization of large matrices required in Refs.[68][69][70].
The integration of Eq. ( 6) does however necessitate the calculation of the inverse matrices of S and B. Since these matrices may be singular, their inverse is often illdefined.In order to progress, regularization schemes are typically employed.A common approach involves adding a small diagonal contribution to the matrices to mitigate the impact of small eigenvalues.Yet, we have found that the effect of such an addition is not well controlled during time evolution, particularly when using an adaptive integrator.Similar observations were reported in Ref. [78].Taking inspiration from this work, we adopt a regularization scheme based on the singular value decomposition of S = U ΣV † , where Σ = diag(σ 2 1 , . . ., σ 2 M ).Specifically, we define the pseudo-inverse of S as S −1 = V Σ + U † where We do the same for B. The value of λ in chosen adaptively at each iteration according to (B3) Throughout the paper we set atol = 10 −6 .Note that the Moore-Penrose pseudo-inverse corresponds to choosing the discontinuous function f (σ 2 ) = θ(σ 2 − λ 2 ).In line with the results from Ref. [78], we find that opting for a smooth functional form for f (σ 2 ) significantly enhances the stability and efficiency of the adaptive timestepping in the integration routine.This improvement is illustrated in Fig. 11, where we plot the number of integration steps required as a function of time.
Appendix C: The low-rank dynamical truncation scheme We review here the ensemble truncation methods (or low-rank dynamical truncation schemes) introduced in Refs.[68][69][70].The starting assumption is that, at time t, the system is described by a LR density matrix in diagonal form ρ(t) = with | φj ⟩ the first M (t + dt) eigenvectors of T T † with the largest eigenvalues.The apparent difficulty in diagonalizing the N H × N H density matrix T T † is avoided by noticing that the M × M matrix T † T has the same nonvanishing eigenvalues as its adjoint with the associated eigenvectors being linearly related through T .
1. First-order equivalence of the two schemes In this section we set out to prove the equivalence, to first order in perturbation theory, between the schemes in Refs.[68][69][70] and the present LR-TDVP method summarized by the EOM (6).To do so we make the following two assumptions.First, we assume that at time t the system density matrix ρ(t) is in its diagonal spectral form (C1), that is: (C8) Second, we assume the rank M of ρ(t) to be unchanged at time t + dt.Within these assumptions, the variational Eq. ( 6) can be cast in the compact form Ḃ = L, (C9) Note that since the dynamical truncation scheme does not enforce trace preservation, in this section neither do we.Indeed, Eq. (C9) does not include the additional term −1 Tr{L}/M enforcing trace preservation in Eq. ( 6).The proof now consists in showing that the diagonal matrix (C7) produced by the dynamical truncation scheme coincides with that generated by Eqs.(C9) and (C10) to leading order in dt.To complete the proof, we express the density operator ρ(t + dt) to leading order in dt as ρ(t + dt) =  C10) with the last term on the rhs of Eq. (C13).The remaining terms in Eq. (C13), together with Eq. (C12), are the perturbation-theory expansion of the M × M matrix in Eq. (C11).This concludes the proof.
The relation between the truncation error ϵ M in Eq. (C6) and the approximated variational error in Eq. ( 10) is also made clear from the first-order perturbative expansion.Indeed, if j p j = 1, then where we made use of the fact that S −1 = 1 for an orthonormal basis.

FIG. 1 .
FIG. 1. Illustrative overview of the LR-TDVP Method.(a) Top: Pictorial representation of the density matrix ρ(t) (depicted in red) at a generic time t during the dynamics.The variational manifold selectively encompasses a specific region of the Hilbert space, capturing the predominant components of the statistical ensemble at time t.Bottom: Illustration of how memory efficiency is achieved by decomposing the NH × NH density matrix into smaller NH × M and M × M matrices.(b)Schematic representation of the dynamic adjustment of the variational manifold's dimension.Both the manifold and the rank of the truncated density matrix are dynamically modified to maintain a consistent level of accuracy, adapting to changes in the system's entropy.This adjustment reflects the principle that highly entropic systems necessitate a larger number of states for accurate representation, whereas low-entropy ensembles can be effectively described with fewer states.
FIG.2.Time evolution of the approximate variational error, χ(t), as defined in Eq.(10).The upper bound for χ is set at ϵmax = 10 −4 , which strictly applies for ∆t > 0. The inset displays the overlap between the true solution and the truncated one as obtained from Eq. (7).The physical model underlying this example is the XYZ Heisenberg model (detailed in Sec.III A) with N = 9 spins and Jy = 1.The values of all remaining physical parameters [see Eqs.(11) and (12)] are the same as those used in Fig.4.

FIG. 3 .
FIG. 3.Variation of the rank M and of the overlap O (between the true and approximated steady-state solutions) with the threshold ϵmax.We present results for all the choices of control parameter χ discussed in the text.The underlying physical model is the XYZ Heisenberg model with N = 6 spins and Jy = 1.The remaining physical parameters take the same values as in Fig.4.

FIG. 5 .
FIG.5.Von Neumann entropy S and its derivative ∂S/∂Jy as a function of the coupling parameter Jy for lattices with different numbers of spins N (legend in Fig.4).The inset displays the maximum value of the derivative as a function of the linear lattice length L = √ N .The solid line is a powerlaw fit of the finite-size scaling.We use the same physical parameters as in Fig.4.
FIG. 7.Memory allocations as a function of number of simulated spins N for the full simulation and the LR-TDVP one.We set the value of Jy = 0.98.For each value of N , we choose the rank M ensuring the relative error on the magnetization Mz to be smaller than 0.1%, namely err(Mz) = |M full

( 1 ) 1 , 2 =
FIG. 8.The first-order correlation function g (1) 1,2 (a) and the von Neumann entropy S (b) as a function of the amplitude of the two-photon driving G for the systems with N = 2, 3 cavities.Dashed grey lines are used to display the expected plateaux values of the different quantities in the limit of large driving.Parameters: γ = η = 1, U = 10, ∆ = J = −10.

FIG. 10 .
FIG. 10.Error probabilities of N -qubit Z gates, with N 1, 2, 3. (a) Phase-flip error probability PZ [Eq.(30)] as a function of the photon number α 2 .We display as a dashed gray line the analytical approximation PZ ≈ κ1T α 2 + ϵ 2 Z T /α 2 (for N = 1) provided in Ref.[117].(b) Bit-flip error probability PX [Eq.(31)] as a function of the photon number α 2 .To improve visualization, as the curves would largely overlap, we multiply each curve by a different constant, effectively shifting them vertically on the plot.The constants are R = 10 2 , 1, 10 −2 for N = 1, 2, 3 respectively.Full lines are used to display the exponential fit of the LR data.Because of computational constraints, the exact simulations for N = 3 were only possible up to α 2 = 5.The inset shows the von Neumann entropy S at time T as a function of α 2 for N = 1.In both panels, the results obtained from the L simulations (squares) are in agreement with the full solution (yellow crosses).Parameters: κ2 = 1, κ1 = 1/1000, ϵZ = ϵZZ = ϵZZZ = 1/20.

2 − λ 2 )f (σ 2 − λ 2 )FIG. 11 .
FIG. 11.Number of integration steps in of the adaptive integrator as a function of time.We simulate the evolution of the XYZ model with N = 9 spins and Jy = 1.The remaining physical parameters are set to the values discussed in Sec.III evolve while adjusting the rank with a threshold of ϵmax = 10 −4 .The regularization parameters are: atol = 10 −6 and rtol = 10 −5

p
j |φ j ⟩⟨φ j | + L[ρ(t)]dt = M j=1 p j |φ j ⟩⟨φ j | + N H j,k=1 L jk |φ j ⟩⟨φ k | dt , (C11)where we have completed the orthogonal basis of the full Hilbert space by introducing a set of orthogonal vectors {|φ j ⟩ , j = M + 1, . . ., N H }. The dynamical truncation scheme is achieved by diagonalizing the matrix ρ(t + dt) in Eq. (C11).The first term in Eq. (C11) is diagonal in the LR subspace H M , and thus has nonzero diagonal elements p j for j = 1, . . ., M only.The second term is non-diagonal and small.The eigenvalue perturbation theory, therefore, achieves the diagonalization to leading order in dt.According to perturbation theory, the new eigenvalues and eigenvectors are expressed as pj = p j + L jj dt , (C12)| φj ⟩ = |φ j ⟩ + k ⟩⟨φ k | L(ρ) |φ j ⟩ p j dt ,(C14)where we have distinguished the terms with j = 1, . . ., M from the other terms in the sum over states orthogonal to the LR subspace.Notice that in the last term of Eq. (C13) the sum over the projectors |φ k ⟩⟨φ k | factors out, as in the unperturbed matrix p k = 0 for k > M .Using that 1 − P = N H k=M +1 |φ k ⟩⟨φ k |, we immediately identify Eq. (