Simulation of Quantum Dynamics Based on the Quantum Stochastic Differential Equation

Abstract

The quantum stochastic differential equation derived from the Lindblad form quantum master equation is investigated. The general formulation in terms of environment operators representing the quantum state diffusion is given. The numerical simulation algorithm of stochastic process of direct photodetection of a driven two-level system for the predictions of the dynamical behavior is proposed. The effectiveness and superiority of the algorithm are verified by the performance analysis of the accuracy and the computational cost in comparison with the classical Runge-Kutta algorithm.

1. Introduction

Since Nelson successfully described the kinematics law of the quantum fluctuations by the Itô equation [1] and the Schrödinger equation was derived from Newtonian mechanics; the stochastic interpretation of quantum mechanics was established, in which a diffusion process was used to analyze the quantum fluctuation instead of the wave function. Then, the stochastic mechanics has gradually drawn much attention with research fields ranging from atomic and optical physics to condensed matter physics and quantum information science [2]. It becomes clear that a deep understanding of the effects of environments on a quantum system such as the mechanisms of decoherence and the dynamics of entanglement in the framework of quantum open systems is both of fundamental interest in quantum foundation issues and of practical importance in quantum information sciences. Many scholars have made thorough research on the quantum diffusion movement based on the basic theory and achieved fruitful results. For example, the normative structure of the dynamics equation of the particle fluctuation and its stability analysis methods were determined, followed by the unified interpretation of the Brown motion and the basic equation of the quantum mechanics [3–5]. The quantum stochastic dynamics elaborated the organic link between the microscopic behavior and macroevolution of the system, by which the details of the system evolution from any initial state to the final state can be analyzed. It has been successfully applied to the proliferation of microscopic particles, the molecular motors, the quantum chaos, and so forth [6].

Though the deterministic differential equations of quantum stochastic mechanics are relatively complex, it can be seen that with the development of the computer technology [7], the evolution of a microsystem can be analyzed using the numerical simulation method [8–10]. The quantum trajectory method is a typical one. It can be used for a wide range of open quantum systems to solve the master equation by unraveling the density operator evolution into individual stochastic trajectories in Hilbert space [11]. Over the last twenty years the theory of quantum trajectories has been developed by many researchers for a variety of purposes, including modeling continuously monitored open systems, improving numerical calculation and investigating the problem of quantum measurement [12, 13].

In this paper, based on the quantum stochastic dynamics, the master equation describing the time evolution law of the quantum state and its reduced density operator are investigated, and the effect of nonunitary operators on the evolution of the system is analyzed. Then the quantum stochastic differential equation is established to describe the microkinetic characteristics of the system, and a numerical iterative algorithm for the simulation of the system evolution is proposed. The practicality and advantages of the algorithm are verified by comparison with the classical Runge-Kutta numerical iterative algorithm, which is followed by further discussions on the convergence of the algorithm.

2. Methods

The quantum state diffusion theory replaces the deterministic evolution of the density operator 𝜌 representing an ensemble of open systems [14] ̇𝑖𝜌=−ℏ𝐻,𝜌+_𝑗𝐿_𝑗𝜌𝐿⁺_𝑗−12𝐿⁺_𝑗𝐿_𝑗1𝜌−2𝜌𝐿⁺_𝑗𝐿_𝑗 (1) by a unique stochastic diffusion of a quantum state, representing an individual system of the ensemble in interaction with its environment. 𝐻 is the Hamiltonian, and 𝐿_𝑗 are a set of environment operators which represent the collective effects of interaction with the environment.

However, for some complicated systems, it can be very difficult to get either the analytical or the numerical solution of (1). In that case, it is often advantageous to take alternative ways considering an unraveling of the master equation into individual quantum trajectories. Quantum state diffusion (QSD) is one of these unraveling techniques. The corresponding quantum state diffusion equation is a stochastic differential equation for the normalized state vector |𝜓⟩ representing the pure state of the system that evolves according to the QSD equation. The differerential Itô form is [15] 𝑖d𝜓=−ℏ𝐻𝜓+_𝑗𝐿⁺_𝑗𝐿_𝑗−12𝐿⁺_𝑗𝐿_𝑗−12𝐿⁺_𝑗𝐿_𝑗𝜓d𝑡+_𝑗𝐿_𝑗−𝐿_𝑗𝜓d𝜉_𝑗, (2) where ⟨𝐿_𝑗⟩ are defined by 𝐿_𝑗=𝜓𝐿_𝑗𝜓. (3)

In the QSD equation, the standard normalized terms d𝜉_𝑗 represent independent complex Wiener processes and satisfy the relations 𝐄d𝜉_𝑗=0,𝐄d𝜉_𝑗d𝜉_𝑗′𝐄=0,d𝜉^∗_𝑗d𝜉_𝑗′=𝛿_𝑗𝑗′d𝑡, (4) where 𝐄(⋅) denotes an ensemble average of the noise. QSD reproduces the master equation in the mean 𝐄𝜓𝜓=𝜌. (5) That is to say, the reduced state of the system, 𝜌, is obtained as an ensemble average. And this is what is meant by an unraveling of the master equation. Expectation values for operators obey a similar relationship: 𝑂_𝜌=Tr𝑂𝜌=𝐄𝑂_𝜓. (6) The use of QSD as a practical algorithm to solve master equations has been widely investigated [16, 17]. This includes calculations of output spectra in quantum optics [18]. As a practical method of computation, QSD gains over the direct solution of the master equation, because of a basis of 𝑁 states, QSD needs a computer store with 𝑁 elements, and the time of computation is also proportional to 𝑁. For the direct solution these are proportional to 𝑁².

In this paper, we take advantage of the system simulation method to simulate the evolutionary behavior of the open quantum systems and thus calculate and analyze various physical properties of the ensemble of open quantum systems. However, noting that we cannot usually get the analytical solution of (2), an alternative way is to find the numerical solution of the system evolution and investigate various control algorithms and control strategies based on the simulation method, which is a powerful tool built on the systems science, system identification, control theory, and computer technology for the analysis and synthesis of complex systems, especially large-scale systems [19].

In the system simulation, we should pay attention to the problem that the physical description of the stochastic process is relied on the master equation in the Lindblad operator form. If we want to get the numerical solution, the Lindblad operator must be explicitly quantified. Fortunately, existing results have summarized various forms of decoherence Lindblad operator in open quantum systems for reference; a general form of (2) can be written as [20] d𝜓(𝑡)=𝐷₁𝜓(𝑡)d𝑡+𝐷₂𝜓(𝑡)d𝑊(𝑡) (7) in which d𝑊(𝑡) is a Wiener incremental process and the random items characterizing the system decoherence are 𝐷₁=𝐿⁺1𝐿−2𝐿⁺1𝐿−2𝐿⁺𝐿𝜓,𝐷⟨⟩₂=(𝐿−⟨𝐿⟩)𝜓. (8)(9)

According to the above system model, in the given interval [0,𝑡_𝑓], a sample of realization can be generated by the following algorithm [20].(1) At the initial time 𝑡=0, the initial state of the process 𝜓^𝑟(0) is determined by the initial distribution.(2)It is assumed that at time 𝑡, the normalized state 𝜓^𝑟(𝑡) is reached through a quantum jump; then we set 𝜓^𝑟(𝑡)=𝜓. (3)Determine a random waiting time 𝜏. This can be done, for example, by drawing a random number 𝜂 which is uniformly distributed over the interval [0,1] and by determining 𝜏 from the equation: 𝜂=1−𝐹𝜓,𝜏=𝑒^−𝑖𝐻_𝜏𝜓². (10) First we define the defect of the waiting time distribution 𝑞 by the identity 𝑞≡lim_𝜏→∞exp−𝑖𝐻𝜓². (11)

 For 𝜂>𝑞, a unique solution can be obtained. If 𝜂≤𝑞, we set 𝜏=∞ in which case there will be no further jumps. Within the time interval [𝑡,𝑡+𝜏] the realization follows the deterministic time evolution: 𝜓^𝑟()=𝑒𝑡+𝑠^−i𝐻𝑠𝜓𝑒^−i𝐻𝑠𝜓,0≤𝑠≤𝜏. (12)(4) At time 𝑡+𝜏 (if 𝜏 is finite and 𝑡+𝜏<𝑡_𝑓), one of the possible jumps labeled by the index 𝑖 occurs according to (𝑡)id𝜓=−i𝐻+2_𝑖𝛾_𝑖𝐿_𝑖𝜓(𝑡)²𝜓(𝑡)+𝑑𝑡_𝑖𝐿_𝑖𝜓(𝑡)𝐿_𝑖𝜓(𝑡(𝑡)𝜓)−𝜓d𝑁_𝑖(𝑡). (13) Then we select a specific jump of type 𝑖 with probability 𝑝_𝑖=𝛾_𝑖𝐿_𝑖𝜓^𝑟(𝑡+𝜏)²∑_𝑖𝛾_𝑖𝐿_𝑖𝜓^𝑟(𝑡+𝜏)². (14) And then we update the state of system as 𝜓^𝑟()=𝐿𝑡+𝜏^𝑖𝜓^𝑟(𝑡+𝜏)𝐿_𝑖𝜓^𝑟(.𝑡+𝜏) (15)(5) Repeat steps (1) to (4) until the desired final time 𝑡_𝑓 is reached, which yields the realization 𝜓^𝑟(𝑡) over the whole time interval [0,𝑡_𝑓]. Once a sample of realizations 𝜓^𝑟(𝑡),𝑟=1,2,3,…,𝑅, has been generated according to this algorithm, any statistical quantity can be estimated through an appropriate ensemble average.

According to the above framework, an iterative algorithm is described as follows: 𝜓_𝑘+1=𝜓_𝑘+12𝐷₁𝜓_𝑘+𝐷₁𝜓_𝑘+1Δ𝑡4𝐷₂𝜓⁺_𝑘+𝐷₂𝜓⁻_𝑘+2𝐷₂𝜓_𝑘Δ𝑊_𝑘+14𝐷₂𝜓⁺_𝑘−𝐷₂𝜓⁻_𝑘Δ𝑊_𝑘²−Δ𝑡Δ𝑡^−1/2, (16) where 𝜓_𝑘=𝜓_𝑘+𝐷₁𝜓_𝑘Δ𝑡+𝐷₂𝜓_𝑘Δ𝑊_𝑘,𝜓^±_𝑘=𝜓_𝑘+𝐷₁𝜓_𝑘Δ𝑡±𝐷₂𝜓_𝑘√Δ𝑡. (17)

3. Results and Discussion

We consider the process describing the direct photodetection of a driven two-level system, and the piecewise deterministic process is given by the following equation [21, 22]: d𝜓(𝑡)=−i𝐻+i𝛾2𝜎₋𝜓(𝑡)²𝜓(𝑡)+𝜎d𝑡₋𝜓(𝑡)𝜎₋𝜓(𝑡(𝑡(𝑡).)−𝜓)d𝑁 (18)

The Poisson increments d𝑁(𝑡) satisfy i𝐻=𝐻−2_𝑖𝛾_𝑖𝐋^†_𝑖𝐋_𝑖,d𝑁(𝑡)²(𝑡),𝐄[(𝑡=𝜎=d𝑁d𝑁)]₋𝜓(𝑡)²d𝑡. (19)

The corresponding stochastic Schrodinger equation takes the formd𝜓(𝑡)=−i𝐻_𝐿𝜓(𝑡)+𝛾d𝑡2𝜎₋+𝜎_+𝜓𝜎₋−𝜎₋𝜎₊−14𝜎₋+𝜎₊²_𝜓𝜓(𝑡)+√d𝑡𝛾𝜎₋−12𝜎₋+𝜎_+𝜓𝜓(𝑡)d𝑊(𝑡), (20) where 𝐻_𝐿Ω=−2𝜎₋+𝜎₊. (21)

In the numerical simulations, it is assumed that the atom is in its ground state |𝑔⟩ and the probability of finding the atom in the excited state |𝑒⟩ can be calculated as follows: 𝜌₁₁(𝑡)=𝑒𝜌(𝑡)𝑒. (22)

From a sample of realizations this probability is estimated by determining the average 𝑀_𝑡=1𝑅^𝑅_𝑟=1𝑒𝜓^𝑟(𝑡)². (23)

An appropriate estimator for the corresponding statistical errors in the determination of 𝑀_𝑡 from a finite sample of size 𝑅 is given by 𝜎²_𝑡=1𝑅(𝑅−1)^𝑅_𝑟=1𝑒𝜓^𝑟(𝑡)²−𝑀_𝑡². (24)

According [21, 22], the analytical solution of the process is 𝜌₁₁(𝑡)=Ω²𝛾²₀+2Ω²1−𝑒^{−3𝛾0𝑡/4}cos𝜇𝑡+3𝛾₀,4𝜇sin𝜇𝑡 (25) where 𝜇=Ω²−𝛾4². (26)

Hence, it is possible to compare the numerical results with the analytical results.

3.1. Simulation Results

A sample of realizations for 𝜌₁₁ calculated from (20) is as the dashed line shows in Figure 1 with the following parameters: step size Δ𝑡=0.01, Ω=0.45, and 𝛾=0.3, while the smooth line gives the analytical solution.

Figure 1 

The comparison between the numerical solution of 𝜌11 and the analytical solution direct photodetection of a driven two-level system with the parameters: Δ𝑡=0.01, Ω=0.45, and 𝛾=0.3.

In order to discuss the performance of the algorithm, we introduced the classic Runge-Kutta iterative algorithm for generating the sample of realizations as follows [23]: 𝜓_𝑘+1=𝜓_𝑘+16𝜓¹_𝑘+2𝜓²_𝑘+2𝜓³_𝑘+𝜓⁴_𝑘Δ𝑡+𝐷₂𝜓_𝑘Δ𝑊_𝑘, (27) where 𝜓¹_𝑘=𝐷₁𝜓(𝑘,𝜓_2𝑘)=𝐷₁𝜓_𝑘+12Δ𝑡𝜓¹_𝑘,𝜓³_𝑘=𝐷₁𝜓_𝑘+12Δ𝑡𝜓²_𝑘,𝜓⁴_𝑘=𝐷₁𝜓_𝑘+Δ𝑡𝜓³_𝑘. (28)

In this paper, the computer simulation platform is the Intel(R) Core (TM)2 Duo CPU E7500 @ 2.93 GHz under the Windows XP operating system with the numerical calculation software Matlab7. By selecting different simulation step sizes Δ𝑡, we can get set of results according to the corresponding Δ𝑡 with different algorithms:(1)the estimated means of |⟨𝑒|𝜓^𝑟(𝑇)⟩|² minus the exact values obtained by the analytical solution of 𝜌₁₁(𝑡) for different step sizes and methods, (|⟨𝑒|𝜓^𝑟(𝑇)⟩|²−𝜌₁₁(𝑡)), as Table 1 shows,(2)the normalized CPU time versus the absolute error for the data points, as Table 2 shows.

Table 1 

The estimated means of |⟨𝑒|𝜓^𝑟(𝑇)⟩|² minus the exact values obtained by the analytical solution of 𝜌₁₁(𝑡), that is, (|⟨𝑒|𝜓^𝑟(𝑇)⟩|²−𝜌₁₁(𝑡)), for different step sizes Δ𝑡 and methods.

Table 2 

The normalized CPU time for different step sizes Δ𝑡 and methods.

As can be seen from Table 1, when the simulation step size is increased, there is little change in the errors in the proposed algorithm, while the errors using the Runge-Kutta method increase linearly. That is, in a certain step length range, the proposed iterative algorithm can generate more accurate approximation numerical solution than the Runge-Kutta algorithm in comparison to the analytical solution.

At the same time, we can draw a conclusion from Table 2: when obtaining a more accurate numerical solution using the proposed algorithm, the computational time is on the same order of magnitude as the Runge-Kutta method consumes. And it is not difficult to find that, when the step length is gradually reduced, accompanied by improving the accuracy, the computational time of the Runge-Kutta algorithm grows faster.

Taking the above two advantages comparing to the classical Runge-Kutta algorithm, the proposed algorithm can generate a more accurate sample of realizations while paying lower computational costs, which reflects the superiority and practicality of the proposed algorithm.

3.2. Convergence Analysis

Generally speaking, the smaller the step size Δ𝑡 is, the closer the numerical solution matches the true solution, and consequently the convergence seems to take place. We denote 𝑋(𝑡_𝑘) as the true analytical solution and 𝑋_𝑘 as the numerical approximation. Noting that 𝑋(𝑡_𝑘) and 𝑋_𝑘 are random variables, we can measure the difference using 𝐄(𝑋_𝑘−𝑋(𝑡_𝑘)), where 𝐄(⋅) represents the expected value. Then a method is said to have strong order of convergence equal to 𝛽 if there exists a constant 𝑀 such that 𝐄𝑋_𝑘−𝑋𝑡_𝑘≤𝑀Δ𝑡^𝛽 (29) for any fixed 𝑡_𝑘≡𝑡₀+𝑘Δ_𝑡∈[0,𝑡_𝑓] and Δ𝑡 sufficiently small.

The strong order of convergence (29) measures the rate at which the “mean of the error” decays as Δ𝑡→0. A less demanding alternative is to measure the rate of decay of the “error of the means.” This leads to the concept of weak convergence. A method is said to have weak order of convergence equal to 𝛽 if there exists a constant 𝑀 such that for all functions 𝑝 in some class 𝐄_𝑝𝑋_𝑘−𝐄_𝑝𝑋𝑡_𝑘≤𝑀Δ𝑡^𝛽 (30) for any fixed 𝑡_𝑘≡𝑡₀+𝑘Δ_𝑡∈[0,𝑡_𝑓] and Δ𝑡 sufficiently small [24].

According the above theory, we consider a numerical algorithm to integrate (7) in the time period [𝑡₀,𝑡₀+𝑇]. Such an algorithm will generate a discrete time approximate realizations 𝜓_𝑘 for the exact process 𝜓(𝑡_𝑘) at the given time 𝑡_𝑘≡𝑡₀+𝑘Δ_𝑡, where 𝑘=0,1,…,𝑛=𝑇/Δ𝑡. Before our discussion, several conventions should be made as follows.(1)In this section, 𝜓_𝑘 always represents a numerical approximation, while 𝜓(𝑡_𝑘) stands for the exact process described by (7). (2)Similarly, we define the discrete time approximation of the density matrix as 𝜌_𝑘=𝐄[|𝜓_𝑘⟩⟨𝜓_𝑘|], whereas 𝜌(𝑡)=𝐄[|𝜓(𝑡)⟩⟨𝜓(𝑡)|] is the exact density matrix satisfying the Lindblad equation  (9).(3)We always set a deterministic initial state 𝜓₀≡𝜓(𝑡₀).

In order to illustrate the numerical convergence of the algorithm, we will compare the Taylor expansion of the exact density matrix 𝜌(𝑡) which is given by (31) with the generated approximation 𝜌_𝑘: 𝜌()(𝑡)(𝑡)1𝑡+Δ𝑡=𝜌+𝐋𝜌Δ𝑡+2𝐋²𝜌(𝑡)Δ𝑡²+𝑂Δ𝑡³. (31)

We compare the numerical simulation approximate obtained by the proposed algorithm to the true evolution 𝜌(𝑡). The single-step error of a certain numerical scheme may then be expressed through the difference 𝜌₁−𝜌𝑡₁=𝑂Δ𝑡^𝛾 (32) which means that the strategy reproduces the Taylor expansion of 𝜌(𝑡) including terms of order 𝛾−1 in Δ𝑡.

Thus, it is direct to prove that the integration over a finite time period [𝑡₀,𝑡₀+𝑇] decreases the convergence order by one since 𝑛=𝑇/Δ𝑡 time steps are needed to compute the density at time 𝑡_𝑛=𝑡₀+𝑇, that is, 𝜌_𝑛−𝜌𝑡_𝑛=𝑂Δ𝑡^𝛽 (33) with 𝛽=𝛾−1. If (33) can be satisfied, the numerical strategy is defined to be a strategy of order 𝛽 [15]. It should be noted that (33) is a special case of weak convergence of order 𝛽 describing the degree of proximity of the probability distributions of 𝜓_𝑛 and 𝜓(𝑡_𝑛) which is a much weaker criterion comparing to strong convergence defined in (29). In actual applications, one tends to care about this weaker form of convergence especially when considering the approximation of functionals of the stochastic variable.

When investigating the stochastic differential equations and numerical simulation solution process, if the issues are related to the numerical simulation of the stochastic process, the evaluation criteria of the solutions convergence are usually defined as the numerical approximation curve. It must be sufficiently close to the real trace of the evolution. That is to say, the higher the convergence order is, the closer the distribution of the numerical solution is with the analytical solution of the distribution. So it requires that the probability distribution 𝜓_𝑛 obtained by the simulation iterative algorithm is close enough to the probability distribution determined by the quantum master equation [24].

Performing the error analysis according to (16) one can find it converging with order 𝛽=2 in contrast to the classical Runge-Kutta method with 𝛽=1. Thus, it is really a higher-order strategy in the weaker convergence sense.

4. Conclusions

The decoherence of open quantum systems usually makes the system evolve from the initial pure state to mixed states (in some cases, may also be mixed state into a pure state). Being a powerful tool for investigating the open quantum systems, the quantum master equation can give a quantitative description of the transition, dissipation, and decoherence caused by the interaction between the closed system and the environment. Taking this as the starting point of our research, in order to obtain the evolution of the open quantum systems according to its dynamic characteristics, we used the system simulation method to get the numerical solution to the reduced density operator of a typical open quantum system. And its effectiveness and superiority were verified in comparison with the classical algorithm. Further research includes the control scheme [25, 26] for quantum manipulation based on the characteristics of quantum dynamics.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (61203062) and the Doctor Start Foundation of Guangdong University of Technology (405105028).