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 displaced Fock states, which reduces to a Gibbs state of effective temperature decided by the defect at zero driving limit. Furthermore, 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 (Ref. 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 toward 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 modeled as two-level systems.^{5} In particular, these two-level defects are in general categorized into two types: fluctuators, with level spacing $ \Delta E \u2009 \u2272 \u2009 k B T$, and coherent two-level systems (cTLS), with level spacing $ \Delta E \u226b 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–21} where the theoretical analyses are modeled on the interactions between Pauli operators (for cTLS) and anharmonic oscillators (for qubits).^{22–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 toward 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–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 $ ( \alpha , \alpha \u2217 )$-plane of coherent states by regarding the driving as a displacement that translates the creation and annihilation 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 ( \alpha , \alpha \u2217 , 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 on both the qubit effective decay mediated by the cTLS and the driving strength. The transient behaviors toward 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. II, 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

*ω*and anharmonicity

_{q}*χ*. 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.

^{35–37}The second term $ H c = \omega c \sigma 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 \u2020 \sigma \u2212 + a \sigma + )$ and for the semiclassical field driving $ W ( t ) = \Omega ( a \u2020 e \u2212 i \omega d t + h . c . )$, where Ω is the driving strength and

*ω*the field frequency.

_{d}*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

*ω*, where Ω and the qubit-drive detuning $ \Delta d = \omega q \u2032 \u2212 \omega d$ determine the amount of displacement $ \alpha d = \Omega / \Delta d$ from the vacuum state, i.e., $ H q + W ( t ) = D ( \alpha d ) H q D \u2020 ( \alpha d )$. Given weak driving Ω, the magnitude of the displacement is small, which renders the qubit to fall into the weak-excitation regime,

_{d}^{38}where the oscillator frequency can be linearized as $ \omega q \u2032 = \omega q \u2212 \chi ( n \xaf q \u2212 1 )$. The number $ n \xaf q$ denotes the average photon number in the qubit, which remains a constant value since the coupling strengths $ g , \Omega \u226a \omega $. By regarding $ D H q D \u2020$ and

*H*as free energies in the coherent-state space, the effective qubit-cTLS coupling under the interaction picture becomes

_{c}*ρ*is derived from the Liouville equation $ \rho \u0307 = \u2212 i [ H eff , \rho ]$ for the coupled system

_{q}*ρ*. 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., $ \rho ( t ) = \rho q ( t ) \u2297 \rho 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

*η*term. With $ \nu = ( exp \u2009 ( \omega c / k B T ) + 1 ) \u2212 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 $ \gamma = 2 g 2 \gamma c / ( \gamma c 2 + \Delta c 2 )$ is the decay rate according to the typical relaxation rate

*γ*of a cTLS.

_{c}^{39}The

*η*term reflects the external driving under the influence of the environment, i.e.,

*α*determined by the driving amplitude. Given the environmental temperature (on the mK range) of superconducting qubits, the inversion

_{d}*v*adopts a negligible value and the existence of the driving field in general accelerates qubit relaxations (contribution of the first term $\u226b$ that of the second term).

## III. SOLVING FOR QUADRATURES

*α*-plane. Under the

*Q*-representation $ Q ( \alpha , \alpha \u2217 , t ) = \u27e8 \alpha | \rho q ( t ) | \alpha \u27e9 / \pi $, the FP equation reads

*ν*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 = ( \alpha e i \Delta d t + h . c . ) / 2$ and $ y = \u2212 i ( \alpha e i \Delta d t \u2212 h . c . ) / 2$ (cf. Appendix B), giving

*x*and

*y*variables, we consider the associated equation of motion for the propagator (viz. the Green's function) $ G ( x , y , t ; x \u2032 , y \u2032 , 0 )$ linking the initial coordinate $ ( x \u2032 , y \u2032 )$ and the current coordinate (

*x*,

*y*) at current time

*t.*

^{40}Since it reduces to the separated $ \delta ( x \u2212 x \u2032 ) \delta ( y \u2212 y \u2032 )$ at the time limit $ t \u2192 0$, the propagator is separable as $ G x ( x , t ; x \u2032 , 0 ) G y ( y , t ; y \u2032 , 0 )$, in general, breaking up its equation of motion into the pair,

^{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,

*α*

_{0}as the displacement of the coherent state. The mean converges to the asymptotic coordinates,

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 $ | \mu ( t ) |$ and the numerical $ | \mu ( t ) |$ is marginal, a difference on the order of $ 10 \u2212 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 $ | \alpha 0 \u27e9$ for which $ Q ( 0 ) \u223c N ( \alpha 0 , 1 / 2 )$, follows a typical trajectory of decay into the quasi-coherent state $ \rho ss$ for which $ Q ( t \u2192 \u221e ) \u223c N ( \mu ss , \nu / 4 \nu \xaf + 1 / 2 )$ independent of the initial state $ | \alpha 0 \u27e9$. The mean $ \mu 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 \u2009 { \u2212 | \alpha | 2 / ( n \xaf + 1 )} / \pi ( n \xaf + 1 )$ in the continuous *α*-plane with $ n \xaf = \nu / 2 \nu \xaf = ( exp \u2009 { \omega c / k B T} \u2212 1 ) \u2212 1$ indicating the equivalent average photon number. In other words, the qubit levels would be populated not according to $ \omega q \u2032$, 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 cTLS-mediated bath at the effective temperature $ T q = T \omega q \u2032 / \omega c$. Given the ultra-low bath temperature

*T*of the dilution fridge, the difference in qubit population inversions under

*T*and

*T*is significant, despite the near-resonant frequencies of $ \omega q \u2032$ and

_{q}*ω*. This temperature difference is recently measured

_{c}^{44}and is useful for determining

*ω*.

_{c}^{39,45}Moreover, the widened variance of $ \sigma ss 2$ by $ \nu / 4 \nu \xaf$ compared to $ \sigma 2 ( 0 )$ also stems from the thermal energy transferred from the cTLS to the qubit.

^{46,47}$ | \mu ss , n \u27e9$, i.e., the corresponding density matrix is

_{d}and the driving strength Ω. At near resonance, the driving becomes more effective in pushing the steady-state $ \mu ss$ away from zero. According to the evolution of the

*Q*-distribution, the energy $ \u27e8 H q ( t ) \u27e9 = \u210f \omega q \u2032 ( | \mu ( t ) | 2 + 2 \sigma 2 ( t ) \u2212 1 )$ of the qubit varies with time in a way that if the driving is weak or far-resonant such that $ | \mu ss | < | \alpha 0 |$, the energy supplied by the driving to the qubit is not sufficient to cover the energy loss toward 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 is 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 \u27e9 \u27e8 e |$ as the qubit initial state. Its *Q*-representation $ | \u27e8 e | \alpha \u27e9 | 2 / \pi $ is the continuous complex function $ e \u2212 | \alpha | 2 | \alpha | 2 / \pi $ on the $ ( \alpha , \alpha * )$ 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 \u2020$, the distribution of *Q*(*x*, *y*) shown in the plot is identical to $ Q ( \alpha , \alpha \u2217 )$ 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 $ \sigma G 2 ( t ) = [ 1 \u2212 exp \u2009 ( \u2212 2 \nu \xaf \gamma t ) ] ( \nu / 4 \nu \xaf + 1 / 2 )$, the solution expressed by Eq. (9) represents directly a propagator $ G ( \alpha , t ; \alpha 0 , 0 )$ of the diffusion process from that point source. Then, the time function written as the integral $ \u222b d 2 \alpha 0 \u2009 | \alpha 0 | 2 \u2009 exp \u2009 { \u2212 | \alpha 0 | 2} G ( \alpha , t ; \alpha 0 , 0 ) / \pi $ describes the entire decoherence process (cf. Appendix E), the illustration of which at moments after *t* = 0 are given in Figs. 3(b)–3(f). The center of the distribution, marked by a red circle, shows a converging trajectory (red solid curve) toward a limit cycle of radius $ | \mu 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 \pi / 4$ contracts into a disk shape of variance $ \nu / 4 \nu \xaf + 1 / 2$. In other words, the decoherence process is reflected as a displacement process on the quadrature plane from the origin to $ | \mu 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 toward it. Shown in Figs. 4(a)–4(b), the trajectories of $ \mu ( 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 $ \alpha 0 = 0$ has a faster approach than that of $ \alpha 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 ( \alpha 0 )$ of the starting coordinate *α*_{0} on the same axis. The faster rate possessed by the starting coordinate $ \alpha 0 = 0$ is associated with the situation that *P*(0) is closer to the fixed point $ P ( \alpha 0 ) = \alpha 0$, i.e., *x* coordinate of the limit cycle, than *P*(1). Moreover, the driving conditions determine a partition point shown as $ \alpha 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 ( \alpha 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 $ ( \alpha , \alpha \u2217 )$-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.

## ACKNOWLEDGMENTS

H.I. thanks the support by FDCT of Macau under Grant Nos. 0130/2019/A3, 0015/2021/AGJ, and 006/2022/ALC.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Yanxiang Wang:** Conceptualization (supporting); Formal analysis (lead); Investigation (equal); Methodology (equal); Validation (supporting); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Ziyang You:** Conceptualization (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Validation (lead); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Hou Ian:** Conceptualization (lead); Formal analysis (supporting); Funding acquisition (lead); Investigation (equal); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### 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 Schrödinger picture. To be exact, we consider the density matrix $ \rho I ( t ) = e i H 0 D t \rho ( t ) e \u2212 i H 0 D t$ in the interaction picture with the free-energy Hamiltonian $ H 0 D = ( 1 / 2 ) \omega c \sigma z + \Delta d a \u2020 a$ given under the rotating frame $ U = exp \u2009 ( \u2212 i a \u2020 a \omega d t )$ and the transformation of displacement $ D ( \alpha d ) = exp \u2009 ( \alpha d a \u2020 \u2212 \alpha d * a )$.

*.,*

^{48,49}

*T*, and

*γ*is relaxation rate of the cTLS. $ \Delta c ( \Delta d )$ is the qubit detuning from the cTLS (the driving). Similarly, the coefficients $ \eta \xb1$ are given by

_{c}*γ*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 $ \rho c I ( t ) = \rho c I ( 0 )$, justifying the application of the approximation.

### APPENDIX B: FOKKER–PLANCK EQUATION

*Q*representation of the density matrix

*ρ*that follows from Eq. (3). This

_{q}*Q*representation of the qubit is the quasi-probability distribution $ Q ( \alpha , \alpha \u2217 , t ) = \u27e8 \alpha | \rho q | \alpha \u27e9 / \pi $ in the complex $ ( \alpha , \alpha \u2217 )$-plane.

^{40}Taking the inner product of

*ρ*according to the definition on both sides of Eq. (3) and using the formulae $ a | \alpha \u27e9 = \alpha | \alpha \u27e9$, $ a \u2020 | \alpha \u27e9 = ( \alpha \u2217 + \u2202 / \u2202 \alpha ) | \alpha \u27e9$, $ | g \u27e9$ and their Hermitian conjugates, the master equation is rendered into the Fokker–Planck (FP) equation,

_{q}*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 $ \alpha = \alpha \u0303 e \u2212 i \Delta d t$ and $ \alpha \u2217 = \alpha \u0303 \u2217 e i \Delta t$ to eliminate the imaginary coefficients in the differential operator on the first line. The FP equation can now be rewritten as

*x*and

*y*alone, e.g., for the excited state, $ Q \u0303 = ( | x | 2 + | y | 2 ) \u2009 exp \u2009 ( \u2212 | x | 2 \u2212 | y | 2 ) / \pi $ is inseparable. The initial value problem constituted by Eq. (6) and an initial condition $ Q \u0303 ( x , y , 0 )$ cannot be 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 \u2032 , y \u2032 , t ; x \u2032 , y \u2032 , 0 )$, the solution of $ Q \u0303$ can then be simply obtained by integrating the product integrand $ G ( x \u2032 , y \u2032 , t ; x \u2032 , y \u2032 , 0 ) Q \u0303 ( x \u2032 , y \u2032 , 0 )$ over $ x \u2032$ and $ y \u2032$.

### APPENDIX C: SOLVING THE FOKKER–PLANCK EQUATION

*t*is

^{40}

^{41}

*x*-quadrature can be written as

*Q*at any arbitrary time

*t*becomes

In particular, if the initial state is a coherent state $ | \alpha 0 \u27e9$, whose *Q*-distribution is $ exp \u2009 { \u2212 | \alpha \u2212 \alpha 0 | 2} / \pi $, the Gaussian Eq. (9) would arise from Eq. (C12). Compared to the variance $ \sigma G$ of the propagator, the variance from this coherent state follows as $ \sigma 2 ( t ) = exp \u2009 { \u2212 2 \gamma \nu \xaf t} / 2 + \sigma G 2$ while the mean $ \mu ( t ) = \mu ss \u2009 exp \u2009 { i \Delta d t} + ( \alpha 0 \u2212 \mu ss ) \u2009 exp \u2009 { \u2212 \nu \xaf \gamma t}$ is the same of that of the propagator except that the arbitrary $ \alpha \u2032$ becomes the coherent state displacement *α*_{0}.

### APPENDIX D: DETERMINING DYNAMICAL PHASES OF EVOLUTION

We analyze the dynamical features reflected in the evolutions of $ Q ( \alpha , \alpha \u2217 , 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 $ \mu ( t )$ and variance $ \sigma ( 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 $ \alpha 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 $ | \mu ( 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 $ | \mu 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 $ \mu ( 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 $ | \mu ( 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 $ \Omega > \Delta 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 $ | \alpha 0 |$. This near-resonant strong driving hence gives rise to the amplifying phase IV.

### APPENDIX E: TRANSIENT DYNAMICS ANALYSIS

*Q*-distribution can be solved analytically by Eq. (C12) as

*t*= 0. The decay rate $ 2 \gamma \nu \xaf$ 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 $ \Omega / 2 \pi = 500 \u2009 kHz$ and $ \Delta d / 2 \pi = \u2212 300 \u2009 kHz$.

If the initial state is rather a coherent state $ | \alpha 0 \u27e9$, 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 $ \alpha 0 = 1$, instead of a ring shape about the origin. The decoherence dynamics then mandates the trajectory of the mean $ \mu ( 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 \u2009 \mu $s [shown in subplot (b)] before converging toward the limit cycle of a radius smaller than $ \alpha 0 = 1$ [shown in subplot (c)]. The variance of the Gaussian distribution is slightly increased throughout the evolution.

### APPENDIX F: COMPARISON OF ANALYTICAL AND NUMERICAL 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 $ | \mu ( t ) |$ with full population at $ | e \u27e9$ 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) shows the magnitude $ | \mu |$ while subplot (b) shows the difference $ \Delta | \mu |$ between the analytical $ | \mu | = | \mu s s \u2009 exp \u2009 ( i \Delta d t ) + ( \alpha 0 \u2212 \mu s s ) \u2009 exp \u2009 ( \u2212 \nu \xaf \gamma t ) |$ and the numerical $ | \mu |$. The numerical $ \mu ( t )$ is identical to $ \u27e8 a ( t ) \u27e9 = tr ( a \rho q ( t ) )$ such that given the numerical solution of $ \rho q ( t )$ at time *t* as some *n* × *n* matrices, 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.

*Q*quasi-probability distribution is given by Eq. (E1) for $ | e \u27e9$ being the initial state. We recall that $ Q ( \alpha , \alpha \u2217 , t )$ has two terms: one ring-shaped distribution and one disk-shaped distribution. The former contributes, up to some constant coefficient, a product of $ \alpha | \alpha \u2212 \mu | 2$ and a Gaussian function in the integrand; the latter contributes a product of

*α*and the same Gaussian. Writing $ \alpha = x + i y$ and $ \mu = A + i B$, $ \alpha | \alpha \u2212 \mu | 2$ contains four terms of the form $ x ( x \u2212 A ) 2 = ( x \u2212 A ) 3 + A ( x \u2212 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 \u2212 A ) 3$ term would vanish after integration due to the symmetry in the integral while the $ A ( x \u2212 A ) 2$ term would contribute $ A \sigma 2$ after integration by parts. The integration over

*y*is about the Gaussian only and contributes only a unity coefficient. Similarly, the term $ i y ( y \u2212 B ) 2$ becomes $ i B \sigma 2$ after integration. The rest two cross terms, such as $ i y ( x \u2212 A ) 2 = i ( y \u2212 B ) ( x \u2212 A ) 2 + i B ( x \u2212 A ) 2$ can be similarly integrated, contributing $ i B \sigma 2$ and $ A \sigma 2$, respectively. To summarize,

*α*and the Gaussian, write $ x = ( x \u2212 A ) + A$ and $ i y = i ( y \u2212 B ) + i B$. Thus, following similar observations about evenness and oddities of the integrands, we have

Using the same physical parameters as those of Fig. 3, i.e., $ \Delta d / 2 \pi = \u2212 300$ kHz and $ \Omega / 2 \pi = 500$ kHz, the plots in 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 discrepancy from the numerical solution is on the order of $ 10 \u2212 4$, negligible compared to the overall magnitude of $ | \mu |$ over the entire evolution duration considered.

## REFERENCES

*Statistical Methods in Quantum Optics 1: Master Equations and Fokker-Planck Equations*

*Handbook of Linear Partial Differential Equations for Engineers and Scientists*