Simulation of quantum dynamics based on the quantum stochastic differential equation.

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] by a unique stochastic diffusion of a quantum state, representing an individual system of the ensemble in interaction with its environment. H is the Hamiltonian, and L 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] where 〈L 〉 are defined by

In the QSD equation, the standard normalized terms dξ represent independent complex Wiener processes and satisfy the relations where E(·) denotes an ensemble average of the noise. QSD reproduces the master equation in the mean 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: 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 N states, QSD needs a computer store with N elements, and the time of computation is also proportional to N. For the direct solution these are proportional to N 2.

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] in which dW(t) is a Wiener incremental process and the random items characterizing the system decoherence are

According to the above system model, in the given interval [0, t ], a sample of realization can be generated by the following algorithm [20].

At the initial time t = 0, the initial state of the process ψ (0) is determined by the initial distribution.

It is assumed that at time t, the normalized state ψ (t) is reached through a quantum jump; then we set .

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:

First we define the defect of the waiting time distribution q by the identity

For η > q, a unique solution can be obtained. If η ≤ q, we set τ = ∞ in which case there will be no further jumps. Within the time interval [t, t + τ] the realization follows the deterministic time evolution:

At time t + τ (if τ is finite and t + τ < t ), one of the possible jumps labeled by the index i occurs according to

Then we select a specific jump of type i with probability

And then we update the state of system as

Repeat steps (1) to (4) until the desired final time t is reached, which yields the realization ψ (t) over the whole time interval [0, t ]. Once a sample of realizations ψ (t),  r = 1,2, 3,…, R, 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: where

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]:

The Poisson increments dN(t) satisfy

The corresponding stochastic Schrodinger equation takes the form where

In the numerical simulations, it is assumed that the atom is in its ground state |g〉 and the probability of finding the atom in the excited state |e〉 can be calculated as follows:

From a sample of realizations this probability is estimated by determining the average

An appropriate estimator for the corresponding statistical errors in the determination of from a finite sample of size R is given by

According [21, 22], the analytical solution of the process is where

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

3.1. Simulation Results

A sample of realizations for ρ 11 calculated from (20) is as the dashed line shows in Figure 1 with the following parameters: step size Δt = 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: Δt = 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]: where

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 Δt, we can get set of results according to the corresponding Δt with different algorithms:

the estimated means of |〈e|ψ (T)〉|2 minus the exact values obtained by the analytical solution of ρ 11(t) for different step sizes and methods, (|〈e|ψ (T)〉|2 − ρ 11(t)), as Table 1 shows,

Table 1

The estimated means of |〈e|ψ (T)〉|2 minus the exact values obtained by the analytical solution of ρ 11(t), that is, (|〈e|ψ (T)〉|2 − ρ 11(t)), for different step sizes Δt and methods.

	0.01	0.02	0.05	0.10	0.15	0.20	0.25
The proposed algorithm	0.00012	0.00018	0.00015	0.00025	0.00021	0.00017	0.00021
Runge-Kutta	0.00026	0.00049	0.00070	0.00121	0.00142	0.00220	0.00311

the normalized CPU time versus the absolute error for the data points, as Table 2 shows.

Table 2

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

	0.01	0.02	0.05	0.10	0.15	0.20	0.25
The proposed algorithm	0.153	0.108	0.081	0.061	0.050	0.038	0.007
Runge-Kutta	1	0.162	0.100	0.083	0.072	0.065	0.030

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 Δt is, the closer the numerical solution matches the true solution, and consequently the convergence seems to take place. We denote X(t ) as the true analytical solution and X as the numerical approximation. Noting that X(t ) and X are random variables, we can measure the difference using E(X − X(t )), where E(·) represents the expected value. Then a method is said to have strong order of convergence equal to β if there exists a constant M such that for any fixed t ≡ t 0 + kΔ ∈ [0, t ] and Δt sufficiently small.

The strong order of convergence (29) measures the rate at which the “mean of the error” decays as Δt → 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 M such that for all functions p in some class for any fixed t ≡ t 0 + kΔ ∈ [0, t ] and Δt sufficiently small [24].

According the above theory, we consider a numerical algorithm to integrate (7) in the time period [t 0, t 0 + T]. Such an algorithm will generate a discrete time approximate realizations ψ for the exact process ψ(t ) at the given time t ≡ t 0 + kΔ, where k = 0,1,…, n = T/Δt. Before our discussion, several conventions should be made as follows.

In this section, ψ always represents a numerical approximation, while ψ(t ) stands for the exact process described by (7).

Similarly, we define the discrete time approximation of the density matrix as ρ = E[|ψ 〉〈ψ |], whereas ρ(t) = E[|ψ(t)〉〈ψ(t)|] is the exact density matrix satisfying the Lindblad equation  (9).

We always set a deterministic initial state ψ 0 ≡ ψ(t 0).

In order to illustrate the numerical convergence of the algorithm, we will compare the Taylor expansion of the exact density matrix ρ(t) which is given by (31) with the generated approximation ρ :

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

Thus, it is direct to prove that the integration over a finite time period [t 0, t 0 + T] decreases the convergence order by one since n = T/Δt time steps are needed to compute the density at time t = t 0 + T, that is, 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 ψ(t ) 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.