# Quantum Synchronization and Entanglement of Dissipative Qubits Coupled to a Resonator
## 2 Model Description
### Rate Equation Perturbation Theory  
The above approach for steady-state computation can 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 the density matrix; it also allows us to investigate if the steady state of model Equation (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 time scales between the Hamiltonian dynamics and the slower dissipative processes. The crossover from quantum master equations to incoherent rate dynamics is a fundamental problem that has been studied previously in different contexts [43; 44; 45; 46]. In the following treatment, we adapt and optimize the rate equation approach for the quantum synchronization of a driven cavity coupled to qubits, thereby increasing the convergence radius of the dissipative series. We describe this approach by starting from the general steady-state Lindblad equation and introducing a formal expansion parameter \(\epsilon\) to keep track of the order of the dissipative parameters in the expansion:  
\[\mathcal{L}(\hat{\rho}) = -i[\hat{H},\hat{\rho}]+\epsilon\sum_{k}\biggl{(}\hat{L}_{k}\hat{ \rho}\hat{L}_{k}^{+}-\frac{1}{2}\hat{L}_{k}^{+}\hat{L}_{k}\hat{\rho}-\frac{1} {2}\hat{\rho}\hat{L}_{k}^{+}\hat{L}_{k}\biggr{)} \tag{15}\] \[= \mathcal{L}_{\mathcal{H}}(\hat{\rho})+\epsilon\mathcal{L}_{ \gamma}(\hat{\rho})\] (16) \[\mathcal{L}_{\mathcal{H}}(\hat{\rho}) = -i[\hat{H},\hat{\rho}]\] (17) \[\mathcal{L}_{\gamma}(\hat{\rho}) = \sum_{k}\biggl{(}\hat{L}_{k}\hat{\rho}\hat{L}_{k}^{+}-\frac{1}{2 }\hat{L}_{k}^{+}\hat{L}_{k}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{L}_{k}^{+} \hat{L}_{k}\biggr{)} \tag{18}\]  
here, \(\hat{L}_{k}\)--as previously--denotes the relaxation operators, and \(\hat{H}\) is a steady-state Hamiltonian, which in our case is given by RWA Equation (3). To find the steady state, we need to solve \(\mathcal{L}(\hat{\rho})=0\), and we proceed to make a formal expansion in \(\epsilon\):  
\[(\mathcal{L}_{\mathcal{H}}+\epsilon\mathcal{L}_{\gamma})(\hat{\rho}_{0}+\epsilon \hat{\rho}_{1}+...+\epsilon^{i}\hat{\rho}_{j}+...)=0 \tag{19}\]  
The equations for the lowest-order approximation \(\hat{\rho}_{0}\) are as follows:  
\[\mathcal{L}_{\mathcal{H}}\hat{\rho}_{0}=0\;,\;\mathcal{P}_{\mathcal{D}} \mathcal{L}_{\gamma}\hat{\rho}_{0}=0\;,\;\text{Tr}\;\hat{\rho}_{0}=1 \tag{20}\]  
where the operator \(\mathcal{P}_{\mathcal{D}}\) acts on a density matrix defined in the eigenbasis of \(\hat{H}\) by keeping only its diagonal part. The first equation \(\mathcal{L}_{\mathcal{H}}\hat{\rho}_{0}=0\) imposes that \(\hat{\rho}_{0}\) is diagonal in the eigenbasis \(|n\rangle\) of the Hamiltonian \(\hat{H}=\sum_{n}e_{n}|n\rangle\langle n|\). Thus, \(\hat{\rho}_{0}\) has the form \(\hat{\rho}_{0}=\sum_{n}P_{n}|n\rangle\langle n|\) where \(P_{n}\) can be interpreted as the probability of eigenstate \(|n\rangle\). The last two equations are, thus, the rate equation determining \(P_{n}\) and the normalization condition \(\sum_{n}P_{n}=1\). The solution to this equation involves building the matrix, as follows:  
\[\mathcal{K}_{mn}=\langle m|\mathcal{L}_{\gamma}(|n\rangle\langle n|)|m\rangle. \tag{21}\]  
We then find non-zero solutions \(\mathcal{K}|\psi\rangle=0\), such solutions always exist since \(\langle 1,1,...,1|\mathcal{K}=0\).