Variational neural network ansatz for steady states in open quantum systems

TL;DR

This paper introduces a neural network-based variational method to efficiently find steady states of open quantum lattice systems, demonstrated on the dissipative quantum transverse Ising model.

Contribution

It develops a neural network ansatz for the steady state density matrix using purification and variational minimization, enabling new computational approaches for open quantum systems.

Findings

Successfully applied to the dissipative quantum transverse Ising model

Demonstrates the effectiveness of neural network ansatz in open quantum systems

Provides a scalable method for steady state computation

Abstract

We present a general variational approach to determine the steady state of open quantum lattice systems via a neural network approach. The steady-state density matrix of the lattice system is constructed via a purified neural network ansatz in an extended Hilbert space with ancillary degrees of freedom. The variational minimization of cost functions associated to the master equation can be performed using a Markov chain Monte Carlo sampling. As a first application and proof-of-principle, we apply the method to the dissipative quantum transverse Ising model.

Equations

$$\dot{\hat{\rho }}=\mathcal{L}\hat{\rho }=-i\text{commutator}\hat{H}\hat{\rho }+\sum _j\frac{{\gamma }_j}{2}\left[2{\hat{L}}_j\hat{\rho }{\hat{L}}_j^{†}-\{{\hat{L}}_j^{†}{\hat{L}}_j,\hat{\rho }\}\right],$$

$$\hat{\rho }(\mathbf{v})=\sum _{\mathbf{\sigma },{\mathbf{\sigma }}^{\prime }}{\rho }_{\mathbf{v}}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })|\mathbf{\sigma }⟩⟨{\mathbf{\sigma }}^{\prime }|.$$

$${\rho }_{\mathbf{v}}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })=\sum _{\mathbf{a}}{\psi }_{\mathbf{v}}(\mathbf{\sigma },\mathbf{a}){\psi }_{\mathbf{v}}^{\star }({\mathbf{\sigma }}^{\prime },\mathbf{a}).$$

$${\psi }_{\mathbf{v}}(\mathbf{\sigma },\mathbf{a})=\sqrt{{\mathcal{P}}_{{\mathbf{v}}_A}(\mathbf{\sigma },\mathbf{a})}\exp [-1/2\log ({\mathcal{P}}_{{\mathbf{v}}_{\theta }}(\mathbf{\sigma },\mathbf{a}))].$$

$$E_{\nu }(\mathbf{\sigma },\mathbf{a},\mathbf{h})=\mathbf{\sigma }\cdot {\mathbf{b}}_{\nu }^{(\sigma )}+\mathbf{a}\cdot {\mathbf{b}}_{\nu }^{(a)}+\mathbf{h}\cdot {\mathbf{b}}_{\nu }^{(h)}+{\mathbf{\sigma }}^T{\mathbf{W}}_{\nu }\mathbf{h}+{\mathbf{\sigma }}^T{\mathbf{U}}_{\nu }\mathbf{a}.$$

$${\rho }_{\mathbf{v}}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })=\exp [{\Gamma }_{\mathbf{v}}^{-}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })+{\Gamma }_{\mathbf{v}}^{+}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })+{\Pi }_{\mathbf{v}}(\mathbf{\sigma },{\mathbf{\sigma }}^{\prime })]$$

$$\mathcal{C}(\mathbf{v})=\frac{\text{norm}{d{\hat{\rho }}_{\mathbf{v}}/dt}_2^2}{\text{norm}{{\hat{\rho }}_{\mathbf{v}}}_2^2}=\frac{\text{Tr}[{\hat{\rho }}_{\mathbf{v}}^{†}{\mathcal{L}}^{†}\mathcal{L}{\hat{\rho }}_{\mathbf{v}}]}{\text{Tr}[{\hat{\rho }}_{\mathbf{v}}^{†}{\hat{\rho }}_{\mathbf{v}}]},$$

$$\mathcal{C}(\mathbf{v})=\sum _{\mathbf{\sigma },\tilde{\mathbf{\sigma }}}{p}_{\mathbf{v}}(\mathbf{\sigma },\tilde{\mathbf{\sigma }})|{\mathcal{C}}^{\text{Loc}}(\mathbf{v},\mathbf{\sigma },\tilde{\mathbf{\sigma }}){|}^2,$$

$${\mathcal{C}}^{\text{Loc}}(\mathbf{v},\mathbf{\sigma },\tilde{\mathbf{\sigma }})=\sum _{{\mathbf{\sigma }}^{\prime },{\tilde{\mathbf{\sigma }}}^{\prime }}\mathcal{L}(\mathbf{\sigma },\tilde{\mathbf{\sigma }};{\mathbf{\sigma }}^{\prime },{\tilde{\mathbf{\sigma }}}^{\prime })\frac{{\rho }_{\mathbf{v}}({\mathbf{\sigma }}^{\prime },{\tilde{\mathbf{\sigma }}}^{\prime })}{{\rho }_{\mathbf{v}}(\mathbf{\sigma },\tilde{\mathbf{\sigma }})}.$$

$$\nabla_{\bm{v}}\ \mathcal{C}({\bm{v}})=\sum_{\bm{\sigma},\bm{\tilde{\sigma}}}p_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})\mathcal{C}^{\text{Loc}}(\bm{v},\bm{\sigma},\bm{\tilde{\sigma}})^{\star}\Big{[}\sum_{\sigma^{\prime},\tilde{\sigma}^{\prime}}\mathcal{L}(\bm{\sigma},\bm{\tilde{\sigma}};\bm{\sigma}^{\prime},\bm{\tilde{\sigma}}^{\prime})\cr\frac{\rho_{\bm{v}}(\bm{\sigma}^{\prime},\bm{\tilde{\sigma}}^{\prime})}{\rho_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})}\mathcal{O}_{\bm{v}}(\bm{\sigma}^{\prime},\bm{\tilde{\sigma}}^{\prime})\Big{]}-\mathcal{C}(\bm{v})\mathcal{O}_{\bm{v}},$$

$$\text{expectationvalue}\ast \hat{\Theta }=\frac{\text{Tr}[\hat{\rho }\hat{\Theta }]}{\text{Tr}[\hat{\rho }]}=\sum _{\mathbf{\sigma }}{p}_{\mathbf{v}}^{\mathrm{obs}}(\mathbf{\sigma })\sum _{\tilde{\mathbf{\sigma }}}\frac{{\rho }_{\mathbf{v}}(\mathbf{\sigma },\tilde{\mathbf{\sigma }})\Theta (\tilde{\mathbf{\sigma }},\mathbf{\sigma })}{{\rho }_{\mathbf{v}}(\mathbf{\sigma },\mathbf{\sigma })},$$

$$H=\frac{V}{4}\sum _{⟨j,l⟩}{\hat{\sigma }}_j^z{\hat{\sigma }}_l^z+\frac{g}{2}\sum _j{\hat{\sigma }}_j^x,$$

Full text

In spite of the tremendous experimental progress in the isolation of quantum systems, a finite coupling to the environment Breuer and Petruccione (2007) is unavoidable and certainly plays a crucial role in the practical implementation of quantum information and quantum simulation protocols Devoret and Schoelkopf (2013). Moreover, through an active control of the environment via the so-called reservoir engineering, an open quantum manybody system can be prepared in non-trivial phases Carusotto and Ciuti (2013); Noh and Angelakis (2016); Hartmann (2016) with also possible quantum applications Verstraete et al. (2009); Barreiro et al. (2011). The theoretical description of open quantum manybody systems is in general out-of-the equilibrium and much less developed than for equilibrium systems. A mixed state with a finite entropy can be described by a density matrix, whose evolution is described by a master equation. Recently, a few theoretical methods have been developed to solve the master equation of open quantum manybody systems, including analytical approaches based on the Keldysh formalism Sieberer et al. (2016); Maghrebi and Gorshkov (2016), numerical algorithms based on matrix product operator and tensor-network techniques Mascarenhas et al. (2015); Cui et al. (2015); Jaschke et al. (2019); Werner et al. (2016); Kshetrimayum et al. (2017), cluster mean-field methods Biella et al. (2018); Jin et al. (2016), corner-space renormalization Finazzi et al. (2015); Rota et al. (2017, 2018), Gutzwiller mean-field Casteels et al. (2018), full configuration-interaction Monte Carlo Nagy and Savona (2018), permutation-invariant solvers Shammah et al. (2018) or efficient stochastic unravelings for disordered systems Vicentini et al. (2018). The research in the field is very active, since the different methods are optimal for different specific regimes. For example, the corner-space renormalization method is best suited for systems with moderate entropy, while matrix product operator techniques to systems with short-range quantum correlations.

In the last decade, the field of artificial neural networks has enjoyed a dramatic expansion and success thanks to remarkable applications in the recognition of complex patterns such as visual images or human speech (for a recent review see, e.g., LeCun et al. (2015)). The optimization (supervised learning) of the network is obtained by tuning the weights quantifying the connections between neural units via a variational minimization of a properly defined cost function. The wavefunction of a manybody system is in general a complex quantity, which is hard to be recognized. Recent works have proposed to exploit artificial neural networks to construct trial wavefunctions, where the connection weights in the network play the role of variational parameters Carleo and Troyer (2017); Glasser et al. (2018). Neural network approaches have already been succesffuly applied to a wide number (see e.g. van Nieuwenburg et al. (2017); Schindler et al. (2017); Choo et al. (2018); Czischek et al. (2018); Sharir et al. (2019)) of close Hamiltonian systems. However, they have not yet been generalized to the important quantum manybody problem of open systems.

In this Letter, we present a theoretical approach based on a variational neural network ansatz in order to determine the steady state of the master equation of open quantum lattice systems. We construct the ansatz for the mixed density matrix starting from a Restricted Boltzman Machine ansatz for a pure many-body wavefunction in an extended Hilbert space. We determine the optimal variational parameters by minimizing a cost function which involves the Liouvillian superoperator associated to the master equation for the density matrix. As a first application, we have considered the dissipative tranverse field quantum Ising model. We present a proof-of-principle demonstration by benchmarking the neural network calculations of the steady state against numerically exact simulations performed by quantum trajectories in the full Hilbert space Daley (2014). Our minimization of the cost function is performed by Markov chain Monte Carlo sampling of the gradient and is thus scalable to a large number of lattice sites. Perspectives of the present approach are discussed in the conclusions.

The general task that we wish to solve is the determination of the steady state of an open quantum system described by the Lindblad master equation Breuer and Petruccione (2007) for the system reduced density matrix $\hat{\rho}$, which reads (setting $\hbar=1$):

[TABLE]

where $\mathcal{L}$ is the so-called Liouvillian superoperator depending on the system Hamiltonian operator $\hat{H}$. The coupling to the environment is represented by interaction channels with the reservoir characterized by dissipation rates $\gamma_{j}$ and jump operators $\hat{L}_{j}$ acting on the system. Here we will focus on situations where the steady state ($\partial_{t}\hat{\rho}_{ss}=0$) is unique. In this case, the steady-state density matrix can be obtained as $\hat{\rho}_{ss}=\lim_{t\rightarrow\infty}\hat{\rho}(t)$ regardless of the initial condition. Although it is possible to engineer peculiar Liouvillians with more than one steady state Albert and Jiang (2014), typical physical systems with a finite Hilbert space dimension have a unique steady-state Spohn (1977); Minganti et al. (2018); Nigro_2019Nigro_2019.

For the many-body problem an analytical expression for $\hat{\rho}_{ss}$ can be found in very few cases Prosen (2014, 2011). In general, because of the exponential growth of the Hilbert space with the number of lattice sites, describing the full density matrix requires exponentially many complex numbers, which in practice can be done exactly only for a small number of sites. If one wants to attack the problem within a variational framework, the density matrix can be represented by an ansatz $\hat{\rho}_{\bm{v}}$ depending on a set of variational parameters ${\bm{v}}$. If $\{\ket{{\bm{\sigma}}}=\ket{\sigma_{1},\sigma_{2},...,\sigma_{N}}\}$ denotes a basis of states for the system Hilbert space, the density matrix can be expressed in the form

[TABLE]

In order to construct our neural network ansatz for the density matrix, we consider an extended Hilbert space $\mathcal{H}=\mathcal{H}_{S}\otimes\mathcal{H}_{A}$ where $\mathcal{H}_{S,A}$ represents respectively the system and ancillary Hilbert spaces. Such extended space is spanned by the basis set $\{\ket{{\bm{\sigma}},{\bm{a}}}\}$ where ${\bm{a}}=({a_{1},a_{2},...,a_{N_{a}}})$ labels the ancillary degrees of freedom. We start by considering a pure state in the extended Hilbert space, represented by the wavefunction $\psi_{\bm{v}}({\bm{\sigma}},{\bm{a}})$. In this framework the reduced density matrix of the system $S$ is obtained by tracing out the ancillary degrees freedom Torlai and Melko (2018), namely

[TABLE]

The next step is to represent $\psi_{\bm{v}}({\bm{\sigma}},{\bm{a}})$ via a neural network ansatz. This purified procedure automatically ensures that $\hat{\rho}_{\bm{v}}$ is Hermitian and positive semi-definite, as required for a density matrix. In a recent paper, Torlai and Melko Torlai and Melko (2018) proposed to describe purified wavefunctions as

[TABLE]

Both the amplitude ${\mathcal{P}}_{{\bm{v}}_{A}}({\bm{\sigma}},{\bm{a}})$ and phase-related function ${\mathcal{P}}_{{\bm{v}}_{\theta}}({\bm{\sigma}},{\bm{a}})$ of the purified wavefunction are given by the Boltzmann-like expression ${\mathcal{P}}_{\nu}({\bm{\sigma}},{\bm{a}})=\sum_{\bm{h}}\exp[-E_{\nu}({\bm{\sigma}},{\bm{a}},{\bm{h}})]$ (with $\nu\in\{\bm{v}_{A},\bm{v}_{\theta}\}$), where the associated dimensionless energy reads

[TABLE]

Note that the ansatz parameters are ${\bm{v}}=({\bm{v}}_{A},{\bm{v}}_{\theta})$ where ${\bm{v}}_{\nu}=({\bf{b}}^{(\sigma)}_{\nu},{\bm{b}}^{(a)}_{\nu},{{\bm{b}}}^{(h)}_{\nu},{{\bm{W}}}_{\nu},{{\bm{U}}}_{\nu})$. The rectangular matrix ${{\bm{W}}}_{\nu}$ weighs the connections between the system variables (visible layer) to the auxiliary variables (hidden layer), while the weight matrix ${{\bm{U}}}_{\nu}$ quantifies the connection between the system variables and the ancillary ones (ancillary layer). Such neural network ansatz is represented by a tri-partite Restricted Boltzmann Machine depicted in fig. 1. In other words, there are two independent artificial neural networks, one for the amplitude ($\nu=A$) and one for the phase ($\nu=\theta$). By substituting those formulas into Eq. 3 and carrying out the sum over the ancillary degrees of freedom one obtains a closed formula for the entries of the density matrix:

[TABLE]

where the expression of $\Gamma^{+/-}$ and $\Pi$ can be found in the Supplemental Material 111For more details, see Supplemental Material at the end of this PDF document. The representation power Roux and Bengio (2008); Younes (1996); Montufar and Ay (2011) of this ansatz can be systematically improved by increasing the density of the hidden ($\alpha=N_{h}/N$) and ancillary layer ($\beta=N_{a}/N$). It is worth pointing out that this scheme is not specific to this network topology, but relies only on the general fact that if two visible layers are connected by a shallow ancillary layer, the ancilla can be traced out analitically and an efficient neural-network description of the density matrix can be obtained.

Having defined a variational ansatz $\hat{\rho}(\bm{v})$, we now wish to define a variational principle to determine the optimal parameters. In particular, we have to recast the search for the steady state into a minimization problem for a real, positive cost function $\mathcal{C}({\bm{v}})$ which has a global minimum when the master equation ${\mathcal{L}}\hat{\rho}_{{\bm{v}}}=0$ is satisfied Weimer (2015). Moreover, in order to be able to deal with large Hilbert spaces, we need a quantity which can be sampled and computed efficiently. These requirements are met by the following cost function expressed in terms of the 2-norm of the time derivative of the density-matrix:

[TABLE]

as (i) $\mathcal{C}({\bm{v}}_{ss})=0\iff\hat{\rho}({\bm{v}}_{ss})=\hat{\rho}_{ss}$ and (ii) $\mathcal{C}({\bm{v}})\geq 0$.

It is useful to rewrite Eq. 7 as a sum over the whole space of bounded operators on the Hilbert space:

[TABLE]

where $p_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})=\absolutevalue{\rho_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})}^{2}/Z$ corresponds to a probability distribution as $Z=\sum_{\bm{\sigma},\bm{\tilde{\sigma}}}\absolutevalue{\rho_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})}^{2}$ 222Being able to rewrite the cost function (7) into this form allows us to sample it efficiently. This is the main advantage with respect to using the trace norm of $\mathcal{L}\hat{\rho}$ Weimer (2015) as a cost function, which requires computing the singular-value-decomposition of $\mathcal{L}\rho$ at each iteration.. The local contribution reads 333 This is not the only possible expression that allows to sample the cost function $\mathcal{C}(\bm{v})$. However, this choice of the local cost function $\mathcal{C}^{\rm Loc}$ is particularly convenient since it respects the zero-variance property, it is more numerically stable and cheaper to compute. For further details see Sec. 3 of the Supplemental Material. :

[TABLE]

In order to find the global minimum of the cost function (7) by means of gradient-based iterative schemes we need to compute its gradient

[TABLE]

where we have defined the log-derivatives of the density matrix $\mathcal{O}_{\bm{v}}=\sum_{\bm{\sigma},\bm{\tilde{\sigma}}}\mathcal{O}_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})$ and $\mathcal{O}_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})=\nabla_{\bm{v}}\log{\rho_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})}$, which can be efficiently computed for the considered neural network.

The computational complexity of evaluating $\nabla\mathcal{C}({\bm{v}})$ exactly grows exponentially with the size of the system. This cost can be considerably reduced if one only uses an estimate of $\nabla\mathcal{C}({\bm{v}})$ obtained by sampling the values $(\bm{\sigma},\bm{\tilde{\sigma}})$ according to the probability $p_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})$. Because the normalisation factor $Z$ is not fixed, we cannot sample the distribution directly and have to resort to a Markov Chain Monte Carlo Becca and Sorella (2017a) with Metropolis update rules 444We point out that a promising direction of research would be to devise particular trial wave functions where $Z$ is fixed or cheaper to compute, so that a direct sampling of the distribution $p_{\sigma,\tilde{\sigma}}$ without a Markov Chain would lead to easier convergence properties. Indeed, it has been recently shown Sharir et al. (2019) that a direct sampling is possible in some types of networks.. At every sampling step, we propose to update the configuration $(\bm{\sigma},\bm{\tilde{\sigma}})\rightarrow(\bm{\sigma}^{\prime},\bm{\tilde{\sigma}}^{\prime})$ by switching a random number of spins and accept the new configuration with probability $\min(\exp[p_{\bm{v}}(\bm{\sigma},\bm{\tilde{\sigma}})/p_{\bm{v}}(\bm{\sigma}^{\prime},\bm{\tilde{\sigma}}^{\prime})],1)$.

Finally, in order to find the global minimum of the cost function we employ a standard Stochastic Gradient Descent algorithm Bottou (2010). In order to improve the performance of the Stochastic Gradient Descent (i.e. to reduce the number of iterations needed to converge to the global minima of the cost function) we update the variational parameters according to to the metric of the space of density matrices exploiting the Stochastic Reconfiguration Approach Becca and Sorella (2017b). During the optimization procedure we sample the physical observables of interest through another Markov chain as

[TABLE]

where $p^{\rm obs}_{\bm{v}}(\bm{\sigma})=\rho_{\bm{v}}(\bm{\sigma},\bm{\sigma})/\Tr[\hat{\rho}]$.

In order to benchmark our neural-network approach for open quantum systems, we consider here the dissipative quantum transverse Ising model, whose Hamiltonian is

[TABLE]

being $\hat{\sigma}_{j}^{\alpha}$ the Pauli matrices ($\alpha\in\{x,y,z\}$) acting on the $j$-th site. The first term represents the nearest-neighbor spin-spin interaction depending only on the $z$-components, being $V$ the coupling strength. The second term accounts for a local and uniform magnetic field along the transverse direction $x$. We consider local dissipative spin-flip processes described by the site-dependent jump operator $\hat{L}_{j}^{(z)}=\hat{\sigma}_{j}^{-}=\frac{1}{2}(\hat{\sigma}^{x}_{j}-i\hat{\sigma}^{y}_{j})$, which fully determine the master equation in Eq. (1).

Numerical results for steady-state observables of the dissipative quantum transverse Ising model on a 1D periodic chain are reported in Fig. 2. In particular, we report the spatial components of the averaged magnetization as a function of the magnetic field $g$ (in units of the dissipation rate $\gamma$) for $V/\gamma=2$. For $N=16$ lattice sites the predictions of the neural-network variational method (circles) are compared to the results obtained with a brute-force exact integration of the master equation in the whole Hilbert space, showing a good agreement over all the parameter range. For $g\lesssim\gamma$ and $g\gtrsim 2.5\gamma$ a remarkable precision is reached for all the local observables with a low density of the hidden and ancillary layer $\alpha=\beta=1$ and $\mathcal{O}(10^{2})$ minimization steps. For $1\lesssim g/\gamma\lesssim 2.5$ an higher number of variational parameters is required. In particular, as shown in Fig. 3, a systematic improvement of the relative error \epsilon_{\rm rel}\big{[}\expectationvalue{\sigma_{x}}\big{]} with respect to the exact solution can be obtained by increasing $\beta$. Interestingly, for $1\lesssim g/\gamma\lesssim 2.5$, we note that the gradient-descent procedure requires more iterations. This region corresponds to the range of $g/\gamma$ where the smallest nonzero eigenvalue of $\mathcal{L^{\dagger}}\mathcal{L}$ decreases significantly Jin et al. (2018). In this range the steady-state density matrix also displays nontrivial correlations and non-thermal mixness properties Jin et al. (2018). Remarkably, the fidelity of the reconstructed local density matrix with respect to the exact one is alway larger then $0.998$ for all the values of $g/\gamma$ considered.

Finally, as an example of convergence, the top panel of Fig. 4 depicts a typical evolution of the cost function in the iterative minimization procedure for a fixed set of parameters ($g/\gamma=1$), showing a good convergence towards the global minimum. In the bottom panel of Fig. 4, the convergence of the $x$-component of the averaged magnetization is also reported.

In conclusion, we have presented a general variational approach for the steady-state density matrix of open quantum manybody systems based on an artificial neural network scheme. Our method is scalable since the cost function associated to the Liouvillian of the master equation can be calculated via Monte Carlo sampling. We have demonstrated a proof-of-principle of the theoretical scheme by a successful benchmarking to brute-force finite-size simulations in the full Hilbert space for arrays of spins described by the dissipative quantum transverse Ising model. We would like to point out that the present approach does not depend on the specific network topology. Indeed, the variational procedure presented in this Letter is general and can be applied to many other neural-networks or physically-inspired variational ansätze. There are many future developments at the horizon, including the study of dynamical properties, the use of deep neural networks and/or alternative cost functions, comparison with other existing techniques as well as the study of disordered systems without translational invariance. The neural network approach has the potential to pave the way to the theoretical study of a wide spectrum of open quantum manybody systems.