# Entangled Frequency-Tunable Microwave Photons in a Superconducting Circuit
## 4 Numerical Results  
### Creation of the Photon Pairs  
The system is initialized in the two-mode vacuum state \(|\psi_{\psi}\rangle=|0,0\rangle\) with \(\langle\hat{a}^{\dagger}\hat{a}\rangle=\langle\hat{b}^{\dagger}\hat{b}\rangle\equiv 0\). Under the excitation of the pump field, the photons of the two modes are emitted out of the circuit; i.e., the input pump photons are split into two different modes when they pass through the superconducting circuit. The state of an emitter is \(|\psi\rangle\equiv|\psi_{a}\rangle\otimes|\psi_{b}\rangle\), where \(|\psi_{a}\rangle\) and \(|\psi_{b}\rangle\) are the cavity Fock states with twenty dimensions \(\mathcal{R}_{a}\) and \(\mathcal{R}_{b}\), respectively. \(\mathcal{R}_{a}\) and \(\mathcal{R}_{b}\) build the truncated Hilbert space \(\mathcal{G}=\mathcal{R}_{a}\otimes\mathcal{R}_{b}\) (note that the increasing dimension of the truncated Hilbert space \(\mathcal{G}\) does not significantly affect the numerical steady-state results of our system [46; 47; 48; 49; 50; 51]). The output photon modes \(a\) and \(b\) can be described by the corresponding operator \(\hat{p}(t)\) and governed by the master equation, which describes the evolution of the quantum system as  
\[\frac{d\hat{p}}{dt}=-\frac{i}{\hbar}[\hat{H},\hat{p}]+\kappa_{a}\mathcal{D}[ \hat{a}]\hat{p}+\kappa_{b}\mathcal{D}[\hat{b}]\hat{p}, \tag{10}\]  
Figure 3: Dependence of the Bessel function of the first kind of order 0 (red) and 1 (blue) on the dimensionless parameter \(\alpha_{p}\) of the pump field, which determines the strength of \(\chi_{0,1}\), where \(\alpha_{p}=2.404\) is chosen to eliminate the swap interaction between the two modes. where \(\mathcal{D}[\hat{e}]\beta=\hat{e}\hat{\rho}\hat{e}^{\dagger}-\hat{e}^{\dagger}\hat{e }\hat{\rho}/2-\hat{\rho}\hat{e}^{\dagger}\hat{e}/2\), \((\hat{e}=\hat{a},\hat{b})\). \(\mathcal{D}[\hat{e}]\) is the dissipative superoperator and \(\kappa_{a}\) (\(\kappa_{b}\)) is the decay rate of the mode \(a\) (\(b\)). Under the rotating wave approximation with the pump excitation, the master equation of our circuit with effective Hamiltonian \(\hat{H}_{I}\) is  
\[\frac{d\hat{\beta}}{dt}=-\frac{i}{\hbar}[\hat{H}_{I},\hat{\rho}]+\kappa_{a} \mathcal{D}[\hat{a}]\hat{\beta}+\kappa_{b}\mathcal{D}[\hat{b}]\hat{\beta}, \tag{11}\]  
where the \(\hat{\beta}\) is the density matrix of the source under the rotating framework. To demonstrate that the proposed circuit pumped by an external field will generate photons, we use QuTip [52; 53] to numerically calculate the master equation of Equation (11) combined with Equation (9) at different working points \(\Phi_{\text{dc}}\) corresponding to different coupling strengths \(g\). As shown in Figure 4a, the mean microwave photons of both modes flee from the system with the same tendency, which means the two modes are generated simultaneously and have an identical dependence on the coupling strength \(g\), which is a function of \(\Phi_{\text{dc}}\). \(\Phi_{\text{dc}}=0\) with a coupling strength \(g/2\pi=22.889\) MHz is selected as an example to demonstrate the detailed behavior in Figure 4b. In the steady state of the two modes, the mean photon number is \(\langle\hat{a}^{\dagger}\hat{a}\rangle\approx 1.59\) and \(\langle\hat{b}^{\dagger}\hat{b}\rangle\approx 0.80\), respectively. Note that we carry out the numerical calculations with the whole Hamiltonian of our system and compare the results with these calculated with the interaction Hamiltonian. No obvious difference is observed. One can measure the bunching effect of the generated two different modes by observing \(\hat{D}_{j}(t)\), \(\hat{D}_{j}(t+\tau)\) (\(j=a,b\)), which represent the intensity of modes with a delay time \(\tau\) (upper panel in Figure 1b) and record the photon statistics in the time domain. To further demonstrate that the photons are emitted in pairs simultaneously, we reconstruct the second-order coherence function \(g_{ab}^{(2)}(t,\tau)\) of the output microwave signals at different dc working points using numerical simulation.