Constructing neural stationary states for open quantum many-body systems

TL;DR

This paper introduces a neural-network based variational method to efficiently compute stationary states of open quantum many-body systems described by Lindblad equations, using restricted Boltzmann machines.

Contribution

It presents a novel neural stationary state ansatz and applies variational Monte Carlo to find stationary states, extending neural quantum state techniques to open systems.

Findings

Successfully simulates 1D and 2D transverse-field Ising models

Efficiently computes stationary states of the XYZ model in 1D

Demonstrates the method's effectiveness for open quantum systems

Abstract

We propose a new variational scheme based on the neural-network quantum states to simulate the stationary states of open quantum many-body systems. Using the high expressive power of the variational ansatz described by the restricted Boltzmann machines, which we dub as the neural stationary state ansatz, we compute the stationary states of quantum dynamics obeying the Lindblad master equations. The mapping of the stationary-state search problem into finding a zero-energy ground state of an appropriate Hermitian operator allows us to apply the conventional variational Monte Carlo method for the optimization. Our method is shown to simulate various spin systems efficiently, i.e., the transverse-field Ising models in both one and two dimensions and the XYZ model in one dimension.

Equations

$$\frac{d\hat{\rho }(t)}{dt}=\mathcal{L}\hat{\rho }(t):=-i[\hat{H},\hat{\rho }(t)]+\sum _i{\gamma }_i\mathcal{D}[{\hat{\Gamma }}_i]\hat{\rho }(t).$$

$$\mathcal{D}[{\hat{\Gamma }}_i]\hat{\rho }(t)$$

$$\mathcal{L}{\hat{\rho }}_{\mathrm{SS}}=0,$$

$$\hat{\rho }=\sum _{\sigma \tau }{\rho }_{\sigma \tau }|\sigma ⟩⟨\tau |\mapsto |\rho ⟩⟩=\frac{1}{C}\sum _{\sigma \tau }{\rho }_{\sigma \tau }|\sigma ,\tau ⟩⟩,$$

$$\hat{A}\hat{\rho }\hat{B}$$

$$\mapsto |A\rho B⟩⟩$$

$$=$$

$$\frac{d|\rho (t)⟩⟩}{dt}$$

$$=\left(-i(\hat{H}\otimes \text{mathbbm}1-\text{mathbbm}1\otimes {\hat{H}}^T)+\sum _i{\gamma }_i\hat{\mathcal{D}}[{\hat{\Gamma }}_i]\right)|\rho (t)⟩⟩$$

$$\hat{\mathcal{D}}[{\hat{\Gamma }}_i]={\hat{\Gamma }}_i\otimes {\hat{\Gamma }}_i^{\ast }-\frac{1}{2}{\hat{\Gamma }}_i^{†}{\hat{\Gamma }}_i\otimes \text{mathbbm}1-\text{mathbbm}1\otimes \frac{1}{2}{\hat{\Gamma }}_i^T{\hat{\Gamma }}_i^{\ast }.$$

$$\hat{\mathcal{L}}|{\rho }_{\mathrm{SS}}⟩⟩=0,$$

$$\hat{\mathcal{L}}|{\rho }_n⟩⟩={\lambda }_n|{\rho }_n⟩⟩,$$

$$e^{\hat{\mathcal{L}}t}|{\rho }_n⟩⟩=e^{{\lambda }_{n}t}|{\rho }_n⟩⟩,$$

$${\hat{\mathcal{L}}}^{†}\hat{\mathcal{L}}|{\rho }_{\mathrm{SS}}⟩⟩=0.$$

$$⟨⟨\sigma ,\tau ⟩{\rho }_{\mathrm{RBM}}⟩=\frac{1}{Z}$$

$$\times \exp \left(\sum _i{a}_i{\sigma }_i+{\bar{a}}_i{\tau }_i+\sum _j{b}_j{h}_j\right),$$

$${\hat{\rho }}_{\mathrm{RBM}}^{\prime }=\frac{{\hat{\rho }}_{\mathrm{RBM}}+{\hat{\rho }}_{\mathrm{RBM}}^{†}}{2},$$

$$p(\sigma ,\tau )=\frac{|⟨⟨\sigma ,\tau ⟩{\rho }_{\mathrm{RBM}}⟩{|}^2}{⟨⟨{\rho }_{\mathrm{RBM}}⟩{\rho }_{\mathrm{RBM}}⟩}.$$

$$\hat{H}$$

$${\hat{\Gamma }}_i$$

$$F({\hat{\rho }}_1,{\hat{\rho }}_2)={\left(\mathrm{Tr}\sqrt{\sqrt{{\hat{\rho }}_1}{\hat{\rho }}_2\sqrt{{\hat{\rho }}_1}}\right)}^2.$$

$$\hat{H}$$

$${\hat{\Gamma }}_i$$

$$\hat{H}$$

$${\hat{\Gamma }}_i$$

$$\hat{\rho }=\frac{{\hat{X}}^2}{\mathrm{Tr}[{\hat{X}}^2]},$$

Full text

I Introducton

The dramatic development of machine learning techniques has inspired physicists to invent new numerical algorithms that further explore the frontier of condensed matter physics Carrasquilla and Melko (2017); Carleo and Troyer (2017). Successful applications include the phase classification using the well-established algorithms such as the deep learning Wang (2016); Carrasquilla and Melko (2017); Nieuwenburg et al. (2017); Hu et al. (2017); Broecker et al. (2017); Schindler et al. (2017); Liu and van Nieuwenburg (2018); Zhang et al. (2018); Yoshioka et al. (2018a), the acceleration of Monte Carlo simulations Torlai and Melko (2016); Liu et al. (2017); Huang and Wang (2017); Wang (2017); Shen et al. (2018); Yoshioka et al. (2018b), and the representation of the quantum many-body states using the high expressive power of the neural networks Carleo and Troyer (2017); Deng et al. (2017a, b); Gao and Duan (2017); Nomura et al. (2017); Saito (2017); Carleo et al. (2018); Saito and Kato (2018); Glasser et al. (2018); Choo et al. (2018); Kaubruegger et al. (2018); Levine et al. (2019); Sharir et al. (2019). In particular, the variational states based on the restricted Boltzmann machine (RBM) architecture have turned out to express the ground states of quantum many-body Hamiltonians composed of large number of spins efficiently, including one-dimensional (1D) and two-dimensional (2D) systems Carleo and Troyer (2017) and highly entangled systems Deng et al. (2017a).

Despite its rapid progress, however, machine learning has yet to be applied to one of the most challenging problems in modern condensed matter physics – open quantum many-body systems. Although the advancement of experiments Barreiro et al. (2011); Barontini et al. (2013); Labouvie et al. (2016); Tomita et al. (2017); Fitzpatrick et al. (2017) motivates an active field of research on open quantum many-body physics, it is notoriously difficult to solve the fundamental equation of motion, which is often well captured by the Lindblad master equation Lindblad (1976). Due to the growth of the number of parameters in proportion to the square of the Hilbert space dimension, description of the quantum states by density matrices requires additional computational resource compared to the closed system. Accordingly, the simulation of the Lindblad equation with the exact diagonalization method is hard even for small system sizes. It is thus important whether the machine learning techniques help us to simulate open quantum many-body physics. Particularly intriguing are nonequilibrium stationary states of dynamics, which can exhibit exotic structure such as entanglement Kraus et al. (2008); Kastoryano et al. (2011), nontrivial topology Diehl et al. (2011); Bardyn et al. (2013), and novel dissipative phases of matter Diehl et al. (2010); Tomadin et al. (2011); Ates et al. (2012); Gong et al. (2018); Gambetta et al. (2019).

In this work, we present a new scheme for simulating the stationary states of open quantum many-body systems by employing the ansatz which we refer to as the neural stationary state (NSS) in the following. As is schematically illustrated in Fig. 1, our NSS method is constituted by the following three steps:

• (a)

Vector representation: Make a copy of the Hilbert space and define a non-Hermitian operator $\mathcal{\hat{L}}$, which is the generator of Lindblad dynamics in the doubled Hilbert space.

• (b)

Definition of the cost function: Consider a Hermitian positive-semidefinite operator $\mathcal{\hat{L}}^{\dagger}\mathcal{\hat{L}}$, which becomes zero for the stationary state Cui et al. (2015).

• (c)

Optimization: Optimize the NSS ansatz using the variational Monte Carlo method (VMC).

We first demonstrate the expressive power of the ansatz by showing that the generic NSS exhibits volume-law entanglement entropy in the vector representation, which is the so-called operator space entanglement entropy Prosen and Pižorn (2007); Pižorn and Prosen (2009). Next, we show that our NSS ansatz is capable of representing the stationary states of the dissipative transverse-field Ising models in 1D and 2D, and XYZ model in 1D.

We remark that there have been many previous proposals for simulating open quantum many-body systems numerically. For example, the Lindblad dynamics is simulated by the density matrix renormalization group Cai and Barthel (2013); Werner et al. (2016); Gangat et al. (2017); Kshetrimayum et al. (2017) under the tensor network representation, which works very well especially in 1D as long as the operator space entanglement entropy of the density matrix is small. In addition, numerous works have focused particularly on the stationary states of the Lindblad dynamics. Cui et al. Cui et al. (2015) presented an elegant variational method to search for the stationary states of the Lindblad dynamics by minimizing the expectation value of $\mathcal{\hat{L}}^{\dagger}\mathcal{\hat{L}}$ using the matrix product operator (MPO) algorithm, which is powerful for 1D systems. Beyond 1D, Ref. Weimer (2015) treated variational quantum states that take low-order correlations around the product states into account. It is also notable that certain approximations beyond the mean-field theory, e.g., the cluster mean-field theory Jin et al. (2016, 2018), were employed. Few methods have been proposed, however, that can efficiently capture full quantum correlations beyond 1D.

The rest of the paper is organized as follows. A brief overview of open quantum systems in the Lindblad form and its vector representation is given in Sec. II. This representation allows us to map the stationary-state search problem into finding a zero-energy ground state of an appropriate Hermitian operator that is composed of the Lindblad operator and its Hermitian adjoint operator. In Sec. III, we introduce the NSS ansatz, which is optimized via the variational Monte Carlo technique. We show in Sec. IV that our ansatz is capable of expressing density matrices with volume-law operator space entanglement and also the stationary states of various spins systems, i.e., the transverse-field Ising models in both 1D and 2D, and the XYZ model in 1D. Finally, the summary of our work and the discussion on the future directions are presented in Sec. V. For completeness, we discuss the result for fitting random density matrices with the NSS in Appendix A, and the comparison of the computational time between the NSS and Lanczos methods is discussed in Appendix B.

II Open quantum systems in the Lindblad form

In this section, we first give a brief overview of the Lindblad form for describing open quantum systems. To simulate the stationary states efficiently, we introduce the vector representation of mixed states. Consequently, the stationary state of the Lindblad dynamics can be obtained by finding the ground state of an appropriate Hermitian matrix that is composed of the Lindblad operator and its Hermitian adjoint, or $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$. We propose that this problem can be solved efficiently via the conventional VMC method with insight into the optimization quality: the expectation value of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ is regarded as the cost function since the target state corresponds to the zero-energy eigenstate of it.

II.1 Lindblad master equation

Open quantum physics consider situations where a system interacts with their environments outside and follows non-unitary time evolution. Such systems with certain conditions, e.g. the Markovianity, are known to be well described by the Lindblad equation Lindblad (1976), which possesses the completely positive and trace-preserving property. Concretely, the time evolution of a mixed state $\hat{\rho}(t)$ is given by

[TABLE]

Here, $\mathcal{L}$ is the Lindblad superoperator, which is a linear map that takes a density matrix to another density matrix. The first term in the right hand side, given by the commutator $[\hat{H},\hat{\rho}]=\hat{H}\hat{\rho}-\hat{\rho}\hat{H}$, describes the unitary dynamics ruled by the Hamiltonian $\hat{H}$. The second term describes the non-unitary dynamics due to the dissipations. The contribution of the $i$-th term, whose strength is given as $\gamma_{i}$, is governed by a superoperator $\mathcal{D}[\hat{\Gamma}_{i}]$ acting on the density matrix $\hat{\rho}(t)$ as

[TABLE]

Here, $\hat{\Gamma}_{i}$, or the $i$-th jump operator, determines the detail of the dissipations.

II.2 Vector representation of the Lindblad equation

It is known that a time-independent Lindblad equation has at least one stationary state satisfying

[TABLE]

where $\hat{\rho}_{\rm SS}$ is a density matrix of the stationary state Rivas and Huelga (2011). To employ well-established numerical calculation schemes, we first map the density “matrix” $\hat{\rho}$ to an element in the so-called “operator space” as $\ket{\rho}\rangle\in\mathcal{H}\otimes\mathcal{H}$. The new representation of the state, which we call the “vector” representation throughout this manuscript, is explicitly given by

[TABLE]

where $\ket{\sigma,\tau}\rangle=\ket{\sigma}\otimes\ket{\tau}\in\mathcal{H}\otimes\mathcal{H}$ is a spin configuration basis that spans $\mathcal{H}\otimes\mathcal{H}$ and $C=\sqrt{\sum_{\sigma\tau}|\rho_{\sigma\tau}|^{2}}$ denotes the normalization factor. In the doubled Hilbert space, we discriminate the spins denoted by $\sigma$ and $\tau$ by referring to them as the physical and fictitious spins, respectively. We note that the normalizations in two representations are different from each other; the trace of the matrix is set to unity, i.e., $\sum_{\sigma}\rho_{\sigma\sigma}=1$, in the matrix representation, whereas the $L^{2}$-norm of $\ket{\rho}\rangle$, or $\langle\braket{\rho}\rangle$, is unity in the vector representation.

Using the mapping in Eq. (4), operators $\hat{A}$ and $\hat{B}$ that respectively act on $\hat{\rho}$ from left and right are mapped as follows,

[TABLE]

where $C^{\prime}=\sum_{\sigma\tau}\left|\sum_{\mu\nu}A_{\sigma\mu}\rho_{\mu\nu}B_{\nu\tau}\right|^{2}$ is a normalization factor. Applying the mapping in Eq. (5) to Eq. (1), we obtain the vector representation of the Lindblad equation as

[TABLE]

with

[TABLE]

Here, the Lindblad operator is denoted by the operator $\hat{\mathcal{L}}$ acting on the operator space $\mathcal{H}\otimes\mathcal{H}$. Hence, in this new representation the problem of finding the stationary state is expressed in terms of the standard linear algebra. Our goal is, concretely, to solve the equation for a non-Hermitian operator $\hat{\mathcal{L}}$ as follows:

[TABLE]

where $\ket{\rho_{\rm SS}}\rangle$ denotes the vector representation of the stationary state. Note that other right eigenvectors with non-zero eigenvalues satisfy

[TABLE]

where $\ket{\rho_{n}}\rangle$ is the with a right eigenvalue $\lambda_{n}$. For eigenmodes that are not stationary states, the real part of the corresponding right eigenvalues satisfy ${\Re}[\lambda_{n}]<0$ Breuer et al. (2002); Rivas and Huelga (2011), which implies that the modes eventually decay.

II.3 Stationary state as a “ground state” of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$

The Lindblad operator in Eq. (10) is a non-Hermitian matrix, whose eigenvalues are in general complex. In contrast, the product with the Hermitian-conjugated Lindblad operator, $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$, is a Hermitian matrix and has a real-valued non-negative spectrum. In this case, the lowest eigenstate(s) with eigenvalue(s) $\lambda=0$ of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ correspond to the stationary state(s). In other words, $\ket{\rho_{\rm SS}}\rangle$ satisfies

[TABLE]

This allows us to apply the well-established ground-state search technique in closed systems such as the variational approaches, in addition to the Lanczos method, if the first excited energy of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ does not vanish Cui et al. (2015). Therefore, the expectation value of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$ is suited for the cost function in the VMC method.

Note that the uniqueness of the stationary states is confirmed in various systems. For example, if the annihilation operator, or the incoherent spin flip along the $z$-axis in the language of spins, is included as the dissipation for each site, the quantum system has a unique stationary state regardless of the Hamiltonian Schirmer and Wang (2010). Unique stationary states also appear for other types of dissipations, as demonstrated in, e.g., Refs. Prosen (2012); Cai and Barthel (2013); Horstmann et al. (2013).

Let us emphasize that the variational approach has advantage in the sense of the cost function $\langle\braket{\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}}\rangle$, where $\langle\braket{\hat{O}}\rangle$ denotes the expectation value of the operator $\hat{O}$ in the vector representation. Since $\langle\braket{\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}}\rangle$ is exactly zero by construction for stationary states, this indicates the quality of the optimization. Note the difference from usual variational problems of finding ground states of Hamiltonians, for which quantification of the optimization quality is difficult without knowing the ground-state energy. In that case, the convergence of the cost function may both yield the desired state or indicate an ill result due to the local minima.

III Neural Stationary States

In this section, we present the method to compute our NSS for a given Lindbladian of the system. First, we describe the complex-valued RBM $\ket{\rho_{\rm RBM}}\rangle$, which introduces auxiliary binary degrees of freedom to extend expressive power of the ansatz, as follows:

[TABLE]

where $W_{ij}$ ($\bar{W}_{ij}$) denotes complex interaction amplitude between the $i$-th physical (fictitious) spin $\sigma_{i}$ ($\tau_{i}$) and $j$-th hidden spin $h_{j}$, $a_{i}$ ($\bar{a}_{i}$) is a complex magnetic field on the $i$-th physical (fictitious) spin, and $b_{j}$ is a complex magnetic field on the $j$-th hidden spin. The normalization factor $Z$ is determined such that $\langle\braket{\rho_{\rm RBM}}{\rho_{\rm RBM}}\rangle=1$. Denoting the number of the physical, fictitious, and hidden spins as $N,\bar{N}(=N),$ and $M$, respectively, we define the number ratio of the spins as $\alpha=M/(N+\bar{N})$ to compare the performance of the NSS ansatz under different system sizes.

Although the NSS obtained from Eq. (III) is not positive-semidefinite or Hermitian in general, sufficient optimization of the cost function is expected to ensure these two conditions in an approximated way Cui et al. (2015). In fact, we have confirmed that absolute values of unphysical negative eigenvalues, if any, and $||\hat{\rho}_{\rm RBM}-\hat{\rho}_{\rm RBM}^{\dagger}||/||\hat{\rho}_{\rm RBM}+\hat{\rho}_{\rm RBM}^{\dagger}||$ are in the order of $10^{-3}$. Both quantities are sufficiently small compared to unity, which indicates that the NSS method works well, and can be further reduced by, for instance, taking larger $\alpha$. In the following, physical observables such as the entropy are computed using the symmetrized density matrix,

[TABLE]

which assures the physical observables to be real-valued.

We update the parameters given in Eq. (III) so as to approximate the stationary state using the VMC sampling in the vector representation over the probability distribution

[TABLE]

In the following, the number of the sampled spin configurations at each step of optimization is denoted as $N_{s}$. The parameters in the NSS ansatz are updated using the stochastic reconfiguration method Sorella (2001), which is also known as the natural gradient method Amari et al. (1992); Amari (1998). This optimization, being equivalent to the sufficiently long imaginary-time evolution in the truncated Hilbert space spanned by variational ansatz Nomura et al. (2017), successfully avoids the local minima and converges to the desired state. Such an update step is repeated for $N_{\rm it}$ times until the cost function reaches the order of $10^{-3}$ or less.

IV Model and Result

In this section, we first demonstrate that our NSS based on the RBM is capable of simulating a state with large complexity in the sense of the operator space entanglement entropy, which is defined as the entanglement entropy of the mixed state in the vector representation. We then verify our NSS method by applying it to three models that are in principle experimentally realizable using cold atoms or trapped ions Lee et al. (2013): the transverse-field Ising models in 1D and 2D as well as the XYZ model in 1D.

IV.1 Random-valued NSS

The RBMs for pure states are known to be capable of expressing quantum states with large entanglement efficiently. Concretely, Ref. Deng et al. (2017a) has shown that the maximally entangled states can be expressed using only $O(L)$ hidden spins, where $L$ is the total number of spins in the system.

Similarly, we argue that the NSS ansatz given as Eq. (III) efficiently expresses density matrices with large operator space entanglement, namely the entanglement entropy of the density matrix in the vector representation Prosen and Pižorn (2007); Pižorn and Prosen (2009). Concrete definition of the operator space entanglement entropy throughout this paper is given as follows. Let a mixed state in the vector representation, $\ket{\rho\rangle}$, be a pure state on the doubled Hilbert space spanned by $L$ physical and $L$ fictitious spins. After choosing $[L/2]$ physical spins and corresponding $[L/2]$ fictitious spins to form a subsystem $\mathcal{S}$, we compute the entanglement entropy of $\mathrm{Tr}_{\mathcal{\bar{\mathcal{S}}}}[\ket{\rho\rangle}\bra{\langle\rho}]$, where $\bar{\mathcal{S}}$ is the complement of $\mathcal{S}$. Here, $[x]$ denotes the largest integer that does not exceed $x$.

To demonstrate our argument, we show that generic NSS ansatz with random parameters exhibits volume-law scaling of the operator space entanglement entropy. Shown in Fig. 2 is the sistem size dependence of the operator space entanglement entropy in random-valued NSS ansatz characterized by several different parameters (see the next paragraph). The quantum entanglement in the operator space seems to increase along the number of spins, which demonstrates the volume-law scaling. We thus argue that the large operator space entanglement entropy is not necessarily an obstacle for reliable simulations for our NSS, in contrast with methods based on the tensor network ansatz such as the MPO algorithm. As a caveat, we note that not all volume-law states can be expressed efficiently by NSS as is discussed in Appendix A.

The detailed calculation for random-valued RBM is done as follows, which is necessarily to justify the positive-semidefiniteness of the state Torlai and Melko (2018). A subset of hidden spins, which are labeled by $j$ with the number ratio of spins denoted as $\alpha_{1}$, are connected to both physical and fictitious spins. The interactions between the $i$-th physical (fictitious) spins, denoted as $W_{ij}(\bar{W}_{ij})$ are required to satisfy $W_{ij}=\bar{W}_{ij}^{*}$ and also the magnetic field to be $b_{j}=0$. The rest of the hidden spins are connected to either only physical or fictitious spins. Denoting the labels for such hidden spins as $k$ and $\bar{k}$, the parameters obey $W_{ik}=\bar{W}_{i\bar{k}}^{*}$ and $b_{k}={b}_{\bar{k}}^{*}$ while the other parameters are zero. The number ratio of spins for such hidden spins are given by $\alpha_{2}$ each, and hence the total is given as $\alpha=\alpha_{1}+2\alpha_{2}$. Under such conditions, both the real and imaginary parts of the parameters are drawn randomly from a section $[-r,r]$.

IV.2 Transverse-field Ising model in one dimension

We now discuss the validity of our NSS ansatz for the concrete open quantum many-body systems. We first consider the stationary state of a 1D transverse-field Ising model with the length $L$ under the periodic boundary condition. The Hamiltonian and the jump operators are given as

[TABLE]

where $\hat{\sigma}_{i}^{a}\ (a=x,y,z)$ is the Pauli matrix that acts on the $i$-th site, $V$ is strength of the nearest-neighbor interaction, $g$ is amplitude of the transverse field along the $x$-axis, and $\gamma$ gives the magnitude of the homogeneous dissipations. To take advantage of the periodic boundary condition, i.e., $\hat{\sigma}_{L}=\hat{\sigma}_{0}$, we impose translation symmetry on the NSS ansatz.

As was introduced in Sec. III, we optimize the expectation value $\langle\braket{\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}}\rangle$ using the stochastic reconfiguration method. Figure. 3 shows the comparison of stationary-state density matrices obtained by the Lanczos method, which efficiently approximates a subset of eigenvectors and eigenvalues of a sparse matrix Lehoucq et al. (1998), and NSS ansatz with the number ratio of the spins taken as $\alpha=1$. Here, the model parameters are taken as $V=0.3,g=1,$ and $\gamma=0.5$, which results in a stationary state with the volume-law entropy. Figure 3(a)(b) visually illustrates that the approximation of the state with the NSS well represents the stationary state calculated by the Lanczos method. The accuracy of the stationary state is also confirmed quantitatively via the calculation of the fidelity. The fidelity between $\hat{\rho}_{1}$ and $\hat{\rho}_{2}$, which is exclusively considered as the stationary-state density matrices obtained by the NSS optimization and Lanczos method in practice, is defined as Jozsa (1994)

[TABLE]

This corresponds to the largest fidelity between any two purifications of the density matrices. For the current case, we find the fidelity to satisfy $F>0.999$. We also observe in Fig. 3(c) that the expectation value $\langle\braket{\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}}\rangle$, which gives measure of the approximation 111We can show that $|\braket{\braket{{\rho}_{\mathrm{SS}}}{{\rho}_{\mathrm{RBM}}}}{|}^{2}\geq 1-\frac{\langle\braket{\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}}\rangle}{\Delta}$, where $\Delta$ is the spectral gap of $\hat{\mathcal{L}}^{\dagger}\hat{\mathcal{L}}$, is nicely optimized and reaches the order of $10^{-3}$. Accordingly, the physical quantities are in good agreement with the exact results. For example, the entropy contribution for each eigenvalue of the density matrix, i.e., $-p_{n}\ln p_{n}$ for the $n$-th eigenvalue $p_{n}$, is remarkably accurate (see Fig. 3(d)), such that the relative error of the total entropy is the order of $10^{-3}$.

As is the case with other VMC calculations, it must be noted that both numerical cost and required memory for optimizing the NSS ansatz is much suppressed compared to methods that deal with the whole Hilbert space. In particular, the wall time for the NSS and the Lanczos methods are compared in Appendix B.

IV.3 Transverse-field Ising model in two dimension

We next optimize the NSS ansatz for the 2D transverse-field Ising model on the square lattice with system size $L_{x}$ and $L_{y}$ along the $x$- and $y$-axes, respectively. We again take the periodic boundary condition. The Hamiltonian and the jump operators are given as Jin et al. (2018)

[TABLE]

where the summation in the first term of $\hat{H}$ is taken over the edges connecting the neighboring sites, which are denoted as $i$ and $i^{\prime}$.

The cost function in Figure 4(a) shows that our optimization works well even for the 2D case. Indeed, as shown in Fig. 4(b), the NSS simulates the entropy contribution for each eigenvalue of the stationary state with high accuracy. This result strengthens the expectation that our NSS ansatz does not suffer from high dimensionality, which can cause problems for the MPO ansatz.

IV.4 XYZ model in one dimension

Finally, we investigate the 1D XYZ model, in which the dissipations are known to invoke dramatic change of the phase diagram compared with the closed system Jin et al. (2016). The model is defined as

[TABLE]

where $J_{a}$ denotes the interaction for $a\ (a=x,y,z)$ component of the spin. We particularly take $J_{x}=0.9,J_{y}=0.4,J_{z}=1,$ and $\gamma=1$ with the periodic boundary condition, at which the finite system shows remnants of the phase transition predicted by the mean-field approximation Weimer (2015).

Shown in Fig. 5 is the comparison of the translationally symmetric NSS ansatz and the Lanczos method regarding the entropy contribution for each eigenvalue. Even though our choice of parameters leads to the non-simple stationary state of our small systems (as indicated from the peak of the structure factor Jin et al. (2016)), the NSS describes the exact results well.

V Conclusions and outlooks

We have proposed that the neural quantum states are suited for expressing the stationary states of open quantum many-body systems. By mapping the original stationary-state search problem of the Lindblad equation to the zero-energy ground state search problem of an appropriate Hermitian operator, we solve it with a variational ansatz based on the RBM, or the neural stationary state (NSS).

We have confirmed that the NSS can even We have then demonstrated that our NSS ansatz is capable of expressing the stationary states of the dissipative one- and two-dimensional transverse-field Ising models and one-dimensional XYZ model.

While the aim of our work is to make a first attempt to show the adequacy of the NSS for open quantum many-body systems including highly entangled states and two-dimensional states, we leave several intriguing questions as future works. One naive question is whether our ansatz can simulate larger system sizes, in which other methods suffer from expensive numerical cost. Another important question is to clarify the versatility of open quantum many-body systems addressable by our method. We expect that our ansatz performs well regardless of the dimensionality, as suggested in our calculations and the bipartite-graph structure of the RBM, which is free from the geometry of the underlying physical lattice. It is also interesting whether our method works for various long-range interacting systems (such as the Haldane-Shastry model Deng et al. (2017a)) with dissipations, whose mixed stationary states can be highly entangled.

Note added.- After completion of our work, we became aware of some related works. Refs. Hartmann and Carleo (2019); Nagy and Savona (2019) discussed the time evolution and stationary states of open quantum many-body systems by using the complex RBM and Ref. Vicentini et al. (2019) studied the approximation of the stationary states by the RBM ansatz.

Acknowledgements.

We are grateful to Masahito Ueda for reading our manuscripts with valuable discussions and helpful advice. We also thank Hosho Katsura, Zongping Gong, and Shunsuke Furukawa for fruitful comments. The numerical calculations were carried out with the help of NetKet Carleo et al. , SciPy Jones et al. (2001), QuTiP Johansson et al. (2013), and ARPACK Lehoucq et al. (1998). This work was supported by KAKENHI Grant No. JP18H01145 and a Grant-in-Aid for Scientific Research on Innovative Areas “Topological Materials Science” (KAKENHI Grant No. JP15H05855) from the Japan Society for the Promotion of Science. N. Y. and R. H. were supported by Advanced Leading Graduate Course for Photon Science (ALPS) of Japan Society for the Promotion of Science (JSPS). N. Y. was supported by JSPS KAKENHI Grant-in-Aid for JSPS fellows Grant No. JP17J00743. R. H. was supported by JSPS KAKENHI Grant-in-Aid for JSPS fellows Grant No. JP17J03189).

Appendix A Approximating Random Density Matrices by NSS

In this appendix, we randomly generate a density matrix and fit it by the NSS to see that the expressive power of the ansatz does not assure efficient representation of all volume-law states. Here, random density matrices are generated as

[TABLE]

where $\hat{X}$ is sampled from the Gaussian unitary ensembles of random Hermitian matrices. We have numerically checked that the operator space entanglement entropy defined as in the main text exhibits a volume-law scaling, i.e., operator space entanglement entropy $\propto L$ for matrices with size $2^{L}\times 2^{L}$ (data not shown).

Figure 6 shows that while a random density matrix generated following Eq. (24) can be approximated better by the NSS with larger $\alpha$, or the number ratio of the spins, the number of parameters and accordingly the numerical cost required to reach some fixed fidelity increase rapidly.

Appendix B Comparison of Computational Cost

In Fig. 7, we show the scaling of the computational time for calculating the stationary states by optimization of the NSS and the Lanczos method. Our variational method exhibits only polynomial scaling which is, clearly, far more efficient than the exponential scaling in the Lanczos method.