Liouville Space Neural Network Representation of Density Matrices
Neural network quantum states as ansatz wavefunctions have shown a lot of
promise for finding the ground state of spin models. Recently, work has been
focused on extending this idea to mixed states for simulating the dynamics of
open systems. Most approaches so far have used a purification ansatz where a
copy of the system Hilbert space is added which when traced out gives the
correct density matrix. Here, we instead present an extension of the Restricted
Boltzmann Machine which directly represents the density matrix in Liouville
space. This allows the compact representation of states which appear in
mean-field theory. We benchmark our approach on two different version of the
dissipative transverse field Ising model which show our ansatz is able to
compete with other state-of-the-art approaches.
Simon Kothe, Peter Kirton
# Liouville Space Neural Network Representation of Density Matrices

###### Abstract

Neural network quantum states as ansatz wavefunctions have shown a lot of promise for finding the ground state of spin models. Recently, work has been focused on extending this idea to mixed states for simulating the dynamics of open systems. Most approaches so far have used a purification ansatz where a copy of the system Hilbert space is added which when traced out gives the correct density matrix. Here, we instead present an extension of the Restricted Boltzmann Machine which directly represents the density matrix in Liouville space. This allows the compact representation of states which appear in mean-field theory. We benchmark our approach on two different version of the dissipative transverse field Ising model which show our ansatz is able to compete with other state-of-the-art approaches.

## I Introduction

Developing techniques to study the many-body physics of open quantum systems opens up the avenue to access behavior beyond that possible in equilibrium. This is relevant for understanding the dynamics of a variety of experimental platforms, but it also allows us to address fundamental questions about what kinds of physics is realizable in situations where incoherent driving and dissipation compete with coherent Hamiltonian dynamics. For example, in arrays of superconducting qubits coherent hopping of excitations can compete with onsite losses [1], or in semiconductor microcavities non-linearities in the Hamiltonian can compete with photon losses [2].

Performing accurate simulations of the models describing this kind of physics can be very challenging. Even in the simplest case, where the dynamics of the system is well captured by a Markovian master equation [3], the effective size of the space required grows exponentially with the number of degrees of freedom included. This makes it difficult to make concrete statements about what phases are stable at large systems sizes, in the thermodynamic limit. Thus, developing new techniques, able to capture all of the required physical processes is a key challenge to be overcome.

For closed systems a variety of techniques are available, allowing calculation of both ground states and dynamics in a variety of situations. These range from density functional theory, which relies on approximations of the electronic density [4], to Monte Carlo methods [5; 6; 7; 8; 9] which draw a small number of samples from the Hilbert space approximating the correct distribution, and tensor network methods [10; 11; 12], which rely on compression of the full many-body wavefunction.

Recently, there has been much progress in using neural networks to accurately represent wavefunctions of many-body systems [13; 14; 15]. These neural network quantum states (NNQS) rely on the ability of neural networks to learn and represent any sufficiently smooth function [16; 17; 18], allowing an efficient description without requiring an exponential number of parameters. The simplest architecture used for NNQS is the Restricted Boltzmann Machine (RBM). They consist of two fully connected layers, one visible and one hidden. This simple structure of the architecture translates to a function which can be easily implemented numerically, allowing for an efficient Markov sampling of the relevant probability distribution [19; 20] and efficient optimization. By increasing the hidden unit density the expressibility of the ansatz increases, allowing them to represent any state in the limit of infinite parameters. This has already shown to be competitive with other state-of-the-art approaches for finding ground [13] and excited [21] states of 1D and 2D spin systems. This ansatz has been applied to a wide range of problems in this area, from simulating topologically complex states [22] to investigating how the eigenvalue spectrum of the quantum Fisher matrix gives information about how networks learn groundstates [23]. RBMs are able to represent some states with a volume-law entropy scaling [24]. However, this is not a generic property [25] and there are certain volume-law states which still do not permit an efficient representation.

Architectures beyond the RBM have also been used to study related problems. These range from recurrent neural networks [26; 27] to convolutional neural networks [28; 29; 30] and deep autoregressive models, which forgo the need for Monte Carlo sampling altogether [31].

These approaches have also begun to be applied to the dynamics and steady-states of open quantum systems. Neural density machines (NDMs) [32; 33; 34; 35] use a pair of copies of a pure state NNQS coupled together by an extra layer of hidden neurons. The density matrix of interest is then obtained by tracing out the degrees of freedom associated with one of the copies. Such purification schemes have been also used in tensor network approaches [36; 37; 11]. Other developments use a POVM representation of the density matrix, achieving accurate results [38; 39; 40] as well as autoregressive Gram-Hadamard density operators which extend the NDM by adding additional layers which allows the capture of more complex correlations [41]. Another ansatz was introduced by Yoshioka and Hamazaki [42] which uses a binary encoding to vectorise their density matrix. So far, the literature has focused mostly on purifying the density matrix in Hilbert space, with the exception of Ref. [42]. Here, we propose a different approach, which writes the den
sity matrix directly in Liouville space. The Liouville density machine (LDM) ansatz which we propose does not require additional visible units and maintains a correspondence of the visible layer and the concrete physical system. We show that this is able to efficiently represent a larger range of physically relevant states than the NDM improving learning and accuracy, while retaining the simplicity of RBMs.

This paper is organized as follows: In Sec. II we give a brief description of the Lindblad master equation which describes Markovian open system dynamics. Sec. III details how finding the steady-state of this equation can be recast in terms minimization of a cost function and how it can be estimated using Monte Carlo sampling. Then, in section IV, we give a brief introduction to NNQS, using the original ansatz of Ref. [13], going on to show how this can be generalized to open systems. Sec. V gives a detailed comparison of the different approaches to solve this problem, using two versions of the dissipative transverse field Ising model as a benchmark. We gain insights into which kinds of states are easy and difficult to represent efficiently with the architectures described. Finally in Sec. VI we give our conclusions and discuss possible ways to build upon these results. The appendices contain a detailed derivation of the stochastic reconfiguration algorithm and an explanation of how we produce the Markov chains for estimating the required expectation values from the neural network.

## II Open system dynamics

The presence of an external environment fundamentally changes the nature of the dynamics of a quantum system. The state is no longer captured by a wavefunction and the non-unitary dynamics cannot be described by the Schrodinger equation. Instead, the natural description is in terms of a reduced density operator for the system, \(\rho\), and a master equation which defines its time evolution. This then allows for energy dissipation, decoherence and, most importantly for the present manuscript, the relaxation into a non-equilibrium steady-state in the long time limit.

If the coupling to the environment is weak and structureless, one can assume that system-induced correlations within the bath decay faster than the effects of the bath on the system. This leads to the particularly simple Lindblad master equation [3] which takes the form

\[\frac{d\rho}{dt}=\mathcal{L}\rho=-i[H,\rho]+\sum_{i}\gamma_{i}\mathcal{D}[A_{ i}]\rho. \tag{1}\]

Here, the sum runs over the dissipation channels \(i\) which are defined by the jump operators \(A_{i}\) and rates \(\gamma_{i}\). The dissipation superoperators are given by the Lindblad form

\[\mathcal{D}[A]\rho=A\rho A^{\dagger}-\frac{1}{2}\left(A^{\dagger}A\rho+\rho A ^{\dagger}A\right). \tag{2}\]

We only consider time-independent master equations.

The time evolution of the density operator is then governed by the formal expression

\[\rho(t)=\mathrm{e}^{\mathcal{L}t}\rho(0), \tag{3}\]

such that in the long time limit the stationary state satisfies

\[\lim_{t\rightarrow\infty}\mathcal{L}\rho(t)=0. \tag{4}\]

To simplify what follows both mathematically and numerically we recast Eq. (1) in Liouville space (sometimes also referred to as Choi's space). In this space the density matrix is reshaped into a vector,

\[\rho=\sum_{m}\sum_{n}\rho_{m,n}\ket{m}\bra{n}\rightarrow\ket{\rho}\rangle= \sum_{m,n}\rho_{m,n}\ket{m,n}\rangle. \tag{5}\]

Here the density matrix is expanded in the basis \(\{\ket{m}\}\) with expansion coefficients \(\rho_{m,n}\). The double-ket notation, \(\ket{\cdot}\rangle\), used above denotes a vectorized operator which lives in a Hilbert space consisting of two copies of the original, \(\ket{\rho}\rangle\in\mathcal{H}\otimes\mathcal{H}\). Then superoperators which act on these operator-kets are elements of Liouville space \(L\in(\mathcal{H}\otimes\mathcal{H})^{*}\otimes(\mathcal{H}\otimes\mathcal{H})\)[42; 43]. Here we often abbreviate the double index \((m,n)\) with a single index \(s\) which runs over the ket- and bra-indices of the density matrix \(\sum_{m,n}\rho_{m,n}\ket{m,n}\rangle\equiv\sum_{s}\rho(s)\ket{s}\rangle\). With this Eq. (1) now reads:

\[\frac{d}{dt}\ket{\rho}\rangle=L\ket{\rho}\rangle. \tag{6}\]

where

\[L= -i\left(H\otimes\mathbbm{1}-\mathbbm{1}\otimes H^{T}\right)\] \[+\sum_{k}\left[A_{k}\otimes A_{k}^{*}-\frac{1}{2}\left(\mathbbm{ 1}\otimes(A_{k}^{\dagger}A_{k})^{T}+(A_{k}^{\dagger}A_{k})\otimes\mathbbm{1} \right)\right]\]

is the matrix form of the superoperator \(\mathcal{L}\) that acts from the left on a density-ket, \(\ket{\rho}\rangle\).

We can then use standard linear algebra results to find the eigenvalues and eigenkets of the matrix \(L\), such that

\[L\ket{\rho_{i}}\rangle=\lambda_{i}\ket{\rho_{i}}\rangle. \tag{7}\]

A stationary state has \(\lambda_{i}=0\) while all other states have \(\mathrm{Re}\,\lambda_{i}<0\).

As we have seen above, the size of the required Liouville space scales much more quickly than the already exponentially growing Hilbert space. For example, for spin-1/2 lattice problems, the state vector of a closed
system grows as \(\mathcal{O}(2^{N})\), with the system size \(N\). The density matrix, on the other hand, grows as \(\mathcal{O}(4^{N})\) and hence the Liouvillian has up to \(\mathcal{O}(16^{N})\) elements which limits the size of accessible systems even further. In the next sections we will detail our approach of tackling this problem.

## III Finding the steady-state

Many routes to accessing larger system sizes with numerical simulations are variational methods [44; 45]. For example, tensor network methods [46] can be seen as variationally finding the steady-state by optimizing over the set of states captured by a given tensor network architecture, while corner space RG [47] optimizes a particular (small) set of basis states which can be used to build up to larger system sizes. As we will see later, neural network techniques make use of the same underlying mathematical structures. Crucial to these methods is the use an ansatz function which specifies the full state with a finite number of parameters, which we denote by \(\alpha\). We may then write \(\rho(s;t)\rightarrow\rho_{\alpha}(s;t)\) and optimize over \(\alpha\), to find the closest representation of the steady-state within this class. For this approach to be an efficient solution, the number of parameters must grow at most polynomially with the system size.

To find the steady-state, \(\lim_{t\rightarrow\infty}\rho(t)\equiv\rho_{ss}\), we require a cost function which measures how close the ansatz is to the true stationary density operator. To do this we make use of the fact that \(L\left|\rho_{ss}\right\rangle\rangle=0\) and that for all other states this quantity is finite. This gives us a cost function to optimize as well as a metric for how well the variational state parametrizes the true steady-state. The cost function we wish to optimize is given by the quantity

\[\left|\mathrm{C}(\alpha)\right|^{2}=\left|\frac{\langle\langle\rho_{\alpha}|L |\rho_{\alpha}\rangle\rangle}{\langle\langle\rho_{\alpha}|\rho_{\alpha}\rangle \rangle}\right|^{2}, \tag{8}\]

this uniquely vanishes in the steady-state. Note that the vectors appearing above are in Liouville space, hence the denominator involves the sum over all elements of the density matrix \(\langle\langle\rho_{\alpha}|\rho_{\alpha}\rangle\rangle=\sum_{s}|\rho_{\alpha} (s)|^{2}\). Expanding this cost function yields a form that can easily evaluated via Markov chain Monte Carlo

\[C_{\alpha} =\frac{\langle\langle\rho_{\alpha}|L|\rho_{\alpha}\rangle\rangle} {\langle\langle\rho_{\alpha}|\rho_{\alpha}\rangle\rangle}\] \[=\frac{1}{\sum_{s^{\prime}}|\rho_{\alpha}(s^{\prime})|^{2}}\sum_{ s,s^{\prime}}\rho_{\alpha}^{*}(s)\rho_{\alpha}(s^{\prime})\langle\langle s|L|s^{ \prime}\rangle\rangle\] \[=\frac{1}{\sum_{s^{\prime}}|\rho_{\alpha}(s^{\prime})|^{2}}\sum_ {s,s^{\prime}}\frac{|\rho_{\alpha}(s)|^{2}\rho_{\alpha}(s^{\prime})}{\rho_{ \alpha}(s)}\langle\langle s|L|s^{\prime}\rangle\rangle\] \[=\sum_{s}p(s;\alpha)C_{\mathrm{loc}}(s;\alpha),\]

where

\[p(s;\alpha)=\frac{|\rho_{\alpha}(s)|^{2}}{\sum_{s^{\prime}}|\rho_{\alpha}(s^{ \prime})|^{2}}, \tag{9}\]

defines a probability distribution over the entire density matrix which we can draw samples from. The local cost associated to each element of this distribution is then

\[C_{\mathrm{loc}}(s;\alpha)=\sum_{s^{\prime}}\langle\langle s|L|s^{\prime} \rangle\rangle\frac{\rho_{\alpha}(s^{\prime})}{\rho_{\alpha}(s)}. \tag{10}\]

All parts of this local cost can be efficiently calculated. The states \(s^{\prime}\) which connect, via the matrix element of the Liouvillian, to the original state \(s\) are generated during the Metropolis-Hastings step and can be accessed at any later stage of the process, see App. B for details on the algorithm used. Since, for a local Liouvillian, the number of non-zero elements grows only linearly in system size, the evaluation of the matrix elements of \(L\) can also be done efficiently. The samples we draw follow the distribution \(p(s;\alpha)\), and so the cost can be calculated as a simple mean over the local costs [20]

\[C_{\alpha}\approx\frac{1}{N_{s}}\sum_{s}C_{\mathrm{loc}}(s;\alpha), \tag{11}\]

where \(N_{s}\) is the number of samples and the sum runs over the Monte Carlo samples.

Along with an estimate for the cost function we also need a way of updating the parameters, \(\alpha\) such that the state we find is optimized. The simplest way to do this is via stochastic gradient ascent (SGA), where the gradients of \(C_{\alpha}\) are also estimated with the Monte Carlo samples and at each step in the simulation the parameters are updated as

\[\alpha\rightarrow\alpha^{\prime}=\alpha+\eta\nabla_{\alpha}C_{\alpha}, \tag{12}\]

with the learning rate \(\eta\). There are several problems with this approach, e.g. it has been shown [23] that SGA has severe problems with steep energy surfaces. The main problem for the present case is that \(L\) also has a left eigenstate with eigenvalue \(0\),

\[\langle\left\lceil\right\rceil L=0. \tag{13}\]

This is the trace-state which is defined as

\[\langle\left\lceil\right\rceil\rho\rangle\rangle=\mathrm{Tr}\left[\rho\right], \tag{14}\]

and since the dynamics of any physical master equation is necessarily trace preserving we find the result above. This means that there is another state which optimizes the cost function. A solution to both of these problems
is to instead use Stochastic Reconfiguration (SR) [20; 23] to update the parameters. SR can be derived by asking which parameter update \(\gamma\),

\[\alpha\rightarrow\alpha^{\prime}=\alpha+\eta\gamma, \tag{15}\]

best approximates a step in real time (see App. A for a derivation). Since the trace state cannot be found by a real-time evolution generated by \(L\) this guarantees that we will find the correct steady-state when optimizing the cost, Eq. (8). Furthermore, SR takes into account the curvature of the energy landscape, speeding up the optimization on flat areas and slowing down in the presence of strong curvature. The result is that the updates are calculated as

\[\gamma=S^{-1}f, \tag{16}\]

where \(S_{i,j}=\langle O_{i}^{*}O_{j}\rangle-\langle O_{i}^{*}\rangle\langle O_{j}\rangle\) is the quantum Fisher matrix and \(f_{i}=\langle O_{i}^{*}L\rangle-\langle O_{i}^{*}\rangle\langle L\rangle\) is the the gradient of the cost function in Eq. (8). The angle brackets, \(\langle\ \cdot\ \rangle\), denote the expectation value over the Monte Carlo samples and the operator \(O_{i}(s)\) is the logarithmic derivative of \(\rho_{\alpha}(s)\) with respect to the \(i\)th parameter

\[O_{i}(s)=\frac{1}{\rho_{\alpha}(s)}\frac{\partial}{\partial\alpha_{i}}\rho_{ \alpha}(s). \tag{17}\]

The matrix \(S\) is defined to be positive, however when estimating it via Monte Carlo sampling it can happen that some eigenvalues vanish and \(S\) becomes singular. One can work around this problem by either calculating the pseudo-inverse or adding a small regularization, \(\lambda\approx 10^{-3}\), to the diagonal. We employ the latter method here.

We next need to specify the form of our ansatz function for the density operator, \(\rho_{\alpha}\). To do this we first give a brief overview of how neural networks can be used to represent many-body wavefunctions.

## IV Neural network quantum states

Neural network quantum states are a variational ansatz whose parameters and functional form are defined by an underlying neural network architecture. This approach has been recently developed to describe the many-body wavefunction of interacting spin systems and was introduced in Ref. [13]. The key insight here is that the heart of the many-body problem is to find the relevant parts of Hilbert space in which e.g. the ground state of a system lives. This is essentially a problem of dimensionality reduction and feature extraction, which are the two strongest points of neural networks.

NNQS use this ability of neural networks to efficiently represent any continuous function [16; 17; 18], e.g. a wavefunction [13; 21; 22] or density matrix [33; 34; 35] by learning the most important features which define that function.

The simplest NNQS architecture is the RBM. These are bilayer neural networks with one visible and one hidden layer of neurons. The network is parametrized by a bias attached to each neuron and a set of weights which connects the two layers. This setup is shown schematically in Fig. 1. In the same spirit as the discussion in the previous section, it is possible to then write a variational ansatz for a wavefunction for a spin-1/2 lattice model. This takes the form

\[\Psi_{\alpha}(s) =\sum_{\{h_{i}\}}\mathrm{e}^{\sum_{j}a_{j}s_{j}+\sum_{i}b_{i}h_{i }+\sum_{ij}W_{ij}h_{i}s_{j}} \tag{18}\] \[=\mathrm{e}^{\sum_{j}a_{j}s_{j}}\prod_{i=1}^{M}2\cosh\left[b_{i} +\sum_{j}^{N}W_{ij}s_{j}\right]. \tag{19}\]

Where in the second line we use the fact that the sum over configurations can be analytically calculated. Here \(a_{i}(b_{j})\) denotes the visible (hidden) bias of the \(i\)th site(\(j\)th hidden unit), and \(W_{ij}\) the weights connecting the \(j\)th visible unit to the \(i\)th hidden unit. \(M(N)\) denotes the number of hidden(visible) units. The ratio between hidden and visible units is \(\beta=M/N\). For \(\beta=0\) the only states which can be described are those that would arise in a mean-field description, i.e. those which are a product of single site wavefunctions. Increasing the value of \(\beta\) allows a systematic increase in the amount of entanglement which can be described [13]. A further strength of the RBM wavefunction is the simplicity of its derivatives with respect to the parameters. This allows for efficient evaluation of all of the derivatives required to update the state.

In its simplest form the RBM ansatz described above is only able to discriminate between two different visible states, \(s=[-1,1]\), which allows at most for the representation of pure states of spin-1/2 systems. However, this is not sufficient to describe a vectorised density operator in Liouville space which, even for a spin-1/2 system has 4 elements. To overcome this issue we use an approach

Figure 1: Schematic representation of a RBM. It shows the visible units (green) connected to the hidden units (blue) via the weight matrix (double-headed arrows). The single-headed arrows labeled by \(a_{i}\) and \(b_{i}\) show the biases applied to the neurons.
based on that introduced in Ref. [19] for the study of pure states of spin-1 lattices. This was achieved by adding an additional set of biases and weights which allows discrimination between the \(s=[-1,0,1]\) states.

We can adapt this idea to provide a basis for a NNQS ansatz for density-kets, \(\ket{\rho}\)). Such an ansatz needs to be able to discriminate between the four states of the local density operator, while maintaining a one-to-one correspondence between the visible layer and the physical system. By adding yet one more set of visible biases and weights we are able to discriminate the four different local states which allows us to represent a density-ket. Furthermore, we do not require any additional visible nodes to achieve this, as in Ref. [42]. The total ansatz function now reads:

\[\rho_{\alpha}(s)=\mathrm{e}^{\sum_{j}^{N}a_{j}^{(1)}s_{j}+a_{j}^ {(2)}s_{j}^{2}+a_{j}^{(3)}s_{j}^{3}}\\ \times\prod_{i=1}^{M}2\cosh\left[b_{i}+\sum_{j}^{N}U_{i,j}s_{j}+V _{i,j}s_{j}^{2}+W_{i,j}s_{j}^{3}\right], \tag{20}\]

where \(N\) and \(M\) are the number of visible and hidden nodes respectively. We have introduced a set of biases, \(a^{(n)}\) to the visible units which allow us to distinguish different visible biases, \(b\) again gives the bias on each hidden neuron and \(U\), \(V\), \(W\) are the weight matrices that connect the two layers. Theses biases and connections are shown schematically in Fig. 2. Similar to the pure-state case of Eq. (18), this ansatz has labeling freedom and very simple derivatives which can be found analytically and efficiently implemented. To label the density matrix elements on each site we use the mapping

\[s=\begin{cases}&2\rightarrow\ket{\uparrow\uparrow}\rangle\\ &1\rightarrow\ket{\uparrow\downarrow}\rangle\\ -1\rightarrow\ket{\downarrow\uparrow}\rangle\\ -2\rightarrow\ket{\downarrow\downarrow}\rangle\end{cases} \tag{21}\]

where \(s=\pm 2\) denote the diagonal density matrix elements, while \(s=\pm 1\) are the off-diagonal elements. If the hidden unit density, \(M=0\), this ansatz is able to exactly capture those states described by mean-field theory, i.e. where the full density matrix can be factorized onto individual sites. This then allows the parameter \(M\) to systematically increase the amount the correlations which can be represented in a similar way to the bond dimension of an MPS.

This ansatz, which we refer to as a _Liouville Density Machine_ (LDM), bears some similarities to that proposed in Ref. [42]. In both cases the problem is transformed into Liouville space but in Ref. [42] this is achieved by adding an additional visible unit on each site. This allows them to represent the four required states. However, using two visible units to represent the state of a single spin leads to difficulties representing certain product states, limiting the expressibility of the ansatz. These approaches should be seen in contrast to those in Refs. [32; 33; 34; 35]. In these papers a purification ansatz is used, an extra layer of hidden units is added to represent the state of an auxiliary system which, when traced over, leaves the appropriate density matrix for the original system. Similar techniques have been applied to MPS simulations [36; 37]. In what follows we will refer to this purification ansatz as a _Neural Density Machine_ (NDM).

## V Results

To benchmark and understand the strengths and limitations of the LDM ansatz we study the stationary state of the 1D dissipative transverse-field Ising (TFI) model [40; 42; 48; 49; 50; 35]. The system consists of a chain of spin-1/2 particles. The Hamiltonian part of the evolution is governed by

\[H=\frac{J}{4}\sum_{i}^{N-1}\sigma_{i}^{z}\sigma_{i+1}^{z}+\frac{h}{2}\sum_{i}^ {N}\sigma_{i}^{x}. \tag{22}\]

Here, \(N\) is the number of sites, \(\sigma\) are the usual Pauli matrices, \(J\) is the interaction strength and \(h\) is the strength of the transverse field. The dissipation is governed by excitation loss on each site, so the jump operators which appear in the Liouvillian are \(A_{i}=\sqrt{\gamma}\sigma_{i}^{-}\). We choose the units of the interaction and field such that the dissipation strength is \(\gamma=1\). In all calculations that follow we set the interaction strength to \(J/\gamma=2\).

The steady-state of this model then has simple solutions in two limiting cases. When \(h\to 0\) the dissipation dominates the dynamics and the stationary state is a pure product-state with all the spins pointing down

\[\lim_{h\to 0}\rho_{ss}=\bigotimes^{N}\ket{\downarrow}\bra{\downarrow}. \tag{23}\]

In the opposite limit where \(h\rightarrow\infty\) there is only competition between the local onsite field and the dissipation

Figure 2: Schematic representation of an LDM. Similar to the RBM, it has visible (green) and hidden (blue) units. LDMs, however have three visible biases \(a_{i}^{(1,2,3)}\) as well as three weight matrices \(U\), \(V\), and \(W\) (double arrows). The weight matrices couple to different powers of the configuration vector \(s\) and encode correlations between different local states.
and so the steady-state ends up again as a product-state, but this time the state on each site is mixed

\[\lim_{h\rightarrow\infty}\rho_{ss}=\bigotimes\limits^{N}\frac{1}{2}\left(\ket{ \uparrow}\bra{\uparrow}+\ket{\downarrow}\bra{\downarrow}\right). \tag{24}\]

At intermediate values of \(h\) the steady-state interpolates between these two, building up complex long range classical and quantum correlations.

We will compare the results of the LDM ansatz to those obtained using the NDM approach as implemented in the NetKet library [51]. For small system sizes we will also be able to compare to results obtained from exact diagonalization, these are calculated using QuTIP [52].

As a first example we look at the case of a small system with \(N=6\). Our findings are summarized in Fig. 3. In each case we randomly initialize the parameter values and at each step we take Monte Carlo samples to approximate the cost function and the best updates to the parameters. At the end of each run we produce a new Markov chain, but this time sample from a probability distribution which follows the diagonal of the density matrix. This allows us to estimate the expectation value of observables of interest. Since there are fewer diagonal states than entries in the full density matrix we usually only need about 500-800 diagonal samples. In panels (a)-(b) of Fig. 3 we show the expectation value of \(\sigma_{i}^{x}\) on the central site and the \(\sigma_{i}^{z}\sigma_{i+1}^{z}\) correlation function on the central pair of sites as a function of the field strength, \(h\). This gives a good indication of how well the various approaches are able to produce single site observables. We see that, in general, there is good agreement between the exact results and those obtained using both the LDM and NDM approaches. In all cases the agreement is worst in the central region where \(1\leq h/\gamma\leq 2.5\), which is in agreement with the results of [35]. For the LDM a hidden unit density of \(\beta=1\) corresponds to 132 parameters for the NDM this is 174 parameters. We see that even with only 3/4 of the parameters the LDM generally gives as good or better results than the NDM. By increasing the value of \(\beta\) in the LDM we also increase the number of parameters so that, at \(\beta=2\), there are 246 parameters. We see that for the LDM ansatz increasing the value of \(\beta\) and hence the number of parameters used is able to significantly decrease the deviation from the exact result, thus we may use \(\beta\) as a way of checking for convergence when exact results are no longer possible. We can also see this effect more clearly in Fig. 3(c) and (d) where we show the Monte Carlo estimated cost function for each ansatz. The value of the cost function is significantly decreased at all values of \(h\) when \(\beta\) is increased. We also see that for the LDM ansatz the cost function is at a maximum in the regions where the convergence to the steady-state is worst, this allows us to use this estimate to again judge the accuracy of our results for system sizes where exact methods are unavailable. This is not true of the NDM approach where the cost function reaches a maximum at intermediate values of \(h\) and does not significantly decrease as \(h\) increases further. This is because the mixed product state, described in Eq. (24), is not so easy to represent in a purification ansatz; this mixed state requires a large amount of entanglement between the real and auxiliary spins. This is not the case for the LDM approach which can represent this state exactly with a hidden unit density \(\beta=0\).

A further comparison is shown in Fig. 4. This figure shows how the cost function of the LDM evolves for two

Figure 3: Comparison of the steady-state of Eq. (22) using both the NDM and LDM approaches. The system size, \(N=6\), is small enough that exact diagonalization is possible. (a) The expectation value of \(\sigma_{i}^{x}\) as a function of \(h\). (b) The expectation value of the \(\sigma_{i}^{z}\sigma_{i+1}^{z}\) correlation function. In both cases the expectation value is taken on the central site(s). The blue lines show the exact result, the NDM with \(\beta=1\) is green, the LDM with \(\beta=1(2)\) is orange (red). (c) The absolute value squared of the cost function in Eq. 8 for the two LDM results. (d) The expectation value \(\langle L^{\dagger}L\rangle\) used as a cost function for the NDM ansatz employed by NetKet. For both cases with \(\beta=1\) we used 4500 samples and optimized for 1000 steps with a learning rate of \(\eta=10^{-2}\) and a regularization of \(\lambda=10^{-2}\), for the expectation values we used 500 diagonal samples. In the \(\beta=2\) case we used 6000 samples and 2000 steps and 800 diagonal samples. The other meta parameters were the same in both cases.

Figure 4: Convergence of the LDM ansatz for two different hidden unit densities. Panel (a) shows the running estimate for the cost function while panel (b) is the variance in the same quantity. Increasing the hidden unit density improves the accuracy of the results. Both calculations were done at \(h/\gamma=1\) using the same parameters as Fig. 3.
different values of \(\beta\) over 6000 steps at one of the most difficult points, \(h=\gamma\). By increasing the number of variational parameters from 132 to 246 we were able to reduce the cost function by an order of magnitude. We also see that simply checking the value of the cost function does not give an accurate stopping condition for the algorithm. After around 1400 steps the cost function for \(\beta=2\) is very small but the variance is quite large. This means that the LDM has not found an eigenstate of the Liouvillian but is still giving a small value for the cost function. We propose that a condition based on a combination of both of these quantities can give a good way to automatically stop the learning process when a good approximation to the steady-state has been reached.

We now go on to examine how the accuracy of these approaches scales to larger system sizes. At \(N=16\) it is impossible to use exact methods to compare against, however this model is straightforward to solve with MPS simulations which we found to be fully converged for a bond-dimension of \(\chi=7\). Results of these calculations are shown in Fig. 5. The number of parameters for the NDM is 1104 while the LDM has 1126. For reference, a bond dimension of \(\chi=7\) corresponds to 1952 matrix elements in the MPS. We see very similar behavior to the \(N=6\) case, both approaches are more difficult to converge in the region of intermediate \(h/\gamma\) and the cost function for the LDM has a peak in this region. In Fig. 6 we show how the convergence can again be improved by increasing the hidden unit density \(\beta\). Here we choose \(h=2\gamma\) as this is the point where the convergence is worst. We see that as \(\beta\) is increased the cost function decreases towards zero and the expectation value moves towards that found in the MPS simulation. The expectation value here is a two-point correlation function which, in general, are harder to converge than single-site operators.

### TFI Model with rotated Hamiltonian

By making a simple change to the model discussed above it is possible to make the convergence of both neural network approaches considerably worse. To do this we change the Hamiltonian to a rotated basis [53, 48, 50]

\[H=\frac{J}{4}\sum_{i=0}^{4}\sigma_{i}^{x}\sigma_{i+1}^{x}+\frac{h}{2}\sum_{i= 0}^{5}\sigma_{i}^{z}, \tag{25}\]

but keep the dissipation processes the same as for the previous model. In this case the dissipation does not explicitly break the \(\mathbb{Z}_{2}\) symmetry of the model as the interaction term is perpendicular to the dissipation. Therefore the competition between the coherent and dissipative dynamics gives rise to complex correlations in the steady-state. This leads to a very rich mean-field phase diagram in high dimensions with possibilities for both first and second order phase transitions between different magnetic orderings [53, 49, 54]. In 1D these phase transitions turn into continuous crossovers, but complex correlations still build up when \(h\sim J\sim\gamma\). For a detailed review of the behavior of this model in 1D see Ref. [50].

We again study the convergence of both the LDM and NDM approaches for finding the steady-state of this model for a system size of \(N=6\). In panels (a) and (b) of Fig. 7 we show how both a single site and two-site observable varies with the applied field \(h\). We see that, even when using a large amount of samples and parameters the NDM ansatz is not able to find a good

Figure 6: Improvement of convergence as a function of \(\beta\) for the \(N=16\) TFI model at \(h=2\gamma\). In (a) we show the \(\sigma_{i}^{x}\sigma_{i+1}^{x}\) correlation function and in panel (b) the estimate for the cost function. All calculations ran for 7000 steps. To accommodate for higher parameter counts we increased the number of samples with \(\beta\) from 9000 at \(\beta=1\) to 17000 at \(\beta=2\). Other parameters are the same as in Fig. 5.

Figure 5: Steady-state of Eq. (22) as a function of field strength, similar to Fig. 3, for a larger system size \(N=16\). The expectation value of \(\sigma_{i}^{x}\) and \(\sigma_{i}^{z}\sigma_{i+1}^{z}\) on the central site(s) are shown in panels (a)–(b) and the relevant cost functions in panels (c)–(d). The LDM and NDM results are compared to those obtained from MPS simulations. We used \(\beta=1.4\) in the case of LDM and \(\beta=1\) for the NDM to ensure that both approaches use a similar number of parameters. The NDM has 1104 parameters, while the LDM has 1126. In both cases we evolved for 7000 steps with a learning rate of \(\eta=10^{-3}\) we took 9000 Monte Carlo samples at each step and a regularization of \(\lambda=3\times 10^{-3}\). We used 800 diagonal samples to estimate the expectation values.
approximation to the exact result, while the LDM is able to get much closer to the expected result, especially at small values of \(h\). We see that for both approaches the cost function estimate is much larger than it was for the simpler model described by Eq. (22). This is because the steady-state in this case has much more complex correlations than in the previous model, simple expressions like those in Eqs. (23)-(24) are not available, except at very large \(h\to\infty\) where the steady-state is the same as given in Eq. (23). We next go on to show how using measures of the entanglement found in the steady-state can give good intuition for when these kinds of difficulties arise.

### Entanglement Properties

To further understand the convergence of these approaches we investigate the entanglement present in the steady-states of both models over a range of parameters. Contrary to pure states, quantifying the amount of correlations that are present in a mixed state isn't as simple as just calculating the entanglement entropy between two halves of the system [55]. For our purposes we find that the negativity [56; 57]

\[\mathcal{N}=\frac{||\rho^{T_{A}}||-1}{2}, \tag{26}\]

provides a useful measure of the correlations which are difficult to represent using the LDM approach described above. Here, \(||\rho^{T_{A}}||\) denotes the trace norm of the partially transposed density matrix with the transpose taken over the degrees of freedom labeled by \(A\). This quantity gives a measure for the separability of a state. If two subsystems are entangled, the partial transpose can lead to negative eigenvalues, which leads to a trace norm greater than one and hence a non-zero negativity.

Panel (a) of Fig. 8 shows the negativity for both models considered in this paper for the three different possible bipartitions of a four-site system. In the case of the simpler model in Eq. (22) with the \(\sigma^{z}\sigma^{z}\) interaction we see a clear peak in the negativity at around \(h\sim 0.9\gamma\) for all partitions and a fast decay to zero at values of h above and below this point. This is because of the two limiting cases described in Eqs. (23) and (24) which both have zero negativity. The peak corresponds well to the range of \(h\)-values which were the the most difficult to find convergence with the neural networks.

In case of the more difficult model described in Eq. (25) with \(\sigma^{\pi}\sigma^{x}\) interactions we see a much higher negativity across the whole range of values of \(h\). This ties in well with our experience that this model is much harder to represent using both the LDM and NDM approaches across the board, with generally slight improvements around \(h=0\) and \(h\gtrsim 4\gamma\). We see that there is no region where the negativity reaches zero. This model does not have a simple product state steady-state anywhere in the observed parameter range.

## VI Conclusion

In summary we have proposed a NNQS ansatz which compactly represents density matrices in Liouville space, allowing us to find the steady-state of lattice models described by a Markovian master equation. This LDM approach was shown to be able to calculate the steady-state of a 1D open transverse field Ising model with 6 and 16 sites. The results were compared to the powerful NDM ansatz as implemented in NetKet. We found that our approach is always able to reach a comparable accuracy and in many cases is better able to find steady-state, especially when it contains a lot of correlations. We were

Figure 7: Optimizing the rotated TFI model as in Eq. 25 for \(N=6\). For both the LDM and NDM the optimization was run for 10000 steps, with hidden-unit-density of \(\beta=2\). At each step we used 15000 samples to estimate the cost function and required gradients. In panel (a) we give the steady-state expectation value of \(\sigma_{i}^{x}\) and in (b) we show the two-point correlation function \(\langle\sigma_{i}^{x}\sigma_{i+1}^{z}\rangle\). Panels (c) and (d) show the relevant cost functions for each ansatz. For both the NDM and LDM we used a learning rate of \(\eta=10^{-2}\), a regularization of \(\lambda=3\times 10^{-3}\) and 800 diagonal samples.

Figure 8: (a) The steady-state entanglement negativity, defined in Eq. (26) and (b) purity, \(\mathcal{P}=\mathrm{Tr}[\rho^{2}]\), for the \(\sigma^{x}\sigma^{x}\)-model (red) and the \(\sigma^{z}\sigma^{z}\)-model (blue). Comparing with the results of Figs. 3 and 7 we observe a correlation between a large negativity and poor accuracy of the neural network.
able to show that the accuracy of this approach is able to be systematically improved by increasing the number of hidden units and hence free parameters in the ansatz.

This permits a clear understanding of the class of states accessible to the LDM ansatz. As we show in Figs. 3 and 5, the most difficult regions to find convergence are very strongly correlated with parameters where the true steady-state has high negativity. This is in contrast with the NDM approach where there is difficulty representing mixed states with no correlations leading to a plateaus in the cost function which do not significantly decrease as the steady-state becomes more separable.

The results shown here are just a starting point for examining the usefulness of neural network approaches to finding the steady-state of open quantum systems. The RBM ansatz we used is the simplest possible architecture and extending the approach here to deep networks such as as deep RBMs [58], Recurrent Neural Networks [59], or transformers [60], which have already proven successful in closed system, provides a route to improving both the accuracy and numerical efficiency. The models studied here are also very simple and are accessible by other methods such as tensor network based techniques. However, the lack of an underlying lattice geometry for these neural networks can be exploited to study models with long range interactions and in higher dimensions which are much more difficult to simulate using other approaches.

###### Acknowledgements.

 We acknowledge useful discussions with Andrew Daley, Damian Hoffman, Daniela Pfannkuche and Filippo Vicentini. SK acknowledges financial support from EPSRC (EP/T517938/1).
