0open Perturbative approach to Markovian
SUBJECT AREAS: QUANTUM PHYSICS THEORETICAL PHYSICS
Received 21 February 2014
Accepted 15 April 2014
Published 8 May 2014
Correspondence and requests for materials should be addressed to A.C.Y.L. (andyli@u. northwestern.edu)
open quantum systems
Andy C. Y. Li1, F. Petruccione2 &Jens Koch1
department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208, USA, 2Quantum Research Group, School of Chemistry and Physics, University of KwaZulu-Natal, Durban 4001, South Africa and National Institute for Theoretical Physics (NITheP), KwaZulu-Natal, South Africa.
The exact treatment of Markovian open quantum systems, when based on numerical diagonalization of the Liouville super-operator or averaging over quantum trajectories, is severely limited by Hilbert space size. Perturbation theory, standard in the investigation of closed quantum systems, has remained much less developed for open quantum systems where a direct application to the Lindblad master equation is desirable. We present such a perturbative treatment which will be useful for an analytical understanding of open quantum systems and for numerical calculation of system observables which would otherwise be impractical.
nderstanding open quantum systems is of fundamental importance in many contexts of physics, ranging from traditional atomic and molecular physics1-3 to recent studies of circuit QED4-7 and optomechanical systems8-10. However, the simulation of open quantum systems is even more demanding than that of closed systems, and analytical and numerical approximation methods are less developed. This makes theoretical studies of large open systems particularly challenging and has spurred immense interest in the development of open-system quantum simulators11 which could also be used to study dissipative phase transitions12-14.
Markovian open quantum systems can be described by the Lindblad master equation15,
P (t) = Mi),
which governs the time evolution of the density matrix p in terms of the generalized Liouville super-operator L Equation (1) allows for the investigation of system dynamics and of the asymptotic steady state by exact diagonalization of L or by simulating quantum trajectories16. However, the exponential growth of the Hilbert space dimension N with system size and the even more dramatic growth of the dimension N2 of the associated space of linear operators severely limit the possibility of obtaining numerically exact solutions. Approaches alternative to the exact diagonalization of L include sophisticated numerics such as matrix product methods17-19 and self-consistent projection operator methods20. Although these techniques are suitable for larger systems, they are only easily applicable to one-dimensional or highly symmetric systems.
For closed quantum systems of large size, perturbation theory is routinely employed as a useful approximation. In the past, specialized perturbative methods have also been applied to open quantum systems in specific contexts including full counting statistics21,22, finite-time evolution23 and critical behavior1314. One example of an open-system perturbative treatment is adiabatic elimination (generalized Schrieffer-Wolff formalism), which gives a simplified Liouville super-operator24-27 whenever the spectrum separates into slow and fast degrees of freedom. Another method consists of a perturbative expansion of dissipative terms23,28, which can yield the approximate finite-time evolution of systems with small dissipation rates. Moreover, a perturbative treatment for open systems
can be formulated within the Keldysh Green's function framework based on bosonic or fermionic field propagators13,14,29.
Here, we present a canonical perturbative method that systematically determines the corrections to eigenstates and eigenvalues of the Liouville super-operator L. Our treatment is not limited to a certain system or system type. Rather, it is applicable to a wide range of perturbations and different open quantum systems. Perturbation theory (PT) similar to the one presented here was previously studied30,31 but results were limited to the steady state and the issue of non-positivity of density matrices due to truncation was not addressed by the authors.
The density-matrix PT we develop yields results both for the steady state as well as for all other eigenstates of L. Further, we propose and derive a new PT based on the amplitude matrix. This amplitude-matrix PT guarantees a properly positive steady-state density matrix. Both kinds of PT can be applied to large open systems with lattice
structure. Understanding the properties of such open systems is of immense interest and would otherwise be impractical due to the system size. The study of such systems is an important motivation for developing the PT presented here. Its application to open lattice systems is beyond the scope of the present paper and will be published elsewhere.
Results
Density-matrix PT. We propose a non-degenerate density-matrix PT based on the quantum master equation (1). The Liouville superoperator L, which serves as the generator of the quantum dynamical semigroup15, is generally not Hermitian, i.e. the adjoint superoperator L{ is not equal to L. The right and left eigenstates um and wm of L associated with the eigenvalue 1meC are defined by
Lum~1mUm, L{ w(—(1()(2)
Here, m is a non-negative integer labeling the eigenstates. Suitably normalized, the left and right eigenstates obey the bi-orthonormality relation (wm, un) = dmn, where (x, y) ; Tr [x{y] is the Hilbert-Schmidt inner product32 for linear operators x and y. Together, the eigenvalues 1m and the associated right eigenstates um contain all information of the steady state (labeled by m = 0) and the dynamics of the system.
In analogy to closed-system PT, the density-matrix PT is based on the series expansion of 1m and um with respect to a small parameter a set by the perturbation. Starting point is, thus, the separation of L into two parts: an unperturbed super-operator L0 and a perturbation a L1, i.e.
L —L0 + a L1. (3)
Like the original Liouville operator, L0 must still be a proper generator of the quantum dynamical semigroup. In addition, L0 should be solvable in the sense that its unperturbed spectrum {im>), wm0), um0)} is known, or, at least a subset of interest. As appropriate for non-degenerate PT, we will assume that this part of the spectrum is non-degenerate.
We next employ a general series expansion of eigenvalues and eigenstates in a,
1m=E um=£ tiuip. (4)
j=0 j=0
Here, the index j counts orders of perturbation theory and ij and u() are hence the j-th-order corrections to the eigenvalue and the right eigenstate. We determine recursion relations for ij and uj by plugging equations (3) and (4) into equation (2) and examining the result order by order in a. For the j-th order in a, we obtain the general recursive expression
(Lo-im0))umj) = -L1um'-1) +a( (5)
with aV-Yi A(f)u(j-k). m / —i m m
So far, our treatment directly mirrors the well-known derivation of
stationary PT for a closed system. Specifically, replace L0, L1, um and
1m by the unperturbed Hamiltonian H0, the perturbation V, the
eigenvectors |ym) and eigenenergies Em of H ; H0 + aV. Then,
equation (5) takes exactly the form of the usual recursion equation
in closed-system PT, namely
(ho-e(0)) |ymj)) — -v|ym'-1) )+d|j), (6)
with d«—Vj E(k) iy|j-k)> . To obtain the recursive relation for m Z_/k—1 m r m '
the eigenenergy correction Ej we multiply equation (6) with I from the left. This yields
Ei — ( ymo)|v |y«-1A EM ym°)|y<rk)). (7) k — 1
Analogously, we take the inner product of equation (5) with the left eigenstate w^-1. This yields the recursion relation for iC^,
imj)—/wmo),L1uc-1A -x; if)( wo^). (8) k— 1
Note that equation (7) is usually simplified further by demanding (ym°) | yj) — 0 for j ? 0. We will instead keep the corresponding term (wC», um'-k)) in equation (8) for reasons that will become clear momentarily.
We next turn to the computation ofthe eigenstate corrections. The Hamiltonian of any closed system is Hermitian. Hence, H0 provides a complete orthonormal eigenbasis {| ^^((0) )}. As a result, |ymj)) in equation (6) can be expanded in this eigenbasis. Solving equation (6) is then straightforward. By contrast, L is not Hermitian and may not even be diagonalizable. As a result, the expansion of umj) in terms of um°) will generally fail. We therefore adopt the different strategy of applying an appropriate generalized inverse to the singular superoperator (L0 — im0)). Several options exist for such a generalized inverse21,22,26,30; here, we choose the Moore-Penrose pseudoinverse which is well-defined for non-invertible matrices A and which we denote by A j1 . A brief review of the Moore-Penrose pseudoinverse is provided in the Supplementary Note. Applying the pseudoinverse to equation (5), we obtain
u<j) — (L0-^(-L^+Ay). (9)
Details of the derivation of equation (9) are given below in the Methods section. We emphasize that this pseudoinverse does not guarantee that (wc"-'^'/-') — 0, explaining our previous remark on keeping this term intact.
The steady-state density matrix ps ; u0 defined by Lps — 0 is of particular interest. As a special case of equations (8) and (9), we can simplify the corrections l"j) and p(j) to
l"j)—0, (10) pM — -LJ1 L1p('-(11)
see Methods for details. Corrections to 10 and ps were previously derived30 without using the Moore-Penrose pseudoinverse. The result for the density-matrix corrections in ref. 30 differs from ours [equation (11)] merely by a shift,
p«?p« + jP0), (12)
where Cj is a constant.
This shift can be interpreted as follows. Note that equation (5) does not have a unique solution since (L0 - iC"') is non-invertible. Indeed,
for any given solution psj) of equation (5) its shifted counterpart from equation (12) is also a solution. The choice of a particular generalized inverse effectively selects a specific set of shift parameters C . The difference between the result in ref. 30 and ours merely reflects the different choices of generalized inverses. Since shifts of the form of equation (12) only affect the overall normalization of ps, our result for the steady state corrections is equivalent to that given in ref. 30. We properly normalize our result by defining
ps —o ajp<j)] (13)
with N[A] : A/Tr [A] effecting normalization for any non-traceless matrix A.
Finite-order PT truncates the series in equation (13) just as in closed-system PT. Let us denote the approximate result up to M-th
order as pas'M = Op^J. We next check whether pas'M
indeed represents a proper density matrix, i.e., whether it is normalized, Hermitian and positive-semidefinite15. By virtue of N, the finite-order result pS;M is explicitly normalized. Hermiticity can be verified by noting that LJ1 and L1 map Hermitian operators to Hermitian operators. However, omission of higher-order terms in the truncation can render pj*;M non-positive. In the examples presented later on, this issue indeed occurs for certain parameter choices. Beyond mere conceptual concerns, non-positivity also prevents direct calculation of quantities such as state fidelity and von-Neumann entropy. This marks a key difference between closed-system PT and density-matrix PT. In closed-system PT, the approximate result is always a proper quantum state. By contrast, for density-matrix PT the approximate result may, strictly speaking, not be a proper density matrix.
Similar issues with non-positivity of approximate density matrices are also encountered in quantum tomography, there caused by measurement errors. In this case, a maximum-likelihood method is typically used to reconstruct a physical density matrix from the nonpositive approximation33,34. However, this method is impractical to apply to the large density matrices we are interested in. We next discuss an alternative strategy which reinstates positivity for large density matrices obtained within PT.
Amplitude-matrix PT. We propose an amplitude-matrix PT to perturbatively construct an approximate steady-state density matrix which is manifestly positive. For this, recall that any Hermitian and positive-semidefinite matrix p can be decomposed in the form35: p = ff{. Following ref. 36, we will refer to f as the amplitude matrix. The decomposition p = ff{ is not unique: there are many choices for f leading to the same matrix p. To eliminate these extra degrees of freedom, we choose f to be lower triangular with real, non-negative diagonal elements. Existence and uniqueness of f are then guaranteed by the Cholesky decomposition. (In the case that p possesses a zero eigenvalue, the Cholesky decomposition is not unique but this issue can be bypassed by known strategies37,38.)
We start from the power series expansion of the steady-state amplitude matrix in a: fs _ of-*. Here, all matrices f j are
again of lower triangular shape. Plugging this expansion into ps _ fsf{ and collecting terms of the same order in a, we obtain
„a;M bi [(-a;M/¡-a;M\t
ps;am=n f; (w )
f0)( f0)) t~rS0),
Now, pi*'am is manife s tly a proper den sity matrix: it is normalized, Hermitian and positive- semidefinite by construction. We note already that expectation values for observables obtained from amplitude-matrix PT, however, are not necessarily more accurate than those obtained from density-matrix PT, even if p*'m is slightly non-positive. The examples in the following section illustrate that the respective accuracies of density-matrix versus amplitude-matrix PT generally depend on the specific perturbation and system parameters.
Amplitude-matrix PT is more involved if one or more eigenvalues of ps0) vanish, e.g. when ps0) represents a pure state. In that case, Z0 in equation (15) is non-invertible (see Methods) and thus a unique solution for equation (15) does not exist. Depending on the specific case, there may be infinitely many solutions or no solution. In the case of infinitely many solutions, we add a small identity-matrix component to pls0), i.e., pls0) ->pls0) + ct with 0 < c<l. The identity matrix then acts as a correction matrix which stabilizes the procedure of solving the linear equation [equation (15)] to provide a unique f). We utilize this correction matrix strategy in the second example discussed below. Whenever equation (15) has no solution, other forms of series expansions would need to be applied. We will not further consider that case in the present paper. The correction-matrix method and the validity of the series expansion are further discussed in the Methods section.
Comparing PT with exact results. To illustrate the use and assess the accuracy of density-matrix and amplitude-matrix PT, we study two example systems, see Fig. 1. These two examples are neither claimed to be original nor of particular intrinsic interest. Our selection is guided by the intention to discuss two systems that are non-trivial, yet sufficiently small to still allow for a numerically exact treatment of the master equation. This way, we can compare steady-state expectation values from second-order PT with those obtained from exact diagonalization of L. We note again that the steady-state result obtained from finite-order density-matrix PT can be non-positive. This is indeed the case for some choices of parameters in the two following examples. Thus, we also apply the amplitude-matrix PT and compare results from the two perturbative treatments.
^^¿r.........^^
Qfc^........ift^
Here, the super-operator Z0 is defined via
The zero-order amplitude matrix f ^ is directly obtained from equation (14) by Cholesky decomposition and fj is determined recursively from the system of linear equations in equation (15). We determine pj in equation (15) by density-matrix PT; in this sense, the amplitude-matrix PT is based on the density-matrix PT.
Once we truncate the amplitude matrix to M-th order,
f^;M : _0 Ofj, we can compute the steady-state density matrix
ps,AM by
Figure 1 | Schematic of the examples. (a) The system consists of a ring of four coupled spins with nearest-neighbor flip-flop interaction of strength t (dashed red lines). The spins are jointly driven with a coherent tone of strength e and frequency (blue arrow). Each spin is further subjected to local dissipation with spin relaxation rate y (green curly arrows). (b) The second system is composed of a three-resonator ring (black rectangles) coupled to a single qubit with Jaynes-Cummings interaction strength g (blue dotted line). The resonators are coupled to each other through photon hopping with rate k (red dashed lines). One of the resonators is driven coherently (blue arrow). Resonators and qubit are subject to local dissipation with photon decay rate yn and qubit relaxation rate yq (green curly arrows).
<°>Г>
Figure 2 | Results for four spins coupled in a ring. Shown are: the unperturbed result with respect to L0 (blue dotted line), the exact result with respect to L (black solid line), the second-order density-matrix-PT result (red dashed line) and the second-order amplitude-matrix-PT result (green thin line).
(a) and (b): s-) | and s- ) are plotted as a function of dm/y for e/y — 0.8 and i/y — 0.4. The exact result is well approximated by the second-order PT result. The amplitude-matrix PT is slightly less accurate in this particular case which is more visible in (s+ s- ). (c) and (d): The same is plotted except for t/y— 0.8. Although the deviation is relatively large, the shape of the exact result is still qualitatively captured by the second-order PT results.
Ring of four coupled spins, driven and damped. We consider four spins coupled in a ring as shown in Fig. 1 (a). The spins are coupled by flip-flop interaction with spin-spin coupling strength t. We assume all spins are driven equally with strength e and drive frequency Within the Rotating Wave Approximation (RWA), the Hamiltonian describing this system is given by
+ ^ (s>„-+i+h.c.).
Here, s+ are the raising or lowering operators for the spin at site n, and ; ю0 _ is the detuning between the spin frequency ю0 and the drive frequency Note that in equation (17), the time dependence of the drive has already been eliminated by working in the rotating frame. All four spins are coupled to a zero-temperature bath, leading to pure spin relaxation with relaxation rate y. Thus, the Liouville super-operator L is given by
Lp=-i[H,p]+y^ D[s-] p, (18)
where D [s- ] p : s- ps z--sz s- p--ps z s- is the usual
LnJ^ n > n ^ n n > ^ n n
dissipator for spin relaxation.
As the perturbation, we now choose the spin-spin coupling and, thus, separate L into two parts, L ~ L0 +1 L1 where
L0p ~= E { - i' s„+ s- + E (s„+ + s-) ,p] + yD [s- ] p}
describes the "atomic limit'' in which the spin-spin coupling is absent, and the perturbation
t = (o„+o„ + i+o„ 0++1),p]
captures the spin-spin coupling. Orders of perturbation theory are counted with respect to t. Even though the system possesses a high
degree of symmetry, the steady states of the unperturbed and perturbed systems are unique. Hence, the zero-eigenvalues of L0 and L are non-degenerate and our non-degenerate PT for the steady state is applicable.
We can thus implement second-order PT to compute the steady-state expectation values for several operators. We choose s1{ and the excitation number s1+s1{ since they fully capture the reduced density matrix of a single spin (which is the same for all spins due to symmetry). In Fig. 2, we compare results from density-matrix PT and amplitude-matrix PT to the exact and the unperturbed results.
We first consider the case shown in Fig. 2(a) and (b) where the coupling strength t represents the smallest energy scale. Consistent with fin dings in ref. 4, we observe two symmetric resonance peaks in | (s- | for the unperturbed result (separate spins). When the coupling is present, the two peaks are shifted in position and become asymmetric, in agreement with the results from ref. 39. For s1+s1{ , we observe a similar shift and asymmetry in the resonance peak due to the spin-spin coupling. Second-order PT nicely captures the above features. Note that the amplitude-matrix-PT result is actually slightly less accurate than the density-matrix-PT result in this particular case.
To illustrate the expected limitations of PT, we next increase the coupling strength so much that it matches the drive strength. Results for this parameter choice are shown in Fig. 2(c) and (d). Qualitatively, the shape of the curves from the exact calculation is still captured by the perturbative results. However, the results from PT show significant deviations from the exact result. Large deviations like this are not surprising since the "perturbation" parameter t now matches both e and y and PT breaks down.
It is an interesting question whether a dimensionless parameter a can be identified that governs the convergence of this PT. To investigate this question, we recall that in simple cases ofclosed-system PT, a may be inferred from the expression for second-order corrections. a is then given by the ratio of the perturbation strength t and the difference between the relevant unperturbed eigenenergies; for instance, for the ground state we have:
0.4 0.2
Figure 3 | Eigenvalue spectrum of a driven-damped spin. The diagram
shows the eigenvalue spectrum for system parameters chosen as
dm/y — 0.2 and varying drive strength from e/y— 0 (red circles) to e/y— 0.25 (purple squares). The fourth eigenvalue (black diamond)
corresponds to the steady state and is always zero. As illustrated by the blue trajectories, the flow of the eigenvalues is non-trivial and the eigenvalue(s) closest to zero switches from the complex pair to the purely real eigenvalue
as e increases.
(closed system) .
For open systems, the unperturbed eigenvalues of the Hamiltonian are replaced by those of the super-operator L0. For the steady state, the relevant eigenvalue difference is that between the steady-state eigenvalue zero and that of the closest non-zero eigenvalue(s),
minm=o|im0)— loi
H~ ^ |d® яПя„+^яПя„ + 1+я„яП + ^| +£^я1+я|
+ dv s+s— +fl2S ^ ■
Here, an (a^ is the annihilation (creation) operator for photons in the resonator at site n and dm ; m — md is the detuning between the bare resonator frequency m and the drive frequency md. The qubit is assumed in resonance with the bare resonator frequency. Qubit and resonators are each coupled to a zero-temperature bath, inducing
qubit relaxation and photon decay with rates y9 and ya, respectively. The super-operator L is thus given by
Lp~—i[H,p]+ya^ D[a„]p + yqD[s—]p.
The qubit-resonator coupling mediates two effects: an indirect coherent drive on the qubit, and correlations between the resonator ring and the qubit. We wish to treat the correlation effects perturbatively. To do so, we separate the two effects by means of a coherent displacement as follows. Note that in the absence of qubit-resonator coupling, the drive places the eigenmodes m = 1,2, 3 of the
resonator ring in coherent states with amplitudes /я^
: âm where
~ e / с 2Pm У a --^ ( ою + 2к cos —---i —
exp( i ^ (24)
and am: рЦ^З я" exp ^i ~ m^j is the annihilation operator for
photons in mode We thus displace ят according to
am ~ am am ' and rewrite the Liouville super-operator as
Lp=—i[H0+gHbp]+ye£) D[äm]p+yqD[s—]p. m
Here, H0 is the unperturbed Hamiltonian
H0~ ^ cos-3^^ (am){am
+ d® s+ s— + (eeff s+ +£*ff s—)
Even for a simple system like a driven-damped spin, the spectrum of complex eigenvalues ij0 depends on the system parameters d®, e and y in a rather non-trivial way, see Fig. 3. The tuning of system parameters can even change the identity of the eigenvalue А^0 that is closest to zero. Hence, it is in general difficult to write the small parameter a in a simple form showing the dependence on the various system parameters.
Qubit coupled to a driven and damped resonator ring. For our
second example, we choose a system composed of a single qubit coupled to a three-resonator ring, see Fig. 1(b). The resonators are coupled among each other by photon hopping with rate k. The first resonator is driven with strength e and drive frequency The qubit is coupled to the second resonator with a coupling strength g. Within RWA and in the frame co-rotating with the drive, the system Hamiltonian H is given by
in which the displaced eigenmodes a( and the qubit remain decoupled. One finds that the effective drive on the qubit is given
by eeff: am exp ^ - i . The perturbation gH1 describes
the remaining coupling between the displaced eigenmodes and the qubit,
gHi~ "TfX! (exn — i~'m'
+h. c.
The perturbative treatment at the level of the master equation simply follows this separation and decomposes L —L0+gL1 into the unperturbed super-operator
Lop—-i[Hosp]+yfl X D[&m]p+yqB[s-]p, (28)
and a perturbation which only captures the remaining coupling,
gL p—-i[gH1 ,p]. (29)
The order of perturbation theory here is counted with respect to g.
The steady state of L0 is a product state composed of a density matrix for the resonator ring and one for the qubit. The resonator ring is in a pure state (all displaced eigenmodes in the vacuum state). As a consequence, the unperturbed density matrix pso) has multiple eigenvalues zero. Therefore, when implementing amplitude-matrix PT, we employ the correction matrix method mentioned above (with parameter C = 10—9).
We will focus on expectation values at site 1, specifically (a1) and (n1) : (a{a1), as a function of the drive frequency, expressed in terms of the detuning dm. In the unperturbed, i.e., decoupled case (g = 0), we expect two resonances at dm = — 2k and dm = k corresponding to the eigenmodes of the resonator ring. Once the qubit is coupled to
// u Amplitude-matrix
Unperturbed J? \\ Perturbation theory 0.003
Density-matrix**" 11 jff
perturbation theory U Jp
Figure 4 | Results for single qubit coupled to a resonator ring. The color scheme of Fig. 2 is used. To apply the amplitude-matrix PT, the correction matrix is employed with the parameter c = 10~9.Weplot |(fl1)|, (n1) and |(s2)| as a function of dm/e for k/e_10, g/e_0.5 and Cq/e_ya/e_ 0.05. The resonance at dm = 0 is well captured by second-order PT. The amplitude-matrix PT and density-matrix PT perform with nearly identical accuracy.
the resonator ring, we expect a resonance at dv = 0 originating from the qubit's response to the drive. We monitor this response by calculating the expectation value of s2.
As shown in Fig. 4, we confirm that the resonance at dm = 0, a key consequence of the coupling, is successfully captured by second-order PT. Specifically, we consider the case of drive and coupling strengths e and g small compared to the hopping rate k but large compared to the relaxation and decay rates cq and ca. Note that the perturbation parameter g is not the smallest energy scale in this case; nonetheless, PT holds. The expectation values of a1, n1 and a2 are shown in Fig. 4(a), (b) and (c) respectively. The results from second-order PT are in good agreement with the exact result. Amplitude-matrix PT and density-matrix PT here give results with nearly identical accuracy. The saturation effect visible in Fig. 4(c) shows that the qubit is in the nonlinear regime. These results also illustrate that the correction matrix method required for the amplitude-matrix PT is reliable and yields results consistent with the exact solution.
Discussion
We have detailed a perturbative approach to Markovian open quantum systems and developed a non-degenerate PT based on the Lindblad master equation. This density-matrix PT recursively determines corrections to the eigenvalues and right eigenstates of the Liouville super-operator L. As a result of finite-order truncation, such a perturbative scheme may lead to non-positive steady-state "density matrices'' which prevent direct calculation of quantities like the state fidelity. The issue of non-physical states, which does not occur in closed-system PT, can be tackled by a modified perturbative scheme based on the amplitude matrix. With two example systems, we have illustrated that the approximate results are in excellent agreement with exact results for representative parameter choices. The expectation values obtained from density-matrix PT showed good agreement in the two examples even if the truncated density matrix was slightly non-positive.
The perturbative treatment presented here is suitable for systems of sizes that cannot be handled by exact solution of the quantum master equation. An interesting future application of this PT thus consists of the study of open quantum systems with lattice structure, such as the open Jaynes-Cummings lattice. Promising experimental progress5 indicates that such open lattices can indeed be realized in the circuit QED architecture, and may serve as valuable open-system quantum simulators11. Openness and relatively large size of such systems make the theoretical investigation challenging. We believe that the developed open-system PT will provide a useful tool in exploring the physics of open lattice systems.
Methods
Perturbative corrections. We wish to prove that the expression of u^'0 in equation (9), rewritten as
»H Lo-lf) j,
is indeed a solution to equation (5), i.e., it satisfies
(Lo —Aj
where the operators j are defined as j Liu^-1) + Am^. Solving equation (31) for um) is a standard linear algebra problem. The necessity for working with a pseudoinverse lies in the fact that ( Lo — ^m^) is singular and non-Hermitian. This prevents us from using the ordinary inverse to solve for Uj. Here, we employ the Moore-Penrose pseudoinverse denoted by Jl. (A brief review of the pseudoinverse is given in the accompanying Supplementary Note.) While this choice is not unique, it is convenient since the Moore-Penrose pseudoinverse can be computed efficiently via singular value decomposition.
After plugging Uj from equation (30) into equation (31), it is clear that the proof amounts to verifying that
(Lo — 40)) (Lo — 40)) fmf
Employing the defining properties of the Moore-Penrose pseudoinverse, one finds that the claim of equation (32) can be written in the equivalent form
-¿0>fm -fim-
where PL _^(o) is the orthogonal projector onto the range of ( Lq — l^)- Since every projector acts as the identity on vectors from its range, equation (33) holds if j belongs to the range < of PL _ ^( o) - Furthermore, since PL _ ^(o) is an orthogonal projector, the range and the null space of IP^ are orthogonal subspaces. This implies that f^ belongs to if and only if/jM is orthogonal to Oî^- To see that note that (w^ Jjp) = 0 as a direct consequence of equation (5). To prove that we now simply make use of the fact that is spanned by the left eigenstate, i.e.,^ = span{w^}.
We verify the last statement as follows. Since wj0 is the left eigenstate of Lo with eigenvalue lj0, it is clear that {wj0, ( Lo — lj°00A) _0 for any matrix A. Thus, w( orthogonal to the range of L0{
range, wj0 is also orthogonal to the range < of PL — ^o). Thus, wj0 is an element of the null space of PL _ . Assuming that is a non-degenerate eigenvalue of D_o, we thus find that is spanned by as stated.
Steady-state corrections. For the steady state ps defined by Lps _ 0, i.e. ps ; uo, the recursion relations (8) and (9) can be simplified significantly. To see this, recall that L and L0 are proper generators of the quantum dynamical semigroup. In order to be trace preserving (a necessary condition for a proper generator), the identity 1 must be the left eigenstate of IL and D_o with eigenvalue zero, i.e. Wq = Wq0'' = 3L- It follows that H_ is also the left eigenstate of L1 with eigenvalue zero since L1_L{L0 and thus TrPLLix] = 0 for any operator x. It is straightforward to show from TrplLix] = 0 and equations (8) and (9) that V/eN,
rs^-—L0J1L1pO'—1)
Series expansion of amplitude matrices. The series expansion of fs is valid if all fsj) can be determined according to
j-1 / ,t f — psj) - X fsk) fs^k)(36) k—1
which is equation (15). Here Z" was defined via Z"A — fso)A{ +A(fso)and fso) is determined through Cholesky decomposition: pso) — fs0)(fs0)){. Equation (36) is a system of linear equations and thus has a unique solution if and only if Z" is invertible.
Whether Z" is invertible depends on the form of fs0) , which in turn depends on the Hermitian and positive-semidefinite pso). If one of the eigenvalues of pso) is zero, there is a corresponding eigenvector |y) such that pso)|y) — 0. Consider the decomposition: pso) — hh where h is a Hermitian matrix. The existence of this decomposition can be simply proven by writing pso) in its eigen-decomposition form. Since |y) is the eigenvector of pso) with eigenvalue zero, |y) is also the eigenvector of h with eigenvalue zero, i.e. h | y) = 0. Moreover, due to the fact that fso) (fso)— pso) — hh, fso) and h are unitarily right equivalent40, i.e. fso) — hS where S is a unitary matrix. Thus, |y) is also the left eigenvector of fso) with eigenvalue zero, i.e. (fso)){|y) — S{h|y) — 0. Now, there must be a right eigenvector of fso), denoted by |W), that corresponds to the same eigenvalue (which is zero), i.e. fso) |W) — 0. We can show by directly substitution that Z"(T|W)(W|) — 0 where T is the matrix corresponding to Gaussian elimination which transforms (|W) (W|) to a lower triangular matrix. Therefore, Z" is not invertible if pso) contains at least one eigenvalue zero.
If Z0 is not invertible, equation (36) either has infinitely many solutions or no solution. The former case, in which there are infinitely many solutions, can be bypassed if we shift the eigenvalues of pso) away from zero. To do so, we consider a shift by an identity-matrix component according to
pS0)-PS0)+C1 (37)
We choose the auxiliary parameter C as close to zero as possible while maintaining numerical stability. In this way, we obtain a procedure to obtain a unique solution to equation (36). This method is similar to the correction matrix method37,38 for the Cholesky decomposition of matrices with eigenvalue(s) zero. There, a small diagonal correction matrix is also added to the original matrix to avoid the eigenvalue(s) zero.
If we encounter the case in which there is no solution, we cannot determine fsj). In fact, a two-level system coupled to finite temperature bath with D [s+ ] term treated as the perturbation belongs to this case. The failure to determine fsj) indicates that the series expansion of fs is invalid. This originates from the fact that if pso) contains any zero eigenvalues, the leading order term of some elements of fs may be of the order of a1/2 (or a3/2, etc.) instead of a0. For real functions, an analogy would be the case y(a) = x2(a) where y can be written as a power series in a. If the leading order term of y is of the order of a, x is a series that only contains half-integer orders of a and thus is not a proper Taylor series. Different expansion types would be needed in this case which we do not discuss further in the present paper.
1. Baumann, K., Guerlin, C., Brennecke, F. & Esslinger, T. Dicke quantum phase transition with a superfluid gas in an optical cavity. Nature 464, 1301 (2010).
2. Barreiro, J. T. et al. An open-system quantum simulator with trapped ions. Nature 470, 486 (2011).
3. Lee, T. E., Häffner, H. & Cross, M. C. Collective quantum jumps of Rydberg atoms. Phys. Rev. Lett. 108, 023602 (2012).
4. Bishop, L. S. et al. Nonlinear response of the vacuum Rabi resonance. Nat Phys 5, 105 (2009).
5. Underwood, D. L., Shanks, W. E., Koch, J. & Houck, A. A. Low-disorder microwave cavity lattices for quantum simulation with photons. Phys. Rev. A 86, 023837 (2012).
6. Jin, J., Rossini, D., Fazio, R., Leib, M. & Hartmann, M. J. Photon solid phases in driven arrays of nonlinearly coupled cavities. Phys. Rev. Lett. 110,163605 (2013).
7. Nissen, F., Fink, J. M., Mlynek, J. A., Wallraff, A. & Keeling, J. Collective suppression of linewidths in circuit qed. Phys. Rev. Lett. 110, 203602 (2013).
8. Teufel, J. D. et al. Sideband cooling of micromechanical motion to the quantum ground state. Nature 475, 359 (2011).
9. Ramos, T., Sudhir, V., Stannigel, K., Zoller, P. & Kippenberg, T. J. Nonlinear quantum optomechanics via individual intrinsic two-level defects. Phys. Rev. Lett. 110, 193602 (2013).
10. Pflanzer, A. C., Romero-Isart, O. & Cirac, J. I. Optomechanics assisted by a qubit: From dissipative state preparation to many-partite systems. Phys. Rev. A 88, 033804 (2013).
11. Houck, A. A., Tureci, H. E. & Koch, J. On-chip quantum simulation with superconducting circuits. Nat Phys 8, 292 (2012).
12. Kessler, E. M. eta/. Dissipative phase transition in a central spin system. Phys. Rev. A 86, 012116 (2012).
13. Sieberer, L. M., Huber, S. D., Altman, E. & Diehl, S. Dynamical critical phenomena in driven-dissipative systems. Phys. Rev. Lett. 110, 195301 (2013).
14. Sieberer, L. M., Huber, S. D., Altman, E. & Diehl, S. Non-equilibrium functional renormalization for driven-dissipative Bose-Einstein condensation. ArXiv e-prints (2013). 1309.7027.
15. Breuer, H. & Petruccione, F. The Theory Of Open Quantum Systems (Oxford University Press, 2002).
16. Dalibard, J., Castin, Y. & Molmer, K. Wave-function approach to dissipative processes in quantum optics. Phys. Rev. Lett. 68, 580-583 (1992).
17. Verstraete, F., Garcia-Ripoll, J. J. & Cirac, J. I. Matrix product density operators: Simulation of finite-temperature and dissipative systems. Phys. Rev. Lett. 93, 207204 (2004).
18. Zwolak, M. & Vidal, G. Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent super-operator renormalization algorithm. Phys. Rev. Lett. 93, 207205 (2004).
19. Hartmann, M. J., Prior, J., Clark, S. R. & Plenio, M. B. Density matrix renormalization group in the Heisenberg picture. Phys. Rev. Lett. 102, 057202 (2009).
20. Degenfeld-Schonburg, P. & Hartmann, M. J. Self-consistent projection operator approach to quantum many-body systems. ArXiv e-prints (2013). 1307.7027.
21. Flindt, C., Novotny, T. c. v., Braggio, A., Sassetti, M. & Jauho, A.-P. Counting statistics of non-markovian quantum stochastic processes. Phys. Rev. Lett. 100, 150601 (2008).
22. Flindt, C., Novotny, T. c. v., Braggio, A. & Jauho, A.-P. Counting statistics of transport through Coulomb blockade nanostructures: High-order cumulants and non-markovian effects. Phys. Rev. B 82,155407 (2010).
23. Yi, X. X., Li, C. & Su, J. C. Perturbative expansion for the master equation and its applications. Phys. Rev. A 62, 013819 (2000).
24. Cirac, J. I., Blatt, R., Zoller, P. & Phillips, W. D. Laser cooling of trapped ions in a standing wave. Phys. Rev. A 46, 2668-2681 (1992).
25. Reiter, F. & Sorensen, A. S. Effective operator formalism for open quantum systems. Phys. Rev. A 85, 032111 (2012).
26. Kessler, E. M. Generalized Schrieffer-Wolff formalism for dissipative systems. Phys. Rev. A 86, 012126 (2012).
27. Cai, Z. & Barthel, T. Algebraic versus exponential decoherence in dissipative many-particle systems. Phys. Rev. Lett. 111, 150403 (2013).
28. Fleming, C. H. & Cummings, N. I. Accuracy of perturbative master equations. Phys. Rev. E 83, 031117 (2011).
29. Kamenev, A. Field Theory Of Non-Equilibrium Systems (Cambridge University Press, 2011).
30. Benatti, F., Nagy, A. & Narnhofer, H. Asymptotic entanglement and Lindblad dynamics: a perturbative approach. J. Phys. A: Math. Theor. 44, 155303 (2011).
31. del Valle, E. & Hartmann, M. J. Correlator expansion approach to stationary states of weakly coupled cavity arrays. J. Phys. B: At. Mol. Opt. Phys. 46, 224023 (2013).
32. Petz, D. Monotone metrics on matrix spaces. Linear Algebra Appl. 244, 81-96 (1996).
33. Rehacek, J., Hradil, Z. & Jezek, M. Iterative algorithm for reconstruction of entangled states. Phys. Rev. A 63, 040303 (2001).
34. Lvovsky, A. Iterative maximum-likelihood reconstruction in quantum homodyne tomography. J. Opt. B: Quantum Semiclass. Opt. 6, S556 (2004).
35. Horn, R. & Johnson, C. Matrix Analysis (Cambridge University Press, 1985).
36. Chruscmski, D. & Kossakowski, A. Feshbach projection formalism for open quantum systems. Phys. Rev. Lett. 111, 050402 (2013).
37. Gill, P. & Murray, W. Newton-type methods for unconstrained and linearly constrained optimization. Math. Program. 7, 311-350 (1974).
38. Fang, H.-r. & O'Leary, D. P. Modified Cholesky algorithms: a catalog with new approaches. Math. Program. 115, 319-349 (2008).
39. Nissen, F. et al. Nonequilibrium dynamics of coupled qubit-cavity arrays. Phys. Rev. Lett. 108, 233603 (2012).
40. Bernstein, D. Matrix Mathematics: Theory, Facts, And Formulas (Second Edition) (Princeton University Press, 2009).
Acknowledgments
We thank Guanyu Zhu and Joshua Dempster for valuable discussions. This research was supported by the NSF under Grant No. PHY-1055993 (A.C.Y.L. and J.K.) and by the South African Research Chair Initiative of the Department of Science and Technology and National Research Foundation (F.P.). J.K. and F.P. thank the National Institute for Theoretical Physics for supporting J.K.s visit through a NITheP visitor program grant.
Author contributions
A.C.Y.L. carried out the calculations and did most of the writing. F.P. and J.K. provided support and supervised the project.
Additional information
Supplementary information accompanies this paper at http://www.nature.com/ scientificreports
Competing financial interests: The authors declare no competing financial interests.
How to cite this article: Li, A.C.Y., Petruccione, F. & Koch, J. Perturbative approach to Markovian open quantum systems. Sci. Rep. 4, 4887; DOI:10.1038/srep04887 (2014).
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License. The images in this article are included in the article's Creative Commons license, unless indicated otherwise in the image credit;
if the image is not included under the Creative Commons license, users will need to obtain permission from the license holder in order to reproduce the image. To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-nd/3.0/
