Dissipative dynamics of a biased qubit coupled to a harmonic oscillator: Analytical results beyond the rotating wave approximation

We study the dissipative dynamics of a biased two-level system (TLS) coupled to a harmonic oscillator (HO), the latter interacting with an Ohmic environment. Using Van-Vleck perturbation theory and going to second order in the coupling between TLS and HO, we show how the Hamiltonian of the TLS-HO system can be diagonalized analytically. Our model represents an improvement to the usually used Jaynes-Cummings Hamiltonian as an initial rotating wave approximation is avoided. By assuming a weak coupling to the thermal bath, analytical expressions for the time evolution of the populations of the TLS are found: the population is characterized by a multiplicity of damped oscillations together with a complex relaxation dynamics towards thermal equilibrium. The long time evolution is characterized by a single relaxation rate, which is largest at resonance and whose expression can be given in closed analytic form.


Introduction
In recent years the spin-boson model [1] - [3] has experienced a strong revival, as it is well suited to describe dissipative and decoherence effects on the dynamics of a two-level system (TLS) or qubit coupled to a bath. Crucial for the effects of the environment on the dynamics of the TLS is the shape of the spectral density of the harmonic bath. It is common to assume an Ohmic spectral density, which is linear in the continuous bath modes. In this work we concentrate on a so-called structured bath, for which the spectral density is Ohmic at low frequencies but exhibits a Lorenztian-shaped peak at a certain frequency Ω. It has been shown in [4] that a spin-boson model with such an effective spectral density can be exactly mapped on the model of a TLS which is coupled to a single harmonic oscillator (HO) of frequency Ω, where the latter feels the influence of an Ohmic bath. Due to its wide applicability the TLS-oscillator system has been object of intense research along the years. So it reflects for example the physics of single atoms with a large electric dipole moment coupled to the microwave photons of a cavity [5], or quantum dots in photonic crystals [6,7]. More recently the model has received quite some attention in the field of quantum computation, where two-level systems are used to implement the two logical states of a qubit. We will especially focus on the solidstate implementation of such systems. Here, two prominent realizations of a qubitoscillator system are the Cooper-pair box (CPB) [8] - [11] coupled to a transmission line resonantor [12] - [17] and the Josephson flux qubit [18] read out by a dc-SQUID [19] - [22]. Inspired by experiments with real atoms interacting with a cavity mode, one speaks for the CPB case of circuit quantum electrodynamics, as now the CPB plays the role of an artificial atom and the waveguide acts as a cavity. From such a setup one expects a huge step towards the realization of a quantum computer, as the transmission line resonator can be used to couple qubits together [9,23], store the information of qubits or to provide non-demolition read-out schemes [12,15]. Concerning the flux qubit, the read-out usually happens through a damped dc-SQUID, which is inductively coupled to the qubit. However, through the SQUID enviromental noise is transferred to the qubit leading to decoherence and dissipation within its dynamics. The effect of this noise on the qubit depends very much on the strength g of the coupling between qubit and SQUID and one faces a conflicting situation. On the one hand one wants a strong coupling for a good read-out resolution. On the other hand the coupling should be minimized to keep the negative effects of the environment as small as possible. In [24,25] it has been shown that the qubit-SQUID system can be described by a spinboson model with an effective spectral density G eff (ω) exhibting a peak at the plasma frequency Ω of the SQUID. Applying the above mentioned mapping an equivalent point of view is to consider the SQUID as an LC-circuit coupled to the Ohmic bath and model it as a harmonic oscillator. A detailed description of a nondestructive read-out scheme is e.g. given in [26]. The spin-boson model can be formally solved using e.g. real-time path integral methods [1,2]. However, in order to get closed-form analytical results, approximations must be invoked. A quite common one is the so-called weak coupling approximation (WCA), which is perturbative in the bath spectral density [2]. However, it has been shown that for strong qubit-HO coupling g and for small detuning δ = Ω − ∆ b , where ∆ b is the qubit energy splitting, such an approximation breaks down [27], as coherent exchange processes between TLS and oscillator are disregared. For an unbiased qubit the noninteracting blip approximation (NIBA) used in [28] - [31] circumvents this problem as it is non-perturbative in the coupling g and therefore takes correctly into account the influence of the oscillator on the TLS. Moreover, it allows an analytic treatment of the dynamics. However, the NIBA is known to break down for a biased qubit at low temperatures [1,2]. Another approach, which treats the system non-perturbatively in the bath is the flow-equation renormalization method [32,33], where the spinboson Hamiltonian is diagonalized using infinitesimal unitary transformations. However, whithin this approach analytical solutions are difficult to find. Recently a polaron transformation was used by Huang et al to obtain analytically the population dynamics and confirm the Shiba's relation for an unbiased TLS [34]. In the case in which the qubit and the HO are considered as the central quantum system being coupled to an Ohmic bath, the numerical, ab-initio quasiadiabatic propagator path-integral (QUAPI) method [35,36] is a nice tool as it enables to cover both the resonant regime, where the oscillator frequency is close to the qubit energy splitting, and the dispersive regime with the oscillator being far detuned from the qubit [27,29,30]. Moreover, it can be applied to a biased as well as to an unbiased TLS and therefore be used as a testbed for analytical results. For qubits being operated at the degeneracy point, which means an unbiased TLS, very often a rotating wave approximation (RWA) is applied [12], which is expected to be valid for small detuning and yields as starting point the Jaynes-Cummings Hamiltonian [37,38]. This model was first used to study a two-state atom interacting with a single, close to resonance cavity mode of the electromagnetic field and predicts e.g. the repeated revival and collapse of Rabi oscillations within the atomic excitation probability. By condsidering the TLS-HO system in the representation of displaced HO states, Brito et al were able to truncate the infinite Hilbert space of this system without loosing the effects of the HO on the TLS dynamics [39]. However, so far none of these works could provide an analytical expression for the dynamics of the dissipative qubit being valid for zero as well as non-zero detuning and for both a biased and unbiased TLS. In this work an analytic expression for the dissipative qubit's dynamics which includes the effects of a finite detuning and of a static bias is derived. Specifically, starting from the qubit-HO perspective, the eigenvalues and eigenfunctions of the non-dissipative TLS-HO system are found approximately using Van-Vleck perturbation theory up to second order in the coupling g. Notice that no rotating wave approximation is required. Dissipation effects are then evaluated by solving a Born-Markov master equation for the reduced density matrix in the system's eigenbasis. The structure of the work is as follows. The dissipative TLS-HO Hamiltonian and the main dynamical quantities are introduced in section 2. Inspired by the work of Goorden et al [29,30], we demonstrate in section 3 how the eigenstates and eigenenergies of the non-dissipative Hamiltonian can be found approximately using Van-Vleck perturbation theory [40,41]. In this way we can provide an analytical formula for the non-dissipative dynamics whichs takes into account the full Hilbert space of the qubit-HO system. After that, we show how for low temperatures (k B T < Ω, ∆ b ) this infinite Hilbert space can be truncated and discuss the relevant contributions of the HO to the dynamics. In section 4 the influence of the environment is investigated, by looking at solutions of the Bloch-Redfield equations. Specifically, analytical expressions fo the TLS dynamics are obtained and compared with numerical solutions. The main physical features of the coupled TLS-HO system are discussed in section 5. To illustrate the effects of counter-rotating terms in the Hamiltonian of the qubit-HO system, which are neglected performing a RWA, we compare in section 6 our calculations to results obtained from the Jaynes-Cummings model.

The model
In this section we introduce the Hamiltonian for a qubit coupled through a harmonic oscillator to a thermal bath. Further, a formula for the population difference between the qubit's two logical states is derived.

The qubit-oscillator-bath system
To set up the model we consider the Hamiltonian of a qubit-HO system, H QHO , which is coupled to an environmental bath, H B , by the interaction Hamiltonian H OB , so that the total Hamiltonian becomes The Hamiltonian, H QHO = H 0 + H Int , consists of the Hamiltonian of the TLS/qubit and the harmonic oscillator, and the interaction term The Hamiltonian of the TLS is given in the subspace {|L , |R }, corresponding to a clockwise or counterclockwise current in the superconducting loop of a three-junction Josephson qubit or more generally to the qubit's two logical states. In the case of a superconducting flux-qubit, the energy bias ε can be tuned by an applied external flux, Φ ext , and is zero at the so-called degeneracy point. The tunnelling amplitude is described by ∆ 0 . For ε ≫ ∆ 0 the states |L and |R are eigenstates of H TLS , whereas at the degeneracy point those eigenstates are a symmetric and antisymmetric superposition of the two logical states. Further, B and B † are the annihilation and creation operator for the HO with frequency Ω, and g characterizes the coupling strength. We also introduce the energy splitting ∆ b ≡ ε 2 + ∆ 2 0 between the groundstate |g and the excited state |e of the TLS. Using the transformation with tan Θ = −∆ 0 /ε and − π 2 ≤ Θ < π 2 , we obtain the Hamiltonian of the TLS in the this basis: 2σ z . The states |R and |L become in the energy basis The Following Caldeira and Leggett [42], we model the environmental influences originating from the circuitry surrounding the qubit and the oscillator as a bath of harmonic oscillators being coupled bilinearly to the HO. Thus, the environment is described by

and the interaction Hamiltonian is
The operators b † k and b k are the creation and destruction operator, respectively, for the k th bath oscillator, ω k is its frequency and ν k gives the coupling strength. The whole bath can be described by its spectral density, which we consider to be Ohmic: In [4] it is shown that the above model is equivalent to that of a TLS being coupled directly to a harmonic bath including the single oscillator of frequency Ω; i.e., a spinboson model [1,2] with a peaked effective spectral density, The relation between α and the coupling parameter g between the qubit and the HO is g = Ω α/(8κ) [24,25]. This second perspective is suitable for calculating the dynamics of the qubit using a path-integral approach, as it was done for example in [31] for the case of an unbiased qubit (ε = 0). The approach in [31], however, being based on the NIBA [2], is not suitable to investigate the low temperature dynamics of a biased TLS. Thus, in this manuscript we will consider the TLS and the single oscillator as central quantum system and solve the Bloch-Redfield master equations for the density matrix of this system, which are valid also for the case of a biased TLS.

The population difference
The main goal of this work is to determine the dynamics P (t) of the qubit. That means, we wish to calculate the population difference between the |R and |L states of the qubit. The reduced density matrix of the TLS, is found after tracing out the oscillator and bath degrees of freedom from the total density matrix W (t) = e − i Ht W (0)e i Ht . In turn ρ(t) = Tr B {W (t)} is the reduced density matrix of the qubit-HO system. How to calculate this density matrix will be shown later. After some algebra, illustrated in more detail in Appendix A, we arrive at an expression for P (t), given in terms of diagonal and off-diagonal elements of ρ(t) in the TLS-HO eigenbasis {|n }. It reads where p nm (t) = 2 j cos Θ jg|n m|jg − je|n m|je with ρ nm (t) = n|ρ(t)|m . How to determine the eigenstates of H QHO is described in the next section.

Energy spectrum and dynamics of the non-dissipative TLS-HO system
In this section we show how to find the eigenvalues of the unperturbed qubit-HO Hamiltonian H QHO approximately by using Van-Vleck perturbation theory [40,41]. The idea is to take advantage of the degenerate or doublet structure of the energy spectrum of the uncoupled (g = 0) TLS-HO system near resonance, e.g. at ∆ b ≈ Ω. Then, as long as the perturbation is small compared to the energy separation of the different doublets, the full Hamiltonian will exhibit a similar spectrum of bundled energy levels.

Energy spectrum
The eigenenergies of the uncoupled TLS-HO system are immediately found by applying the HamiltonianH 0 =H TLS + H HO on the eigenstates in (7): The dashed lines in figure 1 show the energy spectrum corresponding to (15) vs. the oscillator frequency Ω for the five lowest eigenstates. Except for the groundstate, |0g , the states |(j + 1)g and |je are degenerate in the resonant case (Ω = ∆ b ). Close to resonance the spectrum exhibits a doublet structure. With the coupling being switched on, the full Hamiltonian H QHO reads in the basis {|jg ; |je }. In order to diagonalize the HamiltonianH QHO we considerH Int as a small perturbation, which is resonable as long as g ≪ ∆ b , Ω. Applying Van-Vleck perturbation theory we construct an effective Hamiltonian, having the same eigenvalues asH QHO but no matrix elements connecting states which are far off from degeneracy. Thus,H eff will be block-diagonal with all quasi-degenerate energy levels being in one common block. As in our case always two states are nearly degenerate, each block ofH eff builds a two-by-two matrix. This can be easily diagonalized in order to determine the eigenstates. Following [40,41] we calculate the transformation matrix S up to second order in g.
and H eff (2) je;je H eff (2) jg;jg Thus,H eff =H eff has the matrix structurẽ where the section shown corresponds to the basis states |je and |(j + 1)g . From this form it is easy to calculate the eigenstates and eigenenergies. The groundstate |0 eff ≡ |0g , which is an eigenstate ofH eff , has the eigenenergy The other eigenstates ofH eff are, j ≥ 0, corresponding to the eigenenergies By construction these are also eigenenergies ofH QHO . Using the transformation (17) we get the eigenvectors ofH QHO as |0 = e −iS |0 eff , |2j + 1 = e −iS |2j + 1 eff and |2j + 2 = e −iS |2j + 2 eff .
The energy spectrum ofH QHO is shown in figure 1 for the case of an unbiased TLS (ε = 0). We want to emphasize that our findings are also valid for the more general case ε = 0. At resonance, where the spectrum for the uncoupled case is degenerated, avoided crossings can be seen. The gap between two formerly degenerated levels for Ω = ∆ b is which is as predicted by the Jaynes-Cumming model [37,38]. As we will show in section 6, the second order correction W 0 in (21), whichleads to a shift in the resonance frequency, is a result of the counter-rotating terms inH QHO . As such it can be interpreted as a Bloch-Siegert shift [45].

Dynamics of the qubit for the non-dissipative case
With the coupling to the bath being turned off, the time evolution of the density matrix of the qubit-HO system is given by ρ(t) = e − i H QHO t ρ(0)e i H QHO t and consequently with ω nm = 1 (E n − E m ). With that (13) becomes where we defined p 0 ≡ n p nn (0). From (29) we notice that the dynamics of the qubit is characterized by an infinite number of oscillation frequencies rather than showing Rabi oscillations with a single distinct frequency. This is clearly a consequence of the coupling of the HO to the TLS. Further we assume that at t = 0 the qubit starts in the state |R and that the occupation numbers of the HO are Boltzmann distributed, so that where Z = e βΩ/2 /(1 − e −β Ω ) is the partition function of the oscillator and β = (k B T ) −1 denotes the inverse temperature of the system. In the TLS-HO eigenbasis this becomes

Low temperature approximation
With (29) we found a formula which describes using the approximate eigenenergies and eigenstates in (24) and (26) the non-dissipative dynamics up to second order in g, thereby taking into account all oscillator levels. Thus, we still have to deal with an infinite Hilbert space. Typically experiments, see e.g. in [13,21], run in a temperature regime for which β −1 Ω, ∆ b . Considering the exponential function in (31) we assume the higher oscillator levels to be only sparsely populated and the maximum value of the sum in (31) is truncated to j = 1. Nevertheless, states |jg/e with j > 1 still play a role in the dynamics. In fact, due to the Van-Vleck transformation exp(−iS), for example the state yields nonvanishing contributions to the matrix elements n|1g and n|1e occurring in (31) due to the fact that the energy eigenstates (26) of the coupled TLS-HO system are made of linear combinations which involve also these states. Using (14b) together with (28) one finds that coefficients p nm (0) with n ≥ 7 are of higher than second order in g. The same is valid for p 50 , p 60 , p 55 and p 66 . Thus, those terms play no role in our calculation of P (t). Furthermore, e − 3 2 βΩ (g/∆ b Ω) 2 ≪ 1. Neglecting also these contributions we find that p n,m ≪ 1 for n ≥ 5. In the end it will be sufficient to concentrate on eigenstates ofH QHO up to |4 . This trunctation leaves us with ten possible oscillation frequencies ω nm , where n, m = 0, 1, . . . , 4 and n > m.
As an example we calculate the dynamics of an unbiased TLS (ε = 0). Here the coefficients p 0 , p 30 (0), p 40 (0), p 21 (0) and p 43 (0) vanish due to symmetry, so that Additionally as a benchmark we consider the mostly studied resonant case, where Ω = ∆ b = ∆ 0 . In this case we find with (24) the transitions frequencies The dynamical quantity P (t) and its Fourier transform are shown in figure 2. One clearly sees the influence of the coupling to the HO on the dynamics of the TLS. Instead of Rabi oscillations with a single frequency, P (t) oscillates with six different frequencies, which are in the Fourier spectrum symmetrically located around the point ω = ∆ 0 . Among those frequencies ω 10 and ω 20 are dominating. They correspond to transitions between the first or second energy level of the qubit-HO system and its groundstate and their weight is almost equal. To summarize, one notices that due to the coupling with the oscillator additional frequencies are induced into the qubit dynamics. Theoretically, the number of those frequencies is infinite. At low temperatures, however, transitions between the lower energy levels of the system are clearly dominating. Again, for simplicity we have shown here the case of an unbiased TLS. For ε = 0 the behaviour is similar only that in the Fourier spectrum the weight difference of the two dominating peaks will be more pronounced.

The influence of the environment
In the preceding section we neglected the influence of the bath on the qubit-HO system. Yet, in order to model a realistic situation, we have to pay attention to environmental influences, as they lead to decoherence and dissipation in the dynamics of the qubit, which is harmful for quantum computing application. Thus, we will now consider the full Hamiltonian H.

Master equation for the qubit-HO system
As shown in section 2.2, we need for the calculation of the population difference P (t) the density matrix ρ(t) of the qubit-HO system. Starting from the Liouville equation of motion for the full density matrix W (t) of H, where the index stands for the interaction picture and following [43,44], we can provide a Born-Markov master equation for ρ(t) being in the Schrödinger picture and expressed in the basis of the eigenstates ofH QHO : The free dynamics of the system is given by the first term of the right-hand side in the above equation. The rate coefficients are defined as For the derivation of the master equation, besides the Born and Markov approximations, some more assumptions have been made, which we briefly mention. First, we consider our system and the bath to be initially (at t = 0) uncorrelated; i.e., and Z B the partition function of the bath. Further, with the bath consisting of infinite degrees of freedom, we assume the effects of the interaction with the qubit-HO system to dissipate away quickly, such that the bath remains in thermal equilibrium for all times t: W (t) I = ρ(t) I ρ B (0). Besides, an initial slip term which occurs due to the sudden coupling of the system to the bath is neglected [2]. And as last approximation the Lamb-shift of the oscillation frequencies ω nm was not taken into account [44].

Matrix elements
In (37) X nm describes matrix elements of the operator X = (B + B † ) in the qubit-HO eigenbasis. By use of (23a), (23b) and (26) those states were expressed in the basis {|jg ; |je }, and we will also calculate the oscillator matrix elements in this basis. For that purpose the operatorX = e iS B † + B e −iS is defined. Four different situations can be distinguished. There are matrix elements were neither the oscillator nor the qubit state are changed, namely jg|X|jg = −2L 0 and je|X|je = 2L 0 with L 0 = εg/∆ b Ω. We see that those elements are independent of j, the occupation number of the oscillator. Next, we look at the case where a single quantum is emitted or absorbed from the oscillator and get jg|X|(j + 1)g = √ j + 1(1 + L osc ) and For a transition within the qubit we have And finally, if the qubit and the oscillator state are changed simultaneously, one obtains jg|X|(j + 1)e = √ j + 1L + q,osc and je|X|(j Comparing the magnitude of the transition terms, we notice that those consisting in changes of the oscillator occupation only are the dominant ones, as they have a part which is of zeroth order in g. Further, for the case in which the qubit is operated at the degeneracy point L 0 and L +/− q,osc vanish. With those results we can calculate the matrix elements X nm . They are given in Appendix C.

Dynamics in the dissipative case
Like in section 3.3 we assume the system to be operated at low temperatures and thus take as highest qubit-HO state the eigenstate |4 . For determination of P (t) the formulas of section 2.2 can be used. Unlike in the non-dissipative case ρ(t) is not given anymore by the simple expression (28). Rather we have to solve a system of coupled differential equations, namely (36). To do this analytically we will follow three different approaches and compare them finally to the numerical solution of (36). We start by introducing which yields the set of differential equations for σ nm :

Full secular approximation (FSA):
As a first approach we make the full secular approximation; i.e., we neglect fast rotating terms in (41) and keep only contributions where ω nm − ω kl vanishes. In this way the off-diagonal elements of σ nm are decoupled from the diagonal ones so thaṫ The equation for the off-diagonal elements is then which becomes with (40) ρ nm (t) = ρ 0 nm e πLnm,nmt e −iωnmt .
As through the FSA the oscillatory motion of the dynamics is separated from the relaxation one we can divide (13) into two parts, where P relax. (t) = n p nn (t) describes the relaxation and P dephas. (t) = n>m p nm (t) the dephasing parts of the dynamics. With (45) the latter takes the form with the dephasing rates Γ nm ≡ −πL nm,nm . Expressions for the dephasing coefficients L nm,nm can be found in Appendix D and the initial conditions ρ 0 nm = σ 0 nm = ρ nm (0) are given by (31). The diagonal elements are more difficult to obtain, as one has to solve a system of coupled differential equations, (42). Calculating the corresponding rate coefficients of this system for the five lowest eigenstates, we find that there are only eight independent ones, namely L 00,11 , L 00,22 , L 11,22 , L 11,33 , L 11,44 , L 22,33 , L 22,44 and L 33,44 . They are given by where j and k adopt the above values. Furthermore, L 00,33 , L 00,44 , L 33,00 and L 44,00 vanish. The remaining rate coefficients are combinations of the above. We find that (50e) However, the system (42) is still too complicated to be solved analytically. Thus, we invoke a further approximation: we consider the factor N nm + 1 = 1 2 [coth( βω nm /2) + 1] with n < m in (49) and use that lim ω→−∞ coth( βω nm /2) = −1. It depends strongly on the temperature β for which value of ω nm this limit is reached approximately. For the parameters we are working with one usually is in the region where (N nm +1) ≪ 1. Thus, we will neglect in the following terms containing the factor (N nm + 1). Furthermore, one sees from (27) that ω 12 ∽ g and ω 34 ∽ g. With that L 11,22 = O(g 3 ) and L 33,44 = O(g 3 ) can be neglected. Using (50a) -(50e) the matrix of the system (42) becomes The eigenvalues and eigenvectors of this matrix and the associated time evolution of the elements σ nn (t) are given in (5.1) -(5.5) of Appendix E. Unlike for the dephasing part (47), we cannot extract a simple analytical expression for the relaxation rate as P relax. (t) = n p nn (t) now consists of a sum of several exponential functions, cf. (14a) together with (5.1) -(5.5). But still we are able to provide an analytical formula for P (t) using (46).

An ansatz for the long-time dynamics:
In order to obtain a simple expression for the relaxation part, we consider the long-time dynamics of the system. In other words, rather than looking at the many relaxation contributions to the populations σ nn (t), we focus on the smallest eigenvalue of the relaxation coefficients, as it will dominate at long times. Further, we consider only the rate matrix associated to the three lowest qubit-HO eigenstates, |0 , |1 and |2 in (42) and obtain with (50a) -(50c) that Here, we have not neglected the rate coefficients containing the term (N nm + 1) and further took L 11,22 into account despite of being of third order in g as such contribution removes the degeneracy between the two lowest eigenvalues at resonance, cf. inset in With the system being detuned this expression can be simplified further, namely The relaxtion rate Γ r as it is given in (53) drawn against the oscillator frequency Ω (solid line). Used values are ε = 0.5∆ 0 , corresponding to a frequency splitting ∆ b = 1.118∆ 0 , and coupling g = 0.18∆ 0 . Moreover, the damping constant is κ = 0.0154 and β = 10( ∆ 0 ) −1 . At resonance (Ω = ∆ b ) Γ r is maximal. For a comparison also the second smallest eigenvalue is plotted (dashed line). The inset shows the two eigenvalues close to resonance.
In figure 3 the relaxation rate Γ r as obtained from (53) is shown as a function of the oscillator frequency Ω. Clearly, it is maximal at resonance (Ω = ∆ b ), whereas it decays with Ω being detuned from the resonance. This effect has already been predicted by Blais et al [12]. As the qubit is not directly coupled to the bath but rather through the oscillator, the latter being detuned filters out the environmental noise at the qubit transition frequency. Additionally, we show the second smallest eigenvalue of (52). We notice that close to the resonant point (Ω = ∆ b ) there is an avoided crossing. Finally, we find that where like in section 3.2 p 0 ≡ n p nn (0). For getting p ∞ we have in principle to find the steady-state solution of (42). Here, we just assume for t → ∞ a Boltzmann distribution for the qubit-HO system, so that ρ nn (∞) = Z −1 QHO e −βEn with Z QHO = n e −βEn . Thus, The formula for the long-time dynamics is obtained as To get further insight on the dominant frequencies we evaluate the Fourier transform of (57) according to yielding (59)

Partial secular approximation (PSA):
An improvement to the FSA is to take into account certain non-vanishing contributions of ω nm − ω kl . We have to keep in mind, that there are quasi-degenerate levels close to resonance. In our case the first with second energy level and the third with fourth one build a doublet, meaning that they are close together in energy space. The level spacing is approximately proportional to g for the former and √ 2g for the latter. Because of that and as the transitions from level three and four are less probable, we will in the following only consider the first and second level as being almost degenerate. Taking this into account in (41) we arrive for the diagonal elements aṫ σ nn (t) = π k L nn,kk σ kk (t) + πL nn,12 σ 12 (t)e −iω 12 t + πL nn,21 σ 21 (t)e −iω 21 t .
A numerical analysis shows that the effect of the last two terms on the right-hand side of the above equation will in the worst case lead to very small wiggles in σ nn (t) and play no role in P (t). Thus, we finally writė which is the same equation as we got in the FSA approach. However, the off-diagonal contributions σ 01 , σ 02 , σ 13 , σ 23 , σ 14 and σ 24 have to be examined more carefully. From (41) we find that one has to solve the equationṡ Here, the oscillation frequencies and the decay of the off-diagonal elements are given by with The amplitudes of the oscillations are given through the coefficients and Thus, we have again all ingredients to calculate analytically the relaxation and dephasing part of (46). For the PSA we cannot provide a simple expression for the dephasing rates as in the FSA, where we had Γ nm = −πL nm,nm . As one can see from (64)  Similar to our findings for the relaxation rate, also here the smallest eigenvalue will dominate the dephasing behaviour. From the right graph in figure 4 we find that this is ℜ{λ 12 . Comparing it with the dephasing rates we got using the FSA, left graph in figure 4, we see that for negative detuning (Ω < ∆ b ) the rate Γ 02 = −πL 02,02 approximates Γ (+) 12 , whereas for positive detuning (Ω > ∆ b ) this is done by Γ 01 = −πL 01,01 . In the FSA Γ 02 and Γ 01 correspond to the frequencies ω 10 and ω 20 , respectively. In the PSA the frequency ω  figure 4 one notices further that the FSA rate Γ 02 grows linearly with Ω for positive detuning. However, as the weight of the corresponding frequency ω 20 will be almost zero, Γ 02 will give no relevant contribution to P dephas. (t) in this regime but the dephasing will rather be associated to the FSA rate Γ 01 . Hence, out of resonance the FSA will still fairly well describe the dynamics of P (t). Comparing the expressions for L 01,01 and L 02,02 given in Appendix D by (4.1) and (4.2) with the approximative expressions for the relaxation rate at positive and negative detuning (54), we see that for zero bias (ε = 0) the PSA dephasing rate is equal to Γ r /2. For a biased system an additional term is added depending on the spectral density of the bath at ω = 0. In figure 5 we compare the three analytical solutions described  above to the numerical solution of the master equation for the case of an unbiased TLS being at resonance with the oscillator. Concerning both the dynamics of P (t) and its Fourier spectrum we see a good agreement between the different solutions. The one being closest to the numerical solution is the PSA solution. We also want to mention that going to stronger damping κ, the FSA results start to show deviations from the numerical solution. Here, one should use the PSA only. However, for the parameter regime used in the following, we will mainly apply (57) due to its simple, analytical form.

Discussion of the results
Having solved the master equation (36) analytically and numerically we can examine the dynamics of the system and its Fourier transform for different situations. First, we will look at a qubit operated at the degeneracy point (ε = 0) being in and out of resonance with the oscillator. Then, we will concentrate on the biased qubit in the same regime of parameters.

The unbiased qubit
For unbiased qubits we can compare our predictions with the analytical results obtained in [31] by starting from a spin-boson model with the effective spectral density (10). In [31] a so-called weak damping approximation (WDA) based on the non-interacting blip approximation (NIBA) is applied. The WDA allows a non-perturbative treatment of the coupling between the TLS and HO and hence can reproduce the occurence of two dominating frequencies as expected e.g. from exact QUAPI calculations [27]. The NIBA, and hence the WDA, however, become not reliable for a biased TLS. We find that the overall agreement between our approach and the WDA is very good. However, in the WDA solution the frequencies are slightly shifted compared to the ones obtained from our master equation. This may result from the perturbative expansion we have performed with respect to g by applying the Van-Vleck perturbation theory. First, we look at the resonant case shown in figure 6. In agreement with previous works [27,31], we find that the dynamics is dominated by two frequencies corresponding to ω 10 and ω 20 with separation being approximately 2g. The weight of the latter is a bit larger. The reason for the bigger weight is that at resonance (Ω = ∆ b ) the qubit-HO eigenstate |j is not a symmetric superposition of the states |j, e and |j + 1, g unlike it is predicted by the Janyes-Cummings model (cf e.g. [12]). We notice that the two unequal peaks have indeed been experimentally observed in [13] (see Fig. 4b  therein). Considering the states |1 eff and |2 eff in (23a) and (23b), one already sees that for a symmetric superposition of these states we need that δ 0 vanishes or that (78)). Besides, in order to get the qubit-HO eigenstates one still has to perform the Van-Vleck transformation, which adds contributions to |1 and |2 from states corresponding to oscillator levels higher than j = 1. Thus, our system behaves for Ω = ∆ 0 as being negatively detuned, which means that the peak belonging to the higher frequency dominates, as we will show below. Slightly increasing Ω will give a stronger weight to the peak at ω 10 . This effect is not very pronounced for the non-dissipative dynamics of the unbiased qubit (figure 2), as there the two frequencies are still almost equally weighted. Looking however at the Fourier transform of the dissipative dynamics (59), one notices that the relaxation rate also contributes to the weight of the peaks with a prefactor Γ −1 nm . As for a negative detuned system Γ 01 is slightly bigger than Γ 02 , the difference between the two peaks becomes more clear in the dissipative case. For ε = 0 the effect can already be noticed in the non-dissipative case. Next, we consider in figure 7 the case of negative detuning, where Ω < ∆ 0 . No matter which approach one is looking at, clearly the frequency ω 20 is dominating. Furthermore, paying attention to the timescale of the dynamics, one notices that the relaxation time is enhanced compared to the one we found for the resonant system. This behaviour was already explained by the formula (53) for the relaxation rate. Again, the numerical and the solution obtained by using the long-time ansatz in section 4.3.2 agree quite well with each other, whereas the amplitude of the oscillation with frequency ω 20 is stronger in the WDA approach. Also remarkable is the fact that looking at the Fourier transform in figure 7 one sees in the inset already small contributions of the higher oscillator levels.
The transitions corresponding to ω 24 and ω 23 give raise to small additional peaks, while the contributions of ω 13 and ω 14 are negatively weighted and cause dips. The WDA approach does not show this additional contributions. They are, however, confirmed by the numerical QUAPI calculations in [27] (see figure 2 therein). In the case of positive detuning (Ω > ∆ 0 ) shown in figure 8 we find a quite good agreement between all three approaches. Also for postive detuning the relaxation time is enhanced compared to the resonant case. Contrary to the negatively detuned situation the additional peaks have  vanished. Besides, now the frequency ω 10 is dominating the dynamics. This behaviour, namely that for negative detuning ω 20 and for positive detuning ω 10 is dominating, was already found in [27]. We will briefly explain how one can explain this observation physically. For this we look at figure 9. For a detuned system (Ω = ∆ b ) the qubit-HO eigenstates are not symmetric superpositions of the states |jg and |je . They rather asymptotically approach the eigenstates of the uncoupled qubit-HO Hamiltonian. In figure 1 we see that for a negatively detuned system (line a) the qubit-HO eigenstate |2j + 1 approaches the state |(j + 1)g , whereas the main contribution to the state |2j + 2 will come from the  state |je . From the left diagram in figure 9 we see that the state |2 is energetically higher than the state |1 . However, due to the Boltzmann distributed occupation of the oscillator, the state |0e will be more populated than the state |1g and consequently also |2 will exhibit a larger population than |1 , as the latter only feels a small contribution from the state |0e . Thus, transitions from |2 to the groundstate are more likely to occur than those from |1 to the groundstate. This explains the dominance of ω 20 in figure 7 and figure 11. In this case the frequency ω 20 ≈ ∆ b and ω 10 ≈ Ω. As far as not excluded by selection rules, minor peaks from transitions to the levels lying in between can be also seen. For positive detuning (line c in figure 1) |2j + 1 approaches |je , while |2j + 2 is close to |(j + 1)g . From the right graph in figure 9 we see that the state |1 , being lowest in energy apart from the groundstate, is now also more probable to be occupied than |2 . Therefore, as confirmed by figure 8 and figure 12, the frequency ω 10 is dominating whereas ω 20 is represented only by a small peak in the Fourier spectrum. Furthermore, as there are no additional energy levels between the state |1 , which is most probably to be populated, and the ground level, other transitions than those corresponding to ω 10 or ω 20 are very unlikely to occur. In figure 12 the dip corresponding to ω 21 appears only very faintly.

The biased qubit
We will now examine a qubit being operated at finite bias. We consider the case ε > 0. For negative bias-offset the behaviour is analogous. Again three different situations are taken into account: the qubit being in resonance with the oscillator (∆ b = Ω), negative (Ω < ∆ b ) and positive (Ω > ∆ b ) detuning. For the resonant case (Ω = ∆ b ) depicted in figure 10 we see a similiar behaviour as for the unbiased qubit. Again two frequencies, ω 10 and ω 20 , are dominating the dynamics. Left to the peak at ω 10 a small dip can be found in the Fourier spectrum. This corresponds to the transition ω 21 . For infinite time the dynamics relaxes to an equilibrium value which is nonzero in contrast to the unbiased case. This can be seen in the Fourier spectrum through a relaxation peak at ω = 0. The peak arises because of the term in (59). The first part of this sum gives rise to the negative shift of this peak. The reason that for the analytical solution the peak is not as strongly shifted as for the numerical one is technical: in order to plot the delta function in (70) we gave it a finite width, which surpresses the negative contribution of the first term in (70). Like for the unbiased qubit the highest energy level playing a role for the dynamics is E 2 ; i.e., only the ground and first excited level of the oscillator are of importance. In figure 11 the dynamics and its Fourier transform for a negatively detuned qubit-oscillator system with ε = 0 are shown. Like for the unbiased case detuning gives raise to longer relaxation times for the qubit. Also in agreement with the unbiased case is the dominance of the frequency ω 20 . We see that for small t the long-time solution (57) slightly overestimates the maxima of the oscillations and underestimats its minima. Furthermore, we get here the unphysical situation that the maximum of the third oscillation in P (t) exceeds the value of one. The reason for that behaviour is that, by construction, we underestimate with (57) the relaxation at short times. As for certain paramteres the term (p 0 − p ∞ ) in (57) can become negative, it increases too fast towards the equilibrium and gives thus raise to the observed deviations in the short time behaviour. On a longer timescale both graphs agree quite well. For the case of positive detuning (Ω > ∆ b ), which is presented in figure 12, the upward shift of the dynamics obtained from (57) and (59) compared to the numerical graph of P (t) at small times is even stronger. To visualize that it is not a failure of the FSA approach we show in figure 11 and figure 12 additionally the analytical FSA solution of (42) and (43)   To conclude this paragraph we want to mention that all the results found both for the unbiased and the biased qubit confirm the numerical QUAPI results in [27].

Symmetrized correlation function
So far we have always considered the qubit for certain values of ε and finite or zero detuning. In this section, we fix the oscillator frequency at Ω = ∆ 0 . That means that an unbiased qubit will be at resonance with the oscillator. Changing the bias to positive or negative values will always lead to negative detuning, as ∆ b ≥ ∆ 0 . Figure 13 shows a density plot of the Fourier transform of the symmetrized correlation function against the bias of the qubit and the Fourier frequency ω. We consider this correlation function rather than P (t), as it is symmetric in the bias ε. The symmetrized correlation function is defined as follows [2]: where σ z (t) = e iHt/ σ z e −iHt/ . Expressed in terms of the population difference P (t) this becomes, with P s (t) and P a (t) being symmetric and antisymmetric in ε and P (t) = P s (t) + P a (t). The Fourier transform of S(t) is defined as Considering now figure 13 we see that for any bias the spectrum is dominated by two frequencies, namely ω 10 and ω 20 . Detuning the system ω 20 gets more and more important, as we could already observe in the two previous sections for the positively detuned systems. Furthermore, the peaks are shifted to higher frequency values and at ω = 0 the relaxation peak occurs. We want to compare these results to a circuit QED experiment performed by Wallraff et al [13]. There the qubit is realized by a Cooper pair box, which is coupled to a superconducting transmission line resonator. The properties of the system are determined by probing the resonator spectroscopically. The amplitude of a microwave probe beam transmitted through the resonator is measured versus the probe frequency and the gate charge of the Cooper pair box (see figure 4 in [13]). Via the gate charge the qubit can be detuned in situ from the degeneracy point. The frequency of the resonator is chosen in such a way that it is in resonance with a qubit being operated at the degeneracy point. For the resonant case two dominating frequencies, being almost equally weighted and symmetrically positioned around the cavity frequency, are observed. Going away from the degeneracy point the system becomes detuned and the frequency of the cavity dominates. The behaviour we observe in figure 13 is similar. However, as we are looking at the dynamics of the qubit, it corresponds to a spectroscopic measurement on the TLS rather than on the oscillator. As explained above the two lowest excited states of the coupled TLS-HO system, namely |1 and |2 , evolve from an almost symmetric superposition of basis states {|jg , |je } at resonance (ε = 0) to the states |1g and |0e (cf the left graph in figure 9). For ε = 0 the two peaks of the Rabi splitting are observed. For ε = 0, which means negative detuning in this case, the peak with the lower frequency corresponding to ω 10 approaches more and more the frequency Ω of the oscillator, as the state |1 becomes |1g for large detuning and then ω 10 ≈ E |1g − E |0g = Ω. Furthermore, the transition peak at ω 10 gets weaker as also the occupation probability of |1 decreases. At ε ≈ ±0.8∆ 0 the symmetrized correlation function vanishes at ω 10 and increases again for higher values of |ε|. Here, the amplitude p 10 in P (t) changes its sign. In contrast the peak at ω 20 becomes stronger with the detuning and approaches more and more the qubit splitting energy ∆ b , as |2 approaches |0e and then ω 20 ≈ E |0e − E |0g = ∆ b . Additionally, looking at the logarithmic plot one sees around ω = 0.4∆ 0 a peak appearing, which corresponds to the frequency ω 21 and is forbidden at ε = 0. For large detuning it arises from transitions from |0e to |1g and therefore has the value ω 21 ≈ ∆ b − Ω. The amplitude of this peak is very small compared to the peaks at ω 10 and ω 20 and is not resolved in the experiment of Wallraff et al .

Comparison with the Jaynes-Cummings model
Van-Vleck perturbation theory enabled us to find approximately the eigenstates and eigenenergies of the full Hamiltonian of the TLS-HO system without performing a rotating-wave approximation. Using those eigenstates and eigenenergies in a Born-Markov master equation we could calculate the dynamics of such a system under the influence of an environmental bath. In the following we will show how the results change if we neglect counter-rotating terms in the TLS-HO Hamiltonian (16) for ε = 0. For this we rewrite the interaction part in (16) as where we identified withH R Int andH CR Int a rotating and counter-rotating part ofH Int , respectively, and introduced the two-level transition operatorsσ ± = 1 2 (σ x ± iσ y ). Neglecting the counter-rotating partH CR Int inH QHO leads to the Jaynes-Cummings This Hamiltonian can be diagonalized exactly and its eigenstate and eigenvalues can for example be found in [46]. In order to see the effect of not taking into account the counter-rotating terms we diagonalizeH JC using Van-Vleck perturbation theory. Looking at the formula fo the effective Hamiltonian (2.6) in Appendix B and keeping in mind that we set ε to zero, we see that the second order contributions in g vanish neglectingH CR Int ; i.e., W 1 , W 0 are zero in (21). Considering further the transformation matrix S, equations (2.4) and (2.5) in Appendix B show that S = 0 forH CR Int = 0. Thus, with (17) we find that for the Jaynes-Cummings modelH eff is identical toH JC and therefore the eigenstates ofH eff are simultaneously eigenstates ofH JC . Consequently, one can determine from (22) -(24) the eigenstates and eigenenergies ofH JC . The energy of the groundstate |0 JC = |0g JC is E JC 0 = − ∆ b /2. For the higher states we get corresponding to the eigenenergies with δ JC = ∆ b − Ω and tan α JC j = 2 √ j + 1|∆|/δ JC . Comparing these eigenstates and eigenenergies to the ones found forH QHO , (23a), (23b) and (26), we see that the counterrotating terms yield second order corrections in g not present in the Jaynes-Cummings Hamiltonian. These corrections give rise to a very prominent effect concerning the resonance condition between TLS and HO. From δ JC we find the TLS being in resonance with the oscillator for Ω = ∆ b . Considering δ j in (25) this resonance condition is shifted to This second order correction to the resonance frequency due to counter-rotating terms is known as Bloch-Siegert shift [45]. The eigenstates (76a) and (76b) are always a superposition of two basis states of the unperturbed system. This is like for the eigenstates (23a) and (23b) of the effective Hamiltonian (17). However, in order to find the eigenstates |n ofH QHO we had to apply the transformation exp(−iS) on the effective eigenstates so that |n is in the end a superpostion of several states of the basis {|jg ; |je }.
To calculate the reduced density matrix of the qubit-HO system described byH JC and taking into account the influence of an environmental oscillator bath we can again use the Born-Markov master equation (36). We just have to use (76a) and (76b) as basis states and the corresponding eigenenergies. For the population difference P (t) we apply (13), which becomes for ε = 0

Selection rules
Before doing a qualitative comparison between the results obtained from the original HamiltonianH QHO and the simplified formH JC , we want to analyse which transitions between the different eigenstates ofH JC yield contributions to P (t). For this it is helpful to rewrite (76a) and (76b). In the following we neglect the upper index JC denoting an eigenstate ofH JC . For the state |n we have three possibilities: the first one corresponds to n = 0. In this case the only non-vanishing component of |0 is 0g|0 . Second, n can be an even number, which means, expressed in terms of oscillator quanta, that n = 2j + 2. Then And third for an odd state n = 2j + 1 we find |n od = cos(αn−1 2 /2)| n + 1 2 g + sin(αn−1 It is quite easy to see that the part p nn (t) in (79) vanishes for any n. That means that the diagonal elements of the reduced density matrix yield no contributions to P (t) for ε = 0 and the equilibrium value of the dynamics will be also zero. Not so easy to see are the combinations of n and m yielding contributions to the off-diagonal part p nm (t). Considering transitions to the groundstate (m = 0), we find using that 0g|0 = 1 that p nev0 = 2 0e|n ev ℜ{ρ nev0 (t)} and p n od 0 = 2 0e|n od ℜ{ρ n od 0 (t)}.
With (80)  degenerate levels are forbidden. This behaviour we have also found in section 3.3 for the non-dissipative dynamics resulting fromH QHO .

Comparison of the two models
In (Ω < ∆ 0 ), which is shown in figure 16 and where ω 20 dominates, the equilibrium value is reached too fast within the Jaynes-Cummings approach. Furthermore, considering the graph of the Fourier transform we find that small contributions which come from higher level transitions, and which have already been discussed in section 5.1, are not caught be the Jaynes-Cummings approach. For the resonant case (Ω = ∆ 0 ) we find that the Jaynes-Cummings method predicts ω 10 to be slightly dominating whereas the approach starting from H QHO results in ω 20 being dominating. The reason for this discrepancy is that due to the counter-rotating terms we have for H QHO no symmetric or antisymmetric superposition of the unperturbed eigenstates at Ω = ∆ 0 in contrast to the Jaynes-Cummings model. To conclude this section we can say that for an unbiased TLS-HO system the Jaynes-Cummings model gives a good insight in the qualitative behaviour of P (t) both for a slightly detuned and a non-detuned system. However, it under-or overestimates dephasing times for the system. Furthermore, we find taking into account counterrotating terms inH QHO that at Ω = ∆ b the dressed eigenstates are not a symmetric or antisymmetric superposition of the uncoupled states. Moreover, the effects of transitions between states of different manifolds are neglected.

Conclusions
In conclusion, we discussed the dynamics of a biased and unbiased TLS coupled through a harmonic oscillator to an environmental bath described by an Ohmic spectral density. In particular, we examined the regime of weak damping and moderate coupling between oscillator and TLS. An equivalent description of our system is provided by the spinboson model with a structured spectral density. In contrast to many other works in this field, our starting point was not the Jaynes-Cummings Hamiltonian, but a more general one given in (1), where no initial rotating wave approximation has been applied. In section 2.2 we provided with (13) a formal expression for the population difference and showed in section 3.1 how the Hamiltonian of the coupled qubit-HO system can be diagonalized approximately using Van-Vleck perturbation theory. This approach is valid both for the oscillator being in resonance with the qubit (Ω = ∆ b ) and for finite detuning (|δ| = 0). In section 3.2 an analytical expression for the non-dissipative dynamics was provided up to second order in the qubit-HO coupling g taking into account the infinite Hilbert space of the system. For low temperatures (k B T < Ω, ∆ b ) we truncated the Hilbert space and found that transition processes between the groundstate and the two first excited energy levels of the qubit-HO system dominate the dynamics (section 3.3). In section 4 the influence of the bath was taken into account by solving the Born-Markov master equation for the density matrix of the qubit-HO system. To do this analytically we considered two variants of the secular approximation: first, in section 4.3.1, the full secular approximation (FSA), where all fast oscillating terms are neglected, and second, in section 4.3.3, the partial secular approximation, where attention was paid to the fact that the first two excited energy levels are almost degenerate. Using an ansatz for the long time dynamics in section 4.3.2, we could provide a general expression for the relaxation and dephasing rates of the qubit, showing that the relaxation time can be enhanced by detuning the oscillator into the off-resonant regime. It was found that all three approaches agree quite well with the numerical solution of the master equation. The dynamics of both a biased and an unbiased qubit were intesively studied for zero and finite detuning in section 5. The results agree qualitatively with the numerical findings within the ab-initio QUAPI approach [27]. Furthermore, in section 5.1 a good agreement with the results of the weak damping approximation performed in [31] for a symmetric spin-boson model was found. Besides, we saw that at resonance (Ω = ∆ b ) the first two excited qubit oscillator states are not a symmetric or antisymmetric superposition of the states |0e and |1g , as predicted by the Jaynes-Cummings model, and thus could give an explanation for the differently weighted peaks of ω 10 and ω 20 in the Fourier spectrum of the dynamics. We further could explain the dominance of frequency ω 20 in the case of negative detuning (Ω < ∆ b ) and of frequency ω 10 for positive detuning, respectively. Moreover, we showed that for large negative detuning ω 20 approaches the energy splitting ∆ b of the qubit, whereas ω 10 approximates the oscillator frequency Ω. This behaviour agrees nicely with spectroscopic experiments performed on a circuit QED architecture [13]. In section 6 we compared our results for an unbiased system to the ones obtained starting with the Jaynes-Cummings model. We visualize the effects of the counter-rotating terms in the Bloch-Siegert shift of the resonance frequency and in contributions of states with larger oscillator number to the TLS dynamics. Apart from this the Jaynes-Cummings model and our approach agree quite well for the nondissipative case. Also for the dissipative case at resonance an initial RWA represents a good approximation to our starting Hamiltonian and seems to be favourable as it is analytically exactly diagonalizable. For detuned systems, however, we find discrepancies concerning relaxation and dephasing times and a RWA becomes less appropriate to give precise results. Thus, we think that our approach represents an improvement as it is valid in a wider parameter range avoiding an initial rotating wave approximation. To our knowledge it provides for the first time analytical results for the dynamics of an unbiased and biased qubit coupled to a structured environment being valid both in the resonant and off-resonant regime. Furthermore, due to the generality of the qubitoscillator model, we expect our results to be of interest for a wide range of experimental applications.
In the case of the HamiltonianH QHO the first order matrix elements are iS (1) e j−1 e j = j ε ∆ b Ω g, (2.7a)