{'fraction_non_alphanumeric': 0.029836316965029835, 'fraction_numerical': 0.0070007000700070005, 'mean_word_length': 4.868153364632238, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A driven-dissipative quantum Monte Carlo method for open quantum systems

We develop a real-time Full Configuration Interaction Quantum Monte Carlo approach for the modeling of driven-dissipative open quantum systems. The method enables stochastic sampling of the Liouville-von-Neumann time evolution of the density matrix, thanks to a massively parallel algorithm, thus providing estimates of observables on the non-equilibrium steady state. We present the underlying theory, and introduce initiator technique and importance sampling to reduce the statistical error. Finally, we demonstrate the efficiency of our approach by applying it to the driven- dissipative two-dimensional XYZ spin model on lattice.


INTRODUCTION
The study of the nonequilibrium dynamics of manybody open quantum systems has attracted increasing attention in recent years, due to the progress in several experimental areas, including ultracold atomic gases, trapped ions, and superconducting circuits [1][2][3][4]. A common feature of these systems is the coupling to an external environment in the form of coherent or incoherent input and output channels. The time evolution of the system is then governed by the Liouville-von-Neumann equation which -in the case of stationary external conditions -typically drives it into a nonequilibrium steady state (NESS). Here, the competition between the coherent and incoherent dynamics gives rise to a multitude of novel phenomena, including nonequilibrium dissipative phase transitions [5,6].
Generally, the steady state density matrix can be obtained in two different ways. First, one can integrate the time evolution until the stationary state is reached, or second, one can find the solution associated with the null eigenvalue of the Liouvillian super-operator. The theoretical description of open many body systems represents a major challenge, and in spite of the numerous improvements, the numerical modeling can be handled only for small system sizes. While several studies have been restricted to mean-field approximations [7][8][9], in the case of one-dimensional systems highly accurate results were obtained via the density matrix renormalization group (DMRG) technique [10][11][12], and the equivalent variational approach based on the matrix product state (MPS) ansatz [13,14]. Alternatively, a stochastic method, i.e the Monte Carlo Wave Function (MCWF) [15,16], can be used which unravels the system-bath interaction onto a stochastic process that adds to the unitary Hamiltonian dynamics and to the effective damping terms. Furthermore, a spatial renormalization approach -the corner space renormalization method -was recently proposed, which relies on an ad hoc spatial decimation protocol for the density matrix. [17,18].
For closed, Hamiltonian systems, various quantum Monte Carlo approaches have been the election tools to stochastically sample system properties, both at zero and finite temperature [19,20]. A class of methods generally known as projector Monte Carlo (PMC) [21] -such as diffusion Monte Carlo (DMC) [22,23] and Green's function Monte Carlo (GFMC) [24] -enables modeling the ground state properties at zero temperature by stochastically sampling the time evolution of the imaginary-time Schrödinger equation. However, PMC methods may suffer from the sign problem which results in a computational cost growing exponentially with the system size. Recently, a novel approach called Full Configuration Interaction Monte Carlo (FCIQMC) has been introduced for quantum chemistry simulations [25][26][27][28] and for the simulation of other strongly correlated systems [29,30]. As a projector technique for zero temperature, it has features in common with DMC and GFMC, although a radically different sampling protocol was introduced which has proven to significantly alleviate the sign problem [26,30,31]. Nevertheless, the sign problem is NP-hard problem [32],and is not completely solved by FCIQMC.
A mutual feature of the Liouvillian dynamics and the imaginary-time Schrödinger equation is the fact, that in the long-time limit, the eigenstate with the smallest-realpart-eigenvalue will dominate. In the Liouvillian case, where the eigenvalues of the Liouvillian super-operator are complex valued with negative real part [33], this corresponds to the null eigenvalue solution -to the NESS [ Table I]. It would therefore be natural to apply projector Monte Carlo techniques to the simulation of the NESS. In this case however, we expect the sign problem to be highly relevant, as the complex valued density matrix cannot be expected to possess elements of the same sign only.
In this paper, we develop a real-time FCIQMC approach to open quantum systems, which we call drivendissipative quantum Monte Carlo (DDQMC). DDQMC shares many of the features of FCIQMC, but it samples the elements of the complex-valued density matrix instead of the wavefunction. The method does not truncate the Hilbert-space and contrary to tensor network methods, its applicability is not bound to the dimensionality of the system. In order to demonstrate the use of DDQMC, we simulate a two-dimensional spin lattice governed by the Heisenberg XYZ Hamiltonian interacting with a dissipative environment. The study of this model has recently attracted interest as an example of a dissipative phase transition resulting from the non-trivial competition between the coherent and incoherent dynam- Real-timė ρ(t) = Lρ

Long-time limit
Ground state Nonequilibrium steady state ics [18,34]. The single-site Gutzwiller mean-field study of the system predicts a phase transition from a paramagnetic phase to a magnetically ordered one [35,36], while a recent analysis showed that this transition should survive in 2D and disappear in case of one-dimensional lattices [37]. The paper is organized as follows. In Section II we give an overview of the original FCIQMC algorithm, setting the basis to the formulation of DDQMC in Section III. In Section IV the initiator approach and the importance sampling are introduced. The method is then applied to the XYZ Heisenberg lattice in Section V and the results are compared to those obtained by an optimized exact diagonalization method and by quantum trajectories. We finally discuss the effectiveness of the approach and offer some concluding remarks.

II. OVERVIEW OF THE FCIQMC ALGORITHM
We begin by giving a short overview of the FCIQMC method. For a more complete derivation readers are referred to Refs. [25,26,31,38,39].
In general, PMC methods are stochastic implementations of the power method which aims at computing the expectation values of operators on the dominant eigenstate of the projector. They prove to be particularly useful when the Hilbert-space is so large that the storage of matrix and vector representations becomes computationally unfeasible. PMC techniques get around this memory limitation by storing at any instant in time only a random sample of vector and matrix elements. The expectation values are then computed as time averages. In common with PMC, FCIQMC also performs the long-time integration of the imaginary-time Schrödinger equation. However, unlike PMC, this is achieved with a completely different sampling strategy.
Consider the imaginary-time Schrödinger equation (we assume here and in what follows = 1) the general solution of the equation is Once expanded onto the basis spanned by the eigenvectors of the Hamiltonian {|φ i }, the wavefunction results in a sum of exponentially decaying terms. In order to prevent the ground state component from decaying in the infinite time limit, a constant energy shift E 0 can be introduced, where E 0 is the ground state energy. Since the value of E 0 is unknown in advance, one solves the equation with an arbitrary shift Ṡ |ψ = −(Ĥ − S1)|ψ .
During the simulated time evolution the value of the S is slowly adjusted in order to maintain a constant normalization and -at convergence -provides an estimate of the actual ground state energy E 0 . FCIQMC stochastically samples the first order Euler approximation of eq. (3). Furthermore, the algorithm works on a discrete basis set {|φ i }, and the Hamiltonian and the wavefunction are projected onto the space spanned by the basis elements The evolution of the expansion coefficients is then governed by where H ij = φ i |Ĥ|φ j . In order to stochastically represent eq. (5), we introduce a fundamental unit called walker. Each walker has a sign (q = ±1), and contributes to sampling the amplitude of one of the |φ i basis states. Let n + i be the number of walkers with positive sign on a given state and n − i that of walkers with negative sign. Then the amplitude of a basis state in the expansion is proportional to the net walker number residing there: c i ∝ n + i − n − i . Starting from an initial distribution, the walkers evolve following a set of rules designed to sample the time evolution of eq. (5) over a time step. This dynamics is then iterated until convergence is reached.
We point out here, that this approach sets an upper bound to the time step. Since the initial state is driven into the dominant eigenstate of the projector, the solution will converge to the ground state only if |1 − ∆τ (E i − S)| ≤ 1 for all eigenvalues E i . This corresponds to the requirement ∆τ ≤ 1/∆, where ∆ is the full spectral width of the Hamiltonian under consideration.
The rules for evolving the walker population can be summarized as follows. At each time step we loop over the entire walker population and perform the following operations: (i) Spawning: For a walker residing on site i a connected site j is chosen randomly and a spawning event is made possible with a probability p(j|i) ∝ |H ji |∆τ (connected sites are linked by non-zero off-diagonal Hamiltonian elements H ij ). If the attempt is successful, walkers are born at site j with sign q j = sign(H ji )q i . If p(j|i) > 1 then the corresponding integer number of walkers are realized deterministically, and the fractional part stochastically [25].
(ii) Clone/Death: For each walker on a given site, a death event is sampled with probability If p death > 0, the walker is removed from the population. If p death < 0, a walker of opposite sign is created. In case |p death | > 1, the integer part of p death is realized deterministically, and the fractional part stochastically.
(iii) Annihilation: pairs of walkers of opposite sign residing on the same basis state are annihilated. Therefore, at the end of each time step, each state is solely occupied by walkers having the same sign.
One of the most significant advantages of FCIQMC is due to the aforementioned annihilation procedure. This step does not alter the evolution of the ground state, but was shown to be crucial in systems with sign problem [26]. The sign problem in FCIQMC simulation is manifested as the fast growth of an unphysical solution dominating the ground-state result. The annihilation procedure can suppress this growth and allow the simulation to converge to the physical solution, but only if a minimal and system-dependent walker population is present. Henceforth, building on the massively parallel nature of the method, a computationally efficient implementation can offer an insight into the study of systems with severe sign problem.

III. DRIVEN-DISSIPATIVE QUANTUM MONTE CARLO
We describe now how the dynamics of open quantum systems following the Liouville-von-Neumann equation can be cast into a Monte Carlo algorithm.

A. Theory
The general problem we aim to solve is that of quantum system with several degrees of freedom, in the presence of external driving fields and Markovian coupling to the environment. The evolution of the steady matrixρ is then governed by the Liouville-von-Neumann master equation [40] The dissipative part of the dynamics is described by whereF i are the jump operators, characterizing the transitions induced by the environment, and γ i are the corresponding transition rates. Contrary to the Hamiltonian, the Liouville superoperator is not hermitian. Dissipative systems evolve under a one parameter semi-group (e Lt , t > 0), generated by the Liouvillian, resulting in a time evolution which is no longer unitary. Due to its non-hermiticity, L has complex eigenvalues with negative real part. It can be shown that the density matrix, under very general assumptions will evolve into an asymptotic steady state, corresponding to the zero eigenvalue of L [33]. By introducing an additional shift into eq. (7) for a diagonal population control, Eq. (9) can be stochastically sampled similarly to the Hamiltonian case and the NESS is obtained as a Monte Carlo average of the long time limit. As in the case of FCIQMC, we take the first-order Euler approximation, and we introduce a set of walkers which now sample the amplitudes of basis operators |φ i φ j |, from now on referred as "configurations". The stochastic sampling of the unnormalized density matrix gives access to the expectation value of any quantum mechanical observable. The expectation value of observableÔ at a given instant in time is computed by Once the simulation has asymptotically approached the steady state -i.e. the shift S reached the steady state eigenvalue of the Liouvillian, S = 0 -the numerator and the denominator can be averaged separately over a sufficiently large number of iteration steps.

B. Multinomial formalism
The original FCIQMC sampling protocol was described, e.g. in Ref. [25]. Here we developed a variant which optimizes the computational cost of the evolution generated by off-diagonal Hamiltonian elements. In the followings we will refer to this variant as multinomial formalism.
In the original scheme, in order to perform the stochastical evolution induced by the off-diagonal H ij elements in eq. (5), the algorithm requires a loop over the entire walker population at each time step. This method becomes computationally heavy as the walker population increases.
Here we introduce an alternative strategy for the spawning generation. Let p(j|i) be the probability of choosing the j-th child starting from site i Then the number of actual spawning events N i sp occurring for N i walkers residing on site i is determined by a stochastic process following a binomial distribution Therefore, in each time step we perform a loop over the currently populated basis states rather than the whole walker population. In systems with local coupling the Hamiltonian is represented by a highly sparse matrix, and a computationally effective state representation makes the extra memory allocation negligible (for an occupied site i it is necessary to store all the possible connected states with the corresponding probabilities). Efficient algorithms for binomial and multinomial random number generation are also present in the literature [41][42][43].

C. Algorithm
The dynamics of the walker population is determined by a set of rules designed to stochastically sample eq. (10). However, the elements of the density matrix are complex valued. Similarly to [29], we can sample a complex density matrix with two types of walkers, respectively for the real and imaginary parts. If the density matrix is expressed in vectorized form, the shifted Liouvillian superoperator can be expressed im matrix form using Kronecker products as [44] Then eq. (9) can be written in the form of where L lm ij are the matrix elements of the superoperator. Here ρ ij represents the now complex valued population on a given configuration |φ i φ j |.
For the sampling protocol we use the multinomial formalism, introduced in Sec. III B. Let us refer to this scheme as a function which returns randomly drawn number of walkers spawned from configuration ij to lm given the matrix element connecting them. At each time step, we loop over the list of currently occupied configurations and perform the following steps: (i) Spawning: Consider the complex walker population residing on ρ ij and perform spawning to all the connected configurations. The real (ℜ) and imaginary (ℑ) parts of L lm ij are considered in turn and two spawning attempts are realized respectively for real and imaginary parents. The number of walkers spawned to each ρ lm are determined by the multinomial formalism.

For real parents
and for imaginary parent walkers where N ℜ lm and N ℑ lm are the number of real and imaginary walkers being spawned to configuration ρ lm and 'sign' indicates the sign of the progeny.
(ii) Clone/Death: This step is required as a real (imaginary) walker can produce an imaginary (real) walker on the same configuration. The spawning occurs on-site with a population determined by the binomial distribution. The real and imaginary parts of L ij ij are considered in turn and two spawning attempts are realized respectively for the real and imaginary population.
(iii) Annihilation: on a given site the real and imaginary population are considered in turn, and pairs of walkers having opposite signs are removed from the simulation.
Unlike in FCIQMC, here the target value of the diagonal shift S is known. In the infinite time limit the master equation is assumed to drive the density matrix into the steady state, which by definition belongs to the zero eigenvalue of the Lindbladian. Therefore, the time evolution of the shift S will indicate if convergence is reached and one can start gathering statistics on the observables. As before, the shift is slowly adjusted in order to maintain a constant walker population. Since estimators for most operators of interest only receive contributions from walkers on or near the diagonal elements, we chose to control the amount of population distributed along the diagonal of the density matrix. The value of the shift is then adjusted according to the familiar shift-update algorithm implemented in FCIQMC [25] calculations where δ is a damping parameter, and N w is the total weight of real walkers residing on diagonal density matrix elements. The method does not have a built-in constraint on the diagonal elements being real. The value of the imaginary part fluctuates around zero and its expectation value naturally vanishes during the simulation.

IV. INITIATOR APPROACH & IMPORTANCE SAMPLING
The algorithm described in Section III allows the stochastical sampling of the steady state density matrix of open quantum systems whose dynamics follows the Liouville-von-Neumann equation. However, with the increase in the number of configurations, the walker population tends to become dilute, resulting in a situation where the simulation contains only a few walkers per basis operator. This leads to an increased statistical error thus reducing dramatically the accuracy of the sampling. In order to address this issue one needs to increase the walker number in the system which -with large system sizes -becomes computationally unfeasible.
In order to overcome this issue two different methods were introduced: initiator approach and importance sampling. Each of the techniques reduces the minimal required walker population by decreasing the number of simultaneously occupied configurations, however the strategy of selecting the configurations to be sampled is fundamentally different. The initiator approach allows the significant configurations to emerge naturally during the simulation, whilst importance sampling gives the possibility to drive the walker population to a selected subset of presumably relevant configurations.
These methods can improve the sampling quality with great success, however, they have to be applied carefully, as both introduce a bias on the result.

A. Initiator approach
Our initiator approach is based on an FCIQMC adaptation detailed in [28]. For the newly spawned walkers an additional survival criteria is introduced which can dramatically reduce the population required to reach convergence. Some of the basis operators are tagged as 'initiator' which have the ability to spawn progeny onto unoccupied basis states. However, progeny of the noninitiator states can only survive if they spawn to basis operators which were previously occupied or to diagonal elements. All the diagonal basis operators are initiators by definition, and during the simulation a basis element might become initiator, if its population exceeds a preset value (I limit ).
This method results in a series of systematically improvable approximations which will tend to the original algorithm in three limits: (i) with decreasing I limit every basis element will become initiator. All the progeny survives regardless of the parent state, which is equivalent to the original method; (ii) with increasing total population all basis element will acquire walkers, therefore, all spawned children will survive regardless of the flag of the parent state; (iii) extending the initiator space by definition will result again in all the basis operators becoming initiators, consequently all progeny will survive.
Setting an initiator limit, introduces a dynamical truncation on the available basis operators, leading to a biased result. In order to compute the unbiased expectation values, we progressively decrease the initiator limit in different simulations and fit the estimated expectation values, thus extrapolating the value in the limit I limit → 0 ( Figure 1).

B. Importance sampling
We start by identifying the basis elements whose sampling needs to be improved. Then a straightforward way to do so is by reducing the probability of spawning out these configurations. Walkers that do reside on unessential elements are given a correspondingly larger weight, hence the expectation values of the observables will be unchanged. We define the following simple importance sampling procedure.
The evolution of the density matrix in the DDQMC formalism follows eq. (16). In order to associate weights to the matrix elements depending on their importance we multiply them a factor w ij The time evolution of the importance-sampled density matrix then becomes which is fully analogous to eq. (16) and can be simulated by the DDQMC method. The expectation value of an observableÔ is thus In this work we introduce a single importance sampling parameter p > 0, and give all the off-diagonal elements a weight w ij = e −p . Meanwhile the diagonal coefficients are not altered. This strategy focuses on sampling the diagonal density matrix elements and gives an easy access to tune the strength of the importance sampling. Fig. 2 shows the amount of simultaneously occupied density matrix elements before and after using the initiator approach and importance sampling given the XYZ Heisenberg lattice model that we describe later.

V. RESULTS
In order to demonstrate the effectiveness of DDQMC, we simulated the two-dimensional spin-1/2 XYZ Heisenberg lattice in the presence of a dissipating channel which tends to relax each spin into the |s z = −1/2 state.
whereŜ α i are the spin operators matrices acting on the i-th spin, J α are the coupling constants between nearest neigbour spins, γ is the dissipation rate, andŜ ± j = S x j ± iŜ y j . Recently, the system has attracted significant interest since the competition between the coherent Hamiltonian dynamics and the incoherent spin flips leads to a dissipative phase transition. Due to the anisotropy, the Hamiltonian part induces a nonzero spin expectation value on the xy plane, while dissipation would drive each site to the spin down state. This competition leads to a phase transition from a paramagnetic phase (with no magnetization in the xy plane) to a magnetically organized phase (which presents a finite polarization in the xy plane). Both the Gutzwiller mean-field theory [35][36][37] and the corner space renormalization method [18,34] predicts this dissipative phase transition.
In order to study the model we chose three different lattice sizes: 2 × 2, 3 × 3 and 4 × 4. The first two sizes are small enough to derive an exact numerical solution of the master equation in the steady state, thus allowing a direct check of the accuracy of our DDQMC results. In case of the 4×4 lattice, the magnetization is compared to those obtained by Monte Carlo wave function technique.

A. Magnetization in the steady-state
The steady-state magnetization per site is defined as Tr(ρσ z i ), (28) where N is the number of lattice sites. Figure 3(a) shows the magnetization of the 3 × 3 lattice as a function of the Monte Carlo iteration step with a diagonal population of 50 4 walkers. The exact solution obtained by directly solving the linear system is also plotted. Increasing the diagonal population to 20 6 reduces the statistical error as seen in the corresponding result in Fig. 3(b).
In Fig. 4, we present the magnetization per site M z as a function of the normalized coupling constant J y /γ for square lattices of different size. The exact and MCWF solutions are also plotted and are in agreement with the results obtained by the DDQMC algorithm.

B. Angularly-averaged susceptibility
Following the scheme presented in [18], we study the system in the presence of an applied magnetic field in the xy plane The linear response is then summarized in the 2 × 2 susceptibility tensor with α, β = x, y.
For a more complete derivation readers are referred to [18]. In Fig. 5, we present the angularly-averaged susceptibility χ av as a function of the normalized coupling parameter J y /γ for different lattice sizes.
The magnetic susceptibility for the different lattice sizes exhibits a peak of increasing height which qualitatively corresponds to the results obtained in [18].

VI. OUTLOOK
The DDQMC method presented here, constitutes a basic PMC approach to the non-equilibrium steady state of open quantum systems. As such, it contains only the basic building blocks of the PMC method, and its effectiveness may be considerably improved by introducing any of the several tools that are common in other PMC schemes. Here, we describe as an outlook three such possible improvements.
The first possible improvement consists in the implementation of a mixed-estimator scheme, in analogy with the one used in projector and diffusion Monte Carlo to find the ground state of Hamiltonian systems [19][20][21]. Here, a possible mixed estimator strategy may consist in formally carrying out an exact real-time evolution, starting from a DDQMC sampled density matrix. More specifically, let us assume that at time t the current DDQMC sample of the density matrix isρ(t). We can formally apply the exact time evolution for an additional time T and then evaluate the expectation value of an observableÔ at time t + T as Ô = Tr(Ôe LTρ (t)) Tr(e LTρ (t)) = Tr(Ôe LTρ (t)) Tr(ρ(t)) , where the second equality results from the trace preserving character of the time evolution. In the limit T → ∞, Eq. (33) provides the steady-state expectation value independently of the actual value of the density matrixρ(t), when assuming that a unique steady state exists. A mixed estimator strategy would then consist in building a "trial" observableÔ T which can still be efficiently computed element-wise, and such thatÔ T ≃Ô H (T ) =Ôe LT . Here,Ô H (t) represents the Heisenberg picture of the observableÔ and, for time-independent Liouvillian maps, it obeys the adjoint quantum master equation dÔH (t) dt = L †Ô H (t) [45]. Hence, the mixed estimator approach in the present case would require the knowledge of an approximate time dependence forÔ H (t), which may be obtained, for instance, from a time-dependent variational principle [46,47] applied to a separable or short-range-correlated ansatz for the observable.
A second improvement would consist in using a "guiding density matrix" for the importance sampling. A natural choice for such a guiding density matrix would again be a variational ansatz, as the variational principle for the NESS is now well established [6], and some variational approaches have already been developed [48][49][50].
Finally, the present scheme is based on the Euler method (10) for the numerical solution of the timedependent master equation. The Euler method is a first order method in the time step, and is only stable if ∆t is chosen to be smaller than the inverse of the full spectral width of the master equation. In PMC, several approaches have been proposed to sample a higher-order discrete-timestep propagator [23], or even the exact one [51]. While a similar approach would be highly beneficial to FCIQMC and DDQMC, the question is still open, whether higher-order propagators may be efficiently sampled within the spawn-annihilation sampling protocol characterizing these Monte Carlo methods.

VII. CONCLUSIONS
We have introduced a quantum Monte Carlo approach to open many-body quantum systems, called DDQMC. The method is based on the FCIQMC algorithm exploiting the analogy between the long-time dynamics of the Lindbladian master equation and the imaginary-time Schrödinger equation. DDQMC allows direct sampling of the steady state density matrices in any discrete basis set, and in all cases studied it has proven to be accurate.
DDQMC, as FCIQMC, uses an annihilation procedure which helps to alleviate the sign problem. The introduction of the initiator approach and importance sampling can lead to a significant improvement in the statis-tical accuracy and reduce the required walker population. The validity of the method was proven by investigating a dissipative phase transition on the two-dimensional Heisenberg-model. The defining feature of DDQMC is that it samples the whole density matrix and it does not introduce a truncation in Hilbert-space. Experience showed that the applicability of the code does not solely depend on the system size, but also on the correlations characterizing the steady state. The application presented in this work is a proof of principle, demonstrating the possibility to stochastically sample the Liouville-von-Neumann equation in a finite difference approximation. DDQMC holds promise as a powerful tool in the study of open quantum systems.