Variational Quantum Monte Carlo Method with a Neural-Network Ansatz for Open Quantum Systems

TL;DR

This paper introduces a variational Monte Carlo approach using neural networks to efficiently simulate the steady states of large open quantum systems, overcoming exponential complexity in density matrix representation.

Contribution

It develops a novel neural network-based variational method combined with stochastic reconfiguration to simulate non-equilibrium steady states of open quantum systems.

Findings

Successfully modeled the 2D dissipative XYZ spin model.

Demonstrated efficiency in handling large system sizes.

Validated the method against known quantum system behaviors.

Abstract

The possibility to simulate the properties of many-body open quantum systems with a large number of degrees of freedom is the premise to the solution of several outstanding problems in quantum science and quantum information. The challenge posed by this task lies in the complexity of the density matrix increasing exponentially with the system size. Here, we develop a variational method to efficiently simulate the non-equilibrium steady state of Markovian open quantum systems based on variational Monte Carlo and on a neural network representation of the density matrix. Thanks to the stochastic reconfiguration scheme, the application of the variational principle is translated into the actual integration of the quantum master equation. We test the effectiveness of the method by modeling the two-dimensional dissipative XYZ spin model on a lattice.

Figures

Equations

$$\frac{d\hat{\rho }}{dt}=-i[\hat{H},\hat{\rho }]-\sum _i\frac{{\gamma }_i}{2}\left[\left\{{\hat{F}}_i^{†}{\hat{F}}_i,\hat{\rho }\right\}-2{\hat{F}}_i\hat{\rho }{\hat{F}}_i^{†}\right],$$

$$\mathcal{L}({\hat{\rho }}_{ss})=0,$$

$${\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })=\sum _{j=1}^J{p}_j(\chi )\cdot {\psi }_j(\mathbf{\sigma },\chi ){\psi }_j^{\ast }(\mathbf{\eta },\chi )$$

$$\begin{array}{rl}{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })& =8\exp \left(\sum _i{a}_i{\sigma }_i\right)\exp \left(\sum _i{a}_i^{\ast }{\eta }_i\right)\\ & \times \prod _{l=1}^L\cosh \left(c_l+\sum _i{W}_{li}{\sigma }_i+\sum _i{W}_{li}^{\ast }{\eta }_i\right)\\ & \times \prod _{m=1}^M\cosh \left(b_m+\sum _i{X}_{mi}{\sigma }_i\right)\\ & \times \prod _{n=1}^M\cosh \left(b_n^{\ast }+\sum _i{X}_{ni}^{\ast }{\eta }_i\right).\end{array}$$

$$\chi (n+1)=\chi (n)+\nu \cdot S^{-1}(n)F(n),$$

$$\begin{array}{rl}{O}_k(\mathbf{\sigma },\mathbf{\eta })& =\frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial {\chi }_k}\\ F_k(n)& =⟨⟨O_k^{\ast }\mathcal{L}⟩⟩-⟨⟨\mathcal{L}⟩⟩⟨⟨O_k^{\ast }⟩⟩\\ S_{k{k}^{\prime }}(n)& =⟨⟨O_k^{\ast }{O}_{k^{\prime }}⟩⟩-⟨⟨O_k^{\ast }⟩⟩⟨⟨O_{k^{\prime }}⟩⟩,\end{array}$$

$$⟨\hat{\mathcal{O}}⟩=\mathrm{Tr}(\hat{\mathcal{O}}{\hat{\rho }}_{\chi })=\sum _{\mathbf{\sigma },\mathbf{\eta }}|{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta }){|}^2\cdot \frac{\mathcal{O}(\mathbf{\eta },\mathbf{\sigma })}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta }{)}^{\ast }},$$

$$\hat{H}=\sum _{⟨i,j⟩}\left(J_x{\hat{\sigma }}_i^x{\hat{\sigma }}_j^x+J_y{\hat{\sigma }}_i^y{\hat{\sigma }}_j^y+J_z{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z\right)$$

$$\frac{\mathrm{d}\hat{\rho }}{\mathrm{d}t}=-i[\hat{H},\hat{\rho }]-\frac{\gamma }{2}\sum _k\left[\left\{{\hat{\sigma }}_k^{+}{\hat{\sigma }}_k^{-},\hat{\rho }\right\}-2{\hat{\sigma }}_k^{-}\hat{\rho }{\hat{\sigma }}_k^{+}\right]$$

$$M_z=\frac{1}{N}\sum _{i=1}^N\mathrm{Tr}(\hat{\rho }{\hat{\sigma }}_i^z),$$

$$S_{ss}^{xx}(\text{k})=\frac{1}{N(N-1)}\sum _{\text{j}\ne \text{l}}{e}^{-i\text{k}(\text{j}-\text{l})}⟨{\hat{\sigma }}_{\text{j}}^x{\hat{\sigma }}_{\text{l}}^x⟩,$$

$${\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })=\sum _{j=1}^J{p}_j(\chi )\cdot {\psi }_j(\mathbf{\sigma },\chi ){\psi }_j^{\ast }(\mathbf{\eta },\chi )$$

$$\psi (\mathbf{\sigma })=\sum _{\{q\}}\exp \left(\sum _i{a}_i{\sigma }_i+\sum _m{b}_m{q}_m+\sum _{m,i}{q}_m{\sigma }_i{X}_{mi}\right)$$

$$\begin{array}{rl}{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })& =\sum _{\{h\}}\sum _{\{q\}}\sum _{\{r\}}\exp \left(\sum _l{c}_l{h}_l\right)\\ & \times \exp \left(\sum _i{a}_i{\sigma }_i+\sum _m{b}_m{q}_m+\sum _{m,i}{q}_m{\sigma }_i{X}_{mi}+\sum _{l,i}{h}_l{\sigma }_i{W}_{li}\right)\\ & \times \exp \left(\sum _i{a}_i^{\ast }{\eta }_i+\sum _n{b}_n^{\ast }{r}_n+\sum _{n,i}{r}_n{\eta }_i{X}_{ni}^{\ast }+\sum _{l,i}{h}_l{\eta }_i{W}_{li}^{\ast }\right).\end{array}$$

$$\begin{array}{rl}{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })& =8\exp \left(\sum _i{a}_i{\sigma }_i\right)\exp \left(\sum _i{a}_i^{\ast }{\eta }_i\right)\\ & \times \prod _{l=1}^L\cosh \left(c_l+\sum _i{W}_{li}{\sigma }_i+\sum _i{W}_{li}^{\ast }{\eta }_i\right)\\ & \times \prod _{m=1}^M\cosh \left(b_m+\sum _i{X}_{mi}{\sigma }_i\right)\\ & \times \prod _{n=1}^M\cosh \left(b_n^{\ast }+\sum _i{X}_{ni}^{\ast }{\eta }_i\right).\end{array}$$

$$\begin{array}{rl}{\tilde{\theta }}_l(\mathbf{\sigma },\mathbf{\eta })& =c_l+\sum _i{W}_{li}{\sigma }_i+\sum _i{W}_{li}^{\ast }{\eta }_i\\ {\theta }_m(\mathbf{\sigma })& =b_m+\sum _i{X}_{mi}{\sigma }_i,\end{array}$$

$${\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })=8\exp \left(\sum _i{a}_i{\sigma }_i\right)\exp \left(\sum _i{a}_i^{\ast }{\eta }_i\right)\prod _{l=1}^L\prod _{m=1}^M\prod _{n=1}^M\cosh {\tilde{\theta }}_l(\mathbf{\sigma },\mathbf{\eta })\cosh {\theta }_m(\mathbf{\sigma })\cosh {\theta }_n^{\ast }(\mathbf{\eta }).$$

$$\begin{array}{rl}\frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Re (a_k)}& ={\sigma }_k+{\eta }_k\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Im (a_k)}& =i({\sigma }_k-{\eta }_k)\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Re (b_k)}& =\tanh \left({\theta }_k(\mathbf{\sigma })\right)+\tanh \left({\theta }_k^{\ast }(\mathbf{\eta })\right)\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Im (b_k)}& =i\left[\tanh \left({\theta }_k(\mathbf{\sigma })\right)-\tanh \left({\theta }_k^{\ast }(\mathbf{\eta })\right)\right]\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Re (X_{kl})}& =\tanh \left({\theta }_k(\mathbf{\sigma })\right){\sigma }_l+\tanh \left({\theta }_k^{\ast }(\mathbf{\eta })\right){\eta }_l\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Im (X_{kl})}& =i\left[\tanh \left({\theta }_k(\mathbf{\sigma })\right){\sigma }_l-\tanh \left({\theta }_k^{\ast }(\mathbf{\eta })\right){\eta }_l\right]\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial c_k}& =\tanh \left({\tilde{\theta }}_l(\mathbf{\sigma },\mathbf{\eta })\right)\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Re (W_{kl})}& =\tanh \left({\tilde{\theta }}_l(\mathbf{\sigma },\mathbf{\eta })\right)({\sigma }_l+{\eta }_l)\\ \frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial \Im (W_{kl})}& =i\tanh \left({\tilde{\theta }}_l(\mathbf{\sigma },\mathbf{\eta })\right)({\sigma }_l-{\eta }_l)\end{array}$$

$${\chi }_k^{\prime }={\chi }_k+\delta {\chi }_k.$$

$$|{\rho }_{{\chi }^{\prime }}⟩\simeq \sum _k\delta {\chi }_k{O}_k|{\rho }_{\chi }⟩,$$

$$O_k(x)=O_k(\mathbf{\sigma },\mathbf{\eta })=\frac{1}{{\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}\cdot \frac{\partial {\rho }_{\chi }(\mathbf{\sigma },\mathbf{\eta })}{\partial {\chi }_k}.$$

$$|{\bar{\rho }}_{\chi }⟩=e^{\mathcal{L}t}|{\rho }_{\chi }⟩\simeq (1+\mathcal{L})|{\rho }_{\chi }⟩.$$

$$\delta {\chi }_k=\nu \sum _{k^{\prime }}{S}_{k{k}^{\prime }}^{-1}{F}_{k^{\prime }},$$

$$\begin{array}{rl}{F}_k(n)& =⟨⟨O_k^{\ast }\mathcal{L}⟩⟩-⟨⟨\mathcal{L}⟩⟩⟨⟨O_k^{\ast }⟩⟩\\ S_{k{k}^{\prime }}(n)& =⟨⟨O_k^{\ast }{O}_{k^{\prime }}⟩⟩-⟨⟨O_k^{\ast }⟩⟩⟨⟨O_{k^{\prime }}⟩⟩.\end{array}$$

$$S_{k{k}^{\prime }}^{\mathrm{reg}}=S_{k{k}^{\prime }}+\lambda (n){\delta }_{k,k^{\prime }}{S}_{k{k}^{\prime }}.$$

$$\begin{array}{rl}{F}_k(n)& =\sum _x|⟨x|{\rho }_{\chi }⟩{|}^2\cdot {\left(\frac{\partial \ln ⟨x|{\rho }_{\chi }⟩}{\partial {\chi }_k}\right)}^{\ast }\frac{⟨x|\mathcal{L}|{\rho }_{\chi }⟩}{⟨x|{\rho }_{\chi }⟩}\\ & -\sum _x|⟨x|{\rho }_{\chi }⟩{|}^2\cdot \frac{⟨x|\mathcal{L}|{\rho }_{\chi }⟩}{⟨x|{\rho }_{\chi }⟩}\sum _{x^{\prime }}|⟨x^{\prime }|{\rho }_{\chi }⟩{|}^2\cdot {\left(\frac{\partial \ln ⟨x^{\prime }|{\rho }_{\chi }⟩}{\partial {\chi }_k}\right)}^{\ast }\\ S_{k{k}^{\prime }}(n)& =\sum _x|⟨x|{\rho }_{\chi }⟩{|}^2\cdot {\left(\frac{\partial \ln ⟨x|{\rho }_{\chi }⟩}{\partial {\chi }_k}\right)}^{\ast }\left(\frac{\partial \ln ⟨x|{\rho }_{\chi }⟩}{\partial {\chi }_{k^{\prime }}}\right)\\ & -\sum _x|⟨x|{\rho }_{\chi }⟩{|}^2\cdot {\left(\frac{\partial \ln ⟨x|{\rho }_{\chi }⟩}{\partial {\chi }_k}\right)}^{\ast }\sum _{x^{\prime }}|⟨x^{\prime }|{\rho }_{\chi }⟩{|}^2\cdot \left(\frac{\partial \ln ⟨x^{\prime }|{\rho }_{\chi }⟩}{\partial {\chi }_{k^{\prime }}}\right),\end{array}$$

$$\mathcal{A}(|x_y⟩\to |x_z⟩)=\min \left(1,\frac{|⟨x_z|{\rho }_{\chi }⟩{|}^2\cdot \mathcal{P}(y|z)}{|⟨x_y|{\rho }_{\chi }⟩{|}^2\cdot \mathcal{P}(z|y)}\right).$$

Full text

Open quantum systems have evolved into a major field of studies in recent years. Focus of these studies are the characterization of emergent phenomena and dissipative phase transitions Carusotto and Ciuti (2013); Hartmann (2016); Noh and Angelakis (2017); Bartolo et al. (2016); Biella et al. (2017); Biondi et al. (2017); Carmichael (2015); Casteels et al. (2017, 2016); Fink et al. (2017, 2018); Fitzpatrick et al. (2017); Foss-Feig et al. (2017); Kessler et al. (2012); Marino and Diehl (2016); Savona (2017); Sieberer et al. (2013); Vicentini et al. (2018); Jin et al. (2016); Lee et al. (2013); Rota et al. (2018, 2017); Casteels et al. (2018), as well as the ongoing debate about whether quantum computing schemes are still hard to simulate classically – and thus achieve quantum supremacy – when in presence of some degree of noise-induced decoherence Bremner et al. (2016); Gao and Duan (2018); Harrow and Montanaro (2017); Preskill (2018).

Assuming a Markovian interaction with the environment, the dynamics of open quantum systems is governed by the quantum master equation in Lindblad form Breuer and Petruccione (2002). Only few models within this description admit an analytical solution Prosen (2011, 2014). The quest for efficient numerical methods to simulate the dynamics and the asymptotic steady state resulting from the Lindblad master equation, is a research field that is still in its infancy. Many recent tools have been developed following in the footsteps of well established numerical methods for the simulations of closed, Hamiltonian quantum systems. In particular, matrix-product state and tensor network schemes Cui et al. (2015); Kshetrimayum et al. (2017); Mascarenhas et al. (2015); Werner et al. (2016), a real-space renormalization approach Finazzi et al. (2015), cluster mean-field Jin et al. (2016), and other ad-hoc approximation schemes Casteels et al. (2018); Nagy and Savona (2018) have recently emerged.

A groundbreaking progress in the numerical simulation of both the ground state and the dynamics of closed quantum systems has recently been made with the introduction of the neural-network variational ansatz Cai and Liu (2018); Carleo et al. (2018); Carleo and Troyer (2017); Freitas et al. (2018); Glasser et al. (2018); Nomura et al. (2017); Gao and Duan (2017), which efficiently represents highly correlated quantum states and whose parameters are easily optimized by means of the variational Monte Carlo (VMC) method. Recently, a self-adjoint and positive semi-definite parametrization of the density matrix, in terms of a neural network has been introduced Torlai and Melko (2018).

The steady state of an open quantum system can be characterized by a variational principle Cui et al. (2015); Jakob and Stenholm (2003); Mascarenhas et al. (2015); Weimer (2015), whereby the dissipative part of the real-time dynamics, under quite general conditions, drives the system towards a unique steady state, in analogy to the imaginary-time Schrödinger equation that leads to the ground state of Hamiltonian systems.

In this Letter, we present a VMC approach to simulate the non-equilibrium steady state (NESS) of open quantum systems governed by the quantum master equation in Lindblad form. The density matrix is parametrized using a neural network ansatz Torlai and Melko (2018) and parameters are varied using an extension of the stochastic reconfiguration method Sorella et al. (2007), which is shown to approximate the real-time dynamics of the system. We apply the present VMC to study the steady-state properties of the dissipative XYZ spin model Jin et al. (2016); Lee et al. (2013); Rota et al. (2018, 2017); Casteels et al. (2018), that displays a prototypical second-order dissipative phase transition. Thanks to the Monte-Carlo sampling of expectation values, this method holds promise for the efficient simulation of open quantum systems with a large number of degrees of freedom.

Dynamics of open quantum systems – The dynamics of the density matrix $\hat{\rho}$ of an open quantum system is governed by the quantum master equation which – in case of Markovian coupling to the environment – takes the Lindblad form

[TABLE]

where the curly brackets denote the anti-commutator. The unitary part of the dynamics is generated by the term depending on the Hamiltonian $\hat{H}$, while $\hat{F}_{i}$ are the jump operators associated to the dissipative processes induced by the environment. The equation is typically expressed in terms of the Liouvillian superoperator as $d\hat{\rho}/dt=\mathcal{L}(\hat{\rho})$, whose formal solution is $\hat{\rho}(t)=e^{\mathcal{L}t}\hat{\rho}(0)$ ($t>0$). The existence and uniqueness of a NESS – defined as $\hat{\rho}_{ss}=\lim_{t\to\infty}\hat{\rho}(t)$ – satisfying

[TABLE]

has been demonstrated under quite general assumptions Nigro (2018); Minganti et al. (2018), in particular for finite-size spin and boson lattices Nigro (2018).

The steady state can be computed as the long-time limit of the solution of the quantum master equation, or by directly solving the homogeoeous linear system (2) with an additional condition on the trace of the density matrix. In both cases, the size of the problem is quantified by the square of the Hilbert space dimension, thus becoming computationally prohibitive already for a modest number of degrees of freedom. A promising route to the numerical computation of the NESS is provided by the variational principle. In cases where a unique steady state exists Nigro (2018), the NESS corresponds to the eigen-matrix of the Liouvillian super-operator $\mathcal{L}$ with zero eigenvalue Minganti et al. (2018). As all other eigenvalues have strictly negative real part, the NESS can be formally derived as the matrix that maximizes the real part of the expectation value (computed in matrix space) of the Liouvillian.

Neural Network Density Matrix – We assume that the Hilbert space of the system is spanned by the computational basis $|\bm{\sigma}\rangle$, where $\bm{\sigma}=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$ labels the states of $N$ degrees of freedom that compose the system. Here and in what follows we will assume binary local degrees of freedom, with $\sigma_{i}=\{-1,1\}$, which applies to the broad class of interacting spin-$1/2$ or qubit models. The density matrix in this basis is formally expressed as $\rho(\bm{\sigma},\bm{\eta})=\langle\bm{\sigma}|\hat{\rho}|\bm{\eta}\rangle$ in terms of the density operator $\hat{\rho}$. We denote a specific variational ansatz for the density matrix as $\rho_{\chi}(\bm{\sigma},\bm{\eta})$, where $\chi=(\chi_{1},\,\chi_{2},\ldots,\,\chi_{N_{p}})$ is a set of variational parameters.

A neural network ansatz for a self-adjoint, positive semi-definite density matrix was recently introduced Torlai and Melko (2018) in the specific form of a Restricted Boltzmann Machine (RBM). In a variational approach, RBMs present the significant advantage that the sum over the hidden-spin configurations can be carried out analytically, and the logarithmic derivatives with respect to the variational parameters admit simple expressions Carleo and Troyer (2017). Here we briefly describe how this ansatz can be derived from simple considerations on the density matrix. A self-adjoint, positive semi-definite expression for the density matrix is

[TABLE]

The states $\psi_{j}(\bm{\sigma},\chi)$ are not necessarily mutually orthogonal and the sum extends over $J$ states, with $J\leq d$ and $d=2^{N}$ is the dimension of the Hilbert space under study.

We start by introducing a RBM ansatz for each state $\psi_{j}(\bm{\sigma},\chi)$ entering expression (3). A RBM is composed of two layers of binary valued nodes (see Fig. 1): a visible layer for encoding the physical state and a hidden layer. Each node is associated with a bias ($a$- and $b$-parameters) and nodes in the different layers are connected via a set of weighted edges ($X$-parameters). For a large number of hidden nodes, this structure is known to describe quantum correlations efficiently Carleo et al. (2018); Glasser et al. (2018).

In order to express the mixed structure in eq. (3) as a single RBM, we embed an intermediate set of $L$ hidden nodes that are used to express the probabilities $p_{j}(\chi)$ in RBM form as $p_{j}(\chi)=\exp(\sum_{l}c_{l}h_{l})$, with $h_{l}=\pm 1$ and $c_{l}\in\mathbb{R}$. To index the different states in the mixture accordingly, this new set of hidden nodes must also enter the RBM expression of the wave functions. When carrying out the sums over configurations of hidden nodes, the final expression for the RBM density matrix is sup

[TABLE]

The RBM is sketched in Fig. 1, and $\chi=\{a_{i},b_{m},X_{mi},c_{l},W_{li}\}$ is the final set of parameters, which are assumed as complex valued with the exception of $c_{l}$ that must take real values. The representational power of the RBM is determined by the number of hidden nodes Le Roux and Bengio (2008). Here we set the densities of hidden nodes through the parameters $\alpha=M/N$, $\beta=L/N$, which measure the representational power of the RBM ansatz independently of the size of the spin lattice. When separately accounting for the real and imaginary parts of complex-valued parameters, the total number of computational parameters in the RBM ansatz is $N_{p}=N[(\alpha+\beta)(2N+1)+\alpha+2]$. In what follows, we will always assume $\alpha=\beta$ for simplicity.

Optimization – It is convenient to rewrite eq. (2) in a vectorized form by reshaping $\hat{\rho}$ into a column vector $|\rho\rangle$. Following Jakob and Stenholm (2003), $\mathcal{L}$ takes matrix form and the steady state density matrix fulfills $\langle\rho|\mathcal{L}|\rho\rangle=0$. Therefore the expectation value over the variational density matrix $\langle\langle\mathcal{L}_{\chi}\rangle\rangle=\langle\rho_{\chi}|\mathcal{L}|\rho_{\chi}\rangle/\langle\rho_{\chi}|\rho_{\chi}\rangle$ is a function of the variational parameters $\chi$. The parameter values that best approximate $\langle\langle\mathcal{L}_{\chi}\rangle\rangle=0$ can be found by means of various optimization procedures Torlai and Melko (2018); Torlai et al. (2018); Carleo and Troyer (2017); Nomura et al. (2017). In this Letter, we choose to adopt the Stochastic Reconfiguration (SR) scheme by Sorella et al. Sorella et al. (2007) which we extend to open quantum systems. The parameters are initialized to a small random value and, at each iteration, they are updated as

[TABLE]

where the learning rate $\nu$ is small enough to guarantee convergence. It can be shown sup that the SR scheme induces, at each iteration, a variation in the parameters that best approximates the time evolution of the density matrix over a time step $\nu$. Here, we define the covariance matrix $S$, the vector of forces $F$, and the logarithmic derivatives $O$ as

[TABLE]

where $k,k^{\prime}=1,\,2,\ldots,\,N_{p}$. The notation $\langle\langle\cdot\rangle\rangle$ denotes the normalized expectation value taken over the variational density matrix $|\rho_{\chi}\rangle$, and the derivatives $O_{k}(\bm{\sigma},\bm{\eta})$ are taken as diagonal operators in these expectation values. We point out that, while the expression for $S$ in (6) results in the VMC iterations following the real time evolution, minimization can be achieved by using any positive-definite covariance matrix. In particular, setting $S$ as the identity results in the steepest descent procedure. Since $S$ can be non-invertible, we apply an explicit regularization scheme, as introduced in Carleo and Troyer (2017): $S_{kk^{\prime}}^{\mathrm{reg}}=S_{kk^{\prime}}+\lambda(n)\delta_{k,k^{\prime}}S_{kk^{\prime}}$, where $\lambda(n)=\max(\lambda_{0}b^{n},\lambda_{\min})$. For the present calculations, they were set to $\lambda_{0}=100$, $b=0.998$ and $\lambda_{\min}=10^{-2}$.

Sampling – The various expectation values in (6) must be evaluated at each iteration step. We evaluate these quantities stochastically over a Markov-chain of $N_{MH}$ configurations $(\bm{\sigma},\bm{\eta})$ sampling the square modulus of the density matrix $|\rho_{\chi}(\bm{\sigma},\bm{\eta})|^{2}$. For this we adopt the Metropolis-Hastings algorithm Metropolis et al. (1953). In the limit of $N_{MH}\rightarrow\infty$, the statistical error decays as $1/\sqrt{N_{MH}}$. Choosing an appropriate set of rules for the random walk is key to an efficient Monte Carlo sampling. Here we randomly choose each move among those allowed by the Liouvillian superoperator sup .

Observables – Once the optimal parameter values have been determined, the expectation value of any quantum mechanical observable $\hat{\mathcal{O}}$ over the steady state can be expressed as

[TABLE]

which can also be evaluated using the Metropolis-Hastings algorithm. For all the quantities considered here, the expectation values were additionally averaged over 100 sets of parameter values $\chi(n)$ chosen in the asymptotic region of the SR interation, in order to improve the statistical accuracy. The overall error in the sampled observables has, in addition to the contribution from the Metropolis-Hastings algorithms, a contribution from the SR scheme and a systematic contribution related to the representational power of the RBM ansatz, as measured by the $\alpha$ and $\beta$ parameters.

Computational cost – The number of floating point operations to evaluate Eq. (6) scales as $N_{p}^{3}$, if we assume that the number of Metropolis-Hastings steps $N_{MH}$ is set to roughly the number of parameters $N_{p}$, as in Ref. Carleo and Troyer (2017). The Metropolis-Hastings procedure also scales with the number of connected states $N_{c}$, i.e. with the average number of nonzero elements in a column of the Liouvillian matrix. Finally, the efficiency of the whole procedure thus scales as $O(N_{p}^{3}+N_{p}N_{c})$.

Results – To assess the effectiveness of the method, we study a spin-$1/2$ XYZ model on a two-dimensional lattice with periodic boundary condition. Each spin is subject to a dissipation process into the $|\sigma^{z}=-1\rangle$ state. This model has been already widely investigated and is known to display a dissipative phase transition between a paramegnetic and a ferromagnetic phase Jin et al. (2016); Lee et al. (2013); Rota et al. (2018, 2017); Casteels et al. (2018). The Hamiltonian and the quantum master equation read ($\hbar=1$)

[TABLE]

where $\hat{\sigma}_{j}^{x}$, $\hat{\sigma}_{j}^{y}$, $\hat{\sigma}_{j}^{z}$ are the Pauli matrices, $\hat{\sigma}^{\pm}_{j}=(\hat{\sigma}^{x}_{j}\pm i\hat{\sigma}^{y}_{j})/2$, $J_{\alpha}$ are the coupling constants between nearest neighbour spins and $\gamma$ is the dissipation rate. The excitations in the system – induced by the anisotropic spin coupling – compete with the isotropic dissipative process, and this competition is at the origin of the dissipative phase transition Jin et al. (2016); Lee et al. (2013); Rota et al. (2018, 2017); Casteels et al. (2018). The effectiveness of the neural network ansatz is demonstrated by studying the system observables across a phase boundary.

In addition to the expectation value $\langle\langle\mathcal{L}_{\chi}\rangle\rangle$, we study the local magnetization

[TABLE]

and the steady-state structure factor

[TABLE]

computed for the asymptotic steady state.

Fig. 2 shows the convergence of $S_{ss}^{xx}(\mathrm{\textbf{k}}=\textbf{0})$ and $S_{ss}^{xx}(\mathrm{\textbf{k}}=(2\pi/3,0))$ to the exact result for a $3\times 3$ lattice, as $\alpha=\beta$ are increased. The parameters $\chi$ are initialized randomly and updated at each VMC step according to the SR scheme sup . The parameters of the model are chosen to lie in the vicinity of the dissipative phase transition, i.e. $J_{x}/\gamma=0.9$, $J_{y}/\gamma=1.2$, $J_{z}/\gamma=1.0$. A clear convergence towards the exact value upon increasing $\alpha=\beta$ is found. The inset in Fig. 2 shows the SR evolution of $\mathrm{Re}(\langle\langle\mathcal{L}_{\chi}\rangle\rangle)$ over a typical VMC run. The oscillations at early times are a feature of the unitary part of the dynamics in the quantum master equation.

In Fig. 3 we display the magnetization as computed for different lattice sizes and as a function of the coupling parameter $J_{y}/\gamma$. For this choice of parameters, a para-to-ferromagnetic phase transition is expected to occur when increasing the coupling through the value $J_{y}\gtrapprox 1.04$ Rota et al. (2017, 2018), while a second phase boundary between a ferromagnetic and a paramagnetic region has been predicted by cluster mean-field calculations at around $J_{y}\gtrapprox 1.4$. For $2\times 2$ and $3\times 3$ lattices the VMC result agrees well with the exact calculation for a large enough number of variational parameters.

In Fig. 4 we display the spin structure factor $S_{ss}^{xx}(\mathrm{\textbf{k}}=\textbf{0})$ for the same parameters as in Fig. 3. The quantity $S_{ss}^{xx}(\mathrm{\textbf{k}}=\textbf{0})$ vanishes when in a paramagnetic phase, while it takes a finite value in the ferromagnetic region of the phase diagram. This behaviour is displayed both by the exact calculation for small lattices, and by the VMC data, in the vicinity of the phase boundary at $J_{y}\gtrapprox 1.04$. For values $J_{y}>1.4$ the system should become again paramagnetic in the thermodynamic limit of large lattices, but this feature was not displayed by the present data up to the largest lattice under study, in agreement also with recent stochastic Gutzwiller calculations Casteels et al. (2018).

   In all the calculations, special care was devoted to the choice of the SR time step $\nu$. The unitary part of the real-time dynamics generated by Eq. (1) makes the differential equation stiff, thus requiring to scale down $\nu$ appropriately as the system size – and thus the spectral width of the differential operator – is increased. A possible workaround would be to study an effective, purely dissipative dynamics using the super-operator $\mathcal{L}^{\dagger}\mathcal{L}$ as a generator. In the case of a unique steady state, this super-operator is self-adjoint and positive semi-defined, with the only null eigenvalue being associated to the steady-state solution. We argue that this effective dynamics would be more robust to the choice of the time step. The super-operator $\mathcal{L}^{\dagger}\mathcal{L}$ is however less sparse than $\mathcal{L}$ on the computational basis, calling for an efficient sampling scheme.

   Existing numerical approaches to the simulation of the steady state of a Markovian open quantum system either require the full representation of the Hilbert space into memory, or rely on a properly chosen truncation of the Hilbert space to a relevant subspace. The present VMC approach is free of these two limitations, thanks to the stochastic evaluation of expectation values by means of the Metropolis-Hastings algorithm. The neural network ansatz in terms of a RBM is highly representative of quantum correlated statistical mixtures, while being simple to handle numerically. In cases with very strong quantum correlations, this ansatz could be extended to deep network representations, as was recently done in the case of Hamiltonian problems Carleo et al. (2018); Cai and Liu (2018); Gao and Duan (2017); Choo et al. (2018). For some of these networks Choo et al. (2018), the hidden degrees of freedom can still be summed analytically, as for RBMs. Neural network representations are not restricted to spin degrees of freedom and have been successfully adopted to represent bosonic many-body states efficiently Saito (2017). For these reasons, the present VMC approach may emerge as the election tool to numerically model open quantum systems, with considerable impact on the study of fundamental physics and on the modeling of near-term, noisy quantum information platforms Preskill (2018).

I SUPPLEMENTAL MATERIAL

Neural Network Representation of the Density Matrix

In the main text we approximate the density matrix as a mapping $\rho_{\chi}$ with parameters $\chi$ which, given two input states $|\bm{\sigma}\rangle$ and $|\bm{\eta}\rangle$, returns the matrix element $\langle\bm{\sigma}|\hat{\rho}_{\chi}|\bm{\eta}\rangle=\rho_{\chi}(\bm{\sigma},\bm{\eta})$. We describe quantum systems with $N$ degrees of freedom on a computational basis $\bm{\sigma}=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$ assuming two-dimensional local Hilbert-spaces $\sigma_{i}=\{-1,1\}$. This choice generally applies to spin-$1/2$ and quantum bits. A general expression, ensuring the positive semi-definite property of the density matrix, is

[TABLE]

where $\chi$ are the parameters. The states $\psi_{j}(\bm{\sigma},\chi)$ are not necessarily mutually orthogonal and, together with the associated probabilities $1\geq p_{j}(\chi)\geq 0$, are an expression of the statistical mixture described by $\hat{\rho}_{\chi}$. Here, the sum extends over $J$ states, with $J\leq d$ and $d=2^{N}$ is the dimension of the Hilbert-space under study.

Following Ref. Carleo and Troyer (2017), we may represent any wave function $\psi(\bm{\sigma})$ as a Restricted Boltzmann Machine (RBM),

[TABLE]

where $q_{m}=\pm 1$ are a set of $M=\alpha\times N$ hidden spin variables and the leftmost sum runs over all possible hidden spin configurations $\{q\}$. The coupling and bias parameters $\{a_{i},b_{m},X_{mi}\}$ are, in general, complex-valued. In order to express the mixed state (12) in terms of a single RBM, we introduce a second set of $L\beta\times N$ hidden spins that we use to represent the probabilities $p_{j}(\chi)$ in RBM form as $p_{j}(\chi)=\exp(\sum_{l}c_{l}h_{l})$, with $h_{l}=\pm 1$ and $c_{l}\in\mathbb{R}$. This new set of hidden nodes should also enter in the expression of the wave functions, in order to index the different wave functions that make the mixed state. Hence, the full RBM form of the density matrix that we propose is

[TABLE]

Since no intra-layer connection is allowed, the hidden variables can be explicitly traced out and the neural network density matrix reads as

[TABLE]

The variational parameters in this RBM are $\chi=\{a_{i},b_{m},X_{mi},c_{l},W_{li}\}$ and, with the exception of the $c_{l}$’s, are all complex valued. When separately accounting for the real and imaginary parts of each parameter, the total number of real-valued computational parameters in the RBM ansatz is $N_{p}=N[(\alpha+\beta)(2N+1)+\alpha+2]$. We introduce the so called effective angles Carleo and Troyer (2017)

[TABLE]

finally obtaining

[TABLE]

In order to find the best set of parameters to describe the steady state, we choose to use the Stochastic Reconfiguration (SR) algorithm. This requires the calculation of the logarithmic derivatives of $\rho_{\chi}(\bm{\sigma},\bm{\eta})$ with respect to the real and imaginary parts of all variational parameters. Given expression (17), these are written in a rather compact form as

[TABLE]

Stochastic Reconfiguration for Open Quantum Systems

We approximate the density matrix with a neural network ansatz $\rho_{\chi}$, in order to find the best representation of the steady state. In the following, we use the vectorized formalism of the Liouville-von-Neumann master equation by reshaping $\hat{\rho}$ into a column vector $|\rho\rangle$ and the Liouvillian superoperator into a matrix $\mathcal{L}$ of dimension $d^{2}\times d^{2}$ where $d$ is the dimension of the Hilbert-space under study. In order to find the optimal set of parameters, we choose to extend the SR method Sorella et al. (2007) to open quantum systems.

In the SR optimization the variational parameters are changed at each iteration step by

[TABLE]

After a small perturbation in linear approximation the density matrix reads as

[TABLE]

where the diagonal operators $O_{k}$ are defined for any configuration $|\bm{\sigma},\bm{\eta}\rangle=|x\rangle$ as the variational derivative with respect to the $k$-th parameter

[TABLE]

Open quantum systems evolve under a one-parameters semigroup dictated by the Liouvillian superoperator as

[TABLE]

Under general assumptions, this dynamics asymptotically converges to the non-equilibrium steady state. We can enforce the real-time dynamics generated by the Liouvillian within the SR scheme. The method can be then interpreted as a real time evolution in the variational subspace.

To achieve this goal, we equate eq. (20) and (22) in the subspace spanned by the vectors $\{O_{k}|\rho_{\chi}\rangle\}$, thus obtaining

[TABLE]

where $\nu$ is small enough to guarantee convergence and we introduced the generalized forces and the covariance matrix for a given set of variational parameters at iteration step $n$

[TABLE]

Here, $\langle\langle\cdot\rangle\rangle$ denotes the normalized expectational value taken over the variational density matrix $|\rho_{\chi}\rangle$. The covariance matrix can be any positive-definite matrix. Choosing it as the identity will correspond to the steepest descent optimization method. Since $S$ can be non-invertible we apply an explicit regularization scheme, as introduced in Carleo and Troyer (2017)

[TABLE]

Here, the $\lambda(n)$ parameter is a function of the iteration step $n$ and decays with $\lambda(n)=\max(\lambda_{0}b^{n},\lambda_{\min})$. For the present calculations, they were set to $\lambda_{0}=100$, $b=0.998$ and $\lambda_{\min}=10^{-2}$.

Stochastic Sampling

For the SR optimization scheme, one must evaluate (24) at each iteration step. In order to perform a stochastic sampling of these quantities, we write the general forces and the covariance matrix into a suitable form as follows

[TABLE]

where $|x\rangle=|\bm{\sigma},\bm{\eta}\rangle$ are the basis elements spanning the space of the vectorized density matrix. Therefore, for a given $\chi$ we can generate a Markov-chain of $N_{MH}$ many-body configurations $|x_{1}\rangle\rightarrow|x_{2}\rangle\rightarrow\dots\rightarrow|x_{N_{MH}}\rangle$ sampling the square modulus of the density matrix elements $|\langle x|\rho_{\chi}\rangle|^{2}=|\rho_{\chi}(\bm{\sigma},\bm{\eta})|^{2}$ via a standard Metropolis-Hastings algorithm Metropolis et al. (1953). The configurations are generated based on a proposition scheme, and accepted according to the probability

[TABLE]

where $\mathcal{P}(z|y)$ is the conditional probability of proposing a state $z$ given $y$. Our move generation is based on the transitions dictated by the Liouvillian superoperator which occur by applying $\mathcal{L}|x\rangle$. Doing so decreases the required number of thermalization and sampling steps by promoting, in the random walk, the configurations with the highest importance. However, in case of a close-to-pure steady-state, a handful of density matrix elements have magnitude dominating the others. This causes the random walk to perform poorly, as some regions in configuration space become scarcely accessible. Hence, we introduced additional moves which occur with a low probability and reach a configuration by randomly flipping one spin in both $\bm{\sigma}$ and $\bm{\eta}$. The propositions always need to obey detailed balance.

For sampling $|\rho_{\chi}(\bm{\sigma},\bm{\eta})|^{2}$ in the case of a two-dimensional spin lattice model, we introduce the following possible moves on the configuration space spanned by $|\bm{\sigma},\bm{\eta}\rangle$

Column-hopping in $\bm{\sigma}$: a site $\sigma_{j}$ and its right neighbour are flipped. This move accounts for single spin hopping as well as double excitation or loss. 2. 2.

Row-hopping in $\bm{\sigma}$: a site $\sigma_{j}$ and its down neighbour are flipped. 3. 3.

Column-hopping in $\bm{\eta}$: a site $\eta_{j}$ and its right neighbour are flipped. 4. 4.

Row-hopping in $\bm{\eta}$: a site $\eta_{j}$ and its down neighbour are flipped. 5. 5.

Excitation in $\bm{\sigma}$: a site $\sigma_{j}$ is flipped. This move corresponds to the effect of having an external field in the direction perpendicular to the quantization axis. 6. 6.

Excitation in $\bm{\eta}$: a site $\eta_{j}$ is flipped. 7. 7.

Dissipator: sites $\sigma_{j}=\eta_{j}$ are flipped with an asymmetrical transition ratio. The dissipative moves are always proposed, while excitation are only proposed with ten-percent probability. 8. 8.

Jumper: two randomly chosen spins $\sigma_{l},\eta_{m}$ are flipped.

For the move proposition we choose one of the actions with uniform probability.

Detail on numerical approach

The VMC code was written in Python. The Metropolis Hastings algorithm was efficiently parallelized by splitting the Markov chain into several independent chains that were run in parallel using MVAPICH2. Both a CPU and a GPU version of the code were developed. GPU calculations allowed to significantly speed up the update of the $O_{k}$ expectation values, as well as the iterative solution of the linear system in the SR, using the MINRES-QLP algorithm Choi et al. (2011); Liu (2019) which can correctly handle the case of a singular matrix. The GPU version brought considerable advantage over CPU when the number of variational parameters was larger than 1000. In the Metropolis-Hastings sampling we set the number of accepted moves instead of the total number of moves. This is necessary as, for some choices of the physical parameters, the density matrix is highly concentrated on a few matrix elements and consequently MH moves are seldom accepted. The number of accepted moves was set to be roughly twice the number of variational parameters in all VMC runs.