Inferring Markovian quantum master equations of few-body observables in interacting spin chains

Full information about a many-body quantum system is usually out-of-reach due to the exponential growth -- with the size of the system -- of the number of parameters needed to encode its state. Nonetheless, in order to understand the complex phenomenology that can be observed in these systems, it is often sufficient to consider dynamical or stationary properties of local observables or, at most, of few-body correlation functions. These quantities are typically studied by singling out a specific subsystem of interest and regarding the remainder of the many-body system as an effective bath. In the simplest scenario, the subsystem dynamics, which is in fact an open quantum dynamics, can be approximated through Markovian quantum master equations. Here, we formulate the problem of finding the generator of the subsystem dynamics as a variational problem, which we solve using the standard toolbox of machine learning for optimization. This dynamical or ``Lindblad"generator provides the relevant dynamical parameters for the subsystem of interest. Importantly, the algorithm we develop is constructed such that the learned generator implements a physically consistent open quantum time-evolution. We exploit this to learn the generator of the dynamics of a subsystem of a many-body system subject to a unitary quantum dynamics. We explore the capability of our method to recover the time-evolution of a two-body subsystem and exploit the physical consistency of the generator to make predictions on the stationary state of the subsystem dynamics.


Introduction
Artificial neural network methods have established themselves as a versatile tool in many areas of physics [1,2,3]. Just to mention few examples, their application ranges from the classification of phases of matter [4] over scattering and reflectivity analysis [5,6] all the way to the learning of topological states [7]. Moreover, neural networks have been effectively employed to encode the quantum state of both closed [1,8] and open [9,10,11,12]  ! ⇢ S (t) (a) Global quantum state ⇢(t) (b) Unitary dynamics Lindblad dynamics Figure 1. E↵ective Lindblad generator for the dynamics of a subsystem embedded in a many-body quantum system. The full quantum state ⇢(t) of a many-body system subject to a unitary quantum dynamics evolves according to the unitary (Schrödinger) operator U t , as ⇢(t) = U t ⇢(0)U † t . (a) Properties about a subsystem S of such a many-body system are described by the so-called reduced quantum state ⇢ S , which is obtained by tracing out the remainder of the system (or, as we call it here, the bath), ⇢ S (t) = Tr B ⇢(t). (b) The time-evolution of the reduced subsystem state ⇢ S (t) can thus be obtained by tracing out the bath in the evolved many-body quantum state: ⇢ S (t) = Tr B ⇢(t) (see rightmost arrow in the diagram). If the emergent subsystem dynamics is Markovian, the time-evolution of any initial state of the subsystem, ⇢ S (0), can be formulated as an open quantum dynamics, characterised by a time-independent generator L, which e↵ectively accounts for the presence of the bath (see bottom arrow in the diagram).
from the classification of phases of matter [4] over scattering and reflectivity analysis [5,6] all the way to the learning of topological states [7]. Moreover, neural networks have been e↵ectively employed to encode the quantum state of both closed [1,8] and open [9,10,11,12] quantum systems. Recently, machine learning methods have also been applied to obtain e↵ective dynamical generators for local degrees of freedom of a many-body quantum system [13]. The idea is the following: let us assume that one is interested in the time-evolution of a subsystem S of a larger quantum system, as depicted in the sketch in Fig. 1(a). For convenience, we refer to the complement of S, i.e. the remainder of the quantum system, as the "bath" denoted by B. The full information about properties of S is contained in the so-called reduced quantum state ⇢ S , obtained by "integrating out" the bath B. Assuming that the full quantum system is subject to a Hamiltonian quantum dynamics implemented by the unitary operator U t , this simply means that the time-evolved reduced quantum state of S is given by ⇢ S (t) = Tr B (U t ⇢U † t ) [see Fig. 1(b)], where Tr B indicates the trace over the degrees of freedom of B. Under certain physical conditions, e.g. weak S-B coupling and very large bath B, such a subsystem time-evolution is well captured by Markovian open quantum dynamics [14,15,16]. Formally, this entails the existence of a time-independent generator L of a dissipative quantum dynamics such that [see also Fig. 1 Figure 1. Effective Lindblad generator for the dynamics of a subsystem embedded in a many-body quantum system. The full quantum state ρ(t) of a many-body system subject to a unitary quantum dynamics evolves according to the unitary (Schrödinger) operator U t , as ρ(t) = U t ρ(0)U † t . (a) Properties about a subsystem S of such a many-body system are described by the so-called reduced quantum state ρ S , which is obtained by tracing out the remainder of the system (or, as we call it here, the bath), ρ S (t) = Tr B ρ(t). (b) The time-evolution of the reduced subsystem state ρ S (t) can thus be obtained by tracing out the bath in the evolved many-body quantum state: ρ S (t) = Tr B ρ(t) (see rightmost arrow in the diagram). If the emergent subsystem dynamics is Markovian, the time-evolution of any initial state of the subsystem, ρ S (0), can be formulated as an open quantum dynamics, characterised by a time-independent generator L, which effectively accounts for the presence of the bath (see bottom arrow in the diagram).
been applied to obtain effective dynamical generators for local degrees of freedom of a many-body quantum system [13]. The idea is the following: let us assume that one is interested in the time-evolution of a subsystem S of a larger quantum system, as depicted in the sketch in Fig. 1(a). For convenience, we refer to the complement of S, i.e. the remainder of the quantum system, as the "bath" denoted by B. The full information about properties of S is contained in the so-called reduced quantum state ρ S , obtained by "integrating out" the bath B. Assuming that the full quantum system is subject to a Hamiltonian quantum dynamics implemented by the unitary operator U t , this simply means that the time-evolved reduced quantum state of S is given by ρ S (t) = Tr B (U t ρU † t ) [see Fig. 1 [14,15,16]. Formally, this entails the existence of a time-independent generator L of a dissipative quantum dynamics such that [see also Fig. 1 where ρ(0) = ρ S (0) ⊗ ρ B (0) and ρ S/B (0) is the initial state of the subsystem/bath, respectively. Interestingly, by exploiting simple neural network architectures, it was shown that such an approximate generator L can be found also beyond the typically Inferring Markovian quantum master equations 3 considered settings, and that it can satisfactorily describe the subsystem dynamics whenever non-Markovian effects are negligible [13]. These findings are relevant since they allow one to infer, from the time-evolution of local degrees of freedom, the physical processes underlying the dynamics of quantum (sub)systems. This can, for instance, deliver insight into the dynamical effects behind an order-parameter time-evolution in non-equilibrium phase transitions or shed light on relaxation and thermalization effects in closed systems [17]. Moreover, as we show here, since the neural network provides the dynamical generator of a Markovian time-evolution, the latter can be used to directly target stationary properties of the subsystem. However, the generator L learned by the neural network architecture of Ref. [13] is not guaranteed to implement a valid, i.e. physical, dynamics (see below for details). This is a problem since it can lead to physically inconsistent predictions on the dynamics or on stationary properties of the subsystem state [18,19,20]. In this paper we use a matrix parametrization of the Lindblad generator, and optimise it using data generated with exact numerical simulations. One can see this as an encoding of the matrix elements of the generator into the variational parameters of a simple neural network [see Fig. 2 below]. This interpretation allows us to "learn" the optimal variational parameters by means of the standard tools of machine learning. For actual implementation and optimization of the network, we used the automatic differentiation toolbox Pytorch [21]. Importantly, our parametrization is such that the learned dynamical generator L is constrained to be physically consistent. Such a generator is meant to approximate the dynamics of the subsystem state ρ S (t) in the sense of Eq. (1). While the neural network, which is agnostic to the fact that the dynamics of ρ S (t) is a result of tracing out the degrees of freedom of the bath B, always retrieves a time-independent generator L. How well such an approximation can reproduce the dynamics of the subsystem S [see Eq. (1)] clearly depends on the dynamical regime considered, e.g. strong or weak interactions, determining whether the subsystem dynamics is effectively Markovian or not. Here, we show how the learned generator L can be used both to extrapolate the subsystem dynamics to times which have not been used to train the network and to predict stationary state properties of the subsystem. We illustrate our ideas by applying our neural network architecture to the reduced dynamics of a two-spin subsystem embedded in a larger quantum spin chain.

Formulation of the problem
We consider a many-body quantum system partitioned into a subsystem S and the remainder B, which in our setting plays the role of a bath.
The separation implies that the full Hilbert space H is obtained as the tensor product H = H S ⊗ H B of the Hilbert space of the subsystem (H S ) and of the bath (H B ). Here, we work under the assumption that H is finite-dimensional with dimension m. For instance, for the spin-1/2 models considered later [see Fig. 3 and Eqs. (15)-(16)] m = 2 N , where N is the total number of spins. The quantum state of the full system is described Inferring Markovian quantum master equations 4 by means of a density matrix ρ(t), which must be positive semi-definite and must have trace equal to one in order to comply with the probabilistic interpretation of quantum mechanics. As such, the space of all possible states of the many-body system is the (convex) subspace of all positive semi-definite unit-trace matrices, S(H) ⊂ M(m), where M(m) is the algebra of square matrices of dimensions m. Under the assumptions that the many-body system is subject to a unitary dynamics implemented by the full system Hamiltonian H S+B , its state at time t is given by The full information about the time-evolution of the degrees of freedom belonging to subsystem S is contained in the reduced density matrix ρ S (t). As mentioned in the introduction, in certain settings the dynamics of the subsystem state can be approximated through a Markovian open quantum dynamics implemented by a so-called Lindblad generator L [15,16], see Eq. (1). In these cases, ρ S (t) would effectively obey a quantum master equation where we have dropped the explicit time dependence. The map H accounts for the Hamiltonian contributions to the dynamics while dissipative effects, uniquely associated with the interaction between S and B, are encoded in D. The general form of these terms (for finite-dimensional systems) is The operator H is the Hermitian Hamiltonian of the subsystem only and describes its quantum coherent evolution. The operators F i , with i = 1, 2, . . . d 2 and d being the dimension of the subsystem Hilbert space, form an orthonormal basis of the subsystem algebra of operators M(d). Here, without loss of generality, we choose this basis with the property that all its elements are Hermitian F i = F † i and that F d 2 is proportional to the identity, F d 2 = 1 √ d 1. (Note that this latter term is not included in the double sum appearing in D.) The orthonormality condition thus reads Tr(F i F j ) = δ ij , where δ ij is the Kronecker delta, showing that all elements F i but F d 2 are traceless. The matrix c, with elements c ij is called the Kossakowski matrix, and it describes dissipative effects on S due to its interaction with B. Such a matrix must be positive semi-definite in order for the dynamical map implemented by L, T t = e Lt , to be completely positive [15,16]; we recall here that a map T : M(d) → M(d) is said to be completely positive if, for any dimension k, the map T ⊗ 1 k : positive, that is, it preserves the semi-positivity of the matrices it acts upon. Since c is semidefinite positive, it can be diagonalized by means of an unitary transformation h: h † ch = diag(γ 1 , ..., γ d 2 −1 ), with non-negative eigenvalues γ i . The Lindblad generator can then be represented in its "diagonal" form, via the operators Inferring Markovian quantum master equations 5 where the operators J i are called jump operators. The problem we address in this work is that of learning -by means of neural network methods -a time-independent physical dynamical generator L for the dynamics of a subsystem embedded in a unitarily evolving many-body system (see Fig. 1). We are moreover interested in investigating in which parameter regimes the dynamics implemented by such a generator satisfactorily reproduces the time-evolution of the subsystem. From a technical perspective, it is convenient to pass from a representation of L as a map acting on the density matrix ρ S to that of a matrix acting on a vector representation of ρ S itself. We therefore discuss in the following how the information contained in ρ S can be encoded in a vector using the coherence-vector formalism [19], and subsequently derive the ensuing matrixrepresentation of the Lindblad generator L. We note that while there are several ways of such vectorizing a the density matrix, the one adopted here allows us to rewrite the full quantum dynamical problem in terms of real numbers only (see below), which is convenient for the implementation of neural network algorithms.
Employing the (Hermitian) orthonormal basis of M(d) introduced above, we can write any density matrix ρ S ∈ M(d) as a linear combination of the F i . The coherence is then the vector in R d 2 that gathers the corresponding expansion coefficients: We note that the trace-normalization condition, Tr(ρ S ) = 1, implies [v] d 2 = 1/ √ d which has been taken out of the sum. Furthermore, we note that v is a real vector since both ρ S and the F i are Hermitian. Equation (5) shows the one-to-one correspondence between ρ S and the coherence vector v, which we can exploit to conveniently represent the subsystem's density matrix. To make this idea more concrete, we discuss the example of spin-1/2 particle, whose Hilbert space dimension is d = 2. In this case, the basis can be chosen to be proportional to the Pauli matrices such that the density matrix takes the form The coherence vector parametrization is then , which is reminiscent of the usual Bloch vector (a part for the last unimportant element and the different normalization of the Pauli matrices). To obtain the coherence vector elements from the density matrix one computes the traces Inferring Markovian quantum master equations 6 as highlighted in the above equation, knowing the coherence vector basically amounts to knowing the expectation value of all possible observables of the subsystem. We now need to understand how the quantum master equation (2) translates into an evolution equation for the elements of the coherence vector. This can be done by taking a generic density matrix ρ S as in Eq. (5), computing the action of the generator on it, L[ρ S ], and expanding this matrix into the basis formed by the operators F i . This procedure, detailed in Appendix A, provides the matrix representation L of the map L which evolves the coherence vector v through the equation Here, analogously to what done for L, we have decomposed the matrix representation of the generator L into a coherent Hamiltonian part, H, and a dissipative one, D. We note that, since the generator L is Hermiticity-preserving, the matrices H and D are real-valued. In particular, as shown in Appendix A, the matrix H is given by where ω defines the expansion of the Hamiltonian matrix H in Eq.
The dissipative part is instead given by (see detailed derivation in Appendix A) where the so-called symmetric structure constants d ijk are given by , and the matrix c is the Kossakowski matrix appearing in Eq. (3).
As we show below, by exploiting a neural network we can learn the matrix representation L, which propagates the degrees of freedom of a subsystem of a manybody system undergoing unitary quantum dynamics. The learned Lindblad generator L will by construction be physically valid, i.e. it will implement a completely positive and trace-preserving dynamics. However -while it is always possible to find L -how well the learned subsystem dynamics reproduces the exact one depends on how much the latter can be actually approximated by a Markovian open quantum dynamics.

The architecture and the training procedure
In the following we discuss how the matrix L is learned by a neural network through training it with data obtained from simulating the exact unitary many-body quantum dynamics, see sketch in Fig. 2. According to Eq. (9) we have the decomposition L = H + D. Here, H is defined in terms of the d 2 − 1 unconstrained real parameters ω appearing in Eq. (10). The matrix D is parametrised by the complex Hermitian matrix c [see Eq. (11)], which is constrained to be positive semi-definite. To enforce this constraint by construction, we express c as c = Z † Z for some complex matrix Z = X + iY , where X and Y are real matrices. The parameters of our neural network which are to be found via training are thus θ = {ω, X, Y }, and our parametrization automatically ensures the complete positivity of the dynamics generated by L. Note that, since we want to learn a Markovian dynamics, we assume the parameters θ = {ω, X, Y } to be time independent. Moreover, only the symmetric real matrix Re(c) = X t X + Y t Y and the skew symmetric real matrix Im(c) = X t Y − Y t X appear in Eq. (11) [see also Eqs. (A.9) and (A.10)], showing that the network indeed only makes use of real parameters.
Given an input coherence vector at time t, v in = v(t), the network M (see Fig. 2), which is a function of parameters θ, is trained to output , that is the coherence vector after a discrete time step of length dt. In this way, the network learns how to propagate the coherence vector by an infinitesimal time-step and this information fully specifies the matrix L [cf. Eq. (9)].
Training data consist of a set of trajectories of the time-evolution of v, starting from an initial state v 0 up to a fixed time T ‡. Each trajectory is thus a list of snapshots Here, the expectation E is taken over the training dataset D.
To summarise, the parameters of the network are θ = {ω, X, Y }, where ω is a real vector defining H(ω) and X, Y are real matrices defining D(X, Y ). In particular, θ = {ω, X, Y } provides the "weights" of the network. One needs first to use the matrices X and Y to compute the semidefinite positive matrix c = (X − iY ) t (X + iY ). After a choice of the basis for the Hilbert space of the subsystem has been made, one can compute the structure constants and use them to map the variational parameters to the matrices H and D via Eqs. (10) and (11). Then H(ω) and D(X, Y ) are applied to v in as 13) ‡ More specifically, the training data is formed by 50 exactly evolved trajectories, each with the initial conditions as discussed in subsection (4.1). The batch size is 256 with 512 batches per epoch, over 20 epochs. Of the generated data, 80% is used as training set, while 20% is employed as validation, to check that the model does not over fit. § Hyperparameters used in Adam were learning rate Alpha = 10 −3 , Beta 1 = 0.9, Beta 2 = 0.999 and Epsilon = 10 −8 . . The goal is to find the optimal parameters θ defining the matrices H and D by training the network with exact simulation data on the input and output coherence vector. The network is parametrized in a way that the learned generator is automatically guaranteed to be physically consistent. The elements of the matrices H and D are functions of the variational parameters -the weights of the networkθ = {ω, X, Y } through equations (10) and (11) and must be optimized according to the loss function, Eq. (12). and the variational parameters are optimized by minimizing the loss function in Eq. (12). Here dt is the time step used in the exact integration of the dynamics upon which the network is trained. In general, in order to learn all the relevant dynamical features of the subsystems one should make sure that dt is sufficiently small, so that all timescales are captured. We have considered dt = 0.01/Ω, as we observed that with such a time step the details of the exact evolution were well reproduced. Since creating the artificial dataset though exact diagonalization was the most time consuming task. Therefore, it not efficient to use a too small dt. Once the network has been trained, it is possible to retrieve the Hamiltonian part H and the Kossakowski matrix c, we give a concrete example of this in Appendix C, where we also give a brief explanation of their forms.

Many-body models and subsystem of interest
For testing our ideas, we use the above-discussed neural network architecture to learn the generator of the reduced dynamics for two neighbouring spins embedded in a onedimensional chain. The spin chain is thus partitioned into two parts as shown in Fig. 3. The first one is formed by two nearest-neighbouring sites which for convenience are labelled as spin 1 and spin 2. This part acts as the "subsystem of interest". The other part is the remainder of the lattice and is regarded as the bath. We assume that a Hamiltonian H S+B governs the dynamics of the full system quantum state ρ(t). Training data is computed by evolving this state according to the unitary evolution operator U t = e −itH S+B and tracing out the degrees of freedom of the bath. As initial Inferring Markovian quantum master equations 9 state, ρ 0 = ρ(t = 0), of the unitary many-body dynamics we consider product states of the form The reduced density matrix at time t is given by ρ S (t) = Tr B U t ρ 0 U † t . We consider two different different quantum spin chain models, as shown in Fig. 3, and refer to them as model I and II. Model I has closed boundaries, i.e. the spins form a ring, while model II has open boundaries. For Model I, the full subsystem-bath dynamics is governed by a nearest-neighbour interaction Hamiltonian The operator n i = 1+σ z i 2 denotes the projector onto the spin up state of the i-th spin. The constants V and V are the strength of the interaction among neighbouring spins: V is associated with interactions among the bath spins, while V is associated with interactions between the subsystem sites and between the subsystem and the bath sites next to it [see also Fig. 3(a)]. The contribution proportional to Ω, which is the same for all sites, parameterises a transverse "laser" field and drives Rabi oscillations between single spin states. The motivation behind studying the Hamiltonian H I is that the first and the second spin define the subsystem and their interaction with the bath is assumed to be controllable.
Model II, which features open boundary conditions, is governed by a Hamiltonian that features power-law interactions: where α ≥ 0 defines the power with which the interactions decay over distance. Note, that unlike in Model I, there is no specific choice for the coupling constant among bath and system spins. In Appendix C we adopt Hamiltonian II to give an example of the form of the retrieved Hamiltonian part and Kossakowski matrix of the Lindblad generator. There, we deem the restriction of Eq. (16) to a system of two spins as H II , and undertake a brief study of how the bath affects the form of this Hamiltonian. Hamiltonians (15) and (16) are variants of the Ising model with transverse and longitudinal field as well as to the so-called PXP model [23,24]. Experimentally they can be realized, for example, with Rydberg atoms [25,26,27,28,29].
Given that we focus on a system whose reduced density matrix is that of two spins, the dimension d of the reduced Hilbert space H S is four, and a natural choice for the basis {F i } d 2 i=1 is given by the two-spin Pauli group generators Hence, {F 1 , F 2 , ..., F d 2 } = {1 2 ⊗ σ x /2, 1 2 ⊗ σ y /2, ..., 1 4 /2}. If N is the length of the spin chain, the subsystem of interest is identified by the spins at positions 1 and 2, such that each element in v corresponds to either a (17)

The initial conditions
For both models, the initial conditions of the system are chosen as the product of a thermal state for the bath ρ Bath ∝ e −βH I/II , with β the inverse temperature, and a valid density matrix ρ S for the system S. Model II was studied to investigate whether and how the algorithm is able to capture the dynamics in the presence of longer-ranged interactions. To be able to focus on this aspect, we decided to take a fixed initial state for the bath (without changing the temperature) and, for simplicity, we considered the infinite-temperature ones (β = 0). For the study of both models we choose the initial density matrix ρ S in the following way: we consider two random d × d real matrices M and N whose entries are taken from a Gaussian distribution centred at zero and with standard deviation one and we construct the density matrix according to Inferring Markovian quantum master equations

11
This ensures that ρ S (0) is a Hermitian positive semi-definite unit-trace matrix.

Results and discussion
In the standard treatment of open quantum systems, the dynamics of the reduced state of a subsystem can be approximated by means of Markovian open quantum dynamics only when certain conditions are met [14] (see e.g. Ref. [30] and references therein for a different derivation of a Markovian quantum master equations). These include, for instance, a weak subsystem-bath coupling, an infinitely large bath with a continuous dispersion relation and a large separation of time-scales between the subsystem and the bath dynamics. For the Models I and II [see Eq. (15) and (16)], not all of the above conditions are met. For instance, while it is possible to tune the parameters in a way that the interaction between S and B is weak, the bath B will always be, in our setting, a finitedimensional object whose Hamiltonian possesses a discrete spectrum. Nonetheless, given that our network is capable of retrieving a time-independent generator for the subsystem from training data, it is natural to ask whether the dynamics implemented by such a generator can provide an approximate description of the subsystem state also in the considered settings.
To determine whether, and for which parameter range, a Markovian dynamical description is a valid one for the subsystem and, thus, whether the network can be used to predict the dynamical behaviour of its observables, we investigate the accuracy of the Markovian approximation obtained through the network by varying the parameters (β, V ) and (α, V ) for model I and II, respectively. To quantify the error made in the approximation, we define the error measure Here, ρ exact is the (time-evolved) quantum state for S obtained from the exact diagonalization of the full many-body problem, while ρ network is the subsystem dynamics as predicted by the network. The trace norm σ 1 for a d × d complex matrix σ, with eigenvalues λ 1 , . . . , λ d is given by In order to remove the dependence of the error from the specific initial condition considered, we consider the quantity I Err , defined in Eq. (19), averaged over 10 different trajectories, each with an independent initial condition (18) [shown in Figs. 4 and 5].

Subsystem dynamics
We start our discussion by considering Model I, see Eq. (15). As apparent from Fig. 4(a), for relatively low values of the coupling between the subsystem and the bath, the generator learned by the neural network can provide a good approximation of the reduced subsystem dynamics. This is also manifest from the bottom panels in Fig. 4(b); the time-evolution of the expectation values of subsystem observables, as predicted by the network, is in very good agreement with exact numerical data, both for single-site observables and for two-site ones. On the other hand, Fig. 4(a) clearly shows that the neural network approximation becomes less and less reliable upon increasing the coupling between the subsystem and the bath. We stress that this is not witnessing an issue occurring in the training of the neural network for these parameter regimes. The reason for this increase of the error is that, consistently with what expected form the weak-coupling approximation, for large enough coupling the reduced   Fig. 4(a)-(c), no distinction is made between extrapolation and interpolation region. For Model II the error is the highest for the smallest of the system sizes investigated, as one may expect. No dependence on system size is discernible instead for model I.
subsystem dynamics becomes non-Markovian and cannot be approximated by the timeindependent generator learned by the network. We also note that a similar worsening of the Markovian approximation of the dynamics can be observed when lowering the temperature of the bath (increasing β), as becomes transparent from Fig. 4(a). As a first application of the learned generator of the subsystem dynamics, one can use it to extrapolate the dynamics of the subsystem observables to times which have not been seen during the training procedure. This is in particular true for the range of parameters for which the error in approximating the dynamics with the one retrieved by the neural network is relatively small, as shown in Fig. 4(b). The rationale is that whenever the subsystem dynamics is Markovian, the dynamical generator is time-independent and, once this is learned, it can also be used outside the training time-window.
For Model II we find that when the interaction between sites decays fast enough (α ≥ 1) and when the coupling is low (V /Ω ≈ 0.1), the error I Err remains relatively small [see Fig. 4(c)]. In this regime, just as before, the learned generator can be used also for making predictions beyond training times, as shown in Fig. 4(d). On the other hand, when the power α decreases, we observe an increase in the approximation error [see Fig. 4(c)]. This is due to the fact that the concomitant increase of the range of the interactions amounts to increase the interaction between S and B. It is thus reasonable to expect that non-Markovian effects become more pronounced.
In the following we investigate the behaviour of the error, for both models, when changing the length of the spin chain, i.e. modifying the size of the bath, see Fig. 5. For Model I, we observe that the error remains very similar for the different system sizes explored. This may suggest that, in the case of nearest-neighbour interactions and for the model at hand, the considered sizes can be already considered sufficiently large for the remainder of the many-body system to act as a proper bath. For Model II, instead, there appears to be an improvement in the approximation upon increasing the size of the bath. This may be related to the fact that for long-range interactions, revivals in the subsystem dynamics -due to the back-flow of information from the bath into the subsystem -have stronger effects for smaller system sizes.
For both models, our observations suggest that a Markovian approximation may be indeed justified, also for a bath formed by discrete degrees of freedom, provided that the latter is sufficiently large. Moreover, we observe that for all parameter regimes, even in those in which the Markovian approximation is not quite satisfactory, the error tends to stay bounded. This means that the predicted long-time value of local observables seems to overall capture the actual behaviour of the subsystem, at least within the considered time-window [cf. Fig. 4(b-d)]. This seems to suggest that even for large errors in the interpolation, one may find actually a relatively small error in the extrapolation, see Fig. 4. This effect is seen in several parameter regimes, where the expectation values of subsystem observables is almost constant, but features small amplitude oscillations in the extrapolation window, as shown in Fig. 4.
In order to have a fairer comparison of the error between initial transient regime with large oscillations and the long-time one in which oscillations have very small amplitude, we introduce a relative error measure. We consider the so-called fraction of variance unexplained, or FVU. We compute this quantity between the coherence vector predicted by the network v network i and the one retrieved from the exact simulation v exact i , as in the following equation .
The set over which the variance Var is computed is given by the time snapshots with discrete time step dt, v i = (v i (t 0 ), v i (t 0 + dt), ...). The fraction of variance unexplained normalises an error to the variance of the signal (here the exact coherence vector).
The results for Model I and II are plotted in Fig. 6 for system sizes N = 11 and N = 10, respectively, while results for smaller system sizes are reported in Appendix B. Comparing Fig. 6 with Fig. 5, we note that considering the FVU provides a smoother behaviour of the error. In particular, it allows to distinguish between the regions for which the error in the approximation I Err is small simply because the overall signal has small variance, from those where the error is small because the neural network dynamics correctly captures the behaviour of the subsystem. The latter situation occurs for small values of the coupling strength in both Models I and II and for fast decaying interactions in model II [cf. Fig. 6]. In this region the Markovian approximation is correct, and the network is able to reconstruct the dynamics almost exactly.

Stationary behaviour
Once the network has learned the matrix-representation L of the Lindblad generator, it is also possible to investigate stationary properties of the system [31,9]. In our case, we can do this by studying the stationary-state coherence vector, which is nothing but the eigenvector v st associated with the zero eigenvalue of L: Lv st = 0 = dvst dt . Such a vector provides the stationary density matrix of the subsystem S, see Eq. (23), which in our case is by construction ensured to describe a physically consistent quantum state. In principle, there is no reason to expect that the long time coherence vector v exact (t) will converge to v st . This is because i) we are training the network for a finite time-window and ii) the bath is finite, and thus one would expect to observe, for long times, recurrence and revivals in the time-evolution of system observables. These effects are associated with a re-entering of the information scrambled from the subsystem into the bath into the subsystem again, and are associated with non-Markovian behaviour. Nonetheless, what we observe (shown in Fig. 7), is that in some parameter regimes the agreement between v st and the long-time behaviour of the exact dynamics is remarkably good.
To be more quantitative, let us define E gap as the smallest, in modulus, among the real parts of the non-zero eigenvalues of L and let τ be 1/E gap . This latter quantity represents the time-scale of the approach to stationarity in a Markovian open quantum system. Moreover, in order to circumvent the issue that the exact subsystem dynamics always presents oscillations (also due to the finiteness of the bath), we define the where Thus, the error measures the distance between the stationary behaviour predicted from the network results and an averaged long-time behaviour of the exact solution. In particular, in order to be sure to address the stationary behaviour, at least in relation to what predicted by the network, we choose a = 5, b = 10. Within this time window -provided that the exact dynamics is correctly captured by the time-independent generator -the expectation value of local observables should have already converged to their stationary value. Given the finiteness of the bath, these will -as mentioned above -display residual oscillations around such an average stationary value which would be described indeed by ρ * exact . The discrete time step dt for the exact integration and for the training of the network is chosen to be dt = 0.01/Ω as before. Once the training was complete, we could retrieve the Lindblad generator L and its eigenvalues, and in particular E gap and the time scale τ = E gap . In Fig. 7 the error in Eq. (23) is shown. It is computed by averaging over 10 different initial conditions for ρ exact . This error is rather different both from the error of Eq. (19) and more importantly from the FVU of Eq. (21). The FVU is indeed used to compare the oscillations of the predicted and of the exact dynamics. However, by construction, the predicted dynamics is bound to reach a stationary state, and the FVU will thus tend to be large. This signals that the learned dynamics is not able to reproduce the small oscillations around the stationary values. Nonetheless, we can investigate how well the learned dynamics can capture the average value of this oscillations, which should provide their stationary value. To this end, a more appropriate error measure is the one given in Eq. (22). One notices that while there are observables that approach the learned steady state v st , e.g. σ 1 y σ 2 x in model I, others do not. In certain cases, the exact time-evolution of the expectation of observables converges towards a long-time behaviour which depends on the initial conditions. This may be related to the fact that we are considering a finite bath. On the other hand, the learned dynamics always approaches the same stationary behaviour since, due to the existence of a finite gap E gap , the steady state of the Lindblad generator is unique.

Conclusions
We have introduced a simple neural network whose parameters can be exactly mapped onto those of a Lindblad generator. Importantly, such a generator which is learned by the network from exact dynamical data is automatically ensured to be a physically consistent generator of a quantum Markovian dynamics. We have investigated the applicability of such an architecture to two different classes of spin models. Even though the considered physical settings are rather different from those known to give rise to Markovian subsystem dynamics, we find that, in certain parameter regimes, the network model provides a faithful approximation of the subsystem time-evolution. Future developments in the same spirit may include the adaptation to architectures capable of encoding time correlations in time series, such as long short term memory (LSTM) [32] or transformers [33], which would allow for the learning of a time-dependent generator. A different path to achieve time dependence would be that of learning the analytical solutions of the differential equations of motion [34], or by numerically solving them by means of neural networks [35].
We have exploited the time-independent generator learned by the network in order to investigate stationary properties of the reduced subsystem state. This idea looks promising as a path towards capturing relevant long-time features such as thermalization effects [17,36]. Finally, instead of learning the linear evolution of the density matrix, one may think of directly learning the evolution of an order parameter, such as the magnetization or particle density [37]. This would entail machine learning of an effectively non-linear dynamical evolution within the state space of the order parameter. This directly leads to the question whether and how more involved neural network architectures permit an increased performance in determining effective generators and possibly allow an improved quantification of time dependence and non-Markovianity.
The presented results open routes towards the understanding of complex nonequilibrium dynamics through a reduced number of (collective) degrees of freedom. This may ultimately allow to develop simplified descriptions of complex dynamical nonequilibrium processes which is not only of interest in fundamental research but may also be of importance when harnessing many-body phenomena in technological applications.  We here give an example of how one can extract information about the subsystem dynamics from the network and interpret it as the Hamiltonian and dissipative part of the Lindblad equation. In particular we focus here on model II for the parameters V /Ω = 0.1 and α = 0.3. In Fig. C1(a) we plot the matrix elements of the Hamiltonian of model II for a chain of only two spins, already discussed below Eq. (16). In Fig. C1(b) is plotted the one retrieved from the network, that is, the Hamiltonian part H of the Lindblad generator. Their difference is given in Fig. C1(c). As shown, this difference is proportional to a term 1 ⊗ σ z + σ z ⊗ 1.
This contribution can be explained by means of a "mean-field" treatment of the Hamiltonian involving the density-density interactions between the bath and the subsystem. Indeed, it results from taking terms like n i n j , where i is a site of the bath and n i stands for expectation value, while j a site of the two-spin subsystems, once terms proportional to the identity are removed. For the same system it is possible to also retrieve the Kossakowski matrix. We plot the real part of its entries in Fig. C2. The basis of choice for the Hilbert space of subsystem is the same as in the main text, and it is given by {F 1 , F 2 , ..., F d 2 } = {1 2 ⊗σ x /2, 1 2 ⊗σ y /2, ..., 1 4 /2}. This shows that the most relevant entries are associated with dephasing implemented by a "collective" jump operator [cf. Eq. (4)] of the form J ∝ σ z ⊗ 1 + 1 ⊗ σ z . Figure C2. Real part of the Kossakowski matrix. Real part of the Kossakowski matrix retrieved from a six spin system evolved according to Hamiltonian II with parameters V /Ω = 0.1 and α = 0.3. This shows that the largest values are associated with combination of operators σ z ⊗ 1 and 1 ⊗ σ z with equal weights. This means that the relevant dissipative effect is a dephasing noise, associated with a jump operator J ∝ σ z ⊗ 1 + 1 ⊗ σ z . Note that, in this case, the entries of the Kossakowski matrix are small. This is due to the fact that we are considering small values of V /Ω and, in this regime, dissipation is expected to be of second-order in V /Ω.