Optimal identification of non-Markovian environments for spin chains
11 Optimal identification of non-Markovianenvironments for spin chains
Shibei Xue
Member, IEEE,
Jun Zhang
Senior Member, IEEE, and Ian R. Petersen
Fellow, IEEE
Abstract —Correlations of an environment are crucial for thedynamics of non-Markovian quantum systems, which may notbe known in advance. In this paper, we propose a gradientalgorithm for identifying the correlations in terms of time-varying damping rate functions in a time-convolution-less masterequation for spin chains. By measuring time trace observablesof the system, the identification procedure can be formulatedas an optimization problem. The gradient algorithm is designedbased on a calculation of the derivative of an objective functionwith respect to the damping rate functions, whose effectiveness isshown in a comparison to a differential approach for a two-qubitspin chain.
I. I
NTRODUCTION
To exactly process quantum information, accurate models,including parameters, structures, and descriptions of dynamics,for quantum information carriers are required. With these mod-els, sophisticated feedback control strategies can be designed,for example, feedback stabilization of a number state in acavity [1], preservation of quantum coherence and entangle-ment for qubit systems [2] and linear quantum systems [3], orcoherent feedback rejection of quantum colored noise [4], [5].However, in practice, an accurate model may not be obtain-able since parts of the parameters, structures or dynamics ofthe quantum system may not be well understood. This wouldlead to unexpected experimental results or degraded controlperformance of a quantum control system. For example, in arecent experiment on quantum dots [6], a calculation based onthe theoretical model has a discrepancy from the experimentaldata on the broadened resonator line-width under suitableparameters, which means that some unknown dynamics of thesystem are not included in the model. Correspondingly, it isshown that degraded estimation performance can be observeddue to ignorance of a noise model for a quantum system [7].Hence, the problem of how to determine the unknown param-eters, structures or dynamics of a “dark” quantum system is an
This research was supported by the National Natural Science Founda-tion of China (No. 61873162) and Shanghai Pujiang Program (Grant No.18PJ1405500). This research was also supported under Australian ResearchCouncil Discovery Projects and Laureate Fellowships funding schemes(Projects DP140101779, DP180101805, and FL110100020) and the Air ForceOffice of Scientific Research (AFOSR) under agreement FA2386-16-1-4065.JZ thanks the National Natural Science Foundation of China (No. 61673264and 61533012), and State Key Laboratory of Precision Spectroscopy, ECNU,China.S. Xue is with Department of Automation, Shanghai Jiao Tong Universityand Key Laboratory of System Control and Information Processing, Ministryof Education of China, Shanghai 200240, China. (e-mail: [email protected])J. Zhang is with the Joint Institute of UM-SJTU and Key Laboratory ofSystem Control and Information Processing, Ministry of Education of China,Shanghai 200240, China. (e-mail: [email protected])I. R. Petersen is with Research School of Engineering, Australian NationalUniversity, Canberra, ACT 2600, Australia.(e-mail: [email protected]) essential problem in quantum control theory which is referredto quantum system identification.A general identification methodology involves finding theunknown parts of a quantum system from measurement data(e.g., the spectrum of an output field or expectations ofobservables,) extracted from the system under excitation [8].The identification method was firstly explored to extractHamiltonian information for closed quantum systems [9],[10], [11]. Moreover, it is worth considering the identificationproblem for open quantum systems where the quantum systemis disturbed by quantum noises arising from an environment.When the correlation function of the quantum noise is theDirac delta function; i.e., the correlation time of the quantumnoise is very short, the noise and the relevant quantum systemrefer to quantum white noise and a Markovian quantumsystem, respectively [12]. For Markovian quantum systems,a continuous-measurement-based method was proposed toidentify unknown parameters in a cavity-atomic system [13]and unknown structures of a spin network [14]. Also, asystem-realization-theory-based method in [11] was extendedto estimate unknown parameters in the Hamiltonian of aMarkovian spin network from measurement time traces [15],whose identifiability was discussed in [16]. In addition, anidentification method for linear Markovian quantum systemswas systematically discussed in [17].In contrast to the Dirac correlation function of quantumwhite noise, correlation functions of quantum colored noisecan be more complicated due to the memory effect of non-Markovian environments. A quantum system disturbed byquantum colored noise is referred to as a non-Markovianquantum system, whose dynamics is quite different from thatof the Markovian quantum systems [12]. To control non-Markovian quantum systems, it is crucial to acquire the knowl-edge of the correlation functions of the non-Markovian envi-ronment beforehand. Wu et al. proposed a frequency domainapproach to the identification of the environment spectrumfor a superconducting single qubit at an optimal point [8].For the quantum dot system in [6], an augmented systemmodel was presented to explore the structure of the non-Markovian environment [18]. Also, a spectroscopic methodwas presented to explore the spectrum of spin baths boththeoretically [19] and experimentally [20], which involveshigh-energy dynamical decoupling pulses.In this paper, we aim to extract the correlation functions ofnon-Markovian environments for spin chains. The dynamicsof the spin chains are described by a time-convolution-lessmaster equation in which time-varying damping rate functionscharacterize the correlation functions of the non-Markovian a r X i v : . [ qu a n t - ph ] O c t (cid:52)(cid:88)(cid:68)(cid:81)(cid:87)(cid:88)(cid:80)(cid:3)(cid:70)(cid:82)(cid:79)(cid:82)(cid:85)(cid:72)(cid:71)(cid:3)(cid:81)(cid:82)(cid:76)(cid:86)(cid:72) (cid:38)(cid:82)(cid:88)(cid:83)(cid:79)(cid:76)(cid:81)(cid:74)(cid:3)(cid:86)(cid:87)(cid:85)(cid:72)(cid:81)(cid:74)(cid:87)(cid:75)(cid:47)(cid:82)(cid:70)(cid:68)(cid:79)(cid:3)(cid:82)(cid:69)(cid:86)(cid:72)(cid:85)(cid:89)(cid:68)(cid:69)(cid:79)(cid:72)(cid:15)(cid:3)(cid:72)(cid:17)(cid:74)(cid:17)(cid:15)(cid:3) (cid:2252) (cid:2235)(cid:2208) (cid:54)(cid:76)(cid:81)(cid:74)(cid:79)(cid:72)(cid:3)(cid:86)(cid:83)(cid:76)(cid:81) Fig. 1. Schematic plot of spin chains in a non-Markovian environment. environment. To identify the damping rate function, we mea-sure the time trace for one local observable of the spin chainswhich correspond to dynamical variables in a coherence vectorrepresentation. Although this formulation for the coherencevector is similar to that in [15], the dynamical equation forthe coherence vector is time-varying due to the damping ratefunctions. Thus, the method in [15] is invalid in our case.We alternatively formulate the identification problem for thedamping rate function as an optimization problem. Also, wedesign a gradient algorithm to iteratively recover the dampingrate functions for the non-Markovian environment.Our contribution in this paper is that we provide a sys-tematic approach to acquire unknown damping rate functionsin a time-convolution-less master equation for non-Markovianquantum systems. In principle, our gradient algorithm canachieve the damping rate function with a high fidelity, whichwould help to obtain an exact model of the spin chainsin an experiment for quantum information processing. Inaddition, our method only requires to measure commutativeobservables of the spin chains, which is easy to be applied inan experiment. Because it avoids measuring non-commutativeobservables which cannot be accurately measured at the sametime according to the uncertainty principle in quantum me-chanics.This paper is organized as follows. In Section II, weintroduce the time-convolution-less master equation for non-Markovian spin chains whose corresponding coherence vectorrepresentation is also introduced. The identification problemformulation and the gradient algorithm are given in Section III.An example for two qubits in a non-Markovian environmentis given in Section IV. Conclusions are drawn in Section V.II. D
ESCRIPTION OF NON -M ARKOVIAN SPIN CHAINS
A. Time-convolution-less master equation
The non-Markovian spin chain we consider in this paper isplotted in Fig.1 where the nodes and the edges represent spincomponents and their couplings, respectively. The componentsare disturbed by quantum colored noise arising from anunknown non-Markovian environment. Here, the time tracesof local observables are measured and the resulting data is tobe processed for identification.In this paper, the non-Markovian dynamics of spin chainsare described by a time-convolution-less master equation [12] ˙ ρ ( t ) = L ρ ( t ) + L γ ( t ) ρ ( t ) , (1) of the spin chains. The superoperator L ρ ( t ) = − i [ H, ρ ( t )] (2)represents the internal dynamics of the spin chains where theHamiltonian of the spin chains H describes the internal energyof the system and their couplings. The commutator [ · , · ] iscalculated as [ A, B ] = AB − BA for two matrices A and B with suitable dimensions. We let (cid:126) = 1 hereafter.The Lindblad dissipative term induced by the non-Markovian environment is expressed as L γ ( t ) ρ ( t ) = 12 N − (cid:88) j,k =1 γ jk ( t )([ L j , ρ ( t ) L † k ] + [ L j ρ ( t ) , L † k ]) . (3)In contrast to the constant damping rates in Markovian masterequations, the damping rate functions γ jk ( t ) for characterizingthe correlation of the non-Markovian environment are time-varying, which encapsule the coupling strengthes betweenthe system and the environment and the density state of theenvironment [12]. Here, we assume that γ jk ( t ) are complexfunctions, i.e., γ jk ( t ) = Re( γ jk ( t )) + i Im( γ jk ( t )) , where Re( · ) and Im( · ) represent real and imaginary parts of afunction, respectively.In addition, the coupling operators L k belong to a set M = { L k , k = 1 , · · · , M = N − } , which are an orthonormalbasis for the Lie algebra su ( N ) . Their commutation and anti-commutation relations can be calculated as [ L j , L k ] = i M (cid:88) l =1 C jkl L l , (4) { L j , L k } = 2 N δ jk I + M (cid:88) l =1 D jkl L l , (5)respectively. The coefficients C jkl , D jkl ∈ R are the com-pletely antisymmetric and symmetric structure constants of theLie algebra su ( N ) , i.e., with respect to interchange of any pairof indices [21]. Here, the anti-commutator {· , ·} is calculatedas { A, B } = AB + BA for two matrices A and B with suitabledimensions and δ is the Kronecker delta function.Note that although the variation of the density matrix inthe time-convolution-less master equation (1) only dependson the current density matrix mathematically, the time-varyingdamping rate functions enable (1) to describe a non-Markovianbehavior physically; i.e., the energy exchanges between thesystem and the environment [12]. B. Dynamical equations in a coherence vector representation
Under the orthonormal basis in M , we can expand thedensity matrix ρ ( t ) and the Hamiltonian H as ρ ( t ) = 1 N I N + M (cid:88) n =1 x n ( t ) L n ,H = M (cid:88) m =1 h m L m , (6)respectively, where h m = tr( HL m ) ∈ R and the coefficients x n ( t ) are the expectation values for the basis L n ; i.e., x n ( t ) = tr( ρ ( t ) L n ) ∈ R . They constitute the so-called coherencevector x ( t ) = [ x ( t ) , · · · , x n ( t ) , · · · , x M ( t )] T ∈ R M . Here, I N is an N × N identity matrix.By using the relations (4), (5), and (6), the time-convolution-less master equation (1) can be rewritten as a dynamicalequation in the coherence vector representation as ˙ x ( t ) = A ( t ) x ( t ) + b ( t ) , x n (0) = tr( L n ρ (0)) , (7)where the initial value of the coherence vector x (0) is determined by the initial density matrix ρ (0) .The elements of the matrices A ( t ) ∈ R M × M and b ( t ) ∈ R M can be calculated as A np ( t ) = Q np + R np ( t ) with Q np = (cid:80) Mm =1 C mnp h m , R np ( t ) = − (cid:80) Ml,j,k =1 ,l ≤ m (2 − δ jk )Re( γ jk ( t )) × ( C jln C klp + C kln C jlp )+ (cid:80) Ml,j,k =1 ,l ≤ m Im( γ jk ( t )) × ( C kln D jlp − C jln D klp ) and b n ( t ) = − N (cid:80) Mj,k =1 ,j The equation (7) describes the full dynamics of the non-Markovian spin chains. However, it is possible that not all thecomponents in x ( t ) are relevant to the observables (8), i.e.,the matrix A ( t ) would have a block-diagonalizable structureand a proper sub-block corresponds to the observables. Wecan use a filtration procedure in [22] to find the accessible set of the which is related to the observables [15].Thus, a reduced coherence vector ˜ x ( t ) corresponding tothe observables (8) can be defined and it obeys a reduceddynamical equation ˙˜ x ( t ) = ˜ A ( t )˜ x ( t ) + ˜ b ( t ) , y ( t ) = ˜ c ˜ x ( t ) , (10)where ˜ A ( t ) , ˜ b ( t ) , and ˜ c with suitable dimensions are sub-matrices of A ( t ) , b ( t ) , and c .III. G RADIENT ALGORITHM FOR IDENTIFYING THEDAMPING RATE FUNCTION A. Optimization formulation for the identification of the damp-ing rate function The identification problem considered in this paper can bestated as follows. Considering the non-Markovian spin chainsgoverned by the dynamical equation of the coherence vector(7), where all the information of the environment are encodedin the unknown damping rate functions, the environment iden-tification problem is to identify the time-varying damping ratefunctions from the data of the time traces of the observables ˆ y ( t ) (9).Note that due to the interactions between the spin chainsand the environment, in general, there is no decoherence-freesubspace for a common non-Markovian system without controlpulses; i.e., decoherence channels affect all the components inthe coherence vector. Thus, the damping rate functions will beimprinted in the observables. Hence, the identification problemfor the damping rate function should be solvable. In addition,since the time-varying damping rate functions γ jk ( t ) result ina linear time-varying system (7), the identification method forMarkovian quantum systems in [15] cannot be applied in ourcase.Generally, it is difficult to identify analytically a time-varying function in a dynamical system. Alternatively, we cannumerically solve this problem by designing an algorithm. Weconsider that the system (7) evolves in a total time T duringwhich the observables are measured K times at equal timeintervals ∆ t , i.e., ∆ t = TK . Thus, a set of real measurementresults ˆ y = [ˆ y (0) , ˆ y (1) , · · · , ˆ y ( K − can be obtained.Also, the time-varying damping rate functions γ jk ( t ) arediscretized and we assume the γ jk ( t ) in each time interval areconstants γ jk ( κ ) , κ = 0 , , , · · · , ( K − . Thus in each timeinterval, we can consider the linear time varying (LTV) systemequation (7) as a linear time-invariant (LTI) system and writethe dynamical equation in a continuous-time form as ˙ x ( τ ) = A ( γ jk ( τ )) x ( τ ) + b (Im( γ jk ( τ ))) , (11) y ( τ ) = cx ( τ ) , τ ∈ [ κ ∆ t, ( κ + 1)∆ t ] (12) where Eq. (11) can be solved as x ( κ + 1) = e A ( γ jk ( κ ))∆ t x ( κ )+ (cid:90) ( κ +1)∆ tκ ∆ t e A ( γ jk ( κ ))(( κ +1)∆ t − τ ) b (Im( γ jk ( κ )))d τ. (13)Here, we have written the initial coherence vector in thecurrent and next time intervals in an abbreviate form as x ( κ ) and x ( κ + 1) , respectively.For a given set of γ jk ( κ ) , κ = 0 , , , · · · , ( K − , anoutput vector y = [ y (0) , y (1) , · · · , y ( K − can be gener-ated by using (12) and (13). A guessed set of γ jk ( κ ) , κ =0 , , , · · · , ( K − may not be identical to the real dampingrate function. Hence, the generated output y will be differentfrom the real data ˆ y . So we define an objective function J = 12 K − (cid:88) κ =0 ( y ( κ ) − ˆ y ( κ )) (14)to evaluate the distance between the two vectors y and ˆ y .Therefore, the identification problem for the damping ratefunctions γ jk ( t ) can be converted into an optimization problemas follows.For real measurement results ˆ y , the optimization problemis to find a set of γ jk ( κ ) , κ = 0 , , , · · · , ( K − such thatthe objective function (14) is minimized, that is min γ jk ( κ ) J s . t . (7) , (9) . (15) B. Gradient algorithm for the identification problem To design a gradient algorithm for the optimization problem(15), it is crucial to obtain the gradient of the objective J withrespect to γ jk ( κ ) in each time interval. By using the chain rule,the gradient can be calculated as d J d γ jk ( κ ) = d J d y · d y d x · d x d γ jk ( κ )= K − (cid:88) κ (cid:48) =0 ( y ( κ (cid:48) ) − ˆ y ( κ (cid:48) )) · c d x ( κ (cid:48) )d γ jk ( κ ) (16)where d x ( κ (cid:48) )d γ jk ( κ ) = κ (cid:48) < κ + 1 d x ( κ +1)d γ jk ( κ ) κ (cid:48) = κ + 1 . (cid:81) κ (cid:48) − κ (cid:48)(cid:48) = κ +1 e A ( γ jk ( κ (cid:48)(cid:48) ))∆ t d x ( κ +1)d γ jk ( κ ) κ (cid:48) > κ + 1 (17)In the following, we will calculate the derivative of x withrespect to the real and imaginary parts of γ jk ( κ ) , respectively.By using a standard formula for derivative of a matrix expo-nential [23] dd x e ( A + xB ) t | x =0 = e At (cid:90) t e Aτ Be − Aτ d τ (18)with matrices A and B of suitable dimensions, we have d e A ( γ jk ( κ ))∆ t dRe( γ jk ( κ )) = ∆ te A ( γ jk ( κ ))∆ t ˜E R , d e A ( γ jk ( κ ))∆ t dIm( γ jk ( κ )) = ∆ te A ( γ jk ( κ ))∆ t ˜E I , (19) where ˜E R = (cid:82) ( κ +1)∆ tκ ∆ t e A ( γ jk ( κ )) τ E R e − A ( γ jk ( κ )) τ d τ and ˜E I = (cid:82) ( κ +1)∆ tκ ∆ t e A ( γ jk ( κ )) τ E I e − A ( γ jk ( κ )) τ d τ . When the timeinterval ∆ t is small enough such that ∆ t (cid:28) || A || − , we have ˜E R = E R , ˜E I = E I , (20)where the elements of E R and E I can be expressed as E Rnp = − (2 − δ jk ) (cid:80) Ml =1 ( C jln C klp + C kln C jlp ) and E Inp = (cid:80) Ml =1 ( C kln D jlp − C jln D klp ) .Thus, the derivative of x ( κ + 1) with respect to the real andimaginary parts of γ jk ( κ ) can be calculated as d x ( κ + 1)dRe( γ jk ( κ )) = ∆ te A ( γ jk ( κ ))∆ t E R x ( κ )+ (cid:90) ( κ +1)∆ tκ ∆ t (( κ + 1)∆ t − τ ) e A ( γ jk ( κ ))(( κ +1)∆ t − τ ) × E R b (Im γ jk ( κ ))d τ, (21) d x ( κ + 1)dIm( γ jk ( κ )) = ∆ te A ( γ jk ( κ ))∆ t E I x ( κ )+ (cid:90) ( κ +1)∆ tκ ∆ t (( κ + 1)∆ t − τ ) e A ( γ jk ( κ ))(( κ +1)∆ t − τ ) × E I b (Im γ jk ( κ ))d τ + (cid:90) ( κ +1)∆ tκ ∆ t e A ( γ jk ( κ ))(( κ +1)∆ t − τ ) F d τ, (22)where the n th row element of the column vector F is − N C jkn .By combining (16), (17), (21), and (22), we obtain thegradient of the objective J with respect to γ jk . Therefore,if we update γ jk as Re( γ jk ) → Re( γ jk ) − (cid:15) R · d J dRe( γ jk ) , Im( γ jk ) → Im( γ jk ) − (cid:15) I · d J dIm( γ jk ) , (23)with small step sizes (cid:15) R and (cid:15) I , we can minimize the objective J . With the relation (23), our gradient algorithm for theidentification problem can be summarized as follows. Step 1 . Choose and measure a set of observables ˆ y ( t ) ;Initialize the state of the system x (0) , the output y (0) , andthe step sizes (cid:15) R and (cid:15) I ; Guess the initial values of { γ jk ( κ ) } ; Step 2 . Calculate x (1) to x ( K − and y (1) to y ( K − with their initial values x (0) and y (0) ; Step 3 . Compute the gradient of the objective J with respectto γ jk ( κ ) ; Step 4 . Update γ jk ( κ ) by using the relation (23); Step 5 . When a termination condition is satisfied, stop thealgorithm; Otherwise, go to the Step 2 and start a new iteration.Note that this algorithm searches optimal solutions ac-cording to the gradient and hence it may stop at a localminimum point. Moreover, the final performance of the al-gorithm depends on the initial guessed solution { γ jk ( κ ) } .When the initial guessed solution is closer to the optimalone, the algorithm will converge faster. The step size mayalso affect the performance of the algorithm. A small stepsize would result in a slow convergence process and a largeone would lead to an algorithm that saturates in the vicinity of a minimum. An adaptive step size can avoid the aboveproblems. In addition, a proper termination condition canbe the completion of a given number of iterations or theattainment of an accuracy of the objective J .IV. A P HYSICAL E XAMPLE In the example, we consider two coupled qubits immersedin an unknown common environment where the dissipativeprocesses for the two qubits can be considered to be thesame. Its non-Markovian dynamics can be described by a time-convolution-less master equation [12] as ˙ ρ q ( t ) = − i [ H q , ρ q ( t )] + γ ( t )2 (cid:88) l =1 ([ σ − l ρ q ( t ) , σ + l ] + [ σ − l , ρ q ( t ) σ + l ]) , (24)where ρ q ( t ) is the density operator for the two qubits andthe ladder operators are defined as σ + = ( σ x + iσ y ) and σ − = ( σ x − iσ y ) with Pauli matrices σ x = (cid:34) (cid:35) , σ y = (cid:34) − ii (cid:35) , σ z = (cid:34) − (cid:35) . (25)Here, the label l is used to index the qubits. The two qubitsare assumed to be coupled in an XY interaction form and thusthe Hamiltonian of the two qubits is written as H q = (cid:88) α =1 ω α σ zα + g σ x σ x + σ y σ y ) , (26)where g is the coupling strength between two qubits and ω α is the angular frequency for the α -th qubit. In thistime-convolution-less master equation (24), all the propertiesof the unknown environment are combined in the dampingrate function γ ( t ) which are the same with respect to eachdissipative channel for the two qubits. Thus the identificationtask is to recover the unknown damping rate function γ ( t ) .To simulate a real damping rate function induced by the non-Markovian environment, we assume that the real damping ratefunction is ˆ γ ( t ) = 2 hλ sinh( dt/ d cosh( dt/ 2) + λ sinh( dt/ , (27)which results from a non-Markovian environment with aLorentzian spectrum [12]. The values of the parameters h , λ , and d will be given in the following paragraph. This realdamping rate function (27) is utilized to generate the realmeasurement results ˆ y ( t ) .Furthermore, we measure the observable σ z forthe first qubit which induces the accessible set { σ z , σ z , σ x σ x , σ x σ y , σ y σ x , σ y σ y } . Due to the interactionsbetween the two qubits, the accessible set alsoincludes the operators for the second qubit. Hence,the reduced dynamical equation for the state vector Iterations 10000 10010(b) Fig. 2. Evolution of the identified damping rate function at specific iterations. ˜ x ( t ) = [ (cid:104) σ z (cid:105) , (cid:104) σ z (cid:105) , (cid:104) σ x σ x (cid:105) , (cid:104) σ x σ y (cid:105) , (cid:104) σ y σ x (cid:105) , (cid:104) σ y σ y (cid:105) ] T is inthe form of (10) with matrices ˜ A ( t ) = − γ ( t ) 0 0 − g g − γ ( t ) 0 g − g 00 0 − γ ( t ) − ω − ω g − g ω − γ ( t ) 0 − ω − g g ω − γ ( t ) − ω ω ω − γ ( t ) , ˜ b ( t ) = (cid:104) − γ ( t ) − γ ( t ) 0 0 0 0 (cid:105) T . In addition, the parameters of the system are chosen asfollows. The angular frequencies of the two qubits are thesame; i.e., ω = ω = 1 . . The coupling strength betweenthe two qubits is g = 1GHz . To simulate the real γ ( t ) , theparameters in (27) are chosen as h = 0 . , d = 0 . and λ = 0 . . The gradient algorithm starts with a guesseddamping rate function γ ( t ) = 0 . 04 cos(0 . t ) + 0 . , (28)which is plotted as a red line at the first iteration in Fig. 2(a).Both the step size in each iteration and the terminationcondition can affect the performance of the gradient algorithm.We choose the step size for updating the damping rate functionas (cid:15) = 0 . . Because a large step size would result inoscillations in the final stage for searching the minimum J anda small one would lead to a slow convergence process. On theother hand, we choose to iterate the algorithm for 20000 timesas the terminal condition. Sufficient iteration times allow thealgorithm to obtain a minimal objective J .With these parameters, the evolution of the damping ratefunction from the initial guess to the identified one are given in -6 -4 -2 Fig. 3. Reduction of the objective function J . Fig. 4. Comparison between the identified and real damping rate functions. Fig. 2. Initially, the gradient algorithm can significantly updatethe identified damping rate function as shown in Fig. 2(a).After iterations, the convergence rate of the identificationprocess slows down as shown in Fig. 2(b). This phenomenonalso reflects in the reduction of the objective function J inFig. 3. With 20000 iterations, the objective function is downto about − and we terminate the algorithm.The final identified damping rate function is plotted as thedashed red line in Fig. 4. Compared to the real dampingrate function represented as the solid blue line, our algorithmcan identify the damping rate function well except at thefinal stage. This is because the measurement results at a finaltime interval contain limited information of the damping ratefunction and thus the resulting gradient provides very smallupdates. We also make a comparison between our methodand a differential approach in [24] whose identified result isplotted as the dashed dot line in Fig. 4. Under the identicalmeasurement result, our method obtains a better result thanthat by using the differential approach. This is because allthe corresponding measured results are utilized to identify theunknown damping rate function in a sample time interval byusing our method. However, the differential approach onlyuse the measurement results at the current and next timesto estimate the damping rate function at the current time.Its identification result can be improved if we measure theobservable densely. V. C ONCLUSIONS In this paper, we have designed a gradient algorithm foridentifying the damping rate function in the time-convolution- less master equation for non-Markovian spin chains. We haveformulated the identification problem as an optimization prob-lem that can be solved iteratively by calculating the gradient ofthe objective with respect to the damping rate function in eachtime interval. The numerical example on a non-Markoviantwo-qubit system demonstrates the efficacy of our gradientalgorithm. Compared to a differential approach, our methodcan identify the damping rate function with a high fidelity.R EFERENCES[1] 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 andstabilizes photon number states,” Nature , vol. 477, no. 7362, pp. 73–77,2011.[2] J. Zhang, R. B. Wu, C. W. Li, and T. J. Tarn, “Protecting coherenceand entanglement by quantum feedback controls,” IEEE Trans. Automa.Control , vol. 55, no. 3, pp. 619–633, 2010.[3] S. Xue, M. R. Hush, and I. R. Petersen, “Feedback tracking control ofnon-markovian quantum systems,” IEEE Trans. on Contr. Syst. Technol. ,vol. 25, no. 5, pp. 1552–1563, 2017.[4] S. Xue, R. B. Wu, W. M. Zhang, J. Zhang, C. W. Li, and T. J.Tarn, “Decoherence suppression via non-Markovian coherent feedbackcontrol,” Phys. Rev. A , vol. 86, p. 052304, 2012.[5] S. Xue, R. Wu, M. R. Hush, and T.-J. Tarn, “Non-markovian coherentfeedback control of quantum dot systems,” Quantum Sci. Technol. ,vol. 2, no. 1, p. 014002, 2017.[6] 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.[7] S. Xue, M. R. James, V. Ugrinovskii, and I. R. Petersen, “LQGfeedback control of a class of linear non-markovian quantum systems,”in American Control Conference , Boston, 2016, pp. 4775–4778.[8] 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 inquantum systems,” Phys. Rev. A , vol. 87, p. 022324, 2013.[9] J. M. Geremia and H. Rabitz, “Optimal identification of hamiltonianinformation by closed-loop laser control of quantum systems,” Phys.Rev. Lett. , vol. 89, p. 263902, 2002.[10] D. Burgarth and K. Yuasa, “Quantum system identification,” Phys. Rev.Lett. , vol. 108, p. 080502, 2012.[11] J. Zhang and M. Sarovar, “Quantum hamiltonian identification frommeasurement time traces,” Phys. Rev. Lett. , vol. 113, p. 080401, 2014.[12] H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems .Oxford: Oxford University Press, 2002.[13] H. Mabuchi, “Dynamical identification of open quantum systems,” Quantum Semiclass. Opt. , vol. 8, pp. 1103–1108, 1996.[14] Y. Kato and N. Yamamoto, “Structure identification and state initializa-tion of spin networks with limited access,” New J. Phys. , vol. 16, no. 2,p. 023024, 2014.[15] J. Zhang and M. Sarovar, “Identification of open quantum systems fromobservable time traces,” Phys. Rev. A , vol. 91, p. 052121, 2015.[16] A. Sone and P. Cappellaro, “Hamiltonian identifiability assisted by asingle-probe measurement,” Phys. Rev. A , vol. 95, p. 022335, 2017.[17] M. Guta and N. Yamamoto, “System identification for passive linearquantum systems,” IEEE Trans. Automa. Control , vol. 61, no. 4, pp.921–936, 2016.[18] S. Xue, T. Nguyen, M. R. James, A. Shabani, V. Ugrinovskii, andI. R. Petersen, “Modelling and filtering for non-markovian quantumsystems,” arXiv:1704.00986 , 2017.[19] G. A. Paz-Silva, L. M. Norris, and L. Viola, “Multiqubit spectroscopyof gaussian quantum noise,” Phys. Rev. A , vol. 95, p. 022121, 2017.[20] F. K. Malinowski, F. Martins, L. Cywi´nski, M. S. Rudner, P. D.Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, andF. Kuemmeth, “Spectrum of the nuclear environment for gaas spinqubits,” Phys. Rev. Lett. , vol. 118, p. 177702, 2017.[21] R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applica-tions , L. N. in Physics Vol. 717, Ed. Springer, Berlin, 2007.[22] S. Sastry, Nonlinear Systems . Springer, New York, 1999.[23] I. Najfeld and T. F. Havel, “Derivatives of the matrix exponential andtheir computation,” Adv. Appl. Math. , vol. 16, no. 3, pp. 321–375, 1995. [24] B. Bellomo, A. De Pasquale, G. Gualdi, and U. Marzolino, “Reconstruc-tion of time-dependent coefficients: A check of approximation schemesfor non-markovian convolutionless dissipative generators,”