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.


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 manybody 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 lay-ers, 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].
Architectures beyond the RBM have also been used to study related problems.These range from recurrent neural networks [25,26] to convolutional neural networks [27][28][29] and deep autoregressive models, which forgo the need for Monte Carlo sampling altogether [30].
These approaches have also begun to be applied to the dynamics and steady-states of open quantum systems.Neural density machines (NDMs) [31][32][33][34] 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 [11,35,36].Other developments use a POVM representation of the density matrix, achieving accurate results [37][38][39] as well as autoregressive Gram-Hadamard density operators which extend the NDM by adding additional layers which allows the capture of more complex correlations [40].Another ansatz was introduced by Yoshioka and Hamazaki [41] 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. [41].Here, we propose a different approach, which writes the density 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 Schrödinger equation.Instead, the natural description is in terms of a reduced density operator for the system, ρ, 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 steadystate 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 Here, the sum runs over the dissipation channels i which are defined by the jump operators A i and rates γ i .The dissipation superoperators are given by the Lindblad form We only consider time-independent master equations.The time evolution of the density operator is then governed by the formal expression such that in the long time limit the stationary state satisfies 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, Here the density matrix is expanded in the basis {|m⟩} with expansion coefficients ρ m,n .The double-ket notation, |•⟩⟩, used above denotes a vectorized operator which lives in a Hilbert space consisting of two copies of the original, |ρ⟩⟩ ∈ H ⊗ H. Then superoperators which act on these operator-kets are elements of Liouville space L ∈ (H ⊗ H) * ⊗ (H ⊗ H) [41,42].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 m,n ρ m,n |m, n⟩⟩ ≡ s ρ(s) |s⟩⟩.With this Eq.( 1) now reads: where is the matrix form of the superoperator L that acts from the left on a density-ket, |ρ⟩⟩.
We can then use standard linear algebra results to find the eigenvalues and eigenkets of the matrix L, such that A stationary state has λ i = 0 while all other states have Re λ 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 O(2 N ), with the system size N .The density matrix, on the other hand, grows as O(4 N ) and hence the Liouvillian has up to 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 [43,44].For example, tensor network methods [45] 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 [46] 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 α.We may then write ρ(s; t) → ρ α (s; t) and optimize over α, 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→∞ ρ(t) ≡ ρ 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 |ρ ss ⟩⟩ = 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 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 ⟨⟨ρ α |ρ α ⟩⟩ = s |ρ α (s)| 2 .Expanding this cost function yields a form that can easily evaluated via Markov chain Monte Carlo 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 All parts of this local cost can be efficiently calculated.
The states s ′ 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; α), and so the cost can be calculated as a simple mean over the local costs [20] 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, α 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 α are also estimated with the Monte Carlo samples and at each step in the simulation the parameters are updated as with the learning rate η.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, This is the trace-state which is defined as 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,47,48] to update the parameters.SR can be derived by asking which parameter update γ, 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 where is the the gradient of the cost function in Eq. ( 8).The angle brackets, ⟨ • ⟩, denote the expectation value over the Monte Carlo samples and the operator O i (s) is the logarithmic derivative of ρ α (s) with respect to the ith parameter 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, λ ≈ 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, ρ α .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 manybody 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.
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 FIG. 1. Schematic representation of a RBM.It shows the visible units (lower row, green) connected to the hidden units (upper row, blue) via the weight matrix (double-headed arrows).The single-headed arrows labeled by ai and bi show the biases applied to the neurons.
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 ith site(jth hidden unit), and W ij the weights connecting the jth visible unit to the ith hidden unit.M (N ) denotes the number of hidden(visible) units.The ratio between hidden and visible units is β = M/N .For β = 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 β 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 betwen 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 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, |ρ⟩⟩.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.To label the density matrix elements on each site we choose to use the mapping where s = ±2 denote the diagonal density matrix elements, while s = ±1 are the off-diagonal elements.These values are chosen as they give a simple mapping from integers to density matrix elements and allow our neural network ansatz to easily distinguish between all of the states.To discriminate these four different local states we add one more set of visible biases and weights, allowing us to represent a density-ket.Furthermore, we do not require any additional visible nodes to achieve this, as in Ref. [41].The total ansatz function now reads: (1) where N and M are the number of visible and hidden nodes respectively.We have introduced a set of complex 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 complex 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.Calculating the derivatives of a complex valued and complex parametrized function is done using Wirtinger calculus [49].
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. [41].In both cases the problem is transformed into Liouville space but in Ref. [41] this is achieved by adding an additional visible unit on each site, allowing them to represent the four required states.A detailed comparison between the expressibility of these approaches is left as an open question.The LDM ansatz should be seen in contrast to those in Refs.[31][32][33][34].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 Schematic representation of an LDM.Similar to the RBM, it has visible (lower row, green) and hidden (upper row, blue) units.LDMs, however have three visible biases a {1,2,3} i 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.
been applied to MPS simulations [35,36].In what follows we will refer to this purification ansatz as a Neural Density Machine (NDM).
Compared to the purification approaches of Refs.[31][32][33][34], the LDM is not guaranteed by construction to be either Hermitian or positive-definite.However, we find that this is not a significant issue in our simulations because of our choice of SR as the descent algorithm.Note that any matrix can be represented in the eigenbasis of the Liouvillian [42] as where c i are expansion coefficients and |ρ i ⟩⟩ are the vectorized eigenmatrices of L. Under real-time evolution, all eigenvalues of eigenstates except for the steady state, |ρ 0 ⟩⟩, have a real-part Re(λ i̸ =0 ) < 0 and hence exponentially decay over time.Hence, any initial matrix, independent of its physicality, must necessarily converge to the steady-state of the Liouvillian.We stress now, that up to first order, the optimization algorithm we employ was derived as an approximation to real-time evolution under a Liouvillian and hence gives a good approximation to this kind of behavior.

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 [34,39,41,[50][51][52].The system consists of a chain of spin-1/2 particles.The Hamiltonian part of the evolution is governed by Here, N is the number of sites, σ 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 = √ γσ − i .We choose the units of the interaction and field such that the dissipation strength is γ = 1.In all calculations that follow we set the interaction strength to J/γ = 2.
The steady-state of this model then has simple solutions in two limiting cases.When h → 0 the dissipation dominates the dynamics and the stationary state is a pure product-state with all the spins pointing down In the opposite limit where h → ∞ there is only competition between the local onsite field and the dissipation and so the steady-state ends up again as a product-state, but this time the state on each site is mixed 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 [53].For small system sizes we will also be able to compare to results obtained from exact diagonalization, these are calculated using QuTIP [54].
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 σ x i on the central site and the σ z i σ z i+1 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 ≤ h/γ ≤ 2.5, which is in agreement with the results of Ref. [34].For the LDM a hidden unit density of β = 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 number of hidden units  23) 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 σ x i as a function of h.(b) The expectation value of the σ z i σ z i+1 correlation function.In both cases the expectation value is taken on the central site(s).The black (solid) lines show the exact result, the NDM with β = 1 is the green (dashed) lines, while the LDM with β = 1 is orange (doted) and with β = 2 is red (dot-dashed).(c) The absolute value squared of the cost function in Eq. 8 for the two LDM results.(d) The expectation value ⟨L † L⟩ used as a cost function for the NDM ansatz employed by NetKet.For both cases with β = 1 we used 4500 samples and optimized for 1000 steps with a learning rate of η = 10 −2 and a regularization of λ = 10 −2 , for the expectation values we used 500 diagonal samples.In the β = 2 case we used 6000 samples and 2000 steps and 800 diagonal samples.The other meta parameters were the same in both cases.
in the LDM we also increase the number of parameters so that, at β = 2, there are 246 parameters.We see that for the LDM ansatz increasing the value of β and hence the number of parameters used is able to significantly decrease the deviation from the exact result, thus we may use the hidden unit density 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 β 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 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. ( 25), 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 without using hidden units.
A further comparison is shown in Fig. 4.This figure Increasing the hidden unit density improves the accuracy of the results.Both calculations were done at h/γ = 1 using the same parameters as Fig. 3.
shows how the cost function of the LDM evolves for two different values of β over 6000 steps at one of the most difficult points, h = γ.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 β = 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.
To test that the state we obtain is physical we construct the full density matrix from the LDM.In the top panel of Fig. 5 we show how the real part of the smallest eigenvalue of this constructed density matrix evolves.For a randomly initiated state the minimum eigenvalue is negative, which indicates a non-physical density matrix.However, as the optimization goes on the minimum eigenvalues become closer to 0 or becomes positive, indicating that the final density matrix is positive semidefinite.The middle panel of this figure shows the imaginary parts of all of the eigenvalues are quickly suppressed which indicates a hermitian matrix and hence we can be we have produced a physical density matrix.Inf Fig. 5(c) we show the fidelity of the state with one obtained from the exact solution of the master equation.We see that, as expected from the results of Fig. 3, the fidelity is highest in the cases where the transverse field strength is either very large or very small and the fidelity is lowest where the steady state has large entanglement around h = 1γ.
We now go on to examine how the accuracy of these approaches scales to larger system sizes.At N = 16 it becomes difficult to use exact methods to compare against, FIG. 6. Steady-state of Eq. ( 23) as a function of field strength, similar to Fig. 3, for a larger system size N = 16.The expectation value of σ x i and σ z i σ z i+1 on the central site(s) are shown in panels (a)-(b) and the relevant cost functions in panels (c)-(d).The LDM (orange, dotted lines) and NDM (green, dashed lines) results are compared to those obtained from MPS simulations (black, solid lines).We used β = 1.4 in the case of LDM and β = 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 η = −3 we took 9000 Monte Carlo samples at each step and a regularization of λ = 3 × 10 −3 .We used 800 diagonal samples to estimate the expectation values.however this model is straightforward to solve with MPS simulations which we found to be fully converged for a bond-dimension of χ = 7. Results of these calculations are shown in Fig. 6.The number of parameters for the NDM is 1104 while the LDM has 1126.For reference, a bond dimension of χ = 7 corresponds to 1952 matrix elements in the MPS.We see very similar behavior to the All calculations ran for 7000 steps.To accommodate for higher parameter counts we increased the number of samples with β from 9000 at β = 1 to 17000 at β = 2.Other parameters are the same as in Fig. 6.N = 6 case, both approaches are more difficult to converge in the region of intermediate h/γ and the cost function for the LDM has a peak in this region.In Fig. 7 we show how the convergence can again be improved by increasing the hidden unit density.Here we choose h = 2γ as this is the point where the convergence is worst.We see that as β 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.

A. 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 [50,52,55] but keep the dissipation processes the same as for the previous model.In this case the dissipation does not explicitly break the 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 [51,55,56].In 1D these phase transitions turn into continuous crossovers, but complex correlations still build up when h ∼ J ∼ γ.For a detailed review of the behavior of this model in 1D see Ref. [52].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. 8 we show how both a single site and twosite 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 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. ( 23).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. ( 24)-( 25) are not available, except at very large h → ∞ where the steady-state is the same as given in Eq. (24).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.

B. 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  27) and (b) purity, P = Tr[ρ 2 ], for the σ x σ x -model (upper, red curves) and the σ z σ z -model (lower, blue curves).Comparing with the results of Figs. 3 and 8 we observe a correlation between a large negativity and poor accuracy of the neural network.halves of the system [57].For our purposes we find that the negativity [58,59] provides a useful measure of the correlations which are difficult to represent using the LDM approach described above.Here, ||ρ 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. 9 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. ( 23) with the σ z σ z interaction we see a clear peak in the negativity at around h ∼ 0.9γ 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. ( 24) and ( 25) 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. ( 26) with σ x σ 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 > ∼ 4γ.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 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 6, 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 [60], Recurrent Neural Networks [61], or transformers [62], 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.
The idea is now to project both sides into the nonorthogonal basis and ask, under which conditions they become equal [22]: The solution of this system of linear equations γ are the parameter updates which most closely resemble a step in real-time of size δt.
several Markov Chains and within these Chains, arguing that their ratio should approach unity in the limit of perfect convergence or infinite samples.This gives a good estimate of how close the samples are to representing the true probability distribution.This measure is known in the literature as the Gelman-Rubin R value.

FIG. 3 .
FIG.3.Comparison of the steady-state of Eq. (23) 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 σ x i as a function of h.(b) The expectation value of the σ z i σ z i+1 correlation function.In both cases the expectation value is taken on the central site(s).The black (solid) lines show the exact result, the NDM with β = 1 is the green (dashed) lines, while the LDM with β = 1 is orange (doted) and with β = 2 is red (dot-dashed).(c) The absolute value squared of the cost function in Eq. 8 for the two LDM results.(d) The expectation value ⟨L † L⟩ used as a cost function for the NDM ansatz employed by NetKet.For both cases with β = 1 we used 4500 samples and optimized for 1000 steps with a learning rate of η = 10 −2 and a regularization of λ = 10 −2 , for the expectation values we used 500 diagonal samples.In the β = 2 case we used 6000 samples and 2000 steps and 800 diagonal samples.The other meta parameters were the same in both cases.

FIG. 4 .
FIG.4.Convergence of the LDM ansatz for two different hidden unit densities, β = 1 (thin, blue line) and β = 2 (orange, thick line) .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/γ = 1 using the same parameters as Fig.3.

FIG. 5 .
FIG. 5. Properties of the state obtained by using SR at three different values of h and all other parameters the same as the β = 1 results in Fig. 3. Panel (a) shows the real part of the minimum eigenvalue of the full density matrix obtained from the LDM.In (b) we show the sum of the absolute value of the imaginary parts of the eigenvalues of the density matrix.Panel (c) gives the fidelity of the ansatz with the exact density matrix.

FIG. 7 .
FIG.7.Improvement of convergence as a function of β for the N = 16 TFI model at h = 2γ (orange, dotted line) compared to the results of an MPS calculation (blue, solid line).In (a) we show the σ x i σ x i+1 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 β from 9000 at β = 1 to 17000 at β = 2.Other parameters are the same as in Fig.6.

FIG. 8 .
FIG.8.Optimizing the rotated TFI model as in Eq.26 for N = 6.The exact results are in blue, those obtained with the LDM are orange and the NDM in green.Different hidden unit densities are shown by different line-styles.In all cases the optimization was run for 4000 steps, with a learning rate of η = 10 −3 and 2000 diagonal samples to estimate the expectation values.For β = 1 we used 4500 samples and a regularization of λ = 10 −3 .For β = 2 the number of samples was increased to 18000 and the regularization was 10 −2 .In panel (a) we give the steady-state expectation value of σ z i and in (b) we show the two-point correlation function ⟨σ z i σ z i+1 ⟩.Panels (c) and (d) show the relevant cost functions for each ansatz.

FIG. 9 .
FIG. 9. (a)The steady-state entanglement negativity, defined in Eq. (27) and (b) purity, P = Tr[ρ 2 ], for the σ x σ x -model (upper, red curves) and the σ z σ z -model (lower, blue curves).Comparing with the results of Figs.3 and 8we observe a correlation between a large negativity and poor accuracy of the neural network.