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

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.¹ 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,² with experiments demonstrating T₁ and T₂ times on the order of 10 (Ref. 3) and 100 μs,⁴ 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,⁵ 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.⁵ 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.¹⁰ Collectively, cTLS and fluctuators contribute a 1/f background noise spectrum,¹³ and their interaction would decide the power dependence of the quality factor of superconducting resonators.¹⁴ 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.¹² We note that in optical amorphous materials, the defect couplings are described by non-Hermitian Hamiltonians.²⁵ 

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 $( α , α ∗ )$-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 ( α , α ∗ , 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.

To verify the validity of the analytical approach, we employ the matrix exponentiation method to numerically solve the master equation in a superoperator form⁴² 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).

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,⁴³ 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 cTLS-mediated 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⁴⁴ and is useful for determining ω_c.^39,45 Moreover, the widened variance of $σ ss 2$ by $ν / 4 ν ¯$ compared to $σ 2 ( 0 )$ also stems from the thermal energy transferred from the cTLS to the qubit.

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 $σ G 2 ( 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)–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 $| μ 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 toward it. Shown in Figs. 4(a)–4(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 α₀ 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.

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.

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

The authors have no conflicts to disclose.

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).

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

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 $ρ I ( t ) = e i H 0 D t ρ ( t ) e − i H 0 D t$ in the interaction picture with the free-energy Hamiltonian $H 0 D = ( 1 / 2 ) ω c σ z + Δ d a † a$ given under the rotating frame $U = exp ( − i a † a ω d t )$ and the transformation of displacement $D ( α d ) = exp ( α d a † − α d * a )$⁠.

In particular, if the initial state is a coherent state $| α 0 ⟩$⁠, whose Q-distribution is $exp { − | α − α 0 | 2 } / π$⁠, the Gaussian Eq. (9) would arise from Eq. (C12). Compared to the variance $σ G$ of the propagator, the variance from this coherent state follows as $σ 2 ( t ) = exp { − 2 γ ν ¯ t } / 2 + σ G 2$ while the mean $μ ( t ) = μ ss exp { i Δ d t } + ( α 0 − μ ss ) exp { − ν ¯ γ t }$ is the same of that of the propagator except that the arbitrary $α ′$ becomes the coherent state displacement α₀.

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.

If the initial state is rather a coherent state $| α 0 ⟩$⁠, the initial Q-distribution on the quadrature plane would be a Gaussian disk about α₀ 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 toward 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.

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) shows the magnitude $| μ |$ while subplot (b) shows the difference $Δ | μ |$ between the analytical $| μ | = | μ s s exp ( i Δ d t ) + ( α 0 − μ s s ) 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 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.

Using the same physical parameters as those of Fig. 3, i.e., $Δ d / 2 π = − 300$ kHz and $Ω / 2 π = 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 − 4$⁠, negligible compared to the overall magnitude of $| μ |$ over the entire evolution duration considered.