Resonant Perturbation Theory of Decoherence and Relaxation of Quantum Bits

Abstract

We describe our recent results on the resonant perturbation theory of decoherence and relaxation for quantum systems with many qubits. The approach represents a rigorous analysis of the phenomenon of decoherence and relaxation for general N-level systems coupled to reservoirs of bosonic fields. We derive a representation of the reduced dynamics valid for all times 𝑡≥0 and for small but fixed interaction strength. Our approach does not involve master equation approximations and applies to a wide variety of systems which are not explicitly solvable.

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: ℏ𝜔>>𝑘_B𝑇, 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 𝑡≥0. We do not make Markov-, Born- or weak coupling (van Hove) approximations.

2. Presentation of Results

We consider an 𝑁-level quantum system S interacting with a heat bath R. The former is described by a Hilbert space 𝔥_𝑠=ℂ^𝑁 and a Hamiltonian𝐻_S=diag𝐸₁,…,𝐸_𝑁.(2.1) The environment R is modelled by a bosonic thermal reservoir with Hamiltonian𝐻_R=_ℝ3𝑎^∗(𝑘)||𝑘||𝑎(𝑘)d³𝑘,(2.2) acting on the reservoir Hilbert space 𝔥_R, and where 𝑎^∗(𝑘) and 𝑎(𝑘) are the usual bosonic creation and annihilation operators satisfying the canonical commutation relations [𝑎(𝑘),𝑎^∗(𝑙)]=𝛿(𝑘−𝑙). It is understood that we consider R in the thermodynamic limit of infinite volume (ℝ³) and fixed temperature 𝑇=1/𝛽>0 (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 𝑎^∗(𝑓)=∫_ℝ3𝑓(𝑘)𝑎^∗(𝑘)d³𝑘 and 𝑎(𝑓)=∫_ℝ3𝑓(𝑘)𝑎(𝑘)d³𝑘, respectively, and the hermitian field operator is𝜙(𝑓)=1√2𝑎^∗(𝑓)+𝑎(𝑓).(2.3) The total Hamiltonian, acting on 𝔥_S⊗𝔥_R, has the form𝐻=𝐻_S⊗𝟙_R+𝟙_S⊗𝐻_R+𝜆𝑣,(2.4) where 𝜆 is a coupling constant and 𝑣 is an interaction operator linear in field operators. For simplicity of exposition, we consider here initial states where S and R are not entangled, and where R is in thermal equilibrium. (Our method applies also to initially entangled states, and arbitrary initial states of R normal w.r.t. the equilibrium state; see [10]). The initial density matrix is thus 𝜌₀=𝜌₀⊗𝜌_R,𝛽,(2.5) where 𝜌₀ is any state of S and 𝜌_R,𝛽 is the equilibrium state of R at temperature 1/𝛽.

Let 𝐴 be an arbitrary observable of the system (an operator on the system Hilbert space 𝔥_S) and set⟨𝐴⟩_𝑡∶=Tr_S𝜌_𝑡𝐴=Tr_S+R𝜌_𝑡𝐴⊗𝟙_R,(2.6) where 𝜌_𝑡 is the density matrix of S+R at time 𝑡 and 𝜌_𝑡=Tr_R𝜌_𝑡(2.7) is the reduced density matrix of S. 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 (𝜆=0), the evolution of the reduced density matrix elements of S, expressed in the energy basis {𝜑_𝑘}^𝑁𝑘=1 of 𝐻_S, is given by𝜌_𝑡_𝑘𝑙=𝜑_𝑘,e^{−i𝑡𝐻S}𝜌₀e^i𝑡𝐻S𝜑_𝑙=e^{i𝑡𝑒𝑙𝑘}𝜌₀_𝑘𝑙,(2.8) where 𝑒_𝑙𝑘=𝐸_𝑙−𝐸_𝑘. As the interaction with the reservoir is turned on, the dynamics (2.8) undergoes two qualitative changes.(1)The “Bohr frequencies” 𝑒∈𝐸−𝐸^∶𝐸,𝐸^∈spec𝐻_S(2.9) in the exponent of (2.8) become complex, 𝑒↦𝜀_𝑒. It can be shown generally that the resonance energies 𝜀_𝑒 have nonnegative imaginary parts, Im𝜀_𝑒≥0. If Im𝜀_𝑒>0, then the corresponding dynamical process is irreversible.(2)The matrix elements do not evolve independently any more. Indeed, the effective energy of S is changed due to the interaction with the reservoirs, leading to a dynamics that does not leave eigenstates of 𝐻_S invariant. (However, to lowest order in the interaction, the eigenspaces of 𝐻_S 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 𝐴∈ℬ(𝔥_S) be an observable of the system S. We show in [9, 10] that the ergodic averages⟨⟨𝐴⟩⟩_∞∶=lim_𝑇→∞1𝑇^𝑇0⟨𝐴⟩_𝑡d𝑡(2.10) exist, that is, ⟨𝐴⟩_𝑡 converges in the ergodic sense as 𝑡→∞. Furthermore, we show that for any 𝑡≥0 and for any 0<𝜔^<2𝜋/𝛽,⟨𝐴⟩_𝑡−⟨⟨𝐴⟩⟩_∞=_𝜀≠0e^i𝑡𝜀𝑅_𝜀(𝐴)+𝑂𝜆²e^{−[𝜔−𝑂(𝜆)]𝑡},(2.11) where the complex numbers 𝜀 are the eigenvalues of a certain explicitly given operator 𝐾(𝜔′), lying in the strip {𝑧∈ℂ|0≤Im𝑧<𝜔′/2}. They have the expansions𝜀≡𝜀^𝑒(𝑠)=𝑒+𝜆²𝛿^𝑒(𝑠)+𝑂𝜆⁴,(2.12) where 𝑒∈spec(𝐻_S⊗𝟙_S−𝟙_S⊗𝐻_S)=spec(𝐻_S)−spec(𝐻_S) and 𝛿^𝑒(𝑠) are the eigenvalues of a matrix Λ_𝑒, called a level-shift operator, acting on the eigenspace of 𝐻_S⊗𝟙_S−𝟙_S⊗𝐻_S corresponding to the eigenvalue 𝑒 (which is a subspace of 𝔥_S⊗𝔥_S). 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𝑅_𝜀(𝐴)=_{(𝑚,𝑛)∈𝐼𝑒}𝜘_𝑚,𝑛𝐴_𝑚,𝑛+𝑂𝜆²,(2.13) where 𝐼_𝑒 is the collection of all pairs of indices such that 𝑒=𝐸_𝑚−𝐸_𝑛, with 𝐸_𝑘 being the eigenvalues of 𝐻_S. Here, 𝐴_𝑚,𝑛 is the (𝑚,𝑛)-matrix element of the observable 𝐴 in the energy basis of 𝐻_S, and 𝜘_𝑚,𝑛 are coefficients depending on the initial state of the system (and on 𝑒, but not on 𝐴 nor on 𝜆).

2.1.1. Discussion

(i)In the absence of interaction (𝜆=0), we have 𝜀=𝑒∈ℝ; see (2.12). Depending on the interaction, each resonance energy 𝜀 may migrate into the upper complex plane, or it may stay on the real axis, as 𝜆≠0.(ii)The averages ⟨𝐴⟩_𝑡 approach their ergodic means ⟨⟨𝐴⟩⟩_∞ if and only if Im𝜀>0 for all 𝜀≠0. In this case, the convergence takes place on the time scale [Im𝜀]⁻¹. Otherwise; ⟨𝐴⟩_𝑡 oscillates. A sufficient condition for decay is that Im𝛿^𝑒(𝑠)>0 (and 𝜆 small, see (2.12)).(iii)The error term in (2.11) is small in 𝜆, uniformly in 𝑡≥0, and it decays in time quicker than any of the main terms in the sum on the r.h.s.: indeed, Im𝜀=𝑂(𝜆²) while 𝜔^−𝑂(𝜆)>𝜔^/2 independent of small values of 𝜆. However, this means that we are in the regime 𝜆²≪𝜔^<2𝜋/𝛽 (see before (2.11)), which implies that 𝜆² must be much smaller than the temperature 𝑇=1/𝛽. Using a more refined analysis, one can get rid of this condition; see also remarks on page 376 of [10].(iv)Relation (2.13) shows that to lowest order in the perturbation, the group of (energy basis) matrix elements of any observable 𝐴 corresponding to a fixed energy difference 𝐸_𝑚−𝐸_𝑛 evolve jointly, while those of different such groups evolve independently.

It is well known that there are two kinds of processes which drive decay (or irreversibility) of S: energy-exchange processes characterized by [𝑣,𝐻_S]≠0 and energy preserving ones where [𝑣,𝐻_S]=0. 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 S and thus typically inducing thermalization of S. 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 S 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 𝜌_𝑡=e^𝑡ℒ𝜌₀, 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 S 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 S+R will approach equilibrium, and hence the asymptotic state of S alone, being the reduction of the coupled equilibrium state, is the Gibbs state of S 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 𝜆→0, 𝑡→∞, 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 𝑡≥0, provided 𝜆 is small. It has been observed [25, 26] that for certain systems of the type S+R, 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 𝐻_S is degenerate, it happens that the second-order contributions to Im𝜀^𝑒(𝑠) 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 S to be a single spin 1/2 with energy gap Δ=𝐸₂−𝐸₁>0. S is coupled to the heat bath R via the operator𝑣=𝑎𝑐𝑐𝑏⊗𝜙(𝑔),(2.14) where 𝜙(𝑔) is the Bose field operator (2.3), smeared out with a coupling function (form factor) 𝑔(𝑘), 𝑘∈ℝ³, and the 2×2 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 S initially in a coherent superposition in the energy basis,𝜌₀=121111.(2.15) In [10] we derive from representation (2.11) the following expressions for the dynamics of matrix elements, for all 𝑡≥0: 𝜌_𝑡_𝑚,𝑚=e𝑍^{−𝛽𝐸𝑚}_S,𝛽+(−21)^𝑚tanh𝛽2Δe^{i𝑡𝜀0(𝜆)}+𝑅_𝑚,𝑚(𝜆,𝑡),𝑚=1,2,(2.16)𝜌_𝑡_1,2=12e^{i𝑡𝜀−Δ(𝜆)}+𝑅_1,2(𝜆,𝑡),(2.17) where the resonance energies 𝜀 are given by𝜀_Δ𝜀₀(𝜆)=i𝜆²𝜋²|𝑐|²𝜉(Δ)+𝑂𝜆⁴,(𝜆)=Δ+𝜆²𝑅+i2𝜆²𝜀𝜋²_−Δ|𝑐|²𝜉(Δ)+(𝑏−𝑎)²𝜉(0)+𝑂𝜆⁴,(𝜆)=−𝜀_Δ(𝜆),(2.18) with𝜉(𝜂)=lim_𝜖↓01𝜋_ℝ³d³𝑘coth𝛽||𝑘||2||𝑔(𝑘)||²𝜖||𝑘||−𝜂²+𝜖²𝑅=1,2𝑏²+1−𝑎²𝑔,𝜔⁻¹𝑔2|𝑐|²P.V._ℝ×𝑆2𝑢²||𝑔(|𝑢|,𝜎)||²coth𝛽|𝑢|2𝑢−1Δ.(2.19) The remainder terms in (2.17), (2.17) satisfy |𝑅_𝑚,𝑛(𝜆,𝑡)|≤𝐶𝜆², uniformly in 𝑡≥0, and they can be decomposed into a sum of a constant and a decaying part, 𝑅_𝑚,𝑛(𝜆,𝑡)=𝑝_𝑛,𝑚_∞−𝛿_𝑚,𝑛e^{−𝛽𝐸𝑚}𝑍_S,𝛽+𝑅^𝑚,𝑛(𝜆,𝑡),(2.20) where |𝑅′_𝑚,𝑛(𝜆,𝑡)|=𝑂(𝜆²e^−𝛾𝑡), with 𝛾=min{Im𝜀₀,Im𝜀_±Δ}. 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 𝜏_𝑇=[Im𝜀₀(𝜆)]⁻¹ and 𝜏_𝐷=[Im𝜀_Δ(𝜆)]⁻¹, 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 𝜆⁻²[|𝑐|²𝜉(Δ)+(𝑏−𝑎)²𝜉(0)]⁻¹; 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 lim_𝑡→∞𝑅_𝑚,𝑛(𝜆,𝑡)=𝑂(𝜆²). This is clear heuristically too, since typically the entire system S+R approaches its joint equilibrium in which S and R are entangled. The reduction of this state to S is the Gibbs state of S modulo 𝑂(𝜆²) terms representing a shift in the effective energy of S due to the interaction with the bath. In this sense, coherence in the energy basis of S is created by thermalization. We have quantified this in [10, Theorem 3.3].(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 S 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𝐻_S=^𝑁_𝑖_,𝑗=1𝐽_𝑖𝑗𝑆^𝑧𝑖𝑆^𝑧𝑗+^𝑁_𝑗=1𝐵_𝑗𝑆^𝑧𝑗,(2.21) where 𝐽_𝑖𝑗 are pair interaction constants and 𝐵_𝑗 is the value of a magnetic field at the location of spin 𝑗. The Pauli spin operator is𝑆^𝑧=10−01(2.22) and 𝑆^𝑧𝑗 is the matrix 𝑆^𝑧 acting on the 𝑗th spin only.

We consider a collective coupling between the register S and the reservoir R: 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))𝜆₁𝑣₁+𝜆₂𝑣₂=𝜆_1𝑁_𝑗=1𝑆^𝑧𝑗⊗𝜙𝑔₁+𝜆_2𝑁_𝑗=1𝑆^𝑥𝑗⊗𝜙𝑔₂.(2.23) 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 𝑆^𝑥=0110(2.24) 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 𝐻_S. These are vectors in 𝔥_S=ℂ²⊗⋯⊗ℂ²=ℂ^2𝑁 indexed by spin configurations𝜎=𝜎₁,…,𝜎_𝑁∈{+1,−1}^𝑁,𝜑_𝜎=𝜑_𝜎1⊗⋯⊗𝜑_𝜎𝑁,(2.25) where𝜑₊=10,𝜑₋=01,(2.26) so that𝐻_S𝜑_𝜎=𝐸𝜎𝜑_𝜎with𝐸𝜎=^𝑁_𝑖_,𝑗=1𝐽_𝑖𝑗𝜎_𝑖𝜎_𝑗+^𝑁_𝑗=1𝐵_𝑗𝜎_𝑗.(2.27) We denote the reduced density matrix elements as𝜌_𝑡_𝜎,𝜏=𝜑_𝜎,𝜌_𝑡𝜑_𝜏.(2.28) The Bohr frequencies (2.9) are now𝑒𝜎,𝜏=𝐸𝜎−𝐸𝜏=^𝑁_𝑖_,𝑗=1𝐽_𝑖𝑗𝜎_𝑖𝜎_𝑗−𝜏_𝑖𝜏_𝑗+^𝑁_𝑗=1𝐵_𝑗𝜎_𝑗−𝜏_𝑗,(2.29) 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 2.1], that for all 𝑡≥0, 𝜌_𝑡_𝜎,𝜏−𝜌_∞_𝜎,𝜏=_{{𝑒∶𝜀𝑒≠0}}e^{i𝑡𝜀𝑒}⎡⎢⎢⎣_{𝜎′,𝜏′}+𝑂𝜆₂₁𝑤^𝜀𝑒_{𝜎,𝜏;𝜎′,𝜏′}𝜌₀_{𝜎′,𝜏′}+𝑂𝜆²¹+𝜆²²⎤⎥⎥⎦+𝜆²²e^{−[𝜔′+𝑂(𝜆)]𝑡}.(2.30) 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), ⟨⟨[𝜌_∞]_𝜎,𝜏⟩⟩=lim_𝑇→∞(1/𝑇)∫^𝑇0[𝜌_𝑡]_𝜎,𝜏d𝑡. The coefficients 𝑤 are overlaps of resonance eigenstates which vanish unless 𝑒=−𝑒(𝜎,𝜏)=−𝑒(𝜎^,𝜏^) (see point (2) after (2.9)). They represent the dominant contribution to the functionals 𝑅_𝜀 in (2.11); see also (2.13). The 𝜀_𝑒 have the expansion𝜀_𝑒≡𝜀^𝑒(𝑠)=𝑒+𝛿^𝑒(𝑠)+𝑂𝜆⁴¹+𝜆⁴²,(2.31) where the label 𝑠=1,…,𝜈(𝑒) indexes the splitting of the eigenvalue 𝑒 into 𝜈(𝑒) distinct resonance energies. The lowest order corrections 𝛿^𝑒(𝑠) satisfy𝛿^𝑒(𝑠)=𝑂𝜆²¹+𝜆²².(2.32) They are the (complex) eigenvalues of an operator Λ_𝑒, called the level shift operator associated to 𝑒. This operator acts on the eigenspace of 𝐿_S 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 Im𝛿^𝑒(𝑠)≥0.

2.2.1. Discussion

(i)To lowest order in the perturbation, the group of reduced density matrix elements [𝜌_𝑡]_𝜎,𝜏 associated to a fixed 𝑒=𝑒(𝜎,𝜏) evolve in a coupled way, while groups of matrix elements associated to different 𝑒 evolve independently.(ii)The density matrix elements of a given group mix and evolve in time according to the weight functions 𝑤 and the exponentials e^{i𝑡𝜀𝑒(𝑠)}. In the absence of interaction (𝜆₁=𝜆₂=0), all the 𝜀^𝑒(𝑠)=𝑒 are real. As the interaction is switched on, the 𝜀^𝑒(𝑠) typically migrate into the upper complex plane, but they may stay on the real line (due to some symmetry or due to an “inefficient coupling”).(iii)The matrix elements [𝜌_𝑡]_𝜎,𝜏 of a group 𝑒 approach their ergodic means if and only if all the nonzero 𝜀^𝑒(𝑠) have strictly positive imaginary part. In this case, the convergence takes place on a time scale of the order 1/𝛾_𝑒, where 𝛾_𝑒=minIm𝜀^𝑒(𝑠)∶𝑠=1,…,𝜈(𝑒)s.t.𝜀^𝑒(𝑠)≠0(2.33) is the decay rate of the group associated to 𝑒. If an 𝜀^𝑒(𝑠) stays real, then the matrix elements of the corresponding group oscillate in time. A sufficient condition for decay of the group associated to 𝑒 is 𝛾_𝑒>0, that is, Im𝛿^𝑒(𝑠)>0 for all 𝑠, and 𝜆₁, 𝜆₂ small.

2.2.2. Decoherence Rates

We illustrate our results on decoherence rates for a qubit register with 𝐽_𝑖𝑗=0 (the general case is treated in [11]). We consider generic magnetic fields defined as follows. For 𝑛_𝑗∈{0,±1,±2}, 𝑗=1,…,𝑁, we have ^𝑁_𝑗=1𝐵_𝑗𝑛_𝑗=0⟺𝑛_𝑗=0∀𝑗.(2.34) 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 𝐵_𝑗=constant). For instance, if the 𝐵_𝑗 are chosen to be independent, and uniformly random from an interval [𝐵_min,𝐵_max], then (2.34) is satisfied with probability one. We show in [11, Theorem 2.3], that the decoherence rates (2.33) are given by𝛾_𝑒=𝜆²¹𝑦₁(𝑒)+𝜆𝜆²²𝑦₂(𝑒)+𝑦₁₂(𝑒),𝑒≠0₂₂𝑦₀,𝑒=0+𝑂𝜆⁴¹+𝜆⁴².(2.35) 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 2]. For the present purpose, we limit ourselves to discussing the properties of the latter quantities.(i)Properties of 𝑦₁(𝑒): 𝑦₁(𝑒) vanishes if either 𝑒 is such that 𝑒₀∶=∑^𝑛𝑗=1(𝜎_𝑗−𝜏_𝑗)=0 or the infrared behaviour of the coupling function 𝑔₁ is too regular (in three dimensions 𝑔₁∝|𝑘|^𝑝 with 𝑝>−1/2). Otherwise, 𝑦₁(𝑒)>0. Moreover, 𝑦₁(𝑒) is proportional to the temperature 𝑇.(ii)Properties of 𝑦₂(𝑒): 𝑦₂(𝑒)>0 if 𝑔₂(2𝐵_𝑗,Σ)≠0 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 𝑦₁₂(𝑒)=0. Otherwise, 𝑦₁₂(𝑒)>0, in which case 𝑦₁₂(𝑒) approaches constant values for both 𝑇→0,∞.(iv)Full decoherence: if 𝛾_𝑒>0 for all 𝑒≠0, 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 𝜆₂≠0 and 𝑔₂(2𝐵_𝑗,Σ)≠0 for all 𝑗, and provided 𝜆₁,𝜆₂ are small enough (so that the remainder term in (2.35) is small). Note that if 𝜆₂=0, then matrix elements associated to energy differences 𝑒 such that 𝑒₀=0 will not decay on the time scale given by the second order in the perturbation (𝜆²¹). We point out that generically, S+R will reach a joint equilibrium as 𝑡→∞, which means that the final reduced density matrix of S is its Gibbs state modulo a peturbation of the order of the interaction between S and R; see [9, 10]. Hence generically, the density matrix of S 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 𝑔₂(2𝐵_𝑗,Σ)≠0 for all 𝑗, then 𝑦₀>0 (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 𝑒=0), is 𝑂(1) in 𝑁.(ii)To determine the order of magnitude of the decay rates of the off-diagonal density matrix elements (corresponding to 𝑒≠0) relative to the register size 𝑁, we assume the magnetic field to have a certain distribution denoted by ⟨⋅⟩. We show in [11] that ⟨𝑦₁⟩=𝑦₁∝𝑒²⁰,⟨𝑦₂⟩=𝐶_𝐵𝔇𝜎−𝜏,⟨𝑦₁₂⟩=𝑐_𝐵𝜆₁,𝜆₂𝑁₀(𝑒),(2.36) where 𝐶_𝐵 and 𝑐_𝐵=𝑐_𝐵(𝜆₁,𝜆₂) are positive constants (independent of 𝑁), with 𝑐_𝐵(𝜆₁,𝜆₂)=𝑂(𝜆²¹+𝜆²²). Here, 𝑁₀(𝑒) is the number of indices 𝑗 such that 𝜎_𝑗=𝜏_𝑗 for each (𝜎,𝜏) s.t. 𝑒(𝜎,𝜏)=𝑒, and 𝔇𝜎−𝜏∶=^𝑁_𝑗=1||𝜎_𝑗−𝜏_𝑗||(2.37) is the Hamming distance between the spin configurations 𝜎 and 𝜏 (which depends on 𝑒 only).(iii)Consider 𝑒≠0. It follows from (2.35)–(2.37) that for purely energy conserving interactions (𝜆₂=0), 𝛾_𝑒∝𝜆²¹𝑒²⁰=𝜆²¹[∑^𝑁𝑗=1(𝜎_𝑗−𝜏_𝑗)]², which can be as large as 𝑂(𝜆²¹𝑁²). On the other hand, for purely energy exchanging interactions (𝜆₁=0), 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 𝜆₂=0, 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𝑣=^𝑅_𝑟=1𝜆_𝑟𝐺_𝑟⊗𝜙𝑔_𝑟,(3.1) where each 𝐺_𝑟 is a hermitian 𝑁×𝑁 matrix, the 𝑔_𝑟(𝑘) are form factors, and the 𝜆_𝑟∈ℝ are coupling constants. Fix any phase 𝜒∈ℝ and define𝑔_𝑟,𝛽(𝑢,𝜎)∶=𝑢1−e^−𝛽𝑢|𝑢|^1/2𝑔^𝑟−e(𝑢_i𝜒,𝜎)if𝑢≥0,𝑔_𝑟(−𝑢,𝜎)if𝑢<0,(3.2) 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 𝐿²(ℝ³,d³𝑘).

Examples of 𝑔 satisfying (A) are given by 𝑔(𝑟,𝜎)=𝑟^𝑝e^−𝑟𝑚𝑔₁(𝜎), where 𝑝=−1/2+𝑛, 𝑛=0,1,…, 𝑚=1,2, and 𝑔₁(𝜎)=e^i𝜙𝑔₁(𝜎).

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 𝑢=0, 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 𝜔^<2𝜋/𝛽.

It is convenient to introduce the doubled Hilbert space ℋ_S=𝔥_S⊗𝔥_S, whose normalized vectors accommodate any state on the system S (pure or mixed). The trace state, or infinite temperature state, is represented by the vectorΩ_S=1√𝑁^𝑁_𝑗=1𝜑_𝑗⊗𝜑_𝑗(3.3) viaℬ𝔥_S∋𝐴⟼⟨Ω_S,(𝐴⊗𝟙)Ω_S⟩.(3.4) Here 𝜑_𝑗 are the orthonormal eigenvectors of 𝐻_S. This is just the Gelfand-Naimark-Segal construction for the trace state. Similarly, let ℋ_R and Ω_R,𝛽 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 ℋ_R=ℱ⊗ℱ, where ℱ is the bosonic Fock space over the one-particle space 𝐿²(ℝ³,d³𝑘) and Ω_R,𝛽=Ω⊗Ω, with Ω being the Fock vacuum of ℱ (see also [10, 11] for more detail). Let 𝜓₀⊗Ω_R,𝛽 be the vector in ℋ_S⊗ℋ_R representing the density matrix at time 𝑡=0. It is not difficult to construct the unique operator in 𝐵∈𝟙_S⊗𝔥_S satisfying 𝐵Ω_S=𝜓₀.(3.5) (See also [10] for concrete examples.) We define the reference vector Ω_ref∶=Ω_S⊗Ω_R,𝛽(3.6) and set 𝜆=max_{𝑟=1,…,𝑅}||𝜆_𝑟||.(3.7)

Theorem 3.1 (Dynamical resonance theory [9–11]). Assume condition (A) with a fixed 𝜔^ satisfying 0<𝜔^<2𝜋/𝛽. There is a constant 𝑐₀ s.t.; if 𝜆≤𝑐₀/𝛽, then the limit ⟨⟨𝐴⟩⟩_∞, (2.10), exists for all observables 𝐴∈ℬ(𝔥_S). Moreover, for all such 𝐴 and for all 𝑡≥0, we have ⟨𝐴⟩_𝑡−⟨⟨𝐴⟩⟩_∞=_{𝑒,𝑠∶𝜀𝑒(𝑠)≠0𝜈(𝑒)}_𝑠=1+𝑂𝜆₂e^{i𝑡𝜀𝑒(𝑠)}𝐵^∗𝜓₀⊗Ω_R,𝛽,𝑄^𝑒(𝑠)𝐴⊗𝟙_SΩ_refe^{−[𝜔′+𝑂(𝜆)]𝑡}.(3.8) The 𝜀^𝑒(𝑠) are given by (2.12), 1≤𝜈(𝑒)≤mult(𝑒) 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 ⟨(𝐵^∗𝜓₀)⊗Ω_R,𝛽,𝑄^𝑒(𝑠)(𝐴⊗𝟙_S)Ω_ref⟩ 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 𝜆≠0. In this case, the projection 𝑄_𝑒≡𝑄^𝑒(𝑠) has rank one, 𝑄_𝑒=|𝜒_𝑒⟩⟨𝜒_𝑒|, where 𝜒_𝑒 and 𝜒_𝑒 are eigenvectors of 𝐾_𝜆 and its adjoint, with eigenvalue 𝜀_𝑒 and its complex conjugate, respectively, and ⟨𝜒_𝑒,𝜒_𝑒⟩=1. From perturbation theory, we obtain 𝜒_𝑒=𝜒_𝑒=𝜑_𝑘⊗𝜑_𝑙⊗Ω_R,𝛽+𝑂(𝜆), where 𝐻_S𝜑_𝑗=𝐸_𝑗𝜑_𝑗 and 𝐸_𝑘−𝐸_𝑙=𝑒. The overlap in the sum of (3.8) becomes𝐵=^∗𝜓₀𝐵⊗_∗Ω_R,𝛽,𝑄_𝑒𝐴⊗𝟙_SΩ_ref=𝐵_∗𝜓₀⊗Ω_R,𝛽,||𝜑_𝑘⊗𝜑_𝑙⊗Ω_R,𝛽𝜑_𝑘⊗𝜑_𝑙⊗Ω_R,𝛽||𝐴⊗𝟙_SΩ_ref+𝑂𝜆²𝜓₀,||𝜑_𝑘⊗𝜑_𝑙⟩⟨𝜑_𝑘⊗𝜑_𝑙||𝐴⊗𝟙_SΩ_S+𝑂𝜆².(3.9) 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))𝐵^∗𝜓₀,||𝜑_𝑘⊗𝜑_𝑙⟩⟨𝜑_𝑘⊗𝜑_𝑙||𝐴⊗𝟙_SΩ_S=1√𝑁𝛿^𝑘𝑛=𝛿_𝑘𝑛𝛿_𝑙𝑚⟨𝐵^∗𝜓₀,𝜑_𝑛⊗𝜑_𝑚⟩=𝛿_𝑘𝑛𝛿_𝑙𝑚𝐵^∗𝜓₀,||𝜑_𝑛⟩⟨𝜑_𝑚||⊗𝟙_SΩ_S=𝛿_𝑘𝑛𝛿_𝑙𝑚𝜓₀,||𝜑_𝑛⟩⟨𝜑_𝑚||⊗𝟙_S𝐵Ω_S𝛿_𝑙𝑚𝜌₀_𝑚𝑛.(3.10) In the second-last step, we commute 𝐵 to the right through |𝜑_𝑛⟩⟨𝜑_𝑚|⊗𝟙_S, since 𝐵 belongs to the commutant of the algebra of observables of S. In the last step, we use 𝐵Ω_S=𝜓₀.

Combining (3.9) and (3.10) with Theorem 3.1 we obtain, in case 𝑒=𝐸_𝑚−𝐸_𝑛 is a simple eigenvalue, 𝜌_𝑡_𝑚𝑛−𝜌_∞_𝑚𝑛=_{{𝑒,𝑠∶𝜀𝑒(𝑠)≠0}}e^{i𝑡𝜀𝑒(𝑠)}𝛿_𝑘𝑛𝛿_𝑙𝑚𝜌₀_𝑚𝑛+𝑂𝜆²+𝑂𝜆²e^{−[𝜔′+𝑂(𝜆)]𝑡}.(3.11) This explains the form (2.30) for a simple Bohr energy 𝑒. The case of degenerate 𝑒 (i.e., where several different pairs of indices 𝑘,𝑙 satisfy 𝐸_𝑘−𝐸_𝑙=𝑒) is analyzed along the same lines; see [11] for details.

3.1. Mechanism of Dynamical Resonance Theory, Outline of Proof of Theorem 3.1

Consider any observable 𝐴∈𝐵(𝔥_S). We have⟨𝐴⟩_𝑡==TTrr_S_S+R𝜌_𝑡𝐴=𝜓₀𝜌_𝑡𝐴⊗𝟙_R,e^{i𝑡𝐿𝜆}𝐴⊗𝟙_S⊗𝟙_Re^{−i𝑡𝐿𝜆}𝜓₀.(3.12) In the last step, we pass to the representation Hilbert space of the system (the GNS Hilbert space), where the initial density matrix is represented by the vector 𝜓₀ (in particular, the Hilbert space of the small system becomes 𝔥_S⊗𝔥_S); see also before (3.3), (3.4). As mentioned above, in this review we consider initial states where S and R are not entangled. The initial state is represented by the product vector 𝜓₀=Ω_S⊗Ω_R,𝛽, where Ω_S is the trace state of S, (3.4), ⟨Ω_S,(𝐴⊗𝟙_S)Ω_S⟩=(1/𝑁)Tr(𝐴), and where Ω_R,𝛽 is the equilibrium state of R at a fixed inverse temperature 0<𝛽<∞. The dynamics is implemented by the group of automorphisms e^{i𝑡𝐿𝜆}⋅e^{−i𝑡𝐿𝜆}. The self-adjoint generator 𝐿_𝜆 is called the Liouville operator. It is of the form 𝐿_𝜆=𝐿₀+𝜆𝑊, where 𝐿₀=𝐿_S+𝐿_R represents the uncoupled Liouville operator, and 𝜆𝑊 is the interaction (3.1) represented in the GNS Hilbert space. We refer to [10, 11] for the specific form of 𝑊.

We borrow a trick from the analysis of open systems far from equilibrium: there is a (nonself-adjoint) generator 𝐾_𝜆 s.t. e^{i𝑡𝐿𝜆}𝐴e^{−i𝑡𝐿𝜆}=e^{i𝑡𝐾𝜆}𝐴e^{−i𝑡𝐾𝜆}𝐾forallobservables𝐴,𝑡≥0,and_𝜆𝜓₀=0.(3.13)

𝐾_𝜆 can be constructed in a standard way, given 𝐿_𝜆 and the reference vector 𝜓₀. 𝐾_𝜆 is of the form 𝐾_𝜆=𝐿₀+𝜆𝐼, where the interaction term undergoes a certain modification (𝑊→𝐼); see for example [10]. As a consequence, formally, we may replace the propagators in (3.12) by those involving 𝐾. The resulting propagator which is directly applied to 𝜓₀ will then just disappear due to the invariance of 𝜓₀. One can carry out this procedure in a rigorous manner, obtaining the following resolvent representation [10]⟨𝐴⟩_𝑡=−12𝜋i_ℝ−i𝜓⁰,𝐾_𝜆(𝜔)−𝑧⁻¹𝐴⊗𝟙_S⊗𝟙_R𝜓₀e^i𝑡𝑧d𝑧,(3.14) where 𝐾_𝜆(𝜔)=𝐿₀(𝜔)+𝜆𝐼(𝜔), 𝐼 is representing the interaction, and 𝜔↦𝐾_𝜆(𝜔) is a spectral deformation (translation) of 𝐾_𝜆. The latter is constructed as follows. There is a deformation transformation 𝑈(𝜔)=e^−i𝜔𝐷, where 𝐷 is the (explicit) self-adjoint generator of translations [10, 11, 30] transforming the operator 𝐾_𝜆 as𝐾_𝜆(𝜔)=𝑈(𝜔)𝐾_𝜆𝑈(𝜔)⁻¹=𝐿₀+𝜔𝑁+𝜆𝐼(𝜔).(3.15)

Here, 𝑁=𝑁₁⊗𝟙+𝟙⊗𝑁₁ is the total number operator of a product of two bosonic Fock spaces ℱ⊗ℱ (the Gelfand-Naimark-Segal Hilbert space of the reservoir), and where 𝑁₁ is the usual number operator on ℱ. 𝑁 has spectrum ℕ∪{0}, where 0 is a simple eigenvalue (with vacuum eigenvector Ω_R,𝛽=Ω⊗Ω). For real values of 𝜔, 𝑈(𝜔) is a group of unitaries. The spectrum of 𝐾_𝜆(𝜔) depends on Im𝜔 and moves according to the value of Im𝜔, whence the name “spectral deformation”. Even though 𝑈(𝜔) becomes unbounded for complex 𝜔, the r.h.s. of (3.15) is a well-defined closed operator on a dense domain, analytic in 𝜔 at zero. Analyticity is used in the derivation of (3.14) and this is where the analyticity condition (A) after (3.2) comes into play. The operator 𝐼(𝜔) is infinitesimally small with respect to the number operator 𝑁. Hence we use perturbation theory in 𝜆 to examine the spectrum of 𝐾_𝜆(𝜔).

The point of the spectral deformation is that the (important part of the) spectrum of 𝐾_𝜆(𝜔) is much easier to analyze than that of 𝐾_𝜆, because the deformation uncovers the resonances of 𝐾_𝜆. We have (see Figure 1) spec𝐾₀(𝜔)=𝐸_𝑖−𝐸_𝑗_{𝑖,𝑗=1,…,𝑁}_𝑛≥1{𝜔𝑛+ℝ},(3.16) because 𝐾₀(𝜔)=𝐿₀+𝜔𝑁, 𝐿₀ and 𝑁 commute, and the eigenvectors of 𝐿₀=𝐿_S+𝐿_R are 𝜑_𝑖⊗𝜑_𝑗⊗Ω_R,𝛽. Here, we have 𝐻_S𝜑_𝑗=𝐸_𝑗𝜑_𝑗. The continuous spectrum is bounded away from the isolated eigenvalues by a gap of size Im𝜔. For values of the coupling parameter 𝜆 small compared to Im𝜔, we can follow the displacements of the eigenvalues by using analytic perturbation theory. (Note that for Im𝜔=0, the eigenvalues are imbedded into the continuous spectrum, and analytic perturbation theory is not valid! The spectral deformation is indeed very useful!)

Figure 1 

Spectrum of 𝐾0(𝜔).

Theorem 3.2 (see [10] and Figure 2). Fix Im𝜔 s.t. 0<Im𝜔<𝜔^ (where 𝜔^ is as in Condition (A)). There is a constant 𝑐₀>0 s.t. if |𝜆|≤𝑐₀/𝛽 then, for all 𝜔 with I𝑚𝜔>7𝜔^/8, the spectrum of 𝐾_𝜆(𝜔) in the complex half-plane {Im𝑧<𝜔^/2} is independent of 𝜔 and consists purely of the distinct eigenvalues 𝜀^𝑒(𝑠)∶𝑒∈spec𝐿_S,𝑠=1,…,𝜈(𝑒),(3.17) where 1≤𝜈(𝑒)≤mult(𝑒) counts the splitting of the eigenvalue 𝑒. Moreover, lim_𝜆→0||𝜀^𝑒(𝑠)(𝜆)−𝑒||=0(3.18) for all 𝑠, and we have Im𝜀^𝑒(𝑠)≥0. Also, the continuous spectrum of 𝐾_𝜆(𝜔) lies in the region {Im𝑧≥3𝜔^/4}.

Figure 2 

Spectrum of 𝐾𝜆(𝜔). Resonances 𝜀𝑒(𝑠) are uncovered.

Next we separate the contributions to the path integral in (3.14) coming from the singularities at the resonance energies and from the continuous spectrum. We deform the path of integration 𝑧=ℝ−i into the line 𝑧=ℝ+i𝜔^/2, thereby picking up the residues of poles of the integrand at 𝜀^𝑒(𝑠) (all 𝑒, 𝑠). Let 𝒞^𝑒(𝑠) be a small circle around 𝜀^𝑒(𝑠), not enclosing or touching any other spectrum of 𝐾_𝜆(𝜔). We introduce the generally nonorthogonal Riesz spectral projections𝑄^𝑒(𝑠)=𝑄^𝑒(𝑠)(𝜔,𝜆)=−12𝜋i_𝒞^𝑒(𝑠)𝐾_𝜆(𝜔)−𝑧⁻¹d𝑧.(3.19)

It follows from (3.14) that⟨𝐴⟩_𝑡=_{𝑒𝜈(𝑒)}_𝑠=1e^{i𝑡𝜀𝑒(𝑠)}𝜓⁰,𝑄^𝑒(𝑠)𝐴⊗𝟙_S⊗𝟙_R𝜓₀+𝑂𝜆²e^{−𝜔′𝑡/2}.(3.20) Note that the imaginary parts of all resonance energies 𝜀^𝑒(𝑠) are smaller than 𝜔^/2, so that the remainder term in (3.20) is not only small in 𝜆, but it also decays faster than all of the terms in the sum. (See also Figure 3.) We point out also that instead of deforming the path integration contour as explained before (3.19), we could choose 𝑧=ℝ+i[𝜔^−𝑂(𝜆)], hence transforming the error term in (3.20) into the one given in (3.8).

Figure 3 

Contour deformation: ∫ℝ−id𝑧=∑𝑒,𝑠∫𝒞𝑒(𝑠)d𝑧+∫ℝ+i𝜔′/2d𝑧.

Finally, we notice that all terms in (3.20) with 𝜀^𝑒(𝑠)≠0 will vanish in the ergodic mean limit, so ⟨⟨𝐴⟩⟩_∞=lim_𝑇→∞1𝑇^𝑇0⟨𝐴⟩_𝑡d𝑡=_{𝑠∶𝜀0(𝑠)=0}𝜓⁰,𝑄^0(𝑠)𝐴⊗𝟙_R⊗𝟙_R𝜓₀.(3.21) We now see that the linear functionals (2.13) are represented as𝑅_𝜀^𝑒(𝑠)(𝐴)=𝜓₀,𝑄^𝑒(𝑠)𝐴⊗𝟙_S⊗𝟙_R𝜓₀.(3.22) This concludes the outline of the proof of Theorem 3.1.

Acknowledgments

This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract no. DE-AC52-06NA25396 and by Lawrence Livermore National Laboratory under Contract no. DE-AC52-07NA27344. This research was funded by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA). All statements of fact, opinion or conclusions contained herein are those of the authors and should not be constructed as representing the official views or polices of IARPA, the ODNI, or the U.S. Government. I. M. Sigal also acknowledges a support by NSERC under Grant NA 7901. M. Merkli also acknowledges a support by NSERC under Grant 205247.