Learning Quantum Dissipation by Neural Ordinary Differential Equation

Quantum dissipation is closely related to such phenomena as decoherence, spectrum broadening and heating, all of which stand as serious obstacles in research areas ranging from quantum computation Ladd2010 ; Nielsen2010 ; Preskill2018 and simulation Georgescu2014 ; Bloch2012 , quantum information storage Lvovsky2009 , to quantum metrology Giovannetti2011 ; Pezze2018 , and sensing Degen2017 . The microscopic origin of dissipation is the breaking of isolation of the quantum system, i.e., the system inevitably interacts with its surrounding environment such that the information leakage occurs as the ambient degrees of freedom are traced out. Dissipation severely impairs the accuracy of modeling and manipulation of quantum systems.

Considerable efforts have been made to counteract the negative effects of dissipation. For example, quantum-nondemolition-mediated feedback has been used in suppressing the decoherence of a Schrödinger cat state in a cavity Wiseman1993 ; Vitali1997 ; spin-echo Hahn1950 and dynamical decoupling protocols Viola1999 have been widely applied to nitrogen-vacancy centers Lange2010 ; Du2009 ; Abobeih2018 ; Laraoui2013 , cold atoms Almog2010 , trapped ions Wang2017 and super-conducting circuits Guo2018 ; Pokharel2018 . These techniques often work with certain prior knowledge of the system-environment interactions. For example, in the nitrogen-vacancy centers, strong bias fields were introduced to circumvent the transverse coupling Du2009 or to effectively establish the bath correlations of the carbon nucleus Laraoui2013 . Furthermore, the coupling manner may also be recognized in advance, if it is magnetic or electric, linear or quadratical. However, such prior information is generally unachievable, especially for many-body systems interacting with complex environments.

In this paper, we propose a data-driven scheme to reconstruct the Markovian open quantum systems, where the dataset comes from the discrete observations of the relaxation dynamics under certain probes. Based on the dataset, we adopt the neural ordinary differential equation (NODE), a recently developed machine learning algorithm Chen2019 , to learn the Liouvillian of the open system by inversely solving the Lindblad (or Gorini-Kossakowski-Sudarshan-Lindblad) master equation Lindblad1976 ; Gorini1976 ; Gardiner2004 . A number of relevant works have been reported on such a topic. The operator expansion Franco2009 and the eigensystem realization algorithm (ERA) Zhang2014 ; Sone2017A were firstly used in the time-trace of the Hamiltonians; the ERA was later extended to the Markovian dissipative systems Zhang2015 and the estimation of the system size Sone2017B . Several eigenstate- Qi2019 ; Bairey2019 or steady-state-based approaches Bairey2020 were further developed. A generic approach for local Hamiltonian tomography was reported using the initial and the final observations of quench dynamics Li2020 . Also for this topic, several traditional machine-learning architectures, such as the fully connected neural network Xin2019 , long short-term memory Che2021 , and convolutional neural network Ma2021 were recently involved.

The novelty and contributions of our present work mainly lie in the following points. At first, our approach can be applied to either the Hamiltonian tomography (closed system) or the reconstruction of the Markovian dissipations (open system), with no need of prior information on the specific structures of the Hamiltonian or Liouvilian. Second, we adopt a relatively new machine-learning algorithm (namely the NODE Chen2019 ) to deal with the gradient, with such issues encountered by the traditional machine-learning algorithms as the gradient explosion and vanishing Pascanu2012 being avoided. Furthermore, we also study the learning efficiency of datasets to facilitate the data acquisition in realistic experiments. We thus expect our work can play an active role in effective modeling and guiding new system-environment decoupling protocols for various open quantum systems.

We consider a quantum system being coupled to the environment, as is schematically displayed in Fig. 1(a). Under the Born-Markov approximation, the equation of motion of the system is governed by the Lindblad master equation Lindblad1976 ; Gorini1976 ; Gardiner2004

where $\rho$ is the density matrix of the quantum system and $L$ is the Liouvilian super-operator in the form of (setting $\hbar=1$)

Here, $H$ denotes the dissipation-independent Hamiltonian, $J_{k}$ is the dissipative operator (or jumping operator) in the $k$-th channel, and $\{J_{k}^{\dagger}J_{k},\rho\}=J_{k}^{\dagger}J_{k}\rho+\rho J_{k}^{\dagger}J_{k}$ accounts for the normalization of $\rho$ as if no jump occurs.

Conventionally, with certain Liouvilian $L$ (namely $H$ and $J_{k}$), one can obtain the solution $\rho(t)$ by numerically solving the Lindblad equation (Eq.(2)) using ordinary differential equation (ODE) solvers such as the Euler or the Runge-Kutta Scherer2013 . In contrast, the goal of this paper is to determine the dissipation $J_{k}$ and even $H$ when the dynamical behaviors of density matrix $\rho(t)$ or certain observations $\langle O(t)\rangle=\text{Tr}[O\rho(t)]$ are given. In other words, we aim to deal with an inverse problem to reproduce $L$ from the sequential dataset $\mathcal{D}=\{\rho(t_{n})\text{ or }\langle O(t_{n})\rangle\}$, as illustrated in Fig. 1(b). Additionally, for the sake of experimental application, we consider that $\mathcal{D}$ is purely constructed by observations in the following discussion, i.e. $\mathcal{D}=\{\langle O(t_{n})\rangle\}$, since compared to $\rho(t)$, acquiring $\{\langle O(t_{n})\rangle\}$ is experimentally simpler and more feasible.

We adopt two different probes to generate the dataset $\mathcal{D}$. i) Time-dependent probe: for a fixed $\rho(t_{0})$, one imposes a t-dependent control $H_{\text{prob}}(t)$, i.e.,

where $H_{\text{prob}}=\mathcal{P}(t)\cdot\mathbf{A}=\sum_{\mu}p_{\mu}(t)A_{\mu}$ with $p_{\mu}(t)$ being some smooth series and $A_{\mu}$ the according control operators. ii) Time-independent probe: with fixed $H$ and $J_{k}$, one diversifies $\langle O(t)\rangle$ by varying $\rho(t_{0})$, which suits the experimental scenarios in which preparing different initial states is more convenient. For each evolution, we perform the measurement $\langle O(t_{n})\rangle$ uniformly at discrete time within the time range $t_{n}\in[0,t_{N})$, where $t_{n}=nt_{N}/N_{\text{ts}}$ with $N_{\text{ts}}$ being the time steps of measurement, which forms a data batch. The entire dataset is constructed by different batches, i.e.,

where $N_{\text{bs}}$ denotes the batch-size. Therefore, the total number of data points in $\mathcal{D}$ is equal to $N=N_{\text{bs}}N_{\text{ts}}$. Here, we mentioned that the discrete-time observations with equal intervals are considered for experimental convenience, which is however not a necessary condition for the learning algorithm that we will show below.

To learn the Liouvilian from $\mathcal{D}$, we adopt a machine learning algorithm called the NODE Chen2019 . The NODE builds upon an ansatz $\dot{\rho}(\rho(t),t,\theta)$ with $\theta=\{\alpha,\gamma\}$, where $\alpha$ and $\gamma$ denote the parameters to be learnt in $H$ and $L_{k}$, respectively. The learning process is illustrated in Fig. 1(c), and it contains two parts. First, with $\dot{\rho}(\rho(t),t,\theta)$ and the initial state $\rho(t_{0})$, we propagate the Lindblad equation (2) forward by using certain ODE solvers, and then obtain a series of predictive solutions $\{\rho_{\text{pred}}(t_{n})\}$; the loss function $\mathcal{L}$ is defined as a functional of $\langle O(t_{n})\rangle_{\text{pred}}$ and $\langle O(t_{n})\rangle_{\text{real}}\in\mathcal{D}$, which effectively measures the distance between $\mathcal{D}$ and the NODE predictions. The purpose of learning is to minimize $\mathcal{L}$ by adjusting $\theta$. To this end, the NODE introduces an adjoint field defined by $a(t)=\partial\mathcal{L}(t)/\partial\rho(t)$. Through backward propagating $a(t)$ from $t_{N}$ to $t_{0}$, we obtain $d\mathcal{L}/d\theta|_{t=t_{0}}$ which is then be used to update the parameters in the way of $\theta\rightarrow\theta-\lambda(d\mathcal{L}/d\theta|_{t=t_{0}})$ with $\lambda$ being the learning rate. More calculation details about the NODE algorithm can be found in Appendix A.

We apply this method to two concrete examples. In the first example, we consider a one-body spin system with the spin quantum number being $S=3/2$. Large spin systems are active in various research areas of quantum physics ranging from high-spin quantum dots Klochan2011 ; Doherty2013 , multi-component quantum gases Kawaguchi2012 to unconventional superconductors Wang2018 . In the second example, our system is a many-body spin-1/2 chain. Spin chains stand as fundamental models in condensed matter physics and quantum computation, which are closely related to quantum criticality Sachdev2011 , topological phase of matters Chiu2016 ; Wen2017 , etc.

In the above, we illustrated the capability of the algorithm in learning open quantum systems. Now, we turn to the question of learning efficiency — how to collect the data points can make the learning more efficient? We emphasize the importance of this question on data acquisition in realistic experiments, especially when collecting data points is expensive or time-consuming.

To make this question simple and clear, let us focus on the situation that $\gamma$ are the only unknown parameters that we would like to learn. Moreover, since different batches are independent of each other in $\mathcal{D}$, we consider $\mathcal{D}$ only contains one batch such that the total number of data points is $N=N_{\text{ts}}$. As mentioned before, the NODE adjusts parameters $\gamma_{\mu}$ according to $d\mathcal{L}/d\gamma_{\mu}$, which motivates us to define

where $\mathcal{L}_{n}$ means the local loss function with respect to a single data point $\langle O(t_{n})\rangle$. The physical meaning of $\eta$ is quite clear, it characterizes the averaged sensitivity of $\mathcal{L}$ on $\gamma$ such that a large $\eta$ would speed up the learning process. To further simplify the Eq. (11), we replace the average of individual derivatives by the derivative to the mean value of $\gamma_{\mu}$, i.e.,

For weak dissipation $\bar{\gamma}\ll\bar{\alpha}$, $\eta$ is closely related to

which is composed of two parts: $|f(t)|$ is a fast-oscillating term characterized by the energy scale of $\bar{\alpha}$, while $te^{-\Delta_{1}t}$ is a slowly varying envelop depicting the dissipation-induced damping of $\langle O(t)\rangle$, where $\Delta_{1}\propto\bar{\gamma}$ is the real part of the Liouvilian gap. Obviously, the envelop of $\chi(t)$ is a non-monotonic function being maximized at $t=t_{\text{dc}}$ (see $\chi(t)$ and the $\chi^{2}(t)$ in Fig. 4(a)), with $t_{\text{dc}}=1/\Delta_{1}$ being the decoherent time.

The most commonly used two loss functions, mean absolute error (MAE) and MSE, are linearly and quadratically proportional to $\chi(t)$ as $\bar{\gamma}$ approaches its realistic value $\bar{\gamma}_{\text{real}}$, i.e.,

where MAE is defined by $\mathcal{L}_{n}=|\langle O(t_{n})\rangle_{\text{pred}}-\langle O(t_{n})\rangle_{\text{real}}|$, while the definition of MSE is presented in Eq. (7). For uniform measurement $t_{n}=nt_{N}/N_{\text{ts}}$, the summations in Eq. (14) can be analytically calculated in the large $N$ limit, i.e., ($t_{N}$ in the unit of $t_{\text{dc}}$)

where both are non-monotonic with maximal values lying at $t_{N}^{\text{opt}}\approx 1.8t_{\text{dc}}$ and $t_{N}^{\text{opt}}\approx 1.7t_{\text{dc}}$, respectively. This indicates that the optimal strategy of uniform data acquisition is to take $t_{N}=t_{N}^{\text{opt}}$, as is schematically illustrated by dots in Fig. 4(a). We benchmark the above analysis on the spin-3/2 learning task using the t-dependent probe and the MSE loss. This task is with $t_{\text{dc}}\approx 9/\bar{\alpha}$. Practically, we set $N=15$ and plot the variations of $\mathcal{E}_{\gamma}$ and $\eta$ in Figs. 4(b) and (c) respectively, where $\eta$ is numerically calculated using Eq. (11). Clearly, $t_{N}=t_{N}^{\text{opt}}=1.7t_{\text{dc}}$ exhibits the largest $\eta$ and the fastest decay rate of $\mathcal{E}_{\gamma}$ throughout the learning process, which is in good agreement with our previous discussion.

Finally, we briefly discuss the effect of the total data number $N$. Therefore, we fix $t_{N}=t_{N}^{\text{opt}}$ and show the dependence of $\mathcal{E}_{\gamma}$ on $N$ in Fig. 4(d). An apparent feature is that the cases with $N<N_{\gamma}=15$ exhibit poor learning results, which is understandable since we need at least $N_{\gamma}$ data points to uniquely determine all $\gamma$, resembling that one needs at least $N_{\gamma}$ equations to determine $N_{\gamma}$ variables. Furthermore, an increasing $N$ can accelerate the learning process, however, this trend will not continue endlessly (comparing the curves $N=15$, $N=25$ and $N=40$ in Fig. 4(d)). For a fairly large $N$, neighboring data points exhibit tiny distinguishment such that adding more data points would be of little help for the learning efficiency.

We proposed a scheme to learn the quantum dissipation of open systems based on the NODE, a machine learning algorithm being able to reproduce the Liouvilian from dynamical observations. The learning process can be accelerated by optimizing the strategy of data collection. There are many follow-up questions. The Lindblad master equation relies on the Born-Markov approximation, which limits the current method to weakly dissipative systems with short correlation time. The generalization from Markov to non-Markov is not so straightforward, since the memory effect leads to an integro-differential master equation Gardiner2000 ; Vega2017 . More advanced techniques are expected to deal with the gradients in relation to the integral kernel. Whether an effective Markovian description can be found for non-Markovian dynamics is still left open. Furthermore, the full-Liouvilian calculation shown above cannot be easily extended to large-scale many-body systems due to the exponential growth of the Hilbert space. One possible solution is to combine the NODE with quantum-trajectory approaches Gardiner2000 ; Daley2014 such as the truncated Wigner method Schachenmayer2015 ; Huber2022 and the tDMRG-quantum-trajectory method Daley2009 . In this regard, several recently developed machine-learning-based solvers Carleo2017 ; Nagy2019 ; Hartmann2019 ; Vicentini2019 ; Mazza2021 ; Liu2022 ; Luo2022 may provide valuable insights. Additionally, with a dissipative model, can machine learning algorithms help design the corresponding protocols for system-environment decoupling? We expect this work, as well as these questions, to stimulate more interdisciplinary studies in the fields of machine learning and open quantum systems.