Optimal identification of non-Markovian environments for spin chains
Correlations of an environment are crucial for the dynamics of non-Markovian quantum systems, which may not be known in advance. In this paper, we propose a gradient algorithm for identifying the correlations in terms of time-varying damping rate functions in a time-convolution-less master equation for spin chains. By measuring time trace observables of the system, the identification procedure can be formulated as an optimization problem. The gradient algorithm is designed based on a calculation of the derivative of an objective function with respect to the damping rate functions, whose effectiveness is shown in a comparison to a differential approach for a two-qubit spin chain.
To exactly process quantum information, accurate models, including parameters, structures, and descriptions of dynamics, for quantum information carriers are required. With these models, sophisticated feedback control strategies can be designed, for example, feedback stabilization of a number state in a cavity , preservation of quantum coherence and entanglement for qubit systems  and linear quantum systems , or coherent feedback rejection of quantum colored noise [4, 5].
However, in practice, an accurate model may not be obtainable since parts of the parameters, structures or dynamics of the quantum system may not be well understood. This would lead to unexpected experimental results or degraded control performance of a quantum control system. For example, in a recent experiment on quantum dots , a calculation based on the theoretical model has a discrepancy from the experimental data on the broadened resonator line-width under suitable parameters, which means that some unknown dynamics of the system are not included in the model. Correspondingly, it is shown that degraded estimation performance can be observed due to ignorance of a noise model for a quantum system . Hence, the problem of how to determine the unknown parameters, structures or dynamics of a “dark” quantum system is an essential problem in quantum control theory which is referred to quantum system identification.
A general identification methodology involves finding the unknown parts of a quantum system from measurement data (e.g., the spectrum of an output field or expectations of observables,) extracted from the system under excitation . The identification method was firstly explored to extract Hamiltonian information for closed quantum systems [9, 10, 11]. Moreover, it is worth considering the identification problem for open quantum systems where the quantum system is disturbed by quantum noises arising from an environment. When the correlation function of the quantum noise is the Dirac delta function; i.e., the correlation time of the quantum noise is very short, the noise and the relevant quantum system refer to quantum white noise and a Markovian quantum system, respectively . For Markovian quantum systems, a continuous-measurement-based method was proposed to identify unknown parameters in a cavity-atomic system  and unknown structures of a spin network . Also, a system-realization-theory-based method in  was extended to estimate unknown parameters in the Hamiltonian of a Markovian spin network from measurement time traces , whose identifiability was discussed in . In addition, an identification method for linear Markovian quantum systems was systematically discussed in .
In contrast to the Dirac correlation function of quantum white noise, correlation functions of quantum colored noise can be more complicated due to the memory effect of non-Markovian environments. A quantum system disturbed by quantum colored noise is referred to as a non-Markovian quantum system, whose dynamics is quite different from that of the Markovian quantum systems . To control non-Markovian quantum systems, it is crucial to acquire the knowledge of the correlation functions of the non-Markovian environment beforehand. Wu et al. proposed a frequency domain approach to the identification of the environment spectrum for a superconducting single qubit at an optimal point . For the quantum dot system in , an augmented system model was presented to explore the structure of the non-Markovian environment . Also, a spectroscopic method was presented to explore the spectrum of spin baths both theoretically  and experimentally , which involves high-energy dynamical decoupling pulses.
In this paper, we aim to extract the correlation functions of non-Markovian environments for spin chains. The dynamics of the spin chains are described by a time-convolution-less master equation in which time-varying damping rate functions characterize the correlation functions of the non-Markovian environment. To identify the damping rate function, we measure the time trace for one local observable of the spin chains which correspond to dynamical variables in a coherence vector representation. Although this formulation for the coherence vector is similar to that in , the dynamical equation for the coherence vector is time-varying due to the damping rate functions. Thus, the method in  is invalid in our case. We alternatively formulate the identification problem for the damping rate function as an optimization problem. Also, we design a gradient algorithm to iteratively recover the damping rate functions for the non-Markovian environment.
Our contribution in this paper is that we provide a systematic approach to acquire unknown damping rate functions in a time-convolution-less master equation for non-Markovian quantum systems. In principle, our gradient algorithm can achieve the damping rate function with a high fidelity, which would help to obtain an exact model of the spin chains in an experiment for quantum information processing. In addition, our method only requires to measure commutative observables of the spin chains, which is easy to be applied in an experiment. Because it avoids measuring non-commutative observables which cannot be accurately measured at the same time according to the uncertainty principle in quantum mechanics.
This paper is organized as follows. In Section II, we introduce the time-convolution-less master equation for non-Markovian spin chains whose corresponding coherence vector representation is also introduced. The identification problem formulation and the gradient algorithm are given in Section III. An example for two qubits in a non-Markovian environment is given in Section IV. Conclusions are drawn in Section V.
Ii Description of non-Markovian spin chains
Ii-a Time-convolution-less master equation
The non-Markovian spin chain we consider in this paper is plotted in Fig.1 where the nodes and the edges represent spin components and their couplings, respectively. The components are disturbed by quantum colored noise arising from an unknown non-Markovian environment. Here, the time traces of local observables are measured and the resulting data is to be processed for identification.
In this paper, the non-Markovian dynamics of spin chains are described by a time-convolution-less master equation 
of the spin chains. The superoperator
represents the internal dynamics of the spin chains where the Hamiltonian of the spin chains describes the internal energy of the system and their couplings. The commutator is calculated as for two matrices and with suitable dimensions. We let hereafter.
The Lindblad dissipative term induced by the non-Markovian environment is expressed as
In contrast to the constant damping rates in Markovian master equations, the damping rate functions for characterizing the correlation of the non-Markovian environment are time-varying, which encapsule the coupling strengthes between the system and the environment and the density state of the environment . Here, we assume that are complex functions, i.e., , where and represent real and imaginary parts of a function, respectively.
In addition, the coupling operators belong to a set , which are an orthonormal basis for the Lie algebra . Their commutation and anti-commutation relations can be calculated as
respectively. The coefficients are the completely antisymmetric and symmetric structure constants of the Lie algebra , i.e., with respect to interchange of any pair of indices . Here, the anti-commutator is calculated as for two matrices and with suitable dimensions and is the Kronecker delta function.
Note that although the variation of the density matrix in the time-convolution-less master equation (1) only depends on the current density matrix mathematically, the time-varying damping rate functions enable (1) to describe a non-Markovian behavior physically; i.e., the energy exchanges between the system and the environment .
Ii-B Dynamical equations in a coherence vector representation
Under the orthonormal basis in , we can expand the density matrix and the Hamiltonian as
respectively, where and the coefficients are the expectation values for the basis ; i.e., . They constitute the so-called coherence vector . Here, is an identity matrix.
where the initial value of the coherence vector is determined by the initial density matrix . The elements of the matrices and can be calculated as with , and , respectively. Here, the corresponding coefficients and are introduced in the calculation of the commutation and anti-commutation relations (4), (5) when transforming Eq. (1) to Eq. (7). Their indexes correspond to the labels of the orthonormal basis for involved in the commutation and anti-commutation relations.
Note that although Eq. (7) is in a similar form as the linear time-invariant dynamical equation for a Markovian quantum system , it is a linear time-varying dynamical equation due to the time-varying damping rate functions .
On the other hand, to extract information of , one can measure expectation values of observables
where are local observables and the symbol represents the expectation of an observable. Also, we can collect the data of as the output of the non-Markovian spin chains. In addition, these observables can be expanded under the orthonormal basis in as with and thus the output (8) can be reexpressed as
where is the th row and th column element of the matrix .
Note that the measured operators in should be commutative since the uncertainty principle in quantum mechanics requires that two non-commutative operators cannot be measured accurately at the same time. Hence, in general, the number of the measured observables should be less than that of the orthonormal basis in , i.e., .
It should be noted that the measurement of expectations of observables is applicable to non-Markovian quantum systems. For example, in a superconducting non-Markovian single qubit system, the measurement of average charge numbers is equivalent to measuring the -component of the coherence vector for the single qubit. The procedure is to measure the same observable many times and then take the average of the results .
Ii-C Reduced dynamical equation
The equation (7) describes the full dynamics of the non-Markovian spin chains. However, it is possible that not all the components in are relevant to the observables (8), i.e., the matrix would have a block-diagonalizable structure and a proper sub-block corresponds to the observables. We can use a filtration procedure in  to find the accessible set of the which is related to the observables .
Thus, a reduced coherence vector corresponding to the observables (8) can be defined and it obeys a reduced dynamical equation
where , , and with suitable dimensions are sub-matrices of , , and .
Iii Gradient algorithm for identifying the damping rate function
Iii-a Optimization formulation for the identification of the damping rate function
The identification problem considered in this paper can be stated as follows. Considering the non-Markovian spin chains governed by the dynamical equation of the coherence vector (7), where all the information of the environment are encoded in the unknown damping rate functions, the environment identification problem is to identify the time-varying damping rate functions from the data of the time traces of the observables (9).
Note that due to the interactions between the spin chains and the environment, in general, there is no decoherence-free subspace for a common non-Markovian system without control pulses; i.e., decoherence channels affect all the components in the coherence vector. Thus, the damping rate functions will be imprinted in the observables. Hence, the identification problem for the damping rate function should be solvable. In addition, since the time-varying damping rate functions result in a linear time-varying system (7), the identification method for Markovian quantum systems in  cannot be applied in our case.
Generally, it is difficult to identify analytically a time-varying function in a dynamical system. Alternatively, we can numerically solve this problem by designing an algorithm. We consider that the system (7) evolves in a total time during which the observables are measured times at equal time intervals , i.e., . Thus, a set of real measurement results can be obtained.
Also, the time-varying damping rate functions are discretized and we assume the in each time interval are constants . Thus in each time interval, we can consider the linear time varying (LTV) system equation (7) as a linear time-invariant (LTI) system and write the dynamical equation in a continuous-time form as
where Eq. (11) can be solved as
Here, we have written the initial coherence vector in the current and next time intervals in an abbreviate form as and , respectively.
For a given set of , an output vector can be generated by using (12) and (III-A). A guessed set of may not be identical to the real damping rate function. Hence, the generated output will be different from the real data . So we define an objective function
to evaluate the distance between the two vectors and .
Therefore, the identification problem for the damping rate functions can be converted into an optimization problem as follows.
For real measurement results , the optimization problem is to find a set of such that the objective function (14) is minimized, that is
Iii-B Gradient algorithm for the identification problem
To design a gradient algorithm for the optimization problem (III-A), it is crucial to obtain the gradient of the objective with respect to in each time interval. By using the chain rule, the gradient can be calculated as
In the following, we will calculate the derivative of with respect to the real and imaginary parts of , respectively. By using a standard formula for derivative of a matrix exponential 
with matrices and of suitable dimensions, we have
where and . When the time interval is small enough such that , we have
where the elements of and can be expressed as and .
Thus, the derivative of with respect to the real and imaginary parts of can be calculated as
where the th row element of the column vector is .
with small step sizes and , we can minimize the objective . With the relation (III-B), our gradient algorithm for the identification problem can be summarized as follows.
Step 1. Choose and measure a set of observables ; Initialize the state of the system , the output , and the step sizes and ; Guess the initial values of ;
Step 2. Calculate to and to with their initial values and ;
Step 3. Compute the gradient of the objective with respect to ;
Step 4. Update by using the relation (III-B);
Step 5. When a termination condition is satisfied, stop the algorithm; Otherwise, go to the Step 2 and start a new iteration.
Note that this algorithm searches optimal solutions according to the gradient and hence it may stop at a local minimum point. Moreover, the final performance of the algorithm depends on the initial guessed solution . When the initial guessed solution is closer to the optimal one, the algorithm will converge faster. The step size may also affect the performance of the algorithm. A small step size would result in a slow convergence process and a large one would lead to an algorithm that saturates in the vicinity of a minimum. An adaptive step size can avoid the above problems. In addition, a proper termination condition can be the completion of a given number of iterations or the attainment of an accuracy of the objective .
Iv A Physical Example
In the example, we consider two coupled qubits immersed in an unknown common environment where the dissipative processes for the two qubits can be considered to be the same. Its non-Markovian dynamics can be described by a time-convolution-less master equation  as
where is the density operator for the two qubits and the ladder operators are defined as and with Pauli matrices
Here, the label is used to index the qubits. The two qubits are assumed to be coupled in an XY interaction form and thus the Hamiltonian of the two qubits is written as
where is the coupling strength between two qubits and is the angular frequency for the -th qubit. In this time-convolution-less master equation (24), all the properties of the unknown environment are combined in the damping rate function which are the same with respect to each dissipative channel for the two qubits. Thus the identification task is to recover the unknown damping rate function . To simulate a real damping rate function induced by the non-Markovian environment, we assume that the real damping rate function is
which results from a non-Markovian environment with a Lorentzian spectrum . The values of the parameters , , and will be given in the following paragraph. This real damping rate function (27) is utilized to generate the real measurement results .
Furthermore, we measure the observable for the first qubit which induces the accessible set . Due to the interactions between the two qubits, the accessible set also includes the operators for the second qubit. Hence, the reduced dynamical equation for the state vector is in the form of (II-C) with matrices
In addition, the parameters of the system are chosen as follows. The angular frequencies of the two qubits are the same; i.e., . The coupling strength between the two qubits is . To simulate the real , the parameters in (27) are chosen as , and . The gradient algorithm starts with a guessed damping rate function
which is plotted as a red line at the first iteration in Fig. 2(a).
Both the step size in each iteration and the termination condition can affect the performance of the gradient algorithm. We choose the step size for updating the damping rate function as . Because a large step size would result in oscillations in the final stage for searching the minimum and a small one would lead to a slow convergence process. On the other hand, we choose to iterate the algorithm for 20000 times as the terminal condition. Sufficient iteration times allow the algorithm to obtain a minimal objective .
With these parameters, the evolution of the damping rate function from the initial guess to the identified one are given in Fig. 2. Initially, the gradient algorithm can significantly update the identified damping rate function as shown in Fig. 2(a). After iterations, the convergence rate of the identification process slows down as shown in Fig. 2(b). This phenomenon also reflects in the reduction of the objective function in Fig. 3. With 20000 iterations, the objective function is down to about and we terminate the algorithm.
The final identified damping rate function is plotted as the dashed red line in Fig. 4. Compared to the real damping rate function represented as the solid blue line, our algorithm can identify the damping rate function well except at the final stage. This is because the measurement results at a final time interval contain limited information of the damping rate function and thus the resulting gradient provides very small updates. We also make a comparison between our method and a differential approach in  whose identified result is plotted as the dashed dot line in Fig. 4. Under the identical measurement result, our method obtains a better result than that by using the differential approach. This is because all the corresponding measured results are utilized to identify the unknown damping rate function in a sample time interval by using our method. However, the differential approach only use the measurement results at the current and next times to estimate the damping rate function at the current time. Its identification result can be improved if we measure the observable densely.
In this paper, we have designed a gradient algorithm for identifying the damping rate function in the time-convolution-less master equation for non-Markovian spin chains. We have formulated the identification problem as an optimization problem that can be solved iteratively by calculating the gradient of the objective with respect to the damping rate function in each time interval. The numerical example on a non-Markovian two-qubit system demonstrates the efficacy of our gradient algorithm. Compared to a differential approach, our method can identify the damping rate function with a high fidelity.
-  C. Sayrin, I. Dotsenko, X. Zhou, B. Peaudecerf, T. Rybarczyk, S. Gleyzes, P. Rouchon, M. Mirrahimi, H. Amini, M. Brune, J.-M. Raimond, and S. Haroche, “Real-time quantum feedback prepares and stabilizes photon number states,” Nature, vol. 477, no. 7362, pp. 73–77, 2011.
-  J. Zhang, R. B. Wu, C. W. Li, and T. J. Tarn, “Protecting coherence and entanglement by quantum feedback controls,” IEEE Trans. Automa. Control, vol. 55, no. 3, pp. 619–633, 2010.
-  S. Xue, M. R. Hush, and I. R. Petersen, “Feedback tracking control of non-markovian quantum systems,” IEEE Trans. on Contr. Syst. Technol., vol. 25, no. 5, pp. 1552–1563, 2017.
-  S. Xue, R. B. Wu, W. M. Zhang, J. Zhang, C. W. Li, and T. J. Tarn, “Decoherence suppression via non-Markovian coherent feedback control,” Phys. Rev. A, vol. 86, p. 052304, 2012.
-  S. Xue, R. Wu, M. R. Hush, and T.-J. Tarn, “Non-markovian coherent feedback control of quantum dot systems,” Quantum Sci. Technol., vol. 2, no. 1, p. 014002, 2017.
-  T. Frey, P. J. Leek, M. Beck, A. Blais, T. Ihn, K. Ensslin, and A. Wallraff, “Dipole coupling of a double quantum dot to a microwave resonator,” Phys. Rev. Lett., vol. 108, p. 046807, 2012.
-  S. Xue, M. R. James, V. Ugrinovskii, and I. R. Petersen, “LQG feedback control of a class of linear non-markovian quantum systems,” in American Control Conference, Boston, 2016, pp. 4775–4778.
-  R. B. Wu, T. F. Li, A. G. Kofman, J. Zhang, Y.-X. Liu, Y. A. Pashkin, J.-S. Tsai, and F. Nori, “Spectral analysis and identification of noises in quantum systems,” Phys. Rev. A, vol. 87, p. 022324, 2013.
-  J. M. Geremia and H. Rabitz, “Optimal identification of hamiltonian information by closed-loop laser control of quantum systems,” Phys. Rev. Lett., vol. 89, p. 263902, 2002.
-  D. Burgarth and K. Yuasa, “Quantum system identification,” Phys. Rev. Lett., vol. 108, p. 080502, 2012.
-  J. Zhang and M. Sarovar, “Quantum hamiltonian identification from measurement time traces,” Phys. Rev. Lett., vol. 113, p. 080401, 2014.
-  H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems. Oxford: Oxford University Press, 2002.
-  H. Mabuchi, “Dynamical identification of open quantum systems,” Quantum Semiclass. Opt., vol. 8, pp. 1103–1108, 1996.
-  Y. Kato and N. Yamamoto, “Structure identification and state initialization of spin networks with limited access,” New J. Phys., vol. 16, no. 2, p. 023024, 2014.
-  J. Zhang and M. Sarovar, “Identification of open quantum systems from observable time traces,” Phys. Rev. A, vol. 91, p. 052121, 2015.
-  A. Sone and P. Cappellaro, “Hamiltonian identifiability assisted by a single-probe measurement,” Phys. Rev. A, vol. 95, p. 022335, 2017.
-  M. Guta and N. Yamamoto, “System identification for passive linear quantum systems,” IEEE Trans. Automa. Control, vol. 61, no. 4, pp. 921–936, 2016.
-  S. Xue, T. Nguyen, M. R. James, A. Shabani, V. Ugrinovskii, and I. R. Petersen, “Modelling and filtering for non-markovian quantum systems,” arXiv:1704.00986, 2017.
-  G. A. Paz-Silva, L. M. Norris, and L. Viola, “Multiqubit spectroscopy of gaussian quantum noise,” Phys. Rev. A, vol. 95, p. 022121, 2017.
-  F. K. Malinowski, F. Martins, L. Cywiński, M. S. Rudner, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, “Spectrum of the nuclear environment for gaas spin qubits,” Phys. Rev. Lett., vol. 118, p. 177702, 2017.
-  R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, L. N. in Physics Vol. 717, Ed. Springer, Berlin, 2007.
-  S. Sastry, Nonlinear Systems. Springer, New York, 1999.
-  I. Najfeld and T. F. Havel, “Derivatives of the matrix exponential and their computation,” Adv. Appl. Math., vol. 16, no. 3, pp. 321–375, 1995.
-  B. Bellomo, A. De Pasquale, G. Gualdi, and U. Marzolino, “Reconstruction of time-dependent coefficients: A check of approximation schemes for non-markovian convolutionless dissipative generators,” Phys. Rev. A, vol. 82, p. 062104, 2010.