Controllable quantum dynamics of inhomogeneous nitrogen-vacancy center ensembles coupled to superconducting resonators

We explore controllable quantum dynamics of a hybrid system, which consists of an array of mutually coupled superconducting resonators (SRs) with each containing a nitrogen-vacancy center spin ensemble (NVE) in the presence of inhomogeneous broadening. We focus on a three-site model, which compared with the two-site case, shows more complicated and richer dynamical behavior, and displays a series of damped oscillations under various experimental situations, reflecting the intricate balance and competition between the NVE-SR collective coupling and the adjacent-site photon hopping. Particularly, we find that the inhomogeneous broadening of the spin ensemble can suppress the population transfer between the SR and the local NVE. In this context, although the inhomogeneous broadening of the spin ensemble diminishes entanglement among the NVEs, optimal entanglement, characterized by averaging the lower bound of concurrence, could be achieved through accurately adjusting the tunable parameters.

that the multi-NVE dynamics itself is very complicated and could exhibit richer dynamical behavior than the two-NVE case by comparing a three-NVE model with a two-NVE one. Then we show how dephasing induced by the inhomogeneous broadening influences quantum dynamics and entanglement generation in a three-NVE case, where we concentrate on a specific frequency distribution of the inhomogeneous NVE. In our study, we find that the dynamics displays a series of oscillations under various experimental situations, which reflects the intricate competition and balance between the NVE-SR collective coupling and the adjacent-site photon hopping. Influenced by the dissipation of SR and NVE, the oscillations show a slight trend of damping. However, for the dephasing effect caused by the inhomogeneous broadening, a counter-intuitive effect is the suppression of the population transfer between the SR and the local NVE, namely, the dephasing effect enlarges the oscillation periods. In this context, the optimal entanglement among the NVEs could be achieved through accurately adjusting the tunable parameters, such as the Rabi frequency due to the external driving field as well as the hopping rates.
On the other hand, preparation of a high-degree entanglement among three or more NVEs is challenging both experimentally and theoretically. A promising way to overcome this obstacle might be offered by recent development of circuit QED technique in 1D (2D) superconducting resonator array [51][52][53][54] , providing more attractive possibilities to realize distributed entanglement with tunable interactions among multiple separated spin ensembles. Our detailed analysis finds a way to extract proper experimental parameters for the optimal entanglement among the NVEs using existing experimental technologies, even in the presence of large inhomogeneous broadening of the spin ensemble. As such, the present system provides a platform to generate multipartite quantum entanglement of the NVEs embedded in distant sites. In this sense such studies not only provide information of how entanglement evolves with time but also suggest ways toward practical purposes because the entanglement can be controlled by adjusting the tunable parameters 55,56 . Therefore, it is desirable to investigate quantum dynamics of the NVEs in different sites, and develop efficient methods for controlling the entanglement dynamics of several distant NVEs.

Results
System and model. As illustrated in Fig. 1, the system under consideration is a circuit QED array consisting of N distant NVEs coupled to N separate SRs respectively. The energy level configuration of the NV − center is shown in Fig. 1(c). Two electrons of the vacancy make a quasi covalent bond with the lone pair of nitrogen atom, and an extra electron located at the vacancy site forms a spin S = 1 pair with the residual electron of vacancy. The zero-field splitting D of the triplet ground state where N 0 is the number of the NV − centers in the spin ensemble, σ = + i 0 2 2 denotes the effective coupling strength between the i-th NV − center and SR, where g 0 is the coupling strength of a NV − center to the resonator, Ω is the Rabi frequency of the external driving field. As shown in Fig. 1(b), there exist four crystalline orientations for each NV − center 14 , and all the angles of the possible four orientations of the NV − center with respect to external field B stat could be identical if the magnetic field is applied along a special direction [100]. In this case, the Zeeman splitting between the states = ± m 1 with μ B the Bohr magneton, ≈ g 2 e the Lande factor for electron spin. δ i is the random offset from the central frequency δ B for the i-th NV − center, describing the inhomogeneous broadening of the spin ensemble induced by local strain and interactions with neighboring electronic spins or nuclear spins 58 . Engineering such a system requires a quantitative understanding of how the inhomogeneous broadening effects on quantum dynamics. Numerically, to fit the inhomogeneous distribution of the frequencies, we adopt the , where γ s is the full width of the frequency at half maximum and the dimensionless parameter 1 < q < 3 covers the Gaussian distribution for q → 1 and the Lorentzian distribution for q = 2 50 . Here the key idea in our numerical treatment is to divide the spin ensemble into several subensembles 43,45,59 and each one is assumed to be homogeneous, described by i are the collective spin operators. Respectively, δ μ and µ  g are the random offset and the effective collective coupling strength for the μ-th subensemble with Nμ NV − centers.
After mapping collective spin-N μ /2 system in Ĥ subs into effective bosonic operators using the Holstein-Primakoff (HP) transformation 60 the creation (annihilation) operators, one can obtain following Hamiltonian where all the nonlinear items containing µ µ μ † b b N / could be safely neglected.
For the whole system, when considering the nonlocal microwave photons hopping between the adjacent sites, we have the total Hamiltonian as where  is the hopping rate between neighboring resonators and l indicates the l-th site. We may assume an uniform distribution of coupling rates across the 1D array, i.e., G l,μ = G μ . The term + µ μˆˆ † † a b ab in the Hamiltonian Ĥ T denotes the linear mixing interaction. Actually, for the single-excitation case, the present model is equivalent to the Jaynes-Cummings-Hubbard model, which describes a single photonic mode strongly coupled to a two-level system in each cavity and photons can hop between the cavities. However, in contrast to single atom dynamics, the NVEs have more possibilities of increasing the collective coupling to a point at which counter-rotating-wave interactions become relevant, benefiting from the collective enhancement of the magnetic-dipole coupling. Additionally, the inhomogeneous broadening effect of spin ensemble, which is absent in the single-atom model, attracts even more attention as it may bring some special or novel effects, such as the results presented in the next section.
Quantum dynamics. As shown below, rich quantum behaviors can be observed using recent experimental technologies even in the presence of large inhomogeneous broadening of spin ensemble. Our system offers a great degree of freedom to control its evolution since the linear mixing process can be controlled by the effective coupling strengths G μ , which are dependent on the Rabi frequencies Ω of the external driving fields. Additionally, the in-situ tunability of the circuit elements also offers a channel to control its evolution by modulating the adjacent-site photon-hopping rates  . Here we define the collective coupling strength δ and investigate the single-excitation dynamics of the whole system in the large coupling case , and the large hopping case   G c , respectively. The equations describing the dynamics of each boson involved in the system can be written as t l l l l hom where κ is the photon loss rate of the SR, γ hom is the dissipative rate of the NVEs, and the dephasing effect is involved by considering the NVE's inhomogeneous broadening.
We first study the ideal case without considering any dissipation/dephasing. In general, multi-NVE dynamics would be more complicated and exhibit richer dynamical behaviors than the two-NVE case. To verify this conjecture, in Fig. 2, we compare the dynamics of a two-site model with a three-site model under the same condition. For example, the initial excitation is seeded in every NVE with equal possibility, i.e., for the two-site model, and for the three-site model, where L, M and R denote the left, middle and right site respectively. For the two-site model, the dynamics of the NVEs (SRs) in the left site and right site are fully overlapped, oscillating with the amplitude around 0.5 in (a1,b1) and 0.07 in (a2,b2). It implies that the full excitation transfers between the spin ensemble and the local resonator take place simultaneously in both sites due to symmetry. For the three-site model, the middle actor breaks the balance of excitation distribution in the two-site case, presenting an obviously different dynamics from the sideward sites. The total dynamics of the system shows small oscillation periods and large oscillation envelopes with different time-dependent amplitudes. These novel features are absent in the two-site model. Besides, there is no evident difference between the large coupling regime and the large hopping regime in the two-site model. In contrast, the dynamics of these two different coupling regimes can be easily distinguished in the three-site model, exposing more complicated and richer evolution behaviors. Now we focus on the dynamics of the three-site system in different coupling regimes as shown in Fig. 3. Initially, the excitation is assumed to be seeded in the left SR, i.e., = In the large coupling regime  G c  , as shown in Fig. 3(a1,b1), the curves representing the population evolution of the NVE in (b1) follow the oscillations of population in the corresponding SR in (a1), implying a rapid and periodic excitation transfer between the SR and NVE at the same site. However, the relatively slow photon-hopping process facilitates a long-period excitation transfer from the left SR to the right one in turn, and so is the reverse. Additionally, the frequency of fast oscillation related to the coupling process is several times as that of the hopping process. Based on the analytical results shown in Method, and using the Fourier transformation, typical oscillation frequencies appear during the dynamical process, as MHz in (a2,b2). Due to symmetry of the dynamics, the thin blue dotteddashed line and the thin red dashed line are overlapping in each panel.
The concrete process can be described in more details: (i) The excitation is initially located in the left SR mode, and then transfers to its corresponding NVE during a period ∼ 2π/G c due to a strong coupling process; (ii) At the same time, through weak hopping mechanism, the energy gradually flows into the middle site and the right site in sequence during a period π ∼2 / ; (iii) The JC-type oscillation also takes place at both the middle and the right sites; (iv) The reverse process occurs, which ultimately leads to excitation of the left SR and the corresponding NVE; (v) The above-mentioned process repeats for several times, as shown in Fig. 3(a1,b1).
In the competition case ≈ G c  , as plotted in Fig. 3(a2,b2), the populations of NVEs and SRs are found to exhibit characteristic oscillations with different frequencies, which is similar to that in the large coupling regime. In the case of = G c  , the typical oscillation frequencies involved during the quantum dynamics are extracted as µ . .
In the large hopping regime   G c , as illustrated in Fig. 3(a3,b3), unlike the behavior in above two regimes, the quantum dynamics exhibits three remarkable characters: (i) The excitation transfer between the adjacent pairs of SRs with fast oscillations dominates the evolution process due to the large hopping rate, and the weak NVE-SR coupling process leads to the excitation transfer between SR and the local NVE with a slow speed. (ii) The populations in the lateral SRs (NVEs) show oscillations between zero and varied maximum, but the population in the middle site oscillates with the fixed amplitudes. (iii) The maximal population in the middle NVE is much smaller than those in the lateral NVEs, namely, the excitation remains mainly in the lateral NVEs. The dynamics of the NVE is a linear superposition of such two oscillations, as shown in Fig. 3(b3), where the ratio between the slow and fast oscillation periods is about 7.5, and a fast and violent oscillation with a weak amplitude is superimposed on the slow oscillation. Here the fast oscillation governed by the photon hopping represents the process of inter-site transfer of excitation, and the slow oscillation governed by the coupling constant G c is actually a complete transfer of excitation between the nonlocal modes of SR. In the case of  = . G 0 1 c , the typical oscillation period times could be calculated as about 0.4 μs and 3 μs.
We now turn our attention to the situation involving both the dissipation and the inhomogeneous broadening of the spin ensembles. The dynamical evolution of the populations regarding SRs and NVEs could be numerically simulated by calculating the motion equations for all the bosonic species involved in the dynamics, where each NVE has a q-Gaussian distribution with L(δ i ) around a central spin frequency δ B , and each NVE is divided into several homogeneous subensembles containing N μ spins. Although q determines the specific profile such as how tall the peak of the distribution is and how fast the tails of the distribution fall off, it turns out that different value of q will bring similar results about the excitation distribution in each site, through numerical calculation. Therefore, here we only present the results of the Lorentzian distribution with q = 2, which is a typical model to study the inhomogeneous broadening of the NV − centers 40,43,50 . As shown in Fig. 4(a1,b1), the initial excitation is identical to the ideal case and the system displays a series of damped oscillations under various experimental situations, which also reflects the intricate balance and competition between the local NVE-SR couplings and the adjacent-site photon hopping. Physically speaking, the time-dependent superposition of the six harmonic wave functions ultimately determines the time evolution of all the one-excitation states. So the overall dynamics may exhibit some complex interferences or competing effects, and it exhibits an oscillating behavior with the amplitude decaying as time goes by and the enlarged oscillation period. To distinguish these two features brought by the different decoherence channels (the dissipation and the dephasing induced by the inhomogeneous broadening), we have studied their influences separately and found that the amplitude decaying mainly results from the dissipative effects with respect to the NVE and SR, and the presence of inhomogeneous broadening suppresses the population transfer between the SR and the local NVE, namely, the dephasing effect enlarges the oscillation period. From the relation 2 , one can find that the values of G μ decrease with growth of the offset frequency δ μ , which is related to inhomogeneous broadening of the spin ensemble. That is the reason why the period of the excitation transfer between SR and the local NVE is enlarged, when the dephasing effect is considered.

Entanglement dynamics.
Combined with the microwave photon-hopping process with great flexibility, a fully tunable model for entangling multi-NVEs can be expected in the strong collective-coupling regime, through engineering a quantum control on resonator-assisted Raman transitions with a high degree of freedom. The dynamics of the system shown above also suggests that the high-degree entanglement among the three NVEs could be achieved through accurately adjusting the tunable parameters in the present model. In our model, both the energy loss of the system and the homogeneous broadening of the ensemble have been considered. With the dissipative and dephasing effects in simulating the entanglement dynamics, an effective approach is the full phenomenological quantum master equation where , and κ (γ hom ) is the decay rate of SR (NVE). Γ is the dephasing rate of the NVE related to the inhomogeneous width γ s of the spin ensemble with γ s = 2Γ under the weak field approximation 50 .
In what follows, we focus on the dynamics of entanglement among three NVEs using the concept of the lower bound of concurrence (LBC) 55,56,61,62 . The concurrence for a general pure tripartite state ϕ ∈  ⊗ 1   ⊗2 3 has the following form  j and k (i, j, k = 1, 2, 3 and i ≠ j ≠ k). This result can be extended to the tripartite mixed state ρ using the convex roof denotes the collection of all the possible pure states into which the mixed state ρ is decomposed as ρ ϕ ϕ = ∑ p i i i i with the positive p i (normalized already). Here, an exact numerical algorithm is required to find the global minimum in the optimization procedure of  ρ ( ) 3 involved in a large number of free parameters, and the lower bound to this measure is available with the following definition 55,56,61,62,67   In our scheme, entanglement of the separate NVEs can be realized through the following channels: (i) Bosonic mode b l in the l-th NVE is entangled with the bosonic mode â l in the l-th SR through the local collective coupling. (ii) This entangled state transfers from the lst SR to the (l ± 1)-th one through the adjacent-site photon hopping with the transfer rate  , where the transfer process could be mediated by the capacitor acting as a tunable coupler. (iii) Entanglement of the bosonic mode b l in l-th NVE with the bosonic mode + b l 1 in (l + 1)-th NVE is achieved via the linear mixing mechanism. Note that both the collective coupling strength G c and the hopping rate  are tunable independently. As a result, entanglement among multiple NVEs in the 1D array of SR can be realized in a controllable way.
In the presence of both the dissipative and the dephasing effects, Fig. 5 characterizes the time-dependent entanglement of three NVE by the LBC. We find that the inhomogeneous broadening has a detrimental effect on the entanglement of spin ensembles, where the optimal LBC becomes smaller and smaller with growth of the dephasing rate Γ . This implies that, to obtain a high-degree entanglement of NVEs, we have to suppress the inhomogeneous broadening as much as we can using the recent experiment methods, such as the spin echo technology 58 .
Next, two distinct cases are considered. The first case plotted in Fig. 5(a1-a3) is for the fixed value of the hopping rate  and the increasing collective coupling strength G c . One can find that, for the same hopping rate  , the increase of the driving frequency G c could induce an obvious growth of oscillation frequency of entanglement and a slight accretion of the maximal amplitude with respect to LBC, in the presence of dephasing effect. It implies that the collective coupling strength G c controlling the linear mixing interaction process has a positive effect on the entanglement generation among NVEs. The second case presented in Fig. 5(b1-b3), with identical coupling strength G c and various hopping rates  , describes a remarkable result that the maximal amplitude of the LBC increases if the hopping rate  has an increment. Note that this feature is obvious only in the region where the hopping and dephasing rates are small. Additionally, as shown in the small hopping case (Fig. 5(b3)), the endurance period of the entanglement among spin ensembles is extended to be twice of those in Fig. 5(b1), where the hopping rate  is relatively large. Given the same evolution time, the higher photon hopping rate means more probabilities for photons shuttling between different sites, which results in a larger amount of entanglement among the NVEs located at different positions. Besides, slower photon-hopping rate gives less chances for the dissipation through the SRs in different sites. Thus, it is the reason that the lifetime of the entanglement becomes longer when the value of  decreases.
Therefore, a rich entanglement dynamics of the NVEs also reflects the intricate balance and competition between the two different kinds of interaction processes. One is the collective magnetic coupling process between the resonator mode and collective mode of the local NVE, and the other is the photon hopping process between the adjacent sites along the whole 1D array. These two related tunable parameters (Ω and  ) influence the LBC dynamics of the NVEs in different ways such as the maximal amplitude and the duration of the oscillation. To provide a more complete picture about how to achieve an optimized entanglement, the dependence of the average LBC during a certain period of time on the parameter space {Ω,  } is plotted in Fig. 6, where the optimal point can be extracted. The average LBC reaches the maximal value 0.43 at the point with Ω/2π = 50 MHz and π = . /2 7 4  MHz. Around the optimal point, a high-degree entanglement of the spin ensemble can still be obtained. It implies that the LBC evolution can be optimized by appropriately tailoring those key parameters. This provides a useful and effective way to controlling the dynamics of entanglement among the long-distance NVEs distributed in different sites.

Discussion
We now survey the relevant experimental parameters. Firstly, the external static magnetic field B stat could be applied by placing a copper box containing the resonator and diamonds in the center of two pairs of perpendicular Helmholtz coils 2,69 . If we set B stat to be 15.4 mT, then the Zeeman splitting becomes δ B /2π ≈ 500 MHz. For the inhomogeneous broadening, we choose the frequency δ i /2π ranging from − 70 MHz to 70 MHz, which ensures both the convergence property and the adiabatic elimination condition δ  B δ i . Besides, the values of Ω/2π are restricted within the regime {0,100} MHz in our model, which satisfy the adiabatic elimination condition, such as δ Ω  B , and are available experimentally 70 . Also, the photon-hopping rate π /2  can be tuned within the regime {0,10} MHz depending on the size of the capacitor 53,54,71 . For a SR with a quality factor Q ∼ 3 × 10 6 , the photon loss rate can be calculated as κ/2π ≈ 0.001 MHz 2 . On the other hand, the electron spin relaxation time T 1 of the NVE can reach up to 10 s at low temperature under an appropriately chosen magnetic field 2,72,73 . The dephasing time reaches T 2 > 600 μs for a lower nitrogen density sample with natural abundance of 13 C at room temperature 58 .
In our model, we have assumed the homogeneous couplings all over the sites of the 1D array, such as the identical number of NV − centers in each spin ensemble, and the identical hopping strength at each pair of SRs. Needless to say, this assumption would not be well satisfied in realistic experiments due to small deviation of the key parameter values induced by the fabrication errors and manipulation inaccuracies. However, treating those imperfections goes beyond the scope of the present paper. Noticeably, the effects of disorder in arrays of coupled cavities have been studied for both small-and large-scale arrays 74,75 . We emphasize that the initial excitation can be seeded in the SR or the NVE of any site. Due to symmetry, there are analogical dynamics of the system for the initial excitation in the left NVE or the right NVE. Note that the distinct advantages of the present system, such as the individual accessibility and high tunability of the parameters, make this NVE-SR array an ideal platform for investigating controllable quantum dynamics and other applications in quantum information processing, such as the excitation transfer and entanglement generation between distant NVEs in different sites, which are prerequisites for realization of spin-based distributed quantum networks.
In conclusion, we have considered a hybrid system consisting of multiple NVEs coupled to an array of SRs respectively, and the adjacent SRs are connected by capacitors. We have quantitatively simulated the controlled evolution by modulating some key external parameters, and found the way to extract proper experimental parameters for optimal entanglement of the NVEs using existing experimental technologies, even in the presence of large inhomogeneous broadening of the spin ensemble. Although what we have presented above is for a three-site model, the method and results can be easily extended to a longer array in one dimension. Therefore, our scheme provides another route towards building a distributed architecture, and our findings demonstrate that this hybrid quantum system provides a realistic platform for investigating quantum dynamics of spin ensembles and generating entanglement among distant spin ensembles. Fig. 1, the Hamiltonian of each site can be described by (in units of = 1  ) Analytical solutions to the simplified model. Here we present the analytical solutions of the ideal model, where both the dissipative and the dephasing effects are not considered, and the motion equations can be reduced to