# Quantum Synchronization and Entanglement of Dissipative Qubits Coupled to a Resonator
## 2 Model Description
### Rate Equation Perturbation Theory
The recurrence equations for the next orders are as follows:  
\[\mathcal{L}_{\mathcal{H}}(1-\mathcal{P}_{\mathcal{D}})\hat{\rho}_{j+1}+ \mathcal{L}_{\gamma}\hat{\rho}_{j}=0 \tag{22}\]  
which defines the off-diagonal components of \(\hat{\rho}_{j+1}\) (in the eigenbasis of \(\hat{H}\)) but leaves its diagonal part undefined. The diagonal part of \(\hat{\rho}_{n+1}\) can be found, requiring the following:  
\[\mathcal{P}_{\mathcal{D}}\mathcal{L}_{\gamma}\hat{\rho}_{j+1}=0 \tag{23}\]  
this can be viewed as higher-order corrections to the rate equations changing the occupation state probabilities, \(P_{n}\). Finally, the matrix \(\hat{\rho}_{j+1}\) is then made traceless by the transformation  
\[\hat{\rho}_{j+1}\rightarrow\hat{\rho}_{j+1}-\hat{\rho}_{0}\text{Tr}\;\hat{ \rho}_{j+1}\;, \tag{24}\]  
this transformation is allowed because the induced error is of a higher order:  
\[= \mathcal{L}_{\mathcal{H}}(\hat{\rho})+\epsilon\mathcal{L}_{ \gamma}(\hat{\rho})\]0  
This can also be viewed as an order-by-order normalization of the density matrix \(\text{Tr}\rho=1\). We denote as \((\mathbf{R})\) the numerical result of the summation of these series. The convergence of this expansion can be checked by evaluating the residual error in Equation (19) at each step of the summation, checking that the error decreases as the number of terms grows. 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 \((\mathbf{R})\) seems to have a rather small convergence radius requiring the amplitude of the dissipation to be small compared to energy level splittings:  
\[= \mathcal{L}_{\mathcal{H}}(\hat{\rho})+\epsilon\mathcal{L}_{ \gamma}(\hat{\rho})\]1  
In practice, some of these transitions are almost forbidden by selection rules involving several spin flips and transitions between several oscillator quanta. We thus attempt to improve the radius of convergence of the rate equation series by incorporating dissipative terms into \(\mathcal{L}_{\mathcal{H}}\). This is done by observing that \(\mathcal{L}_{\mathcal{H}}\) acts on the density matrix \(|n\rangle\langle m|\) as follows:  
\[= \mathcal{L}_{\mathcal{H}}(\hat{\rho})+\epsilon\mathcal{L}_{ \gamma}(\hat{\rho})\]2  
We can, thus, formally include in \(\mathcal{L}_{\mathcal{H}}\) all the dissipative terms, which are diagonal on the same basis. This is done by introducing a new superoperator \(\mathcal{L}_{\gamma\mathcal{K}}\), which in the eigenbasis of \(\mathcal{H}\) is defined as follows:  
\[\mathcal{L}_{\gamma\mathcal{K}}(\hat{\rho})_{nm}=\left\{ \begin{array}{cc}[ \langle n|\mathcal{L}_{\gamma}(|n\rangle\langle m|)|m\rangle]\hat{\rho}_{nm}&( m\
eq n)\\ 0&m=n\end{array} \right. \tag{28}\]  
The expansion \((\mathbf{R})\) is then performed with the updated operators \(\mathcal{L}_{\gamma}\) and \(\mathcal{L}_{\mathcal{H}}\) given by the following:  
\[\mathcal{L}_{\gamma} =\mathcal{L}_{\gamma}-\mathcal{L}_{\gamma\mathcal{K}} \tag{29}\] \[\mathcal{L}_{\mathcal{H}} =\mathcal{L}_{\mathcal{H}}+\mathcal{L}_{\gamma\mathcal{K}} \tag{30}\]  
This approach has some similarities with a Jacobi preconditioning step (see, e.g., [47]), 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, \((\mathbf{R})^{*}\), improves the radius of convergence of the perturbation series, and for low-friction divergence, seems to occur only close to the exact resonance \(\omega=\omega_{0}\), when the RWA Hamiltonian Equation (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.