Synchronization of two cavity-coupled qubits measured by entanglement

Some nonlinear radiations such as superfluorescence can be understood as cooperative effects between atoms. We regard cooperative radiations as a manifested effect secondary to the intrinsic synchronization among the two-level atoms and propose the entanglement measure, concurrence, as a time-resolved measure of synchronization. Modeled on two cavity-coupled qubits, the evolved concurrence monotonically increases to a saturated level. The finite duration required for the rising to saturation coincides with the time delay characteristic to the initiation of superfluorescence, showing the role of synchronization in establishing the cooperation among the qubits. We verify concurrence to be a good measure of synchronization by comparing it with asynchronicity computed from the difference between the density matrices of the qubits. We find that the feature of time delay agrees in both measures and is determined by the coupling regimes of the cavity-qubit interaction. Specifically, synchronization is impossible in the weak coupling regime.

Centuries ago, Huygens studied the correlation among the motions of pendulums and discovered the synchronization pattern of these individual oscillators under the influence of a common oscillator they are coupled to 1 . How synchronizations arise in different situations has since remained a problem of interest 2 . In recent years, the study of the classical phenomenon has been revived under the quantum regime. Synchronization is observed between a pair of nanomechanical oscillators 3 and on the motions of the Cooper pairs among Josephson-insulated superconducting islands 4 . It is also ubiquitously predicted between a qubit and an oscillator 5 , between a cavity field and an oscillator 6 , between two oscillators 7 , among a trapped group of cold atoms 8 , and even between a quantum Van der Pol oscillator and an external drive 9 .
Here we study the synchronization of two qubits coupled indirectly to each other through a cavity field, in specific relevance to the quantum optical phenomenon of superfluorescence 10,11 . This nonlinear fluorescent effect embodies Dicke's formulation of superradiance 12 , where the nonlinearity exemplifies in the N 2 -dependence of the emitted radiation intensity from an ensemble of N atoms. Furthermore, the emitted photonic pulse peaks after a positive time delay 13,14 that bears no direct dependence on the atom-field coupling strength, showing the collective interatomic motion to be a nonlinear process. Experimentally recorded on hydrofluoric gas 15 , cesium 16 , and most recently rubidium vapor 17 , this delay shows the necessity of a finite time during which the atoms establish cooperation before the radiation is initiated 18 .
Such a delay is also characteristic of synchronization: it takes finite time for the classical Huygens pendulums to reach in-phase oscillations. If one regards the correlated atoms as quantized Huygens oscillators concentrated on the lowest two levels at the weak-energy limit, it is not hard to find the resemblance between cooperated radiation and synchronization. In this paper, we formalize the study of this resemblance by analyzing the entanglement dynamics [19][20][21] of two cavity-coupled qubits, which is the minimum-dimensional system that emulates the structure of an atomic ensemble undergoing cooperated radiation.
Many different metrics have been proposed over the years to measure synchronization of different quantum systems, such as a set of two-level systems in a shared bath 22 , quantum dots 23 , and resonator modes 24 . We note that synchronization as a quantum correlation can be studied using quantum discord 25 , as it was for two oscillators 26 , and non-locality. The advantage of entanglement approach is avoiding the unnecessary constraint to two-party correlation. We use specifically, in contrast to quantum discord, concurrence 27  www.nature.com/scientificreports/ is broadly applicable, especially for multipartite system. After it was introduced for two-qubit systems 27 , it was soon extended to the cases of three qubits 28 and N-qubits 29 . Concurrence has then been generalized to two multi-level systems 30,31 and finally to multipartite multi-level systems 32 . Given its generality, concurrence permits synchronization to be exemplified not only between two qubits, but among three parties that include the cavity field. In addition, synchronization is present in a Hilbert space of finite dimensions 33 , so non-locality is not an essential factor to initiate synchronization. Further, concurrence is proved to be a good measure for static analysis of collective phenomena such as radiation 34,35 . Therefore, we extend the employment of multipartite concurrence to the dynamical analysis of the finite-dimensional cavity-qubit system. It is found that the entanglement among the components universally begins at a zero level and rises to a saturated value after a finite delay. The length of the delay and the saturated level depend nonlinearly on the cavity-qubit coupling strength. We use numerical simulations to show that the nonlinear dependence is divisible into weak, strong, and ultrastrong coupling regimes, verifying the role of operation regimes 36 in the dynamics of superconducting circuits. Under all coupling strengths, the saturated entangled state is sustained thereafter, where oscillations in concurrence only exists locally about the saturation offset with an amplitude much smaller than the offset. These distinctive features of a synchronized state stand in constrast to the sudden death or death and revival 37 one usually sees in entanglement evolution.
To verify the coincidence between the maximization of concurrence and the dynamic process of synchronization, we introduce a synchronization measure computed from the density matrix, which is modified from the synchronization measure introduced on the (x, p)-quadratures of quantum oscillators 24 . The transition point in time produced from the synchronization measures matches exactly the delay time found in the concurrence evolution, proving that the qubit-to-qubit synchronization is well registered in the entanglement. Moreover, since two-qubit systems on a superconducting circuit can produce superfluorescent pulses 38 , the synchronization delay is associated with the initiation of superradiance, thereby establishing the dynamic correlation between entanglement and the atomic cooperation for collective phenomena. Results the tripartite system. The derivation for the dynamics discussed above is modeled on a generic cavity quantum electrodynamic (QED) system where each qubit is coupled to the cavity through dipole-field interaction under the rotating wave approximation. The parameters adopted for the numerical analysis are sourced from the superconducting circuit implementation 39,40 of the cavity QED system. Hence, the tripartite system can be illustrated from the model Fig. 1, where the cavity is indicated by the stripped waveguide and the qubits are located at the anti-nodes of the cavity field to ensure maximum coupling.
The total system Hamiltonian H = H 0 + H int + H ext is composed of three parts: the free Hamiltonian, the interactions among the system components, and the external driving, which reads, respectively, ( = 1) In H 0 , ω c denotes the frequency of the cavity mode and A ( B ) denotes the transition frequency of the left (right) qubit, associated with the Pauli matrix σ A,z ( σ B,z ). In H int , η A (η B ) denotes the coupling strength to the left (right) qubit. In H ext , the external driving field has frequency ω D and driving strength ε D .
The combined system of a cavity and two qubits has its bare states described by the tensor product state {|e A �, |g A �} ⊗ {|n�} ⊗ {|e B �, |g B �} , where |e A � ( |g A � ) denotes the excited (ground) state of the left qubit; |n� denotes the Fock number states of the cavity mode; |e B � ( |g B � ) denotes the excited (ground) state of the right qubit. To simplifying the notation, we omit the subscripts A and B when writing the product states and let the first letter www.nature.com/scientificreports/ denote the state of the left qubit, the middle letter that of the cavity mode, and the last letter that of the right qubit (e.g. |e, n, g� = |e A � ⊗ |n� ⊗ |g B �).
The free Hamiltonian H 0 and the interaction Hamiltonian H int constitute a closed subsystem, for which there exist dressed states that diagonalize H 0 + H int . To find an analytical expression for the dressed states, we consider the sets of energy-conserving states |e, n, g� , |g, n + 1, g� , and |g, n, e� , which are resonant within single-photon processes, to contribute to a dressed state for each n. In other words, the state |e, n, e� which is resonant with |g, n + 2, g� through a double-photon process is omitted, allowing the derivation below to be purely analytical. Single-photon processes among the collective states of an ensemble of two-level systems are responsible for the majority of qubit operations useful for quantum information processing, which includes generating the entangled W-state 34 and the multistable attractor states 41 . Such collective states are also experimentally proven on superconducting qubits 42 .
These single-photon resonant states form an invariant subspace, for which the closed Hamiltonian consists of 3 × 3 symmetric block matrices. Therefore, we have the eigen-equation where the eigenvectors |u |g, n, e��g, n + 1, e| + |e, n, g��e, n + 1, g| + √ 2|g, n + 1, g��g, n + 2, g| www.nature.com/scientificreports/ The permitted dressed level transitions induced by the external driving can be found by substituting Eq. (14) into Eq. (3). In the following, we consider the weak-energy limit where the transitions are confined to the lowest two clusters of states ( n = 0 and n = 1 ). This consideration greatly simplify the numerical analysis of the evolution by confining the matrix representation of the equations of motion to 6 dimensions. It is also realistic since higher clusters require the much less likely occurrence of multi-photon resonance driven by the external ε D . Though obtaining states of high photon number are possible with the assistance of auxiliary qubits 43 , they are much shorter-lived. We therefore write the total Hamiltonian as in the dressed space. Introducing the time-dependent state vector in the confined state space and applying it to the Hamiltonian above, one has the Schrödinger equations of the time coefficients where l denotes the weight of the summation over the index l for the system component A (the left qubit), B (the right qubit), or C (the cavity), i.e. A = B = 1 and C = √ 2 . In the equations, we observe the Einstein summation convention.
In the rotating frame j t} , the coupled equations can be written as the linear homogeneous system of differential where the square brackets [·] indicate 3 × 3 submatrices in the matrix A with j and k being the row and the column indices, respectively. We denote for the detuning between the driving and the dressed states. Since A is integrable, then solving the linear system for {c j , d j } and expanding the dressed states by using the bare states in Eq. (16), one can find the expansion coefficients γ of the state vector back in the bare state space. evolution of the state vector. To see that the evolution of the state vector can initiate the synchronization of the two delocalized qubits, we assume the cavity mode is initially driven by the external field to reach a partial population inversion while setting the qubits initially at the ground. In other words, the expansion coefficients at the initial moment are: γ We plot out the evolutions of these coefficients in Fig. 2, using the experimentally accessible parameters of superconducting charge-phase qubits 40 : � A /2π = � B /2π = 6.1 GHz, η A /2π = η B /2π = 500 MHz, and ω c /2π = 6.32 GHz. The frequency of the external field is maintained at ω D /2π = 5.3 GHz. The lower set of states with n = 0 is given in Fig. 2a whereas the upper set with n = 1 is given in Fig. 2b. Note that, for the implementation of cavity QED system on a superconducting circuit, the cavity field refers to the standing microwave field that is sustained in the stripline (Cf. Fig. 1) and coupled capacitively to the Josephson-junction qubit. Hence, by distancing qubits and placing them at different positions relative to the cavity nodes, it is possible to realize arbitrary cavity-qubit interactions. For instance, maximal coupling can be obtained by placing the qubits at antinodes 42 .
We observe that for both the lower set and the upper set of states, there exists a transition point of the oscillations of the coefficients, which is located at about 3.7µs in the plots. In particular, γ B , and γ (1) C are transited from an amplifying region to a region of saturated oscillation envelope. The contrasting behavior of the two sets of coefficients demonstrates that the energy excitation that exists in the cavity mode is transferred to the left and the right qubits whose complementary oscillations imply the build-up of the entanglement between them. www.nature.com/scientificreports/ Bipartite and tripartite concurrences. To fully capture the evolution characteristics of the two cavitycoupled qubits from a holistic point of view, we apply two entanglement measures-bipartite concurrence and tripartite concurrence-to the state vector of the total system. The bipartite concurrence quantifies the inseparability of the joint pure state of two coupled systems of arbitrary dimensions by inverting the density matrix. For our case here, the joint state is the product state |ψ AB � of the indirectly coupled left and right qubits. Thus the inversion is conducted through the superoperator S D 1 ⊗ S D 2 w he re t he d i me ns i ons D 1 = D 2 = 2 and t he bip ar t ite c onc u r re nc e is d e f i ne d as C 2 (ψ AB ) = √ �ψ AB |S 2 ⊗ S 2 (|ψ AB ��ψ AB |)|ψ AB � . Given the consideration of pure states, for which trρ 2 = 1 and trρ 2 A = trρ 2 B , the definition reduces to C 2 (ψ AB ) = 2 1 − tr ρ 2 A where ρ A = tr B (tr C (|ψ��ψ|)) is the reduced density matrix of the left qubit.
Applying |ψ(t)� in Eq. (20) to the formula, we derive the evolution of the bipartite concurrence, shown as the blue curve in Fig. 3. It becomes apparent that the transition point that manifests in Fig. 2a,b signifies the concurrence reaching a maximum after a gradual monotonic increase in the oscillation envelope. This maximum concurrence is retained thereafter. The finite delay time τ D = 3.7µs that the concurrence spends to reach its maximal value reflects the time the two qubits use to reach a maximal synchronization through their mutual couplings to the cavity mode.
The cavity mode plays an active part in initiating the entanglement between the two qubits. From the entanglement-theoretic point of view, the concurrence is distributed among the qubits as well as the cavity. Taking away the pairwise entanglements between any two parties in the tripartite system, one obtains the residual concurrence that remains as an equally distributed entanglement among all three parties 28 . Extending the original formulation on three-qubit systems, we generalize the inversion operations for two arbitrary-dimensional systems given above to three arbitrary-dimensional system. That is, we introduce the superoperator for our D 1 × D 2 × D 3 dimension tripartite system, where D 1 = D 3 = 2 for the qubits and D 2 = n for the cavity mode. In Eq. (21), I denotes the identity matrix while ρ A , ρ C , ρ B , ρ AB , ρ CB , and ρ AC denote the reduced density matrices of the components and the two-component subsystems. Applying the inversion, we thus derive a tripartite residual concurrence Again, using Eq. (20), we plot the tripartite concurrence as the red curve in Fig. 3. One can verify from the plot that the residual concurrence evolves in a similar fashion, which contains a signifying transition point at the exactly same location τ D as that of the bipartite concurrence. Before τ D , it arises from a zero value under a similarly increasing envelope whereas, after τ D , it retains a non-zero saturated value. The identical delay time again demonstrates the duration that the system components spend on cooperation before maximal synchronization is reached. www.nature.com/scientificreports/ Comparing Figs. 2 and 3, one sees that the energy quantum first dwells on the cavity mode ( |g, 1, g� and |g, 2, g� ) without being emitted and absorbed by the qubits. Only when the two qubits start to establish a cooperated motion does the qubit-cavity-qubit resonance become effective such that the qubits be excited to their respective excited states |e, 1, g� and |g, 1, e� . The entanglement is also established among the three components when the excitation commences.
The concurrences plotted in Fig. 3 are computed upon a symmetric setting of system parameters: the level spacings and the coupling strengths of the qubits are assumed identical. The tripartite concurrence of an asymmetric scenario with the right qubit level spacing raised to � B /2π = 7.1 GHz is shown as the green curve in Fig. 4 while the rest of the parameters remain unchanged. For comparison, the symmetric case is plotted as the red curve in the background. We observe that the delay to saturated synchronization is inversely correlated with the larger eigenfrequency out of the two qubits. For the case in Fig. 4, increasing B reduces delay time τ D . On the other hand, symmetric settings lead to maximal synchronization at saturation. The asymmetric case given by the green curve has the saturated synchronization reduced to a lower level. Simulation under parameters set to varied values (not shown in figures) verify these observations. We note that the dynamics originated from the equations of motion (17)- (18) does not take into consideration the relaxation effects which would require a master equation approach under the presence of an environment. Their omission greatly simplifies our derivations and assists the illustration of the temporal features during the synchronization. These features would not be obscured even when environmental effects are observed since their time scales are much shorter than accessible coherence times on superconducting circuits. For instance, the synchronization delays in both settings of Fig. (4) appear within the first 4 µs after initiation. Whereas, superconducting qubits can achieve coherence time up to 100 µs using fabrication techniques 44 and even up to hundreds of µs when an error-correcting feedback mechanism is implemented 45 .
The synchronization between the qubits is also affected by how strong they are driven by the cavity mode, i.e. the coupling strengths η A and η B . Shown in Fig. 5 for the symmetric scenario η A = η B in a semilog plot, the greater is the coupling, the lesser is the delay τ D . The numerical fit shows that the delay obeys a quadratic relation over the exponential increase in coupling strength. After the delay, a shorter delay time is associated with a higher saturated level of concurrence, showing a stronger synchronization between the qubits are reflected in both short delay and higher entanglement measure.
Asynchronicity. Multi-partite concurrence as a measure of synchronization reveals a gradual increase between two qubits in the time domain, explaining the existence of a delay in the superfluorescent pulse of cooperated radiation from a system-intrinsic point of view. This synchronization is affected by many factors, among which the symmetry of the system parameters plays an important part. Tuning the system from a symmetric setting to an asymmetric setting is accompanied by tuning the transition rates of the qubits from a synchronous setting to an asynchronous setting. For the latter, we refer to the scenario where the population of the left qubit oscillates at a Rabi frequency not synchronous to that of the right qubit. www.nature.com/scientificreports/ However, since two oscillators sharing a common oscillating platform are able to synchronize after certain time duration according to classical mechanics, we expect the qubits sharing the cavity resonator would behave similarly. To precisely describe the transition process from asynchronous regime to synchronous regime, we extend the quantum synchronization measure introduced in Ref. 24 for continuous variable systems to discrete systems. We consider instead the measure of asynchronicity that compares the difference between two density matrices for two two-level systems.
To consider the scenario where the sychronization is initiated from the cavity field as the medium, we let the initial state be only populated in the first and second excited state of the cavity mode, i.e.   www.nature.com/scientificreports/ |ψ(0)� = 3|g, 1, g�/ √ 10 + |g, 2, g�/ √ 10 . In a superconducting circuit where the cavity mode corresponds to the photonic number state of the stripline resonator, the population is controlled through a resonant microwave pulse, which has a time duration commensurate with the cavity Rabi period, to be fed into the waveguide input 46 . For example, a π-pulse sends half of the population from |0� to |1� . If imaginary phases in the coefficients of the eigenstates are desired, they can be fine-tuned with the assistance of an auxiliary qubit 43 . |2� is here partially populated to emulate the resonance process where the equally spaced levels of the cavity mode has |1� resonant with |2� simultaneously with |0� . Given the initial condition, we find that a symmetric setting A = B would always lead to a zero asynchronicity throughout independent of the coupling strengths η A and η B . When A = B , the asymmetry leads to coupling-dependent asynchronicity, as shown by the plots given in Fig. 6 for five settings of coupling strengths.
No matter the coupling strength, there exists a transition point after which the asynchronicity remains at a stable value. This transition point is identical to the transition point shown in Fig. 4 (green curves) where the tripartite concurrence reaches a maximal value. The coincidence verifies our expectation that the synchronization is maximized when the asynchronicity is minimized. Therefore, synchronization between two qubits reflects the dynamic identity of the two qubits.
Before reaching the stable minimal value, the asynchronicity increases from a non-zero value for a certain duration, which are spent on the cooperation by the qubits. When the coupling is sufficiently weak (below η ≈ 200 MHz), the minimal stable value is almost vanishing (below 10 −5 ). When the coupling becomes stronger, the feedback from the cavity mode to each of the qubits becomes adverse to the synchronizing motion. However, the feedback effect is not linear. Plotted as a function of the dimensionless coupling strength η/� A at � A /2π = 6.1 GHz and � B /2π = 7.1 GHz in Fig. 7, the stable minimal value Ā of asynchronicity first retains a negligible value in the weak coupling regime. At about η/� A = 0.04 , Ā starts to increase slowly until it reaches a turning point at η/� A = 0.2 . The region between η/� A = 0.04 and η/� A = 0.2 can be regarded as a strong coupling regime for synchronization. After that, Ā enters the ultra-strong coupling regime and increases with the coupling again until η/� A ≈ 0.33 , where stable minimal value of A is no longer discernible.

Discussion
In conclusion, we have studied the synchronization between the two cavity-coupled qubits using multiple concurrence measures and asynchronicity. These real-valued measures are computed as functionals of dressed state vectors that evolve in time as they are driven by an external field. In all these measures acting as time functions, we obtain consistent features of transitions from an arbitrary initial system state to a final synchronized state. The transition in time reveals a synchronization delay that the qubits use to initiate superfluorescent pulse radiation, which explains the cooperation origin of the collective effect of superradiance.
The characteristics of the synchronization process, including the delay and the value of the stabilized asynchronicity, are highly dependent on the coupling strengths of the qubits relative to the their level spacings. They demonstrate from the entanglement perspective the different behaviors that the circuit QED systems adopt when operating in weak, strong, and ultra-strong coupling regimes. In general, synchronization occurs only www.nature.com/scientificreports/ in the strong-and ultrastrong-coupling regimes whereas its level of synchronization at the final state increases exponentially with the qubit-cavity coupling strength.

Methods
For each 3 × 3 block submatrix of the LHS of Eq. (4), we first split it into the sum of two matrices, one being the identity matrix nω c I . The other matrix is the one to be diagonalized, the determinant equation of which is where denotes the eigenvalue to be determined, i.e. E (n) k = + nω c . This determinant equation is equivalent to the cubic equation whose roots can be derived by absorbing the quadratic term through the transform = x + δ/3 . The transformed equation becomes x 3 + px + q = 0 , where In the close cavity-qubit resonance region δ ≈ 0 , the discriminant D is simplified to To let the cubic equation admit three non-degenerate real roots, we consider the range that makes D < 0 . Applying the Vieta's formula, the roots x can be found with a parametric angle θ given by Eq. (6).
Using the eigenvalues given in Eq. (5) and expanding the eigenvector u (n) k in the bare state space given by Eq. (7), we have the column matrix equation www.nature.com/scientificreports/ Letting α C be the proportional constant, we find α A = −η A α C (� − ) and α B = η B α C (� + ) . Then normalizing the coefficients, their expressions are given by Eqs. (8)-(10) for the k-th energy level within the n-th cluster.