Quantum Synchronization and Entanglement of Dissipative Qubits Coupled to a Resonator

In a dissipative regime, we study the properties of several qubits coupled to a driven resonator in the framework of a Jaynes–Cummings model. The time evolution and the steady state of the system are numerically analyzed within the Lindblad master equation, with up to several million components. Two semi-analytical approaches, at weak and strong (semiclassical) dissipations, are developed to describe the steady state of this system and determine its validity by comparing it with the Lindblad equation results. We show that the synchronization of several qubits with the driving phase can be obtained due to their coupling to the resonator. We establish the existence of two different qubit synchronization regimes: In the first one, the semiclassical approach describes well the dynamics of qubits and, thus, their quantum features and entanglement are suppressed by dissipation and the synchronization is essentially classical. In the second one, the entangled steady state of a pair of qubits remains synchronized in the presence of dissipation and decoherence, corresponding to the regime non-existent in classical synchronization.


I. INTRODUCTION
The synchronization of two maritime pendulum clocks, discovered by Christian Huygens in 1665 1 , is at the basis of synchronization phenomena abundantly appearing in various systems ranging from clocks to fireflies, cardiac pacemakers, lasers and Josephson junction (JJ) arrays (see historical survey and overview in [2][3][4] ).
With the modern development of quantum computation and quantum information (see e.g. 5 ) the investigation of quantum synchronization phenomenon in a quantum world becomes of significant importance.The phenomenon of quantum synchronization earns even higher relevance for JJ arrays and superconducting qubits where dissipative and quantum effects play an important role (see e.g. 4 and 6 ).A strong coupling regime of several superconducting qubits with a microwave resonator has been realized experimentally [7][8][9][10][11] that requires to investigate also effects of strong interactions in quantum dissipative driven systems.
The essential element of synchronization in classical mechanics is dissipation that leads to a phase synchrony of coupled autonomous systems with energy supply 2 .Thus the analysis of quantum synchronization requires a use of formalism of dissipative quantum mechanics.This formalism is known and based on the Lindblad master equation for the density matrix ρ of a whole system being developed in 12,13 and reviewed in 14 .It has fundamental grounds and allows to study such complex phenomena as e.g.quantum strange attractors 15,16 .However, its numerical simulation requires integration of time propagation of the whole density matrix with N × N components that is numerically rather heavy.
Another approach to dissipative quantum evolution is based on the method of quantum trajectories 17,18 which stores only a stochastically evolving state vector of size N .An averaging over several realizations allows to characterize the probabilities from density matrix with statistical fluctuations.The application of quantum trajectories to the problem of quantum synchronization has been reported in the early study 19 where it was shown that at small dimensionless values of Planck constant ℏ the synchronization remains robust with respect to quantum fluctuations.These model studies of quantum synchronization 19 were shown to have close similarities with Shapiro steps in Josephson junctions 4 .The method of quantum trajectories was also applied to investigations of quantum synchronization of one and two qubits coupled to a driven dissipative resonator 20,21 .It was shown that the driving phase can impose quantum synchronization of qubit phases and mutual synchronization and entanglement of two qubits.In the frame of synchronization of classical systems this corresponds to the regime of phase synchronization of oscillators by external driving 22 .At the same time the description with quantum trajectories does not provide a complete information about evolution of density matrix.Thus it is important to study the phenomenon of quantum synchronization of qubits using the Lindblad master equation for density matrix.In this work we describe the results obtained with the Lindlad description.
The performed extensive numerical simulations of tens of thousands of Lindblad equation components allow to establish various nontrivial regimes of quantum synchronization and entanglements of one, two, three or four qubits strongly coupled with a driven resonator in presence of dissipation.It is shown that in the regimes of weak or strong dissipation the system behavior can be well described by semi-analytical methods of weak perturbation theory or semiclassical theory respectively.However, certain nontrivial dissipative regimes with preserved entanglement of dissipative qubits remain non-accessible for such a semi-analytical description.
The article is composed as follows: Section II describes the model, numerical methods and semi-analytical approximations, Section III presents obtained numerical results and RWA validity tests of RWA, the results and validity of the weak damping rate equation approach are given in Section IV and those of strong damping semicalssical regime in Section V, the regimes beyond the validity of semi-analytical approaches are analyzed in Section VI, synchronization of several qubits is described in Section VII and discussion is given in Section VIII.

II. MODEL DESCRIPTION
We consider a model of several qubits (spin half systems) interacting with one harmonic cavity driven by an external monochromatic field.The Hamiltonian of this system is as follows: where ω 0 is the frequency of the cavity, ω is the frequency of the driving field and F is the driving strength.The operators â+ and â are the rising/lowering operators for the cavity photons.The qubits in this model are indexed by l (we consider various numbers of qubits, from 1 to 4), their contributions to the Hamiltonian are given by the Pauli matrix operators where Ω l is the energy splitting (Zeeman splitting for spin systems) and λ l are their coupling strengths with the cavity.For F = 0 and one qubit the Hamiltonian ( 1) is reduced to the Jaynes-Cummings model 33 appearing in many systems of quantum optics and other physical systems 34,35 .It was experimentally realized with Rydberg atoms in a resonance cavity 36 .The unitary behavior of the Jaynes-Cummings model with monochromatic driven cavity was studied in 37 .However, here we consider the model with one or several qubits in presence of driving and dissipation that corresponds to a real situation of superconducting qubits 6 .
In the rotating wave approximation (RWA) the Hamiltonian in Eq. (1) becomes: where we kept only the resonant terms in both the driving field and spin-cavity interactions.The Pauli matrices σ+ l and σ− l are the rising/lowering spin operators for spin/qubit l.An advantage of the RWA is that the time dependent Hamiltonian Eq. ( 2) can be made stationary by moving to the rotating frame using the unitary transformation: Û (t) = exp iωâ + ât + iω l σ+ l σ− l t which acts on the density matrix ρ moving it into the rotating frame through ρR = Û (t)ρ Û + (t).It can be checked that this leads to a stationary rotating frame Hamiltonian 34,35 : for simplicity we assumed the same coupling strength between qubit and cavity λ for all qubits.We describe dissipative effects in the framework of the time dependent Lindblad equation 14 , with zero temperature dissipative terms for both qubits and cavity dynamics: where L d gives the dissipative part of the Lindblad dynamics: with γ being the dissipation rate of the cavity and γ s the dissipation rate for qubits.Under the RWA the Lindblad equation Eq. ( 4) becomes: since the Hamiltonian ĤR is stationary the steady-state value of ρR can be found by setting the left hand side of Eq. ( 6) to zero.The advantage is that finding the steady-state does not require numerical integration of the equations of motion Eq. ( 4) and can be found more directly by solving the matrix equation: A. Numerical methods In general Eq. ( 7) is a system of linear equations, with super-operator acting on the unknown density matrix.If the density matrix has size N × N , the superoperator will be a matrix of size N 2 × N 2 and exploiting its sparse structure of the Hamiltonian ĤR and of the dissipative Lindblad terms, finding the solution of Eq. ( 7) becomes numerically prohibitive.Already building the sparse super-operator matrix can be numerically demanding in terms of memory.
Thus we used an alternative approach, exploiting as much as possible the similarity between the Lindblad equation and Sylvester equations which are matrix equations of the form: where Â, B, Ĉ are known matrices and X is the unknown matrix to be found.For this class of matrix equations more efficient solution algorithms are available 38 , in our case we used the Bartels-Stewart algorithm 38 which has a N 3 computation cost where N × N is the size of the density matrix.
The general steady-state Lindblad equation: is not of Sylvester type (here Ĥ is the Hamiltonian and Lk are dissipative operators).So a direct application of the Bartels-Stewart algorithm is not possible.We thus split Lindblad super-operator into a first term which looks like a Sylvester equation and a remainder: The steady-state density matrix L(ρ) = 0 will be a fixed point of: It is thus possible to iterate from an initial guess ρm solving a series of Sylvester equations: We have found, that usually a few or at most tens of iterations are enough to find the fixed point with very high accuracy.The initial guess ρ0 for the iteration can be chosen from one of the two approximate methods presented in the following Sections, we typically used ρ0 obtained from summation of the rate equation perturbation theory.Due to the availability of a good initial ρ0 in our simulations the number of Sylvester equation iterations required for convergence was small.Algorithms for solving generalized forms of the Sylvester equation have also been reported 39 , they can maybe offer a faster converging iterative solution in cases where a good initial guess is not available.It is also possible to find the steady-state by numerical integration of time dependent Lindblad Eqs.(4), (5) or Eq.(6).We used a standard odeint integration library 40 for the integration of equations of motions using a sparse matrix representation for all operators acting on the density matrix which is represented as a dense matrix.Main integrator was an adaptive Runge-Kutta-Dorpi fifth order stepping algorithm, although other integration schemes were also tested with similar results.In the time dependent case it is possible to check the validity of the RWA approximation since we have access to the full time dependent evolution of the density matrix including its oscillations around the RWA steady-state, while integration of the stationary Eq. ( 6) converges to the solution of Eq. (7).

B. Rate equation perturbation theory
The above approach for steady-state computation can still be numerically costly if an important number of Sylvester equation inversions is required.It is thus important to have a good initial guess for the density matrix to start the iteration.In the limit where the damping terms in the Lindblad dynamics are weak compared to characteristic frequencies of the Hamiltonian dynamics it is also natural to use an approximate solution based on a formal expansion of the density matrix in powers of the dissipation operators.In such an approach the dominant terms of the density matrix correspond to weighted eigenstates of the Hamiltonian dynamics and it can be convenient to think of this expansion as a series of rate equations describing the population of the RWA Hamiltonian eigenstates.We use this approach to obtain good initial guesses for density matrix, it also allows us to investigate if the steady-state of model Eq. ( 6) can indeed be described by a perturbative approach based on Hamiltonian eigenstates.Heuristically, this is justified only when there is a clear separation of times scales between the Hamiltonian dynamics and the slower dissipative processes.Our calculations allow to justify this point quantitatively.
We describe this approach starting from the general steady-state Lindblad equation introducing a formal expansion parameter ϵ to keep track of the order of the dissipative parameters in the expansion: here Lk are as previously the relaxation operators and Ĥ is a steady-state Hamiltonian which in our case is given by RWA Eq. (3).
To find the steady-state we need to solve L(ρ) = 0, and we proceed making a formal expansion in ϵ: The equations for the lowest order approximation ρ0 are the following: where the operator P D acts on a density matrix defined in the eigenbasis of Ĥ by keeping only its diagonal part.The first equation L H ρ0 = 0 imposes that ρ0 is diagonal in the eigenbasis |n⟩ of the Hamiltonian Ĥ = n ϵ n |n⟩⟨n|.Thus ρ0 has the form ρ0 = n P n |n⟩⟨n| where P n can be interpreted as the probability of eigenstate |n⟩.The last two equations are thus a rate equation determining P n and a normalization condition n P n = 1.The solution of this equation involves building the matrix: We then find non zero solutions K|ψ⟩ = 0, such solutions always exist since ⟨1, 1, ..., 1|K = 0.The recurrence equations for the next orders are: which defines the off-diagonal components of ρj+1 (in the eigenbasis of Ĥ) but leaves its diagonal part undefined.The diagonal part of ρn+1 is can be found requiring: this can be viewed as higher order corrections to the rate equations changing the occupation state probabilities P n .
Finally the matrix ρj+1 is then made traceless by the transformation this transformation is allowed because the induced error is of higher order: This can also be viewed as an order by order normalization of the density matrix Trρ = 1.
We denote as (R) the numerical result of the summation of these series.Convergence of this expansion can be checked by evaluating the residual error in Eq. ( 19) at each step of the summation checking that the error decreases as the number of terms grow.If the series is convergent we continue summation until close to machine precision in the residual error which is typically reached from the first ten terms, if the series is divergent we stop summation when the residual error starts to increase.
We find that the direct summation of (R) seems to have a rather small convergence radius requiring the amplitude of the dissipation being small compared to energy level splittings: In practice some of these transitions are almost forbidden by selection rules involving several spin flips and transitions between several oscillator quanta.We thus try to improve the radius of convergence of rate equation series by incorporating dissipative terms into L H .This is done by observing that L H acts on the density matrix |n⟩⟨m| as: We can thus formally include into L H all the dissipative terms which are diagonal in the same basis.This is done by introducing a new (super)operator L γK , which in the eigenbasis of H is defined as follows: The expansion (R) is then performed with the updated operators Lγ and LH given by: This approach has some similarity with a Jacobi preconditioning step (see e.g. 41) which can be used to improve the radius of convergence of iterative algorithms to solve linear equations.More physically it corresponds to absorbing some terms of the Lindblad equation into a non-self adjoint dissipate Hamiltonian.We will show that this updated approach (R) * improves the radius of convergence of the perturbation series and for low friction divergence seems to occur only close to the exact resonance ω = ω 0 when the RWA Hamiltonian Eq. ( 3) becomes extremely degenerate.In this case dissipative dynamics can no longer be described satisfactorily based only on populations of the eigenstates of the Hamiltonian.Another point of view on the finite convergence radius of the rate equation series is that formally the series can also be summed for negative values of the dissipation rates which then correspond to gain.The series on the dissipative and gain sides have the same radius of convergence, but if the gain is larger than some threshold the oscillator energy will diverge.This argument provides a maybe more physical explanation for finite convergence radius near resonance.

C. Semiclassical approximation
While the perturbative rate equation approach works at weak dissipation its convergence fails in the opposite limit of strong damping.In this regime, an approximate description of the system steady-state can be obtained from a semiclassical trial density matrix.The advantage of this approach is that the functional to be minimized can be explicitly computed and the parameters have a clear interpretation in terms of cavity coherent states and mean spin orientations.
We start by presenting this functional in the simple model case of a driven cavity, for which the RWA steady-state is a pure coherent state and thus its minimization gives the steady-state exactly and then generalize this approach to a cavity coupled to qubits.
Our approach is to minimize: for some trial form of the density matrix.
We first consider the simple model of a driven dissipative cavity in RWA: where ω r = ω 0 − ω is the detuning between the cavity frequency and the driving field F at frequency ω.
We consider a trial density matrix given by a pure coherent state parameterized by complex α = α x + iα y (α x and α y are the real and imaginary parts of α): where â|α⟩ = α|α⟩.The functional Eq. ( 31) can then be evaluated as: The first bracketed term comes from the Hamiltonian dynamics while the other terms include damping effects.We can check that the minimum S 0 = 0 is achieved at: that is the classical response function of an oscillator at frequency ω 0 excited at a frequency ω with a detuning ω r = ω 0 − ω and a relaxation time τ = 2γ −1 .Since S 0 = 0 this solution is exact.
We now generalize this approach to the case of a cavity coupled to one and then more qubits.In the case of one-qubit with RWA Hamiltonian Eq. ( 3) we generalize the trial form of the density matrix to a product state: which is now parameterized by three real numbers b x,y,z given the Bloch sphere coordinates of the spin density matrix, in addition to the complex number α parametrizing the cavity coherent state.Evaluating Eq. ( 31) we find: where the first term corresponds to the single cavity described previously.In the second term the frequency Ω r = Ω−ω gives the detuning between the qubit level spacing Ω and the driving frequency.Here we omit the qubit index Ω = Ω l and λ is the qubit cavity interaction strength.This term just describes the qubit eigenstates without qubit-cavity coupling corresponding to b x,y = 0.The third term S λ describes the Hamiltonian part of the qubit-cavity coupling and all terms are proportional to λ.
The final term S γs contains all terms which are induced by the qubit damping: where the first term is minimized by the qubit ground state b x , b y = 0 and b z = −1 while the second term mixes qubit-cavity coupling and qubit relaxation.
The cavity-qubit steady-state is then obtained minimizing S as function of five parameters α x,y , b x,y,z .For the case of two-qubits we generalize the trial form of the density matrix to: the expression for S for this trial function is given in the Appendix.
Our trial form for the density matrix does not include any entanglement between qubits and cavity since the trial density matrix is taken as a tensor product.The cavity steady-state is given by a coherent state and the qubits by their Bloch sphere components.We thus call this the semiclassical approximation to the cavity-qubit steady-state.While minimization of the semiclassical function cannot no longer be done analytically for the general case, this minimization is very light computationally compared to full quantum calculations.It is thus important to identify the regimes of the Lindblad equation where the above semi-analytical approaches are able to correctly reproduce its full quantum dissipative solution.

III. NUMERICAL AND SEMI-ANALYTICAL RESULTS
In the previous Section, we described how the time dependent dissipative quantum problem of a driven cavity coupled to qubits can be transformed into a stationary problem by moving to the rotating frame.Finding the steady-state density matrix of the system is then reduced to find the kernel of the RWA Lindblad operator L(ρ).We then presented several approaches to finding this steady-state.An exact numerical approach based on an iteration of Sylvester equation solutions, a weak damping perturbation theory series which can be interpreted as population dynamics of the quantum eigenstates of the system and a semicalssical approximation which we expect to hold at stronger damping here quantum coherence are quickly destroyed.
Here, we first present some numerical results confirming the convergence of the system to its RWA steady-state and then investigate the different regimes of this model where our two semi-analytical approximations can apply.We find that for most parameters either one of these approximations can be used but we also identify a regime where they are not able to describe the entangled dissipative state of the system.

A. RWA validity analysis
To probe the validity of the RWA steady-state we consider an example of Lindblad dynamics for a single non dissipative qubit coupled to a dissipative cavity.This model has been analyzed previously in 20 using the method of quantum trajectories.This method showed a bistability of the quantum trajectory wave-function corresponding to two qubit orientations even for relatively large detuning between the qubit and cavity eigenfrequencies compared to their coupling strength λ.This nontrivial regime seems to be a good test-case to investigate RWA validity.
The bistability discussed in 20 manifests as two branches of the mean cavity occupation ⟨a + a⟩ for the two possible qubit eigenstates.This can be seen as two distinct peaks in the probability distribution of the cavity occupation number for spin down P − (n) = ⟨n, −|ρ|n, −⟩ and spin up respectively P + (n) = ⟨n, +|ρ|n, +⟩, with n the cavity eigenstate index.This bistability is reproduced by integration of time dependent Lindblad evolution in Fig. 1 in the top two panels showing P − (n) and P + (n) as function of the detuning ω r = ω 0 − ω between the cavity and the driving field.Other parameters are fixed to: ∆ = 2λ, where we introduce ∆ = Ω − ω 0 as the detuning between qubit and the cavity frequencies, F = λ and γ = 0.3λ.The cavity frequency is fixed to ω 0 /λ = 10.For convenience here, we consider units where ℏ = 1 and use the cavity qubit coupling strength λ as the energy scale in dimensionless quantities.This choice of units is convenient to analyze RWA validity because only relative energies to the excitation frequency appear in Eq. ( 3) and thus energy differences rather than the absolute energy become physically relevant.Bottom panels in Fig. 1 show the same quantities for the RWA steady-state with good agreement between both methods.
While RWA qualitatively reproduces the bistability behavior, some deviations are also visible.To be more quantitative on the agreement between RWA and exact dynamics we consider how the mean qubit/cavity quantities dependence on ω r = ω 0 − ω are changed with increasing RWA parameter ω 0 /λ which is expected to determine the RWA validity.
Results are shown in Fig. 2 for similar parameters as Fig. 1.While some deviations between integration of time dependent Eq. ( 1) and RWA are visible for ω/λ = 10, the deviation rapidly drops as ω/λ is increased.For ⟨ Ŝx ⟩, ⟨ Ŝz ⟩ and ⟨â + â⟩ the convergence is rapid with identical results for ω/∆ = 10 already, while the qubit projection Ŝy converges more slowly.In the numerical simulations we use 100 cavity levels which was found sufficient being twice larger than the maximum excitation number in Fig. 1.Slow relaxation of one of the spin components in this model is shown in the Appendix Fig. 13.RWA directly computes the steady-state of the system and thus does not require to integrate many oscillation periods before reaching steady-state.The deviations between the RWA steady-state and fully relaxed time integration are small in limit ω/λ ≫ 1 and thus the RWA can be a more reliable way to explore steady-state properties of this regime.Part of the deviations between RWA and full time dependent dynamics can be attributed to the slow relaxation of the spin degrees of freedom in this model, and thus integration over long times are needed to reach the system steady-state.Thus, overall we find that the RWA reproduces the bistability effect quite accurately.In the framework of quantum trajectories description, this bistability becomes visible as long time jumps between quantum states, these jumps on average should reproduce the populations of each bistable state.

IV. WEAK DAMPING RATE EQUATION APPROACH
We then present a typical example of comparison between the RWA steady-state and the summation of the rate equation series.We show the comparison for the case of a single qubit interacting with a cavity in the case where both spin and cavity are dissipative.Fig. 3 shows the mean qubit spin-projection ⟨S x ⟩ in the steady-state as function of the detuning ω r .Similar results are obtained for other spin components ⟨S y,z ⟩.The RWA steady-state found using a Sylvester equation iterations shows series of oscillations when ω r is in the interval [0, ∆].These oscillations correspond to multiphoton resonances in the rotating frame kω r = ∆ (where k is a positive integer).In the laboratory frame these resonances correspond to transitions (k + 1)ω = kω 0 + Ω where k + 1 photons are absorbed to excite k cavity levels and the qubit at energy Ω.The direct rate equation series (R) converge only far from the resonance and fail to reproduce the multiphoton transitions but the corrected series (R) * reproduce this behavior correctly and fail to converge only in a close vicinity to the cavity-driving field resonance when ω r ∼ γ.
The comparison in Fig. 3 is performed at relatively high frictions from the point of view of the rate equation perturbation theory and we see that the multi-photon resonances make the direct rate equation series (R) unstable in a wide range of ω r highlighting the requirement γ ≪ min i,j |ϵ i − ϵ j |.For weaker damping the rate equation series converge in a much wider range of detuning ω r as shown on Fig. 4. In this case even (R) converges almost everywhere except at very narrow resonances where some deviation from the exact steady-state is still visible.

V. STRONG DAMPING SEMICLASSICAL REGIME
We find that in the opposite limit of strong damping the semiclassical approximation is quite successful.To illustrate this we consider a cavity coupled to two qubits with all dissipation terms in Eqs.(4,5) both for cavity and qubits with equal dissipations γ = γ s .In this case the dissipation is chosen to be of the same order of magnitude as the detuning between qubits and cavity.This is thus a regime where we expect the semicalssical theory to give a very good description of the steady-state.This is confirmed by numerical simulations presented in Fig. 5 which shows that both  1).Different traces correspond to quantum dynamics for increasing RWA parameter ω0/λ = 10, 30, 50, 100.The agreement with RWA improves as ω0/λ increases but worsens a bit for the largest value.This behavior is explained as a transient effect in Appendix Fig. 13, where we show that relaxation to steady-state was not completed even after τp = 2 × 10 4 microwave periods.Indeed since the qubit is not dissipative in this simulation (γs = 0) relaxation time scale can be much longer than γ −1 .
the mean spin orientation and the mean cavity susceptibility variables α x = Re⟨â⟩ and α y = Im⟨â⟩ are reproduced almost perfectly by the minimization of the semiclassical functional S. We see that in addition to the main resonance at ω r = 0 the cavity susceptibility shows also features at resonance with the qubits ω = Ω 1,2 in agreement with exact RWA solution.For these strong values of damping the rate equations series are no longer converging for the range of detunings explored here.At weaker frictions the semiclassical results become less accurate but can reproduce still qualitative trends of the qubit response to cavity excitation.To illustrate this we show the semiclassical results for the case from Fig. 6 for one qubit.We see that although this method fails to capture the multiphoton resonance the broad shape of the response is reproduced correctly.

VI. BEYOND SEMICLASSICAL ANALYSIS AND RATE EQUATIONS
We found that the rate equation and semiclassical approaches capture complementary regimes of the dissipative quantum dynamics with respectively low and high damping regimes.It is thus interesting if we can find a regime not accessible to these approaches where exact solution of dissipative quantum dynamics is required.This regime would have to be at a weak damping near resonance drive since this is the only regime which is outside the range of validity of both approaches.However, a strong cavity excitation at a resonance can lead to an effectively semiclassical overheated regime, without special quantum coherence properties.
We thus consider the following setting where two qubits are coupled to the cavity with opposite detunings from the cavity frequency ∆ 1 = −∆ 2 (where ∆ 1 = Ω 1 − ω 0 and ∆ 2 − (Ω 2 − ω 0 ) are the qubit-cavity detunings).The cavity excitation at frequency ω = ω 0 does not break the symmetry between the qubits in RWA, and this seems a good regime to induce entanglements between qubits even in presence of damping.In the following we analyze this model in detail with RWA approach (6) and exact time dependent Lindblad dynamics (4), (5).We also consider the validity The trace (R) corresponds to summation of the series from the direct rate expansion Eq. ( 19) while (R) * exhibiting a larger radius of convergence corresponds to Eq. ( 30).The series (R) qualitatively reproduce the position of the multiphoton resonances but with excessive amplitude and fails to converge.The series (R) * reproduce multiphoton resonace accurately but still fail to converge close to the cavity resonance (ω − ω0)/λ ∼ 1 (further studies are needed to know if divergence occurs on energy scale λ or γ around the resonance). -0.4 FIG. 4: Comparison between RWA and rate equation series for weak damping γ = γS = 0.005λ.As expected in this regime, the rate equation series (R), from Eq. ( 19), converges.Excitation was reduced to avoid overheating at resonance with F = 0.1λ and 0.2λ.As in the previous figures Ω − ω0 = 2λ.
of our two approximate approaches in this regime.
To monitor the entangled state of the two detuned two qubits we consider their mean total spin operator Ŝ2 whose expectation value S(S + 1) = 2 for the triplet S = 1 state of the two qubits and zero for the singlet configuration.The dependence of ⟨ Ŝ2 ⟩ on ω − ω 0 is shown on Fig. 7 for two cases of qubit detuning from the cavity, comparing the antisymmetric detuning introduced above and and more generic values with (Ω 1 − ω 0 )/λ = 2, (Ω 2 − ω 0 )/λ = 3.For antisymmetric detuning a strong reduction of ⟨ Ŝ2 ⟩ is observed when the cavity is excited at resonance, this corresponds to a dynamics of the two qubits where they rotate together, synchronizing in opposite directions and canceling (at For the lowest dissipation γs = 0, γ = 0.3λ the semiclassical theory predicts a reversal of ⟨Sx⟩ in the range (ω − ω0)/λ ∈ (0, 2) which is not present in RWA which gives the exact quantum result but agreement is good outside this range.Adding some dissipation to the spin γs = γ = 0.3λ is enough to observe the polarization reversal in RWA also.Agreement becomes very good for γs = γ = 2λ.As expected at a higher friction the system behaves in a more semiclassical way.
least partially) their total spin.This reduction is not observed for the generic detuning values for which ⟨ Ŝ2 ⟩ stays close to 2, which is its equilibrium value.Fig. 7, corresponds to equal damping rates for cavity and qubits which are below but comparable to the coupling strength γ/λ = 0.2.The entanglement between the two qubits in the asymmetric detuning case can be enhanced when damping rate of the qubits are reduced.Such a situation is shown on Fig. 8 for γ/λ = 0.3 and γ s /λ = 3 × 10 −3 .The total spin ⟨ Ŝ2 ⟩ then reduces to ≃ 0.5 at resonance indicating a higher degree of anti-locking and compensation between the two qubits.
To confirm that this quantum synchronization and entanglement of the two qubits is not an artifact of the RWA we performed direct simulations of the time dependent Lindblad dynamics of Eqs.(4), (5).Fig. 8 presents a comparison between RWA and dynamical equations (4), (5, showing almost perfect agreement for ω/∆ 1 = 30.Interestingly, even if larger deviations from RWA appear when ω/∆ 1 = 30 is lowered, the minimum ⟨ Ŝ2 ⟩ remains constant indicating that the anti-locking of the qubits is robust to non-RWA effects which break the (anti)symmetry between the qubits 7: RWA calculation of the total spin of the qubit pair as function of (ω −ω0)/λ for two values of the qubit-cavity detuning ∆1,2 = Ω1,2 − ω0.The top curve corresponds to ∆1/λ = 2 and ∆2/λ = 3 with only a weak deviation from the equilibrium total spin triplet ⟨S 2 ⟩ = S(S +1) = 2.When the two Zeeman splittings are antisymmetric with respect to the cavity ∆1 = −∆2 = 2λ a significantly stronger reduction of ⟨S 2 ⟩ is observed.Dissipative rates are set to γ = γS = 0.2λ and excitation is F = 1.5λ.As previously 100 oscillator levels are used in the simulation.
shifting the frequency at which minimum ⟨ Ŝ2 ⟩ is achieved from exact resonance.Even if the degree of cancellation between the qubits spins is not perfect it is sufficient to generate steady-states violating Bell inequalities 5 , this is shown in Fig. 9.In the calculation of the Bell inequality the spin projection has to be measured, in two directions by two observers giving the freedom to choose four possible spin projection measurement directions for the computation of the correlators.Since the entangled steady-state is not a pure singlet state we found that violation of the Bell inequality was maximized with a ∼ 20 • rotation of the projection directions compared to the rotation of 45 • which is used for Bell inequality tests on pure Bell (EPR) states.To avoid any dependence on the projection angles, we also show the quantum negativity (see e.g. 42) of the reduced qubit pair density matrix (obtained from the total density matrix by tracing over the cavity).The negativity at resonance in Fig. 9 is 0.36, that is not far from the maximal possible negativity of 0.5 for a two qubit pair.We emphasize that as opposed to entangled states generated by pulse sequences which have finite lifetimes, this stationary entangled state is preserved in time in presence of dissipation and decoherence within our model.The preservation of entangled state of two qubits was also seen in numerical studies with quantum trajectories description but there the dissipation was present only for the cavity and not for qubits 21 .
To conclude this part we investigate if this anti-locked entangled state of the qubits could be reproduced in the terms of the approximate rate equation or semiclassical method.The rate equation series for ∆ 1 = −∆ 2 fails to converge in rather wide range ω r /∆ 1 ∈ (−1, 1) probably due to resonant energy levels which are detrimental for convergence.Treating rate equations as an asymptotic series and summing only the first terms which improve convergence allows to qualitatively reproduce ⟨S x ⟩ in all the range ω r .However the prediction for the total spin ⟨ Ŝ2 ⟩ is then misleading at resonance suggesting a maximum ⟨ Ŝ2 ⟩ in contradiction with exact results, as it is illustrated in Fig. 10.Not surprisingly the semiclassical approximation only succeeds in reproducing qualitative features but completely misses the singlet formation.

VII. SYNCHRONIZATION OF SEVERAL QUBITS
It was shown in 20 , using the method of quantum trajectories, that the phase of rotating qubits can be synchronized with the phase of microwave driving.We confirmed this by integration of the quantum Lindblad evolution Eqs.(4), (5)  (see Appendix Fig. 15).We study this synchronization effect for up to four qubits.Our results are shown on Figs.11,12 where the RWA steady-state is compared with semiclassical (Fig. 11) and rate equation theories (Fig. 12).All the qubits in this simulation are detuned from each other and without external driving their in-plane magnetization will process at different frequencies leading to an average cancellation of the total in plane magnetization.When the cavity is driven two synchronization peaks appear where the total spin projection in the rotating frame ⟨S Here, the driving strength is set to F/λ = 0.25, the dissipative rates are γ = 0.3λ with a weak qubit dissipation γS = 10 −3 γ (for these parameters, at resonance, ⟨n⟩ ≃ 4F 2 /γ 2 ≃ 3).To confirm that the singlet formation is not an artifact of the RWA (black curve), we performed direct integration of the time dependent Lindblad dynamics up to total simulation time 3γ −1 s for increasing RWA parameter ω0/λ (color curves with symbols).The singlet formation is robust to non RWA effects with the minimum ⟨S 2 ⟩ remaining unchanged as ω0/λ is varied by an order of magnitude.Only weak Non RWA effects are visible as a small shift of the minimum from ω = ω0 and an asymmetric ⟨S 2 ⟩ dependence, since non RWA effects break the symmetry between two anti-symmetrically detuned qubits.Maximal negativity for two qubits is 1/2 42 and thus this steady-state shows a high degree of stationary entanglement despite dissipative decoherence of both qubits and cavity.
grows linearly with the number of detuned qubits.The first peak is when the cavity is excited near resonance, the resonant cavity vibrations excite the qubits which then all presses at the same phase and frequency.Surprisingly, a second synchronization peak appears also at a higher frequency detuned from both cavity and qubit resonances.The semiclassical approximation reproduces correctly the two synchronization peaks however the agreement with RWA results is only qualitative.This is because this is a weak dissipation regime with strong cavity qubit coupling ∆ 1 = λ and γ/λ = 0.2.In Appendix Fig. 16 we show that the agreement becomes almost perfect for the simpler weak coupling strong dissipation limit.We note that in this simple regime the linear growth of ⟨S x ⟩ with the number of qubits is observed only at cavity resonance.Fig. 11 suggests that the agreement between RWA and semiclassical approximation tends to improve with the number of qubits, this maybe due to the fact that our semiclassical approximation can also be viewed as a mean-field theory whose accuracy improves with more interacting qubits.The summation of the rate equation series reproduces the RWA result exactly away from cavity resonance on Fig. 12, however it seems that the diverging region around resonance grows slowly with the number of qubits.
To summarize we find that this regime of quantum synchronization, where the total rotating in plane spin grows with the number of qubits, is well described by both semiclassical and rate equation approaches.As discussed in the previous Section, it is also possible to synchronize entangled qubits for antisymmetric qubit-cavity detunings in a vicinity of resonance between microwave driving and cavity frequency.

VIII. DISCUSSION
In this work we used the Lindblad equation density matrix formalism for a description of dissipative monochromatically driven resonator (oscillator cavity) interacting with one or several qubits with very long or moderate dissipative life time.This system extends the seminal Jaynes-Cummings model 33 to a dissipative regime which is typical for superconducting qubits coupled to a driven cavity 6 .The performed numerical simulations show that the results for the full time-dependent Lindblad equation ( 4), (5) are well reproduced by the solution of stationary RWA Lindblad equation (6).Efficient exact methods are developed to find the RWA steady-stae numerically with several qubits and driven frequency being close to the cavity resonance when many cavity eigenstates are excited.For one qubit case a regime of qubit bistability is established confirming previous results obtained with the method of quantum trajectories 20 .
We also developed and tested two semi-analytical approaches which allow to obtain the approximate accurate solutions of exact RWA steady-state density matrix in the regimes of relatively weak and strong dissipation corresponding respectively to the rate equation and semiclassical approximations for the steady-state density matrix.The analytical steps of these two approaches significantly simplify the equations for density matrix with consecutive numerical solution and allow to investigate the system behavior in a vicinity of cavity resonance with many excited states.The numerical verification of these semi-analytical approaches is done by the comparison with the exact RWA solution and confirms their validity for a broad range of system parameters.At the same time we demonstrate the existence of system behavior which cannot be described by these semi-analytical approaches.Thus we find that for two dissipative qubits and dissipative driven cavity there exist a regime when qubits remain entangled, forming a singlet, in the steady-state density matrix.The existence of such a stationary regime of entangled synchronized qubits, created FIG.11: RWA calculation of total spin projection ⟨Sx⟩ of an increasing number of qubits coupled to one cavity (full curves) as function of cavity-excitation detuning and comparison to the semiclassical theory (dashed curves).The qubit cavity detunings are set to ∆1 = λ, ∆2 = 1.5λ, ∆3 = 2λ, ∆4 = 2.5λ (only the first qubits are kept when the number of qubits is smaller than four).Dissipative rates are γ = γS = 0.2λ and excitation is F = 0.77λ.Even if relaxation rates are all small, the semiclassical theory still captures many properties of ⟨Sx⟩, reproducing the increase of ⟨Sx⟩ with the number of qubits which corresponds to a synchronization of qubit rotation by the external drive.In both RWA and semiclassical data synchronization occurs at resonance ω = ω0 but perhaps less expected at a higher frequency detuned from both cavity and qubits ((ω − ω0)/λ increasing from 2 to 3 with the number of qubits).Interestingly the accuracy of the semicalssical approximation seems to improve with the number of qubits, that may be due to its mean-field character.by a monochromatically driven cavity, is really surprising since both qubits and cavity are dissipative.Due to the preservation of entanglement of qubits such a regime can be considered as a real entangled quantum synchronization.We also find regimes (see Fig. 11) when up to 4 qubits are synchronized with a phase of monochromatic cavity driving so that the total spin of the system grows proportionally to the number of system qubits (spin halves).At the same time we show that this synchronized regime is well described by the semiclassical approach so that one can discuss if such synchronized regime of several qubits can be considered as purely quantum or more as a regime of semiclassical synchronization in presence of strong dissipation and noise induced by quantum fluctuations.Indeed, it is known that the classical synchronization is preserved in presence of moderate noise 2 .Even if this synchronization of several qubits can be described in the frame of semiclassical approach and that there is no entanglement one can also argue that spin halves are purely quantum two-level systems and hence their synchronization is also quantum.The regime of entangled quantum synchronization of qubits is purely quantum and cannot be obtained in the frame of semiclassical description.
The obtained results provide a better understanding of nontrivial behavior regimes of several dissipative qubits interacting with dissipative driven cavity.We hope that the developed methods can be useful in other contexts.where ω r = ω 0 − ω and Ω r1 = Ω 1 − ω, Ω r2 = Ω 2 − ω.

FIG. 1 :
FIG.1: Distribution of the oscillator occupation number n as function of the detuning from cavity resonance (ω − ω0)/λ and z-axis spin projection.The color scale shows the amplitude |ρ(n, s; n, s)| where s =↑, ↓ is the spin up/down state.We remind that λ is the strength of the spin cavity interaction, we use it as the energy scale to facilate comparison between direct time integration and RWA which depends only on the difference ∆ = Ω−ω0 between the Zeeman splitting Ω and the cavity frequency ω0, which is set to ∆ = 2λ in this figure.Here, the driving strength is set to F = λ, cavity dissipation is γ = 0.3λ and qubit is dissipationless γs = 0. Top row figures are obtained by direct time integration for ω0/λ = 10, while bottom row are results from RWA.For the Lindblad time integration, numerical data shows the density matrix after τp = 2 × 10 4 excitation periods starting from the system ground-state for F = 0 at t = 0 (we use the same τp for results presentation of other cases of time dependent Lindblad equation).The oscillator phase space in the simulations was truncated to the first 100 oscillator levels (usually used for other cases).

FIG. 2 :
FIG.2: Mean spin projections ⟨Sx,y,z⟩ and oscillator quantum number ⟨n⟩ = ⟨a + a⟩ as function of the detuning (ω − ω0)/λ for Ω − ω0 = 2λ, F = λ, γ = 0.03λ (same values as in Fig.1).Different traces correspond to quantum dynamics for increasing RWA parameter ω0/λ = 10, 30, 50, 100.The agreement with RWA improves as ω0/λ increases but worsens a bit for the largest value.This behavior is explained as a transient effect in Appendix Fig.13, where we show that relaxation to steady-state was not completed even after τp = 2 × 10 4 microwave periods.Indeed since the qubit is not dissipative in this simulation (γs = 0) relaxation time scale can be much longer than γ −1 .

FIG. 3 :
FIG.3: Comparison of RWA simulation with the summation of rate equation series for F = λ, Ω − ω0 = 2λ and γ = γS = 0.3λ.The trace (R) corresponds to summation of the series from the direct rate expansion Eq. (19) while (R) * exhibiting a larger radius of convergence corresponds to Eq. (30).The series (R) qualitatively reproduce the position of the multiphoton resonances but with excessive amplitude and fails to converge.The series (R) * reproduce multiphoton resonace accurately but still fail to converge close to the cavity resonance (ω − ω0)/λ ∼ 1 (further studies are needed to know if divergence occurs on energy scale λ or γ around the resonance).

FIG. 8 :
FIG.8: When dissipative rates are reduced compared to Fig.7, the reduction of ⟨S 2 ⟩ for antisymmetric qubit-cavity detunings ∆1 = −∆2 = λ becomes stronger and the singlet state of the qubit pair becomes the most probable state (75% singlet probability at the minimum of ⟨S 2 ⟩).Here, the driving strength is set to F/λ = 0.25, the dissipative rates are γ = 0.3λ with a weak qubit dissipation γS = 10 −3 γ (for these parameters, at resonance, ⟨n⟩ ≃ 4F 2 /γ 2 ≃ 3).To confirm that the singlet formation is not an artifact of the RWA (black curve), we performed direct integration of the time dependent Lindblad dynamics up to total simulation time 3γ−1

FIG. 9 :
FIG.9: Bell inequality violation and negativity of the steady-state RWA qubit pair with the reduced density matrix (trace done over the cavity) for the parameters of Fig.8.Since the qubit pair is in a mixture of singlet and triplet states, the polarization choice for the Bell inequality has to be adjusted to observe a Bell inequality violation (see Appendix and also Fig.14 there).Maximal negativity for two qubits is 1/2 42 and thus this steady-state shows a high degree of stationary entanglement despite dissipative decoherence of both qubits and cavity.

FIG. 10 :
FIG.10:We test here if our two semi-analytic approaches can reproduce singlet formation for antisymmetric detuning presented in Figs.8,9.The left panel shows the RWA spin projection ⟨S1x⟩ for the first qubit compared with the rate equation series and the variational approximation, the right panel shows the mean total spin ⟨S 2 ⟩.While both approaches reproduce some qualitative features, they both fail to describe the singlet formation at ω − ω0 = 0.

FIG. 12 :
FIG.12:The RWA calculation of total spin projection ⟨Sx⟩ from Fig.11is compared to the result of the summation of the rate equation series (R) * (doted curves).The rate equation series give a result being almost exact away from resonance but they suffer from instabilities near ω = ω0.It seems that the convergence range decreases with the number of qubits with various multi-qubit resonances making the series unstable even at ω − ω0 = 1.5λ for four qubits.

FIG. 15 :
FIG.15: Synchronization between a driven cavity and a qubit in the system steady-state for ∆1 = Ω1 − ω0 = λ, F = λ, γ = λ/3 and RWA parameter ω0/λ = 10.The qubit is non dissipative γs = 0 and data is obtained by integration of Lindblad dynamics.Due to the moderate value of the RWA parameter vibrations around the mean RWA values of the spin projections are clearly visible on the left panel.Right hand panel shows the synchronlisation between the angle θ = arg⟨Sx + iSy⟩ of the qubit in plane spin projection and the phase of the cavity driving field.