Dynamic phases induced by two-level system defects on driven qubits

Recent experimental evidences point to two-level defects, located in the oxides and on the interfaces of the Josephson junctions, as the major constituents of decoherence in superconducting qubits. How these defects affect the qubit evolution with the presence of external driving is less well understood since the semiclassical qubit-field coupling renders the Jaynes-Cummings model for qubit-defect coupling undiagonalizable. We analyze the decoherence dynamics in the continuous coherent state space induced by the driving and solve the master equation endowed with an extra decay-cladded driving term via a Fokker-Planck equation. The solutions for diffusion propagators as Gaussian distributions show four distinct dynamic phases: four types of convergence paths to limit cycles of varying radius by the distribution mean, which are determined by the competing external driving and the defect decays. The qubit trajectory resulted from these solutions is a super-Poissonian over displac ed Fock states, which reduces to a Gibbs state of effective temperature decided by the defect at zero driving limit. Further, the Poincare map shows the dependence of the rate of convergence on the initial state. In other words, the qubit evolution can serve as an indicator of the defect coupling strength through the variation of the driving strength as a parameter.


I. INTRODUCTION
Environmental coupling causes decoherent evolution of qubits.For superconducting qubits, decoherence is especially acute because of their inevitable couplings to the surrounding materials in the solid-state circuits in addition to the vacuum coupling [1].Since finite times are required for performing quantum logical operations, improving decoherence times has been a major field of study for more than a decade [2], with experiments demonstrating T 1 and T 2 times on the order of 10µs [3] and 100µs [4], respectively.Nevertheless, characteristics of decoherence are often precarious and sometimes highly random in fabricated superconducting circuits and an analytical understanding of the sources for decoherence is still required.
Many studies are directed towards sourcing the material origin of qubit decoherence [5], which point to the defects [6,7] located broadly in the interior and on the edges of the sandwiched oxides of the tunnel barriers as well as chip surfaces [8,9].Spectroscopic analyses show that the defects, growing out of the amorphous structure from natural or slow oxidation process, can be modelled as two-level systems [5].In particular, these two-level defects are in general categorized into two types: fluctuators, with level spacing ∆E k B T , and coherent two-level systems (cTLS), with level spacing ∆E ≫ k B T ; the former (latter) are thus strongly (weakly) coupled to the environment at temperature T [10].Collectively, cTLS and fluctuators contribute a 1/f background noise spectrum [13] and their interaction would decide the power dependence of the quality factor of superconducting resonators [14].Yet the greater energy gap places cTLS spectroscopically near the qubit, making the qubit effectively more susceptible to the motions of the cTLS than those of the fluctuators.These coherent two-level defects are thus recognized as the dominant source for relaxing a qubit's coherent dynamics [11,12] and methods to alleviate such relaxations have been proposed [15,16].
Recent experimental progress has characterized the qubit relaxations induced by cTLS [17][18][19][20][21], where the theoretical analyses are modeled on the interactions between Pauli operators (for cTLS) and anharmonic oscillators (for qubits) [22][23][24].Superconducting qubits including the prevalent transmon variants are indeed anharmonic multi-level systems whose higher excited states are often ignored to simplify the discussion of their roles in quantum computations.However, when relaxations towards defects are taken into account, the contributions from the higher levels become significant and must be considered [12].We note that in optical amorphous materials, the defect couplings are described by non-Hermitian Hamiltonians [25].
Since relaxations occur concurrently with driving and their rates are highly dependent on the driving strength [26][27][28][29][30][31][32][33], characterizing the decoherence dynamics under driving becomes necessary and nontrivial because the semiclassical qubit-field coupling removes the discrete diagonalizability of the qubit-defect coupling.The effective dressed relaxation rates are first computed in Ref. [34].
Here, we resort to a solution of the evolution on the continuous complex (α, α * )-plane of coherent states by regarding the driving as a displacement that translates the creation and annhilation operator pair of the anharmonic qubit.Since spectral analysis shows the noise additional to the 1/f background contributed by the fluctuators stems from a single two-level system, we register a single perturbative term provided by the cTLS in the master equation.The qubit decoherence evolution is then fully solved by writing its density matrix in the Q-representation, whence the master equation is replaced by a Fokker-Planck (FP) equation.The solution Q(α, α * , t), obtained by solving two diffusion-type coupled equation of real variables derived from the FP equation, obeys a Gaussian distribution of moving mean and expanding variance.
We find that the steady-state of the moving mean follows a limit cycle determined by both the driving field and the decays induced by the cTLS.The steady state regarded as a thermal state has the qubit thermal inversion decided by the cTLS level spacing in addition to the environmental temperature T .Given an initial point state on the quadrature plane, the radii of the asymptotic limit cycles can be either greater or less than the distance to origin from that initial point, depending cTLS in bath

Driving Field
Figure 1.The system consists of an anharmonic oscillator of base level spacing ωq and anharmonicity χ that models a superconducting qubit, a two-level system of spacing ωc that models a near-resonant defect in the qubit junction, and a microwave driving of frequency ω d to the qubit.The two-level defect is considered interacting with a thermal bath of temperature T at a decay rate of γc and simultaneously interacting with the qubit at strength g.The latter interaction provides a two-way channel of energy flow during the qubit thermalization process with the bath.on both the qubit effective decay mediated by the cTLS and the driving strength.The transient behaviors towards these limit cycles is even more sensitive to the driving condition: it can generally be classified into four dynamic phases of distinct converging behaviors according to the competing driving and decay.This sensitivity translates into an indicator of the decay strength induced by the cTLS under the varying driving strength as a parameter.To understand how these dynamic phases arise, we begin with the discussion of the tripartite interaction model below in Sec. 1, in which the master equation is introduced, and follow it by the process of solution in Sec.III.The classification of the dynamic phases is presented in Sec.IV.The transient dynamics is analyzed in Sec.V before the conclusions are given in Sec.VI.

II. QUBIT-CTLS COUPLING UNDER DISPLACEMENTS
The total Hamiltonian ( = 1) has four terms, accounting for the three comprising parts of the system and the environmental interactions as illustrated in Fig. 1.The first term H q = ω q a † a − χa †2 a 2 denotes the free energy of the qubit with resonance frequency ω q and anharmonicity χ.The anharmonic oscillator model arises from the double-well potential that are quantized from the junction and capacitor energy of the current and voltage in the solid-state based system [36,37].The second term H c = ω c σ z /2 accounts for the free energy of the cTLS.The last two terms denote the interactions, for the qubit-cTLS coupling H int = g a † σ − + aσ + and for the semiclassical field driving W (t) = Ω a † e −iω d t + h.c., where Ω is the driving strength and ω d the field frequency.
Considering the typical scenario where Ω ≫ g (Ω in the range of MHz while g in the range of hundreds of kHz for typical transmon qubits), the driving is relatively strong, compared to which the coupling to cTLS is perturbative.The discrete qubit state space can then be transformed into the continuous coherent-state space under the rotating frame of the driving ω d , where Ω and the qubit-drive detuning ∆ d = ω ′ q − ω d determine the amount of displacement α d = Ω/∆ d from the vacuum state, i.e.
Given weak driving Ω, the magnitude of the displacement is small, which renders the qubit to fall into the weak-excitation regime [38], where the oscillator frequency can be linearized as ω ′ q = ω q − χ(n q − 1).The number nq denotes the average photon number in the qubit, which remains a constant value since the coupling strengths g, Ω ≪ ω.By regarding DH q D † and H c as free energies in the coherent-state space, the effective qubit-cTLS coupling under the interaction picture becomes where ∆ c = ω ′ q − ω c denotes cTLS detuning from the linearized qubit.
Given Eq. ( 2), the master equation of the qubit density matrix ρ q is derived from the Liouville equation ρ = −i [H eff , ρ] for the coupled system ρ.We apply the Born-Oppenheimer-Markov approximation since the motion of the cTLS is relatively much slower than the fast motion of the qubit.Experimental studies [11,17] have shown that cTLS-induced variations in frequency and dephasing rate of the qubit occur over the scale of hours to days; the cTLS evolution is therefore not perturbed by the weak qubit-cTLS coupling over microsecond time scale considered here.Exactly, the density matrix is assumed separable, i.e. ρ(t) = ρ q (t) ⊗ ρ c , at initial moment before the time integration is carried out to the first perturbative order and the cTLS subsystem is traced out.As shown in Appendix A, this leads to in the displaced (α, α * )-space.Besides the two Lindbladians contributed by direct coupling to the cTLS, the indirect coupling skewed by the driving field contributes the extra η term.With ν = (exp(ω c /k B T ) + 1) −1 indicating the thermal inversion of the cTLS, the Lindbladians accounts for the spontaneous qubit emission (the third term) and absorption (the last term), where γ = 2g 2 γ c /(γ 2 c + ∆ 2 c ) is the decay rate according to the typical relaxation rate γ c of a cTLS [39].The η term reflects the external driving under the influence of the environment, i.e.
) which comprises two terms.The first term of Eq. ( 4) is contributed by a two-photon process where the driving field photon and the cTLS radiation photon concurrently mix with the qubit.Whereas, the second term corresponds to a one-photon process where the driving field photon bypasses the qubit and mixes with the cTLS.The resonance at the latter blocks the decay of the qubit that channels into the cTLS.It therefore adversely affects the decay depicted in Eq. ( 3) and has the negative coefficient.On the other hand, both terms share the common factor α d determined by the driving amplitude.Given the environmental temperature (on the mK range) of superconducting qubits, the inversion v adopts a negligible value and the existence of the driving field in general accelerates qubit relaxations (contribution of the first term ≫ that of the second term).

III. SOLVING FOR QUADRATURES
To solve the master equation ( 3) analytically, we first reduce it to a Fokker-Planck (FP) equation concerning the density distribution of the qubit on the complex α-plane.Under the Q-representation Q(α, α * , t) = α| ρ q (t) |α /π, the FP equation reads where the square bracket indicates an Hermitian differential operator and ν = 1/2 − ν indicates the inversion difference between 1/2 at infinite temperate and ν at a given temperature [40].To arrive at an analytical solution to Eq. ( 5), we first rewrite the differential equation in the real quadrature plane under a rotating frame: x = (αe i∆ d t + h.c.)/2 and y = −i(αe Since the initial distribution Q(x, y, 0) is not necessarily separable along the x and y variables, we consider the associated equation of motion for the propagator (viz.the Green's function) G(x, y, t; x ′ , y ′ , 0) linking the initial coordinate (x ′ , y ′ ) and the current coordinate (x, y) at current time t [40].Since it reduces to the separated δ(x − x ′ )δ(y − y ′ ) at the time limit t → 0, the propagator is separable as G x (x, t; x ′ , 0)G y (y, t; y ′ , 0) in general, breaking up its equation of motion into the pair Each equation above is essentially equivalent to a heat diffusion equation, where the first-order derivatives and the linear terms on the right hand side can be eliminated under a new space variable with the time varible transformed accordingly [41].The diffusion equations are analytically solvable using standard Fourier transforms (Cf.Appendix C).Eventually, reversing all the transforms and changes of variables, one obtains the Q-distribution at time t from an initial coherent state distribution under the joint propagator G x G y in the complex plane: which is Gaussian with a moving mean µ(t) = µ ss exp(i∆ d t) + (α 0 − µ ss ) exp(−νγt), assuming α 0 as the displacement of the coherent state.The mean converges to the asymptotic coordinates at steady state, which depends on the external driving as well as the qubit decay.The variance σ 2 = [1 − exp(−2νγt)] ν/4ν + 1/2 is expanding and convergent with time.
To verify the validity of the analytical approach, we employ the matrix exponentiation method to numerically solve the master equation in a superoperator form [42] and find that the discrepancy between the analytical |µ(t)| and the numerical |µ(t)| is marginal, a difference on the order of 10 −4 over unital order of magnitude (Cf.Appendix F for comparison details).

IV. CLASSIFYING EVOLUTION
The evolution depicted by Eq. ( 9) dictates that the qubit, starting from an arbitrary (pure) coherent state |α 0 for which Q(0) ∼ N (α 0 , 1/2), follows a typical trajectory of decay into the quasi-coherent state ρ ss for which Q(t → ∞) ∼ N (µ ss , ν/4ν + 1/2) independent of the initial state |α 0 .The mean µ ss depends only on the parameters of the driving signal and the decay parameters of the cTLS.In fact, if the driving vanishes, the equilibrium distribution Q ss reduces to that of an equivalent thermal Gibbs state [43], which reads exp −|α| 2 /(n + 1) /π(n + 1) in the continuous α-plane with n = ν/2ν = (exp{ω c /k B T } − 1) −1 indicating the equivalent average photon number.In other words, the qubit levels would be populated not according to ω ′ q , but rather to the level spacing ω c of cTLS regarded as an infinite-level harmonic oscillator.One therefore can follow the model in Fig. 1 and regard the qubit-cTLS interaction as a thermalization channel, for which the qubit is thermalized with a cTLSmediated bath at the effective temperature T q = T ω ′ q /ω c .Given the ultra-low bath temperature T of the dilution fridge, the difference in qubit population inversions under T and T q is significant, despite the near-resonant frequencies of ω ′ q and ω c .This temperature difference is recently measured [44]   is useful for determining ω c [39,45].Moreover, the widened variance of σ 2 ss by ν/4ν compared to σ 2 (0) also stems from the thermal energy transferred from the cTLS to the qubit.
For the scenario of a finite driving strength, the steady-state Q ss has a mean µ ss distant from the origin, whose distribution is super-Poissonian over the displaced Fock states [46,47] |µ ss , n , i.e. the corresponding density matrix is which would fall back to the Gibbs state described above when the driving vanishes.When the driving is present, the transient evolutions of Q(α, α * , t) under different driving strengths and frequencies towards Q ss are not unified but classifiable into four dynamic categories or phases (Cf.Appendix D): (I) oscillating, (II) collapse and revive, (III) decay and revive, and (IV) amplifying.These phases demonstrate the competition of the external driving as a restoring energy flow against the qubit decay into the cTLS during the decohering evolution.The typical plots for these four phases of evolutions are shown in Fig. 2, along with the phase partitions over the qubit-driving detuning ∆ d and the driving strength Ω.At near resonance, the driving becomes more effective in pushing the steady-state µ ss away from zero.According to the evolution of the Q-distribution, the energy of the qubit varies with time in a way that if the driving is weak or far-resonant such that |µ ss | < |α 0 |, the energy supplied by the driving to the qubit is not sufficient to cover the energy loss towards the cTLS (phase I to III).At near resonance with sufficiently large driving, however, the qubit can accumulate energy in addition to dissipating into the cTLS, exhibiting the scenario of phase IV.Further analysis of the phase partitions are given in Appendix D.

V. TRANSIENT DYNAMICS
To appreciate the transient evolution of the qubit under the influence of the decohering cTLS, we consider the fully inverted |e e| as the qubit initial state.Its Q-representation | e|α | 2 /π is the continuous complex function e −|α| 2 |α| 2 /π on the (α, α * ) plane, which exhibits a circular-symmetric ring shape in its contour plot over the observable XY -quadrature plane shown in Fig. 3(a).The quadrature operators x and y being linear combinations of a and a † , the distribution of Q(x, y) shown in the plot is identical to Q(α, α * ) up to a proportional constant stemmed from the Jacobian matrix.By shrinking the initial coherent-state distribution of Eq. 9 to a point source, i.e. letting σ 2 G (t) = [1 − exp(−2νγt)] (ν/4ν + 1/2), the solution expressed by Eq. 9 represents directly a propagator G(α, t; α 0 , 0) of the diffusion process from that point source.Then, the time function written as the integral d 2 α 0 |α 0 | 2 exp{−|α 0 | 2 }G(α, t; α 0 , 0)/π describes the entire decoherence process (Cf.Appendix E), the illustration of which at moments after t = 0 are given in Figs.3(b)-(f).The center of the distribution, marked by a red circle, shows a converging trajectory (red solid curve) towards a limit cycle of radius |µ ss |, marked as a white dashed circle.This shows that after the dynamic balance between the driving and the energy decay to the cTLS is reached, the qubit coordinates follow the circular path at the frequency determined by the detuning ∆ d .During the convergence, the ring-shaped distribution of variance 3 √ π/4 contracts into a disk shape of variance ν/4ν + 1/2.In other words, the decoherence process is reflected as a displacement process on the quadrature plane from the origin to |µ ss |, where the excited state distribution decays into the super-Poissonian displaced Fock state of Eq. (11).Not only does the driving strength η under relaxation and the detuning ∆ d determine the radius of the limit cycle orbit, they also decide the rate at which the qubit state gravitates towards it.Shown in Figs.4(a)-(b), the trajectories of µ(t) converge to the same limit cycle under identical driving conditions, but approach it under different rates when commencing from initial points of differing distance to the limit cycle.The illustrated case of α 0 = 0 has a faster approach than that of α 0 = 1.The generalized description of the distinct approaching rates is given conventionally by the Poincare map, shown in Fig. 4(c), where the first recurrence on the positive x-axis is recorded as a function P (α 0 ) of the starting coordinate α 0 on the same axis.The faster rate possessed by the starting coordinate α 0 = 0 is associated with the situation that P (0) is closer  to the fixed point P (α 0 ) = α 0 , i.e. x coordinate of the limit cycle, than P (1).Moreover, the driving conditions determine a partition point, shown as α 0 = 0.89 in the case illustrated in Fig. 4(c), separating a zone of first recurrence occurring at a full circle from another of first occurrence at a half circle.The former has the P (α 0 ) function increase at a slope 0.67 while the latter 0.82.

VI. CONCLUSIONS
We obtain an analytic description of the evolution of a superconducting qubit under the simultaneous influence of a coherent two-level system (cTLS) reflecting the role of a junction defect in a thermal bath and of an external microwave driving field.By solving a master equation in the continuous (α, α * )-plane of coherent state, we has described the thermalization process undergone by the qubit according to an effective temperature determined by the level spacing of the cTLS relative to that of the qubit.The strength of the driving competes with the couplings of the cTLS to the qubit and the environment to determine the effective qubit relaxation, which are classifiable into four phases of distinct dynamic behaviors.The competition has the qubit converge to a limit cycle of finite radius on the quadrature plane, where the rate of convergence is shown through a Poincare map.We also note that a more general scenario would arise when the relaxation is induced by multiple cTLS.The dynamics of the qubit becomes more intricate and it would be a subject of investigation in future research endeavors.

APPENDIX A: MASTER EQUATION OF QUBIT-CTLS EVOLUTION
We derive here the master equation ( 3) from the Liouville equation in the Schrödinger picture with the effective Hamiltonian given in Eq. ( 2).It is first derived in the interaction picture before switching back to the present form under the Schroedinger picture.To be exact, we consider the density matrix ρ I (t) = e iH D 0 t ρ(t)e −iH D 0 t in the interaction picture with the free-energy Hamiltonian H D 0 = 1 2 ω c σ z + ∆ d a † a given under the rotating frame U = exp(−ia † aω d t) and the transformation of displacement D(α d ) = exp(α d a † − α * d a).The evolution of the density matrix follows the Liouville equation ∂ρ I (t)/∂t = −i H eff (t), ρ I (t) in the interaction picture.Because the coupling in the qubit-cTLS system is weak, the RHS of the equation is expanded up to the secondorder terms under the Born approximation.Besides, we assume the initial state between the system and the environment be uncorrelated, i.e. ρ(t) = ρ q (t) ⊗ ρ c (0), so that the firstorder terms becomes zero.Further, under the Markov approximation, the short memory assumption ensures that it is valid to extend the upper limit of the integral in the iterative term to infinity.Therefore, we arrive at the dynamics equation of the qubit-reduced density operator This equation can be further simplified by expanding the commutator and tracing out the degree of freedom in cTLS, viz.
where the rate coefficients are given by According to the Fermi's golden rule, the coefficients γ ± are characterized from the cTLS power spectral density S(ω) at ω = ω ′ q by [48,49] where ν + = 1 − ν and ν − = ν with ν = 1/ [1 + exp (ω c /T )] being the thermal inversion at bath temperature T , and γ c is relaxation rate of the cTLS.∆ c (∆ d ) is the qubit detuning from the cTLS (the driving).Similarly, the coefficients η ± are given by Substituting these rate coefficients into Eq.(A2), we have where The first parenthesis of Eq. (A7) describes the process where an energy quantum is emitted from the qubit into the cTLS, while the second describes the reverse process.The parenthesis on the second line reveals the effect of the driving.
Transforming the master equation back to the Schrödinger picture, one can get the master equation in the rotating displaced frame: If Born-Oppenheimer-Markov approximation is not considered, the evolution of the density matrix ρ I (t) follows the Liouville equation ∂ρ I (t)/∂t = −i H eff (t), ρ I (t) where ρ I (0) = ρ I q (0) ⊗ ρ I c (0) is separable at the initial moment.Substituting Eq. ( 2) and applying σ − e i∆ct → σ − to simplify the expression, the equation reads Compared to Eq. (A7) with the approximation, the first rerm on RHS corresponds to the Lindbladians with coefficient γ when the partial trace over the cTLS is taken.The second term corresponds to the external driving with coefficient η.Due to the cTLS state remaining approximately unchanged during the evolution of the qubit [17], we can approximate ρ I c (t) = ρ I c (0), justifying the application of the approximation.

APPENDIX B: FOKKER-PLANCK EQUATION
Equations ( 7) and ( 8) are the equations of motion for the real and imaginary quadratures of the Q representation of the density matrix ρ q that follows from Eq. ( 3).This Q representation of the qubit is the quasi-probability distribution Q(α, α * , t) = α| ρ q |α /π in the complex (α, α * )plane [40].Taking the inner product of ρ q according to the definition on both sides of Eq. ( 3) and using the formulae a |α = α |α , a † |α = (α * + ∂/∂α) |α , and their Hermitean conjugates, the master equation is rendered into the Fokker-Planck (FP) equation about the time evolution of the probability density function Q.This type of equations is usually solved by separating it into a system of two linear equations for the real and the imaginary quadratures.To do so, we transform the complex variables as α = αe −i∆ d t and α * = α * e i∆t to eliminate the imaginary coefficients in the differential operator on the first line.The FP equation can now be rewritten as about Q(α, α * , t), where the complex coefficients only appear as the oscillating terms on the second line.Recognizing the real quadratures x = ℜ(α) and y = ℑ(α) and consequently rewriting the differential operators such as we obtain Eq. ( 6) about Q(x, y, t), where all imaginary coefficients have been eliminated.In general, the Q function in Eq. ( 6) is not separable into a product of functions of x and y alone, e.g. for the excited state, The initial value problem constituted by Eq. ( 6) and an initial condiction Q(x, y, 0) can not solved by separation of variables.Nevertheless, the equation of the propagator associated with Eq. ( 6) can be solved by separation of variables, regarding the initial condition as a point source delta function.With the solution of the propagator, i.e. a Green's function G(x ′ , y ′ , t; x ′ , y ′ , 0), the solution of Q can then be simply obtained by integrating the product integrand G(x ′ , y ′ , t; x ′ , y ′ , 0) Q(x ′ , y ′ , 0) over x ′ and y ′ .

APPENDIX C: SOLVING THE FOKKER-PLANCK EQUATION
The equation associating with Eq. ( 6) for the propagator during the time duration of 0 to t is [40] where the initial condition is the separated in the complex plane.When assuming G(x, y, t; x ′ , y ′ , 0) = G x (x, t; x ′ , 0)G y (y, t; y ′ , 0), one separates Eq. (C1) into With the help of a transformation of variables [41], with τ ′ = τ | t=0 and ξ ′ = ξ| x=x ′ ,t=0 , the Green's function in the x-quadrature can be written as for which Eq. (C2) becomes the standard diffusion equation The equation admits the typical Gaussian function solution After substituting Eqs.(C4) and (C5) into it, one obtains where Similarly, for Eq.(C3), one obtains Multiplying Eq. (C9) by Eq. (C10) results in the Green's function solution to Eq. (C1): which retains the form of Gaussian distribution with moving mean µ G = µ ss exp {i∆ d t} + (α ′ − µ ss ) exp {−γ νt} in the complex plane, where µ ss = η/(νγ +i∆ d )−α d is the asymptotic mean at steady state.Given Eq. (C11), the distribution Q at any arbitrary time t becomes as explained in the last section.

APPENDIX D: DETERMINING DYNAMICAL PHASES OF EVOLUTION
We analyze the dynamical features reflected in the evolutions of Q(α, α * , t), which are used to obtain the phase diagram Fig. 2(e).This Q function described by Eq. ( 9) is a Gaussian distribution of time-dependent mean µ(t) and variance σ(t).Its general features such as its envelope rate of convergence and asymptotic value are characterized by the moving mean.Starting from the same initial point α 0 = 1, these features become distinct according to two key parameters: (i) the qubit-driving detuning ∆ d and (ii) the driving strength Ω.We categorize the distinction as four dynamic phases shown in Fig. 2.
Phase I occurs when the driving strength is so weak that the energy pumped into the qubit by the driving is not sufficient to cover the energy lost to the cTLS and its environment.Given the finite detuning, the qubit shows an oscillating |µ(t)| from t = 0 to the end of simulation about the reference curve (plotted with no driving).When the driving increases to a sufficient strength, the rate of energy pumped by the driving overcomes the energy loss rate so that asymptotically |µ ss | converges to a finite value at steady state, giving rise to either the dynamic Phase II or the dynamic Phase III.Whether the qubit follows the collapse and revival path (II) or the decay path (III) depends on the detuning ∆ d .At close resonance with small ∆ d , the driving is efficient enough to overcome the energy loss early on so that the qubit trajectory converges directly into the limit cycle without letting the mean µ(t) ever reaching the center of the quadrature plane.In contrast, at large detuning of ∆ d and thus less effective driving, the qubit is halfway between the oscillating Phase I and the decay and reviving Phase III, permitting the |µ(t)| to hit the zero level when the trajectory oscillates below the reference damping curve during the first oscillating period.
Figure 2 shows the cases of trajectories at negative detuning.When the detuning is positive, the behaviors are similar, except that the switched sign inverts the dynamic phase of the trajectories.Consequently, the trajectories undergo an extra π-phase worth of oscillation before exhibiting the typical scenarios of Phases II and III, as illustrated in Fig. 5 below.For instance, shown in Fig. 5(a), the collapse appears later in time and, shown Fig. 5(b), there is a bump before the asymptotic value is reached, compared to the negative detuning cases.
Finally, when the driving strength is strong (typically Ω > ∆ d ), the qubit accumulates energy faster than its dissipation into the cTLS, resulting in a dynamic balance between driving and relaxation that has a limit cycle radius even larger than the initial radius |α 0 |.This near-resonant strong driving hence gives rise to the amplifying phase IV. where The first term inside the square bracket corresponds to a ring-shaped distribution which vanishes at long time limit t → ∞.Whereas, the second term corresponds to a disk-shaped distribution which vanishes at t = 0.The decay rate 2γ ν of the ring-shaped distribution matches the growing rate of the disk distribution.We plot Fig. 3 to indicate these characteristics with the driving parameters Ω/2π = 500 kHz and ∆ d /2π = −300 kHz.
If the initial state is rather a coherent state |α 0 , the initial Q-distribution on the quadrature plane would be a Gaussian disk about α 0 away from the origin, as illustrated in Fig. 6(a) with α 0 = 1, instead of a ring shape about the origin.The decoherence dynamics then mandates the trajectory of the mean µ(t) to converge to the limit cycle, while the distribution always remain Gaussian.Given the operating phase to be type II, the mean reaches the center first at t = 1.142µs (shown in subplot (b)) before converging towards the limit cycle of a radius smaller than α 0 = 1 (shown in subplot (c)).The variance of the Gaussian distribution is slightly increased throughout the evolution.

APPENDIX F: COMPARISON OF ANALYTICAL AND NUMBERICAL SOLUTIONS
To verify the accuracy of our analytical solution, we use matrix exponentiation with Euler's equal step size to numerically solve Eq. (3) in superoperator form and compare it with the analytical solution.Figure 7 shows the evolution of |µ(t)| with full population at |e as the initial state, computed using both the analytical results of Sec.III and the numerical method, which adopts the cases where up to the lowest 4, 8, and 19 excited levels are considered as the computational Hilbert spaces.Subplot (a) show the magnitude |µ| while subplot (b) shows the difference ∆|µ| between the analytical |µ| = |µ ss exp(i∆ d t) + (α 0 − µ ss ) exp(−νγt)| and the numerical |µ|.The numerical µ(t) is identical to a(t) = tr(aρ q (t)) such that given the numerical solution of ρ q (t) at time t as some n × n matrice, the operator a is also expanded under the same Fock basis of n vectors, over which the trace is then taken to compute a numerical value.
The identity between a(t) and µ(t) is recognized if we take the same trace but over the continuous coherent-state space, i.e. a = d 2 α α| ρ q a |α /π = d 2 α αQ(α, α * , t) where the Q quasi-probability distribution is given by Eq. (E1) for |e being the initial state.We recall that Q(α, α * , t) has two terms: one ring-shaped distribution and one disk-shaped distribution.The former contributes, up to some constant coefficent, a product of α|α − µ| 2 and a Gaussian function in the integrand; the latter contributes a product of α and the same Gaussian.Writing α = x+iy and µ = A+iB, α|α−µ| 2 contains four terms of the form x(x−A) 2 = (x−A) 3 +A(x−A) 2 , i.e. an odd power plus an even power about the axis x = A. Since the Gaussian is itself even about x = A, when multiplied by this cubic function, the (x − A) 3 term would vanish after integration due to the symmetry in the integral while the A(x − A) 2 term would contribute Aσ 2 after integration by parts.The integration over y is about the Gaussian only and contributes only a unity coefficient.Similarly, the term iy(y−B) 2 becomes iBσ 2 after integration.The rest two cross terms such as iy(x − A) 2 = i(y − B)(x − A) 2 + iB(x − A) 2 can be similarly integraded, contributing iBσ Summing the two results give rise to A + iB = µ, hence proving µ = a .
Using the same physical parameters as those of Fig 3, i.e. ∆ d /2π = −300 kHz and Ω/2π = 500kHz, the plots of Fig. 7 demonstrate the accuracy of the analytical approach, which is apparently more accurate than the numerical approach when limited number of levels in ρ q is considered.Even when the number of discrete levels is expanded to 20, the discrepany from the numerical solution is on the order of 10 −4 , negligible compared to the overall magnitude of |µ| over the entire evolution duration considered. and

Figure 2 .
Figure 2. Dynamic phases of the decohering evolution of a superconducting qubit under driving, starting from the same initial point α0 = 1 in the complex (α, α * )-plane.The common parameters include the qubit frequency ω ′ q /2π = 4.5 GHz, the qubit-cTLS detuning ∆c/2π = 110 kHz, the interaction strength g/2π = 400 kHz, the cTLS relaxation rate γc/2π = 1 MHz, and the environmental temperature T = 10 mK.The dynamic evolution exhibits distinct behaviors that can be categorized into four phases under distinct driving detunings ∆ d /2π and strengths Ω/2π: (a) −1.5 MHz and 100 kHz for phase I (green); (b) −370 kHz and 100 kHz for phase II (deep blue); (c) 20 kHz and 100 kHz for phase III (light purple); and (d) 100 kHz and 250 kHz for phase IV (pink).The dashed curves in (a)-(d) shows the evolution of µ(t) without driving as a reference.(e) shows the partitions of these dynamic phases over a range of detunings and strengths.

Figure 3 .
Figure 3.The evolution of the qubit Q-distribution plotted as contours on XY-quadrature plane, assuming a fully inverted excited state|e as the initial state.The captures of six time moments from t = 0 to t = 10µs are shown in the subplots (a)-(f), where the red curves show the trajectory of the mean µ(t) and the overlay white dashed curve indicates the limit cycle.The driving field parameters are set to Ω/2π = 500 kHz and ∆ d = −300 kHz, where the other parameters remain the same as those in Fig. 2.

Figure 4 .
Figure 4. Trajectories of µ(t) drawn as red dotted curves on the XYquadrature plane, starting from two initial points: (a) α0 = 0 and (b) α0 = 1, but under the same coupling strength Ω/2π = 1MHz and detuning ∆ = −2.5 MHz.The solid circles indicate the same asymptotic limit cycle of radius |µss| from the distinct initial points.The magenta squares indicate the first recurrence on the positive xaxis.(c) Poincare map showing the coordinate of the first recurrence as a function of initial x-coordinate α0.The black dashed cross shows the fixed point α0 = 0.40 of the recurrence function P (α0) = α0 = |µss|.The function has two rate zones, one before and one after α0 = 0.89.

Figure 6 .
Figure 6.The evolution of the qubit Q-distribution plotted as contours on XY -quadrature plane, with the coherent state |1 as the initial state, shown as subplots at three typical time moments (a) t = 0, (b) t = 1.142µs, and (c) t = 5µs.The red curves show the trajectory of the mean µ(t) and the small red circles mark the locations of µ(t) at the these three moments.The overlay white dashed curve indicates the limit cycle.The driving field parameters are equal to those in Fig. 2(b).

2 andαFigure 7 .
Figure 7. Evolutions of the distribution center µ(t) of the qubit density matrix ρq, expressed by its Q-representation, solved using the analytical as well as numerical approaches.The first excited state |e is assumed as the initial state and physical parameters used are identical to those of Fig 3.(a) shows the magnitude |µ(t)| where the black solid curve corresponds to the analytical solution.The blue dashed curve (eclipsing the black curve) corresponds to the numerical solution with 8 excited levels are simulated while the orange dashed curve corresponds to that with 4 excited levels simulated.(b) shows the difference ∆|µ(t)| of the numerical solutions from the analytical one.The blue curve represents the scenario with 8 excited levels while the yellow represents that with 19 excited levels.