Relaxation of Rabi Dynamics in a Superconducting Multiple-Qubit Circuit

We investigate a superconducting circuit consisting of multiple capacitively-coupled charge qubits. The collective Rabi oscillation of qubits is numerically studied in detail by imitating environmental fluctuations according to the experimental measurement. For the quantum circuit composed of identical qubits, the energy relaxation of the system strongly depends on the interqubit coupling strength. As the qubit-qubit interaction is increased, the system's relaxation rate is enhanced firstly and then significantly reduced. In contrast, the inevitable inhomogeneity caused by the nonideal fabrication always accelerates the collective energy relaxation of the system and weakens the interqubit correlation. However, such an inhomogeneous quantum circuit is an interesting test bed for studying the effect of the system inhomogeneity in quantum many-body simulation.


Introduction
Owing to the fascinating properties such as flexibility, tunability, scalability, and strong interaction with electromagnetic fields, superconducting Josephson-junction circuits provide an outstanding platform for quantum information processing (QIP) [1,2], quantum simulation of many-body physics [3,4,5], and exploring the fundamentals of quantum electrodynamics in and beyond the ultrastrong-coupling regime [6,7]. Additionally, hybridizing these solid-state devices with the atoms may enable the information transfer between macroscopic and microscopic quantum systems [8,9,10,11,12,13], where the superconducting circuits play the role of rapid processor while the atoms act as the long-term memory. Nonetheless, the strong coupling to the environmental noise significantly limits the energy-relaxation (T 1 ) and dephasing (T 2 ) times of superconducting circuits [14,15,16].
Recently, there has been focus on large-scale QIP [17,18]. Several network schemes have already been demonstrated in experiments: one-dimensional spin chain with the nearest-neighbor interaction [19,20], two-dimensional lattice with quantum-bus-linked qubits [21,22], multiple artificial atoms interacting with the same resonator [17], and many cavities coupled to a superconducting qubit [23]. However, only little attention has been paid to the quantum circuit composed of nearly identical superconducting qubits, where an arbitrary individual directly interacts with others and, in addition, all qubits are biased by the same voltage, current, or magnetic bias and are exposed to the same fluctuation source. Studying such a multiple-qubit architecture is of importance to superconducting QIP network. For one thing, the nearest-or nextnearest-neighbour-coupling approximation, which is commonly employed in scalable superconducting schemes [24,25,26], does not always hold the truth in realistic systems. For another, transferring the quantum information from one processor to another over a long distance, which relies on the long-range interqubit coupling, is essential for cluster quantum computing [27] and quantum algorithms [28]. Moreover, in a large-scale network composed of strong or ultrastrong interacting qubits, the energy relaxation and dephasing of one qubit strongly influence the dynamics of others and this fluctuation can be rapidly boosted and affect more qubits. Such a collective dissipation may significantly degrade the information transfer fidelity between two remotely separated qubits. However, to our best knowledge, this collective dissipation of a large ensemble of directly-coupled superconducting qubits has not been explored before.
Here we investigate the collective relaxation effect in a multiple-qubit circuit, where all charge qubits are biased by the same voltage source, capacitively coupled via linking all islands together, and influenced by the same noise source. The system's Rabi oscillation is numerically studied in detail. It is shown that for the homogeneous system with identical qubits, the interqubit coupling strongly modifies the collective relaxation of quantum circuit. In contrast, the nonideal-fabrication-induced extra inhomogeneity always enhances the system's relaxation rate and also provides a platform for studying many-body localization.

Physical model
We consider N capacitively-coupled single-Josephson-junction charge qubits as shown in figure 1. A voltage source V g biases the array of superconducting islands via the identical gate capacitors C g . All Cooper-pair boxes are linked together via the identical coupling capacitors C c , where we have applied the fact that the practical inhomogeneities of capacitors C g and capacitors C c are much smaller than that of Josephson junctions.
For the k-th (k = 1, . . . , N ) charge qubit, we define the self-capacitance of Josephson junction, the total capacitance, the charging energy, and the Josephson energy as C j,k , 2C Σ,k , and E J,k , respectively. The total excess charge Q k of the k-th box distributes on capacitor plates of C g , C j,k , and C c that are involved in the Cooper-pair box. The corresponding charges are defined as Q g,k , Q j,k , and Q c,k , respectively, and we have According to Kirchhoff's circuit laws and charge conservation, one obtains The total charging energy of the system is then given by where N k = − Q k (2e) denotes the number of excess Cooper pairs in the k-th box and N g = CgVg (2e) is the gate-charge bias. We have also defined which is positive and denotes the difference between the electrostatic energy of a Cooper pair in C c and the mean value of charging energies of different Cooper-pair boxes. Adding the tunneling energies of Cooper pairs into E ch [29], the system's Hamiltonian is derived as H = H 0 + H 1 , where in the two-state approximation. The x-and z-components of the Pauli operator are given by σ |0 k and |1 k represent the absence and presence of a single excess Cooper pair in the k-th island.
The multiple-qubit system operates in the charging limit of E C,k E J,k . H 0 gives the total energy of free charge qubits while H 1 corresponds to the interqubit interaction energy. It is seen that besides a constant, the qubit-qubit coupling leads to the linear terms of (1 − 2N g )σ (k) z , which may be involved into H 0 , and the quadratic zz-interaction terms of σ (k 1 ) z σ (k 2 ) z that are diagonal in the charge number basis. In the limit of C c ∼ 0, we have V → ∞, for which H 1 ∼ 0 and charge qubits become independent with each other. As C c is increased, E C,k=1,...,N go down. In this work, we restrict C c within the range satisfying the charging limit condition. Due to the assumption of identical gate and coupling capacitors, the system's inhomogeneity completely comes from the nonidentical Josephson junctions.

Homogeneous Circuit
We first consider the homogeneous system, where all Josephson junctions are identical, i.e., C j,k = C j , E C,k = E C , and E J,k = E J . It is easy to obtain V E C = Cg+C j Cc and H can be simplified as by defining the collective operators x and the function Equation (9) illustrates that multiple qubits behave in the same way if they are initialized in the same state. We define the dimensionless parameter η = E 2 C E J V to measure the interqubit coupling strength.

Ground state of nondissipative system
We are interested in the ground state of the nondissipative system in the limit of N → ∞. The Holstein-Primakoff transformation [30], is employed to map the multiple-qubit system onto a bosonic mode. The annihilation b and creation b † operators fulfill the bosonic commutation relation

to this bosonic mode and obtain a displaced Hamiltonian
H is further expanded as ordered power series in b and b † . In the second-order approximation, choosing β such that the linear terms associated with b and b † vanish leads to the ground state |β and the corresponding energy Finally, we find that β is real and determined by from which multiple solutions may be derived and the accepted one should minimize E. We focus on the observable operator J z . Figure 2(a) shows the ground-state expectation value β|Jz|β N = (2β 2 −1) vs. η and N g . One can see that the large interqubit coupling (η 1) opens a wide intermediate regime, where β|Jz|β N strongly relies on the gate-charge bias N g . In contrast, for η 1 the system returns to an ensemble of independent qubits, and β|Jz|β N ≈ 1 (−1) when N g > 1 2 (< 1 2 ). We should note that, unlike the Dicke model of a single cavity mode interacting with an ensemble of two-level atoms [31,32], both first and second derivatives of E with respect to β are continuous, indicating the absence of phase transitions in our system. In the limit of E J ∼ 0, one obtains the simple solution, β|Jz|β N = −Ξ(V, N g ) and E = 0. Due to the condition of −1 ≤ β|Jz|β N ≤ 1, the ground-state diagram may be divided into three regimes: Rabi and Ramsey oscillations with the system parameters of C g , C j , and E J same to figure (1) and E C = 78E J for C c = 0. For the Rabi oscillation, the qubit is initialized in |1 with N g0 = 1 2 . For performing Ramsey fringes, the qubit is initially prepared in |1 and N g0 is set at 1 2 during two π 2 -pulses while N g0 = 0 in the free evolution period. Solid curves: numerical results. Dash lines: decay-envelope fittings. In the Ramsey oscillation, the detailed behavior around the characteristic time sale T 2 , which is surrounded by the rectangular frame, is zoomed in, and the oscillation frequency is given by For all curves of ensemble average j z (t) E , the size of trajectory ensemble is 10 5 . determines the relaxation time T 1 of the charge qubit, i.e., the characteristic time scale of the damped qubit Rabi oscillation, while the latter primarily affects the dephasing time T 2 , i.e., the characteristic time scale of the damped qubit Ramsey/spin echo oscillation. According to [34,35], the whole noise may be mapped onto the gate-charge bias N g , i. e., the voltage source, in the mathematical treatment. This is still valid in this multiqubit scheme. Since all islands are electrostatically biased by the same voltage source, different qubits are subjected to the same Ohmic noise whose spectrum is proportional to the noise frequency f . In addition, all islands are strongly coupled in the system. When environmental fluctuations interrupt the dynamic of one Cooper-pair box, the local voltage of the corresponding gate capacitor is disturbed. This voltage fluctuation may be mapped onto the voltage source and it further influences the dynamics of other boxes in the same way, leading to the energy relaxation of other charge qubits. Thus, one may map all local voltage fluctuations occurring in different Cooper-pair boxes on to the single voltage source. Moreover, in this work we mainly focus on the qubit Rabi oscillation, whose damping time is T 1 , and the low-frequency 1/f noise in the circuit hardly influences T 1 . Therefore, for these reasons it is valid to model environmental fluctuations by the fluctuation of gate-charge bias (voltage source), i.e., all charge qubits suffer the same noise. This differs from the many-body system composed of excited-state atoms that spontaneously decay independently.
We rewrite the gate-charge bias as i.e., a constant value N g0 plus a fluctuating term δN g (t). The noise spectral density is given by where the typical impedance of the voltage-source circuit is R = 50 Ω and α = 5.0 × 10 −7 [36,37,38]. Accordingly to S δNg (f ), one may numerically generate δN g (t) [see figure 3(a)]. We choose {|n 1 ⊗ · · · ⊗ |n N ; n k=1,...,N = 0, 1} to span the Hilbert space. Using the Schrödinger equation, one can simulate the system's state ψ(t) for a given initial state ψ(0), resulting in the trajectory Repeating the simulation with the same initial condition leads to the ensemble mean observation j z (t) E . Here we use . . . E to denote the ensemble average. As an example, figure 3(b) depicts the Rabi and Ramsey oscillations of a single charge qubit, from which the decoherence times T 1,2 may be extracted.
In the following we only focus on the collective Rabi oscillation of multiple-qubit circuit, where the system is initialized in |1 1 ⊗. . .⊗|1 N and we set N g0 = 1 2 . Figure 4(a) shows the damped j z (t) E for several different η. It is seen that for η 1, i.e., the very weak qubit-qubit interaction, j z (t) E is similar to that of single qubit, meaning that multiple qubits act almost independently. When η is increased, i.e., the interqubit coupling becomes strong, the relaxation time T 1 of j z (t) E is dramatically reduced because the qubits get correlated and their dynamics are affected with each other. When one qubit relaxes, the gate-charge bias is disturbed and this influences the dynamics of other qubits, leading to an enhanced collective relaxation. As η is further increased, the interqubit interactions become very strong.
Surprisingly, j z (t) E relaxes much more slowly than that of single qubit and the oscillation behavior of j z (t) E also disappears. It is understandable from the aspect of interaction-induced detunings [39]. Diagonalizing the Hamiltonian (9) leads to the eigenstates of the system. For the noninteracting system (C c = 0), the eigenstates with the same total number of the qubits in |1 are degenerate. The interqubit coupling removes this degeneracy and gives rise to the energy-level shifts of eigenstates. The much strong interactions among qubits enhance the qubit-state shifts to the values comparable or even larger than the Josephson energy E J . Thus, the system operating point moves far away from the sweet spot N g = 1 2 , resulting in a large qubit detuning and the suppression of oscillation behavior of j z (t) E . In addition, according to the relaxation-time formula derived from the Fermi's golden rule [40], T 1 is extended when the system moves away from the optimal working point. The dependence of T 1 on η is displayed in figure 4(b). It is shown that T 1 starts to rise after η ≈ 1 where, as illustrated by equation (9), the energy-level shift (detuning) E 2 C V is equal to the driving strength E J . We should point out that although T 1 is enhanced in the large interqubit-coupling regime, the qubit dephasing time T 2 is strongly suppressed. At the sweet spot N g = 1 2 , the qubit is minimally sensitive to the 1/f fluctuation in the quantum circuit. When the qubit working point departs from N g = 1 2 , the 1/f noise significantly reduces T 2 [14,15,16]. Nevertheless, the influence of the 1/f noise may be weakened via stabilizing the superconducting qubit to a high-Q resonator [41] and the feedback control method [42,43], which potentially extends the dephasing time T 2 .
After a long enough time (t T 1 ), the multiple-qubit system loses the memory of its initial condition and the ensemble expectation value j z (t) E approaches a steadystate value [44,45]. We use the symbol t → ∞ to denote a large time scale after which the distribution of j z (t) becomes steady, i.e., the ensemble-averaged observables reach steady-state values. As shown in figure 4(a), j z (t → ∞) E is equal to zero and independent of the coupling parameter η. Nevertheless, the steady-state histogram of the ensemble of j z (t → ∞) relies on η. As depicted in figure 4(c), the trajectory ensemble of j z (t → ∞) is distributed over the entire range from -1 to 1 in the weak-coupling limit (η 1). In contrast, as the qubit-qubit interaction is increased, although the qubits are still equally populated on |0 and |1 , the spread of the statistical distribution of j z (t → ∞) is narrowed, indicating the reduction of measurement uncertainty and the enhanced interqubit correlation. Generally, the expectation value of a product of two operators (O 1 and O 2 ) is different from the product of the expectation values of individual operators, i.e., ψ|O 1 O 2 |ψ = ψ|O 1 |ψ ψ|O 2 |ψ , which is attributed to the correlation between two quantities. Accordingly, since the qubits are coupled via the zz-interaction [see equation (8)], we define to measure the total interqubit correlation in the system [46]. For independent qubits one has C zz (t → ∞) = 0. As exhibited in figure 4(b), C zz (t → ∞) approximates zero in the limit of η 1, indicating that the qubits behave independently. As η is increased, C zz (t → ∞) goes up strongly and is saturated eventually.
We should note that C zz (t) differs from the variance ∆j 2 z (t) of the ensemble of trajectories j z (t), the definition of which is ∆j 2 z (t) = j 2 z (t) E − j z (t) 2 E . ∆j z (t) weights the ensemble spread of j z (t). Figure 4(c) shows the time-dependent ensemble distribution of j z (t) and steady-state histogram of j z (t → ∞) for several different η. The distribution of j z (t) diffuses rapidly when η is increased from zero, corresponding to the strong reduction of T 1 . However, the width of steady-state distribution of j z (t → ∞) becomes narrow, i.e., ∆j z (t → ∞) is suppressed, meaning that the interqubit correlation is enhanced. For η 1, the ensemble of j z (t) spreads slowly [see j z (t) with η = 0.77 and η = 7 in figure 4(c)] and the histogram of j z (t → ∞) barely changes compared with that of η ≈ 1, denoting the extension of T 1 and the saturation of C zz (t → ∞).
Actually, the spreading rate of the histogram of j z (t) corresponds to the relaxation rate of the multi-qubit system. As shown in figure 4(c), the system is initially prepared in the ground state |1 1 ⊗. . .⊗|1 N . The histogram width of the trajectory ensemble j z (t) at t = 0 is zero in principle. As the time t is increased, the histogram of j z (t) becomes broader due to the fluctuation in the quantum circuit. For the system with small (large) T 1 , the histogram of j z (t) reaches the steady-state distribution fast (slowly).

Inhomogeneous Circuit
So far, we have only focused on the homogeneous system with identical qubits. However, in the practical circumstance it is extremely difficult to fabricate an array of identical Josephson junctions. Thus, we have to return to the inhomogeneous system described by equations (7) and (8). To take into account this nonideal-fabrication-induced inhomogeneity, we assume a normal distribution with a standard deviation λ for different junction areas whose mean value is normalized to be unity [47]. We use C j and E J to respectively denote the mean values of self-capacitances C j,k=1,...,N and Josephson energies E J,k=1,...,N , i.e., C j = 1 N k C j,k and E J = 1 N k E J,k . The corresponding standard deviations are given by λC j and λE J since both C j,k and E J,k are proportional to the area of the k-th junction. For the charging energies E C,k=1,...,N , the mean value is E C = (2e) 2 2C Σ with C Σ = C g + C j + C c and the standard deviation is derived as λ C j C Σ E C in the limit of C Σ C j . In addition, the ensemble average . . . E needs to involve this extra inhomogeneity, for which we assume that a group of N -qubit systems are prepared in the same initial condition and E C,k=1,...,N and E J,k=1,...,N of each system fulfill the corresponding normal distributions.
As shown in figure 5(a), the inhomogeneity caused by nonidentical Josephson junctions strongly accelerates the relaxation of the collective Rabi oscillation even in the weak-coupling limit (η 1). This corresponds to the inhomogeneous broadening, where the unsynchronized Rabi dynamics of individual qubits destructively interfere with each other. Moreover, the extra inhomogeneity reduces the interqubit correlation and broadens the trajectory distribution [ figure 5(b)]. Therefore, maximally minimizing the inhomogeneity arising from the defective fabrication technology becomes indispensable for the application of multiple-qubit circuits in QIP.
Nevertheless, this inhomogeneous system may be still potentially applied to study the quantum many-body localization. We choose two spin states of the k-th qubit as | ↑ k = 1 √ 2 (|0 k + |1 k ) and | ↓ k = 1 √ 2 (|0 k − |1 k ) and define the new x-and ycomponents of Pauli matrix asσ x,k ≡ σ z,k andσ z,k ≡ σ x,k . Then, the inhomogeneous Hamiltonian H can be mapped onto the Ising model [48] corresponds to the long-range spin-spin interaction, B = E J plays the role of the external field, and the inhomogeneity D k = E J,k − E J acts as the site-dependent disordered potential. Unlike the superconducting circuit in [49], the nearest-neighbor J k,k±1 do not dominate. The last term in H Ising with vanishes at the sweet point N g = 1 2 . However, the nonzero δN g (t) deviates H Ising from the standard disordered Ising model. In addition, the disorder term D k is sampled from a normal distribution, rather than a uniform random variable [49].
According to [48], in the weak-coupling limit (η 1), we have λE J max(J k,k ) and the system may emerge the multiple-spin localization, which is quantified by the normalized Hamming distance [50]  with the initial state ψ(0) being the Néel state. D(t → ∞) arrives at 1 2 for a thermalizing state and remains 0 at a fully localized state. For the disordered dissipative system with (λ = 0, δN g (t) = 0), the Hamming distance D(t → ∞) always approaches 1 2 , reaching the thermal equilibrium [see figure 5(c)]. This is attributed to the effect of environmental fluctuations δN g (t). Thus, besides QIP, suppressing the environmental noise is also the key issue for the application of superconducting circuits in the manybody simulation. Figure 5(d), where we have artificially set δN g (t) = 0, exhibits that the disordered nondissipative system with (λ = 0, δN g (t) = 0) stays nearly fully localized product states in the weak-coupling limit (η 1) while D(t → ∞) rises gradually up to 1 2 as the qubit-qubit interactions become strong. Using the diagonalization method, we have also checked the statistics of the ratio parameter r of adjacent energy-level gaps for the disordered nondissipative system [51,52]. It is seen that r follows the Poisson distribution when η 1 [the inset of figure 5(d)], manifesting the many-body localization, and apparently violates the Poisson distribution for a larger η, suppressing the localization effect.
As experimentally demonstrated in [49], the many-body localized state may be still attainable within a short time scale when the thermalization rate of the system with (λ = 0, δN g (t) = 0) coupling to the environment is much slower than the localization rate of the system with (λ = 0, δN g (t) = 0) (i.e., the localization rate measures how fast a disordered nondissipative system reaches the localized state from an initiallyprepared state). The localization rate γ of the system with (λ = 0, δN g (t) = 0) may be roughly estimated by assuming that the envelop of D(t) follows the exponential law, i.e., D(t → ∞) E (1 − e −γt ) [see the dashed lines in figure 5(d)]. The time scale that measures the thermalization rate of the system with (λ = 0, δN g (t) = 0) coupling to the environment is given by T 1 (λ = 0). The many-body localized state is potentially observed in the disordered dissipative system when γ T −1 1 (λ = 0) which may be fulfilled for a large η.
In comparison between figure 5(c) and 5(d), we find that for η 1 the Hamming distance D(t) of the disordered dissipative system rapidly reaches a metastable value smaller than 1 2 (i.e., the system arrives at the localized state) and then slowly grows up to 1 2 (i.e., the system approaches the thermalizing state). Thus, the many-body localization phenomenon can be observed within a short time scale for η 1. Increasing the disorder strength D k , i.e., increasing λ, may enhance the localization rate of the inhomogeneous system, leading to an easy access to the many-body localized phase. However, a large λ strongly reduces E C,k E J,k , i.e., some qubits may not operate in the charging limit ( E C,k E J,k 1).

Conclusion
We have studied a multiple-charge-qubit system, where the Cooper-pair boxes are all capacitively linked. The qubit dynamics are interfered with environmental fluctuations which are mapped onto the gate voltage source that biases all qubits. The collective Rabi oscillation has been numerically simulated for both homogeneous and inhomogeneous systems. We find the interqubit coupling strongly varies the energy-relaxation rate of the quantum circuit consisting of identical Josephson junctions. For the weak coupling system, the qubit relaxation time T 1 is significantly reduced since the dynamics of one qubit is unavoidably influenced by the fluctuations of other qubits. In contrast, T 1 of the homogeneous system in the strong coupling regime can be enhanced to a value much larger than that of a free qubit. This result is caused by the strong energy-level shifts of qubit states (i.e., the large interaction-induced detunings) and consistent with the expectation of Fermi's golden rule. In QIP, transferring quantum information between two qubits in a homogeneous multi-qubit network is an essential process. The resulting fidelity is limited by the system's decoherence time. Thus, the enhancement of T 1 is crucially relevant to QIP.
The nonideal-fabrication-induced inhomogeneity always expedites the multi-qubit system's collective decay. In addition, we mapped the inhomogeneous system onto the disordered Ising model to probe the many-body localization effect. This enables us to investigate the role of the environmental noise in the quantum simulation. For the inhomogeneous system with the localization rate faster than the thermalization rate, the many-body-localization regime is still accessible.