1. Introduction

Quantum computers (QCs) with large numbers of quantum bits (qubits) promise to solve important problems such as factorization of larger integer numbers, searching large unsorted databases, and simulations of physical processes exponentially faster than digital computers. Recently, many efforts were made for designing scalable (in the number of qubits) QC architectures based on solid-state implementations. One of the most promising designs of a solid-state QC is based on superconducting devices with Josephson junctions and solid-state quantum interference devices (SQUIDs) serving as qubits (effective spins), which operate in the quantum regime: , where is the temperature and is the qubit transition frequency. This condition is widely used in superconducting quantum computation and quantum measurement, when   10–20 mK and   100–150 mK (in temperature units) [1–8] (see also references therein). The main advantages of a QC with superconducting qubits are: (i) the two basic states of a qubit are represented by the states of a superconducting charge or current in the macroscopic (several m size) device. The relatively large scale of this device facilitates the manufacturing, and potential controlling and measuring of the states of qubits. (ii) The states of charge- and current-based qubits can be measured using rapidly developing technologies, such as a single electron transistors, effective resonant oscillators and microcavities with RF amplifiers, and quantum tunneling effects. (iii) The quantum logic operations can be implemented exclusively by switching locally on and off voltages on controlling microcontacts and magnetic fluxes. (iv) The devices based on superconducting qubits can potentially achieve large quantum dephasing and relaxation times of milliseconds and more at low temperatures, allowing quantum coherent computation for long enough times. In spite of significant progress, current devices with superconducting qubits only have one or two qubits operating with low fidelity even for simplest operations.

One of the main problems which must be resolved in order to build a scalable QC is to develop novel approaches for suppression of unwanted effects such as decoherence and noise. This also requires to develop rigorous mathematical tools for analyzing the dynamics of decoherence, entanglement, and thermalization in order to control the quantum protocols with needed fidelity. These theoretical approaches must work for long enough times and be applicable to both solvable and not explicitly solvable (nonintegrable) systems.

Here we present a review of our results [9–11] on the rigorous analysis of the phenomenon of decoherence and relaxation for general -level systems coupled to reservoirs. The latter are described by the bosonic fields. We suggest a new approach which applies to a wide variety of systems which are not explicitly solvable. We analyze in detail the dynamics of an -qubit quantum register collectively coupled to a thermal environment. Each spin experiences the same environment interaction, consisting of an energy conserving and an energy exchange part. We find the decay rates of the reduced density matrix elements in the energy basis. We show that the fastest decay rates of off-diagonal matrix elements induced by the energy conserving interaction are of order , while the one induced by the energy exchange interaction is of the order only. Moreover, the diagonal matrix elements approach their limiting values at a rate independent of . Our method is based on a dynamical quantum resonance theory valid for small, fixed values of the couplings, and uniformly in time for . We do not make Markov-, Born- or weak coupling (van Hove) approximations.

2. Presentation of Results

We consider an -level quantum system interacting with a heat bath . The former is described by a Hilbert space and a Hamiltonian The environment is modelled by a bosonic thermal reservoir with Hamiltonian acting on the reservoir Hilbert space , and where and are the usual bosonic creation and annihilation operators satisfying the canonical commutation relations . It is understood that we consider in the thermodynamic limit of infinite volume () and fixed temperature (in a phase without condensate). Given a form factor , a square integrable function of (momentum representation), the smoothed-out creation and annihilation operators are defined as and , respectively, and the hermitian field operator is The total Hamiltonian, acting on , has the form where is a coupling constant and is an interaction operator linear in field operators. For simplicity of exposition, we consider here initial states where and are not entangled, and where is in thermal equilibrium. (Our method applies also to initially entangled states, and arbitrary initial states of normal w.r.t. the equilibrium state; see [10]). The initial density matrix is thus where is any state of and is the equilibrium state of at temperature .

Let be an arbitrary observable of the system (an operator on the system Hilbert space ) and set where is the density matrix of at time and is the reduced density matrix of . In our approach, the dynamics of the reduced density matrix is expressed in terms of the resonance structure of the system. Under the noninteracting dynamics (), the evolution of the reduced density matrix elements of , expressed in the energy basis of , is given by where . As the interaction with the reservoir is turned on, the dynamics (2.8) undergoes two qualitative changes.(1)The “Bohr frequencies” in the exponent of (2.8) become complex, . It can be shown generally that the resonance energies have nonnegative imaginary parts, . If , then the corresponding dynamical process is irreversible.(2)The matrix elements do not evolve independently any more. Indeed, the effective energy of is changed due to the interaction with the reservoirs, leading to a dynamics that does not leave eigenstates of invariant. (However, to lowest order in the interaction, the eigenspaces of are left invariant and therefore matrix elements with belonging to a fixed energy difference will evolve in a coupled manner.)

Our goal is to derive these two effects from the microscopic (hamiltonian) model and to quantify them. Our analysis yields the thermalization and decoherence times of quantum registers.

2.1. Evolution of Reduced Dynamics of an -Level System

Let be an observable of the system . We show in [9, 10] that the ergodic averages exist, that is, converges in the ergodic sense as . Furthermore, we show that for any and for any , where the complex numbers are the eigenvalues of a certain explicitly given operator , lying in the strip . They have the expansions where and are the eigenvalues of a matrix , called a level-shift operator, acting on the eigenspace of corresponding to the eigenvalue (which is a subspace of ). The in (2.11) are linear functionals of and are given in terms of the initial state, , and certain operators depending on the Hamiltonian . They have the expansion where is the collection of all pairs of indices such that , with being the eigenvalues of . Here, is the -matrix element of the observable in the energy basis of , and are coefficients depending on the initial state of the system (and on , but not on nor on ).

2.1.1. Discussion

It is well known that there are two kinds of processes which drive decay (or irreversibility) of : energy-exchange processes characterized by and energy preserving ones where . The former are induced by interactions having nonvanishing probabilities for processes of absorption and emission of field quanta with energies corresponding to the Bohr frequencies of and thus typically inducing thermalization of . Energy preserving interactions suppress such processes, allowing only for a phase change of the system during the evolution (“phase damping”, [12–18]).

To our knowledge, energy-exchange systems have so far been treated using Born and Markov master equation approximations (Lindblad form of dynamics) or they have been studied numerically, while for energy conserving systems, one often can find an exact solution. The present representation (2.11) gives a detailed picture of the dynamics of averages of observables for interactions with and without energy exchange. The resonance energies and the functionals can be calculated for concrete models, as illustrated in the next section. We mention that the resonance dynamics representation can be used to study the dynamics of entanglement of qubits coupled to local and collective reservoirs, see [19].

The dynamical resonance method can be generalized to time-dependent Hamiltonians. See [20, 21] for time-periodic Hamiltonians.

2.1.2. Contrast with Weak Coupling Approximation

Our representation (2.11) of the true dynamics of relies only on the smallness of the coupling parameter , and no approximation is made. In the absence of an exact solution, it is common to make a weak coupling Lindblad master equation approximation of the dynamics, in which the reduced density matrix evolves according to , where is the Lindblad generator, [22–24]. This approximation can lead to results that differ qualitatively from the true behaviour. For instance, the Lindblad master equation predicts that the system approaches its Gibbs state at the temperature of the reservoir in the limit of large times. However, it is clear that in reality, the coupled system will approach equilibrium, and hence the asymptotic state of alone, being the reduction of the coupled equilibrium state, is the Gibbs state of only to first approximation in the coupling (see also illustration below, and [9, 10]). In particular, the system's asymptotic density matrix is not diagonal in the original energy basis, but it has off-diagonal matrix elements of . Features of this kind cannot be captured by the Lindblad approximation, but are captured in our approach.

It has been shown (see, e.g., [23–26]) that the weak coupling limit dynamics generated by the Lindblad operator is obtained in the regime , , with fixed. One of the strengths of our approach is that we do not impose any relation between and , and our results are valid for all times , provided is small. It has been observed [25, 26] that for certain systems of the type , the second-order contribution of the exponents in (2.12) correspond to eigenvalues of the Lindblad generator. Our resonance method gives the true exponents, that is, we do not lose the contributions of any order in the interaction. If the energy spectrum of is degenerate, it happens that the second-order contributions to vanish. This corresponds to a Lindblad generator having several real eigenvalues. In this situation, the correct dynamics (approach to a final state) can be captured only by taking into account higher-order contributions to the exponents ; see [27]. To our knowledge, so far this can only be done with the method presented in this paper, and is beyond the reach of the weak coupling method.

2.1.3. Illustration: Single Qubit

Consider to be a single spin with energy gap . is coupled to the heat bath via the operator where is the Bose field operator (2.3), smeared out with a coupling function (form factor) , , and the coupling matrix (representing the coupling operator in the energy eigenbasis) is hermitian. The operator (2.14)—or a sum of such terms, for which our technique works equally well—is the most general coupling which is linear in field operators. We refer to [10] for a discussion of the link between (2.14) and the spin-boson model. We take initially in a coherent superposition in the energy basis, In [10] we derive from representation (2.11) the following expressions for the dynamics of matrix elements, for all : where the resonance energies are given by with The remainder terms in (2.17), (2.17) satisfy , uniformly in , and they can be decomposed into a sum of a constant and a decaying part, where , with . These relations show the following.(i)To second order in , convergence of the populations to the equilibrium values (Gibbs law), and decoherence occur exponentially fast, with rates and , respectively. (If either of these imaginary parts vanishes then the corresponding process does not take place, of course.) In particular, coherence of the initial state stays preserved on time scales of the order ; compare for example (2.18).(ii)The final density matrix of the spin is not the Gibbs state of the qubit, and it is not diagonal in the energy basis. The deviation of the final state from the Gibbs state is given by . This is clear heuristically too, since typically the entire system approaches its joint equilibrium in which and are entangled. The reduction of this state to is the Gibbs state of modulo terms representing a shift in the effective energy of due to the interaction with the bath. In this sense, coherence in the energy basis of is created by thermalization. We have quantified this in [10, Theorem ].(iii)In a markovian master equation approach, the above phenomenon (i.e., variations of in the time-asymptotic limit) cannot be detected. Indeed in that approach one would conclude that approaches its Gibbs state as .

2.2. Evolution of Reduced Dynamics of an -Level System

In the sequel, we analyze in more detail the evolution of a qubit register of size . The Hamiltonian is where are pair interaction constants and is the value of a magnetic field at the location of spin . The Pauli spin operator is and is the matrix acting on the th spin only.

We consider a collective coupling between the register and the reservoir : the distance between the qubits is much smaller than the correlation length of the reservoir and as a consequence, each qubit feels the same interaction with the reservoir. The corresponding interaction operator is (compare with (2.4)) Here and are form factors and the coupling constants and measure the strengths of the energy conserving (position-position) coupling, and the energy exchange (spin flip) coupling, respectively. Spin-flips are implemented by the in (2.23), representing the Pauli matrix acting on the th spin. The total Hamiltonian takes the form (2.4) with replaced by (2.23). It is convenient to represent as a matrix in the energy basis, consisting of eigenvectors of . These are vectors in indexed by spin configurations where so that We denote the reduced density matrix elements as The Bohr frequencies (2.9) are now and they become complex resonance energies under perturbation.

Assumption of Nonoverlapping Resonances
The Bohr frequencies (2.29) represent “unperturbed” energy levels and we follow their motion under perturbation (). In this work, we consider the regime of nonoverlapping resonances, meaning that the interaction is small relative to the spacing of the Bohr frequencies.

We show in [10, Theorem ], that for all , This result is obtained by specializing (2.11) to the specific system at hand and considering observables . In (2.30), we have in accordance with (2.10), . The coefficients are overlaps of resonance eigenstates which vanish unless (see point () after (2.9)). They represent the dominant contribution to the functionals in (2.11); see also (2.13). The have the expansion where the label indexes the splitting of the eigenvalue into distinct resonance energies. The lowest order corrections satisfy They are the (complex) eigenvalues of an operator , called the level shift operator associated to . This operator acts on the eigenspace of associated to the eigenvalue (a subspace of the qubit register Hilbert space; see [10, 11] for the formal definition of ). It governs the lowest order shift of eigenvalues under perturbation. One can see by direct calculation that .

2.2.1. Discussion

2.2.2. Decoherence Rates

We illustrate our results on decoherence rates for a qubit register with (the general case is treated in [11]). We consider generic magnetic fields defined as follows. For , , we have Condition (2.34) is satisfied generically in the sense that it does not hold only for very special choices of (one such special choice is ). For instance, if the are chosen to be independent, and uniformly random from an interval , then (2.34) is satisfied with probability one. We show in [11, Theorem ], that the decoherence rates (2.33) are given by Here, is contributions coming from the energy conserving interaction; and are due to the spin flip interaction. The term is due to both interactions and is of . We give explicit expressions for , , , and in [11, Section ]. For the present purpose, we limit ourselves to discussing the properties of the latter quantities.(i)Properties of : vanishes if either is such that or the infrared behaviour of the coupling function is too regular (in three dimensions with ). Otherwise, . Moreover, is proportional to the temperature .(ii)Properties of : if for all (form factor in spherical coordinates). For low temperatures, , , for high temperatures approaches a constant.(iii)Properties of : if either of , or vanish, or if is infrared regular as mentioned above, then . Otherwise, , in which case approaches constant values for both .(iv)Full decoherence: if for all , then all off-diagonal matrix elements approach their limiting values exponentially fast. In this case, we say that full decoherence occurs. It follows from the above points that we have full decoherence if and for all , and provided are small enough (so that the remainder term in (2.35) is small). Note that if , then matrix elements associated to energy differences such that will not decay on the time scale given by the second order in the perturbation (). We point out that generically, will reach a joint equilibrium as , which means that the final reduced density matrix of is its Gibbs state modulo a peturbation of the order of the interaction between and ; see [9, 10]. Hence generically, the density matrix of does not become diagonal in the energy basis as .(v)Properties of : depends on the energy exchange interaction only. This reflects the fact that for a purely energy conserving interaction, the populations are conserved [9, 10, 17]. If for all , then (this is sometimes called the “Fermi Golden Rule Condition”). For small temperatures , , while approaches a finite limit as .

In terms of complexity analysis, it is important to discuss the dependence of on the register size .(i)We show in [11] that is independent of . This means that the thermalization time, or relaxation time of the diagonal matrix elements (corresponding to ), is in .(ii)To determine the order of magnitude of the decay rates of the off-diagonal density matrix elements (corresponding to ) relative to the register size , we assume the magnetic field to have a certain distribution denoted by . We show in [11] that where and are positive constants (independent of ), with . Here, is the number of indices such that for each s.t. , and is the Hamming distance between the spin configurations and (which depends on only).(iii)Consider . It follows from (2.35)–(2.37) that for purely energy conserving interactions (), , which can be as large as . On the other hand, for purely energy exchanging interactions (), we have , which cannot exceed . If both interactions are acting, then we have the additional term , which is of order . This shows the following: The fastest decay rate of reduced off-diagonal density matrix elements due to the energy conserving interaction alone is of order , while the fastest decay rate due to the energy exchange interaction alone is of the order . Moreover, the decay of the diagonal matrix elements is of order , that is, independent of .

Remark 2.2.2 s. (1) For , the model can be solved explicitly [17], and one shows that the fastest decaying matrix elements have decay rate proportional to . Furthermore, the model with a noncollective, energy-conserving interaction, where each qubit is coupled to an independent reservoir, can also be solved explicitly [17]. The fastest decay rate in this case is shown to be proportional to .
(2) As mentioned at the beginning of this section, we take the coupling constants , so small that the resonances do not overlap. Consequently, and are bounded above by a constant proportional to the gradient of the magnetic field in the present situation; see also [11]. Thus the decay rates do not increase indefinitely with increasing in the regime considered here. Rather, are attenuated by small coupling constants for large . They are of the order . We have shown that modulo an overall, common (-dependent) prefactor, the decay rates originating from the energy conserving and exchanging interactions differ by a factor .
(3) Collective decoherence has been studied extensively in the literature. Among the many theoretical, numerical, and experimental works, we mention here only [12, 14, 17, 28, 29], which are closest to the present work. We are not aware of any prior work giving explicit decoherence rates of a register for not explicitly solvable models, and without making master equation technique approximations.

3. Resonance Representation of Reduced Dynamics

The goal of this section is to give a precise statement of the core representation (2.11), and to outline the main ideas behind the proof of it.

The -level system is coupled to the reservoir (see also (2.1), (2.2)) through the operator where each is a hermitian matrix, the are form factors, and the are coupling constants. Fix any phase and define where and . The phase is a parameter which can be chosen appropriately as to satisfy the following condition.(A) The map has an analytic extension to a complex neighbourhood of the origin, as a map from to .

Examples of satisfying (A) are given by , where , , , and .

This condition ensures that the technically simplest version of the dynamical resonance theory, based on complex spectral translations, can be implemented. The technical simplicity comes at a price: on one hand, it limits the class of admissible functions , which have to behave appropriately in the infrared regime so that the parts of (3.2) fit nicely together at , to allow for an analytic continuation. On the other hand, the square root in (3.2) must be analytic as well, which implies the condition .

It is convenient to introduce the doubled Hilbert space , whose normalized vectors accommodate any state on the system (pure or mixed). The trace state, or infinite temperature state, is represented by the vector via Here are the orthonormal eigenvectors of . This is just the Gelfand-Naimark-Segal construction for the trace state. Similarly, let and be the Hilbert space and the vector representing the equilibrium state of the reservoirs at inverse temperature . In the Araki-Woods representation of the field, we have , where is the bosonic Fock space over the one-particle space and , with being the Fock vacuum of (see also [10, 11] for more detail). Let be the vector in representing the density matrix at time . It is not difficult to construct the unique operator in satisfying (See also [10] for concrete examples.) We define the reference vector and set

Theorem 3.1 (Dynamical resonance theory [9–11]). Assume condition (A) with a fixed satisfying . There is a constant s.t.; if , then the limit , (2.10), exists for all observables . Moreover, for all such and for all , we have The are given by (2.12), counts the splitting of the eigenvalue into distinct resonance energies , and the are (nonorthogonal) finite-rank projections.

This result is the basis for a detailed analysis of the reduced dynamics of concrete systems, like the -qubit register introduced in Section 2.2. We obtain (2.30) (in particular, the overlap functions ) from (3.8) by analyzing the projections in more detail. Let us explain how to link the overlap to its initial value for a nondegenerate Bohr energy , and where . (The latter observables used in (2.11) give the matrix elements of the reduced density matrix in the energy basis.)

The is the spectral (Riesz) projection of an operator associated with the eigenvalue ; see (3.19) (In reality, we consider a spectral deformation , where is a complex parameter. This is a technical trick to perform our analysis. Physical quantities do not depend on and therefore, we do not display this parameter here). If a Bohr energy , (2.9), is simple, then there is a single resonance energy bifurcating out of , as . In this case, the projection has rank one, , where and are eigenvectors of and its adjoint, with eigenvalue and its complex conjugate, respectively, and . From perturbation theory, we obtain , where and . The overlap in the sum of (3.8) becomes The choice in (2.6) gives , the reduced density matrix element. With this choice of , the main term in (3.9) becomes (see also (3.3)) In the second-last step, we commute to the right through , since belongs to the commutant of the algebra of observables of . In the last step, we use .