1. Introduction

Over past decades, Schrödinger cat [1] has been the unique example of a paradoxical but fundamental property of quantum mechanics, reflecting that quantum superposition is at the core of quantum theory. Despite experimental evidences in individual atoms [2], photons [3], and molecules [4], this physicists’ favorite feline remains purely hypothetical if its size becomes larger, because objects beyond microscopic regime generally interact strongly with their surroundings, which force them away from superposition.

Great efforts have been devoted to exploring superposition states in mesoscopical systems, e.g., Schrödinger-cat states with six to ten qubits [5, 6], and to generating Schrödinger-cat states in an ultrafast fashion [7]. It is still strongly expected to observe quantum superposition in macroscopic systems for further understanding the peculiar laws in quantum world. Recently, control of vibrational degrees of freedom has been achieved in optomechanics [8] and nanomechanics [9]. Nevertheless, generation of macroscopically distinct superposition states, which are extremely sensitive to decoherence, is still experimentally challenging [10, 11].

In this work, we propose a practical scheme to create macroscopically distinct superposition states, i.e., macroscopic Schrödinger-cat states, in a diamond nitrogen-vacancy ensemble (NVE) coupled to two interacting flux-qubits. The NVE is a promising candidate for quantum information processing, which owns the capability from scalable quantum computing to long-distance quantum information transfer [12–14]. As the other component of our model, the flux-qubit is constructed by a compound Josephson-junction (CJJ) rf-superconducting quantum interference device (rf-SQUID) and we consider two such flux-qubits to be coupled by a coupler [15, 16]. There have been some schemes proposed so far to couple the superconducting flux-qubits to the spin ensembles by a magnetic field [17–20], and also some experimental observations for large couplings between the flux-qubits and the NVEs [21–23]. In contrast, our scheme works based on a giant Kerr nonlinearity appearing in the strong coupling regime of the system under a weak driving field, which is relevant to an N-type level configuration (i.e. the transitions between any two of the levels in the shape of the letter N) [24–27] of the coupled flux-qubits. Compared with the previous method for generating the macroscopic entanglement of NVE [28], the initialization of the quantum data bus into macroscopic cat state is not required in our scheme. Due to the giant Kerr nonlinearity, distinct superposition states could be generated in the NVE within a very short time. The idea can also be applied to generation of large quadrature squeezing states in the NVE and extended to other ensemble systems.

2. Model and Hamiltonian

Consider a superconducting diamond hybrid system as shown in Fig. 1(a) where an NVE, involving M NV-centers, is placed above the two flux-qubits. The latter is of Ising spin glass architecture [15, 16] with coupled CJJ rf-SQUIDs.The system experiences an external magnetic field $B_q^j$ generated by the circularly persistent current of the jth flux-qubit that is in parallel with the direction of S_x of the NVE (defined below). Under the co-resonance conditions [29], the system is described in units of $\hslash =1$ by the Hamiltonian as below,

The first two terms in Eq. (1) denoting the flux-qubits energy, as plotted in Fig. 1(b), can be expressed by the N-type level configuration with the eigenenergies $\left[E_4,E_3,E_2,E_1]=[\sqrt{J^2+{\mathrm{\Omega }}_q^2},J,-J,-\sqrt{J^2+{\mathrm{\Omega }}_q^2}\right]$. Transforming H into a rotating frame through $H^{\text{rot}}=e^{-iRt}H{e}^{iRt}+R$ with $R=D\left|1〉〈1\right|+\left(D-\overline{\omega }\right)\left|2〉〈2\right|-\overline{\omega }\left|4〉〈4\right|+D{\displaystyle \sum }_{i=1}^M|0{〉}_{ii}〈0|$, the effective Hamiltonian under rotating-wave approximation is

 

Fig. 1 (a) Schematic diagram of the hybrid quantum system composed of two coupled CJJ rf-SQUID qubits and an NVE. Inset: the level structure of an NV-center [30] under a weak classical field with an amplitude ξ_p coupling the states $|\pm 〉$ to $|0〉$. (b) The N-type level structure for the coupled flux-qubits with $\omega =J+\sqrt{J^2+{\mathbf{\Omega }}_q^2}$ and $\overline{\omega }=\sqrt{J^2+{\mathbf{\Omega }}_q^2}-J$. A control field with Rabi frequency $\mathbf{\Omega }$ and frequency ${\omega }_c=2J$ is in resonance with the states $|2〉$ and $|3〉$. (c) Three-dimensional plot of $\underset{10}{\log }[g^2\left(0\right)]$ as functions of $\mathbf{\Omega }$ and ϵ_p. The inset denotes the nonlinearity factor χ varying with the amplitude ϵ_p of the driving field along the gradient direction (red curve). Here the coupling strengths $\overline{\lambda }/\kappa =400$ and $\tilde{\lambda }/\kappa =200$. The decay paths are set as ${\gamma }_{41}={\gamma }_{32}=0.1\gamma$ and ${\gamma }_{31}={\gamma }_{21}={\gamma }_{42}={\gamma }_{43}=\gamma$. Where, ${\gamma }_{ij}$represent the decay rate from the states $|i〉$ to $|j〉$. Other parameters are $\mathbf{\Delta }=\delta =\kappa$ and $\gamma /\kappa =0.4$.

Download Full Size | PDF

3. Kerr nonlinearity

Due to no selection rule in the flux-qubits, all the decay paths are possible in the N-type level structure. As such, the Hamiltonian in Eq. (3) can be reduced to an effective self-Kerr nonlinearity interaction form $H_{\text{eff}}=\chi a^{†2}{a}^2$ [27, 29]. By adiabatically eliminating the flux-qubit operator [27, 32], the nonlinearity factor χ is given by $\chi ={\overline{\lambda }}^2/{\mathrm{\Omega }}^2[{\tilde{\lambda }}^2\mathrm{\Delta }/\left({\gamma }_{42}^2+{\mathrm{\Delta }}^2\right)-{\overline{\lambda }}^2\delta /\left[\left({\gamma }_{31}+{\gamma }_{32}{)}^2+{\delta }^2\right]\right]$, which generally holds in the weak coupling limit ${\overline{\lambda }}^2/{\mathrm{\Omega }}^2\ll 1$ [29] (the detail provided in Appendix A). However, when the system has all operational decay paths, a very large nonlinearity is possible [29] and observable by measuring the second-order correlation function $g^{\left(2\right)}\left(\tau \right)=〈a^{†}\left(t\right)a^{†}\left(t+\tau \right)a\left(t+\tau \right)a\left(t\right)〉/|〈a^{†}\left(t\right)a\left(t\right)〉{|}^2$. Under the weak driving field, the effective self-Kerr nonlinearity χ can be obtained from the analytical result of $g^{\left(2\right)}\left(0\right)$ [29, 33] by solving $\dot{\rho }=0$ from the master equation [33],

 

Fig. 2 (a) Time evolution of the fidelity $\mathcal{F}$. (b) Wigner function $W\left(\beta \right)$ at five different time points corresponding to t₁, t₂, t₃, t₄ and t₅ in (a). (c) Probability distribution of the rotated quadrature operator P_X at time points t₃ (blue solid curve) and t₄ (black dashed curve). Here we choose $\kappa =1$, $\gamma =0.4$, $\mathbf{\Delta }=\delta =1$, $\alpha =2e^{i\pi /4}$, ${ϵ}_p=2$, $\tilde{\lambda }=300$, $\mathbf{\Omega }=50$ and decay paths are chosen same as in Fig. 1(c).

Download Full Size | PDF

4. Macroscopic superposition states

Using the self-Kerr nonlinearity, we focus on generating coherent state superpositions of the NVE, such as the Schrödinger-cat state $|\text{cat}〉=\left({\mathrm{e}}^{-\mathrm{i}\pi /4}\left|\mathrm{i}\alpha 〉+{\mathrm{e}}^{\mathrm{i}\pi /4}\right|-\mathrm{i}\alpha 〉\right)/\sqrt{2}$ with the amplitude α [36]. To this end, we start from preparing a coherent state $|\alpha 〉$ of the NVE by the displacement operator $D\left(\alpha \right)=e^{\alpha a^{†}-{\alpha }^{*}a}$, and then the superposition of coherent states is created under the government of the effective self-Kerr nonlinearity $H_{\text{eff}}$ (see Appendix C for detail). Using Eq. (4), we have evaluated the produced state in the presence of dissipations by the Uhlmann fidelity $\mathcal{F}=\sqrt{〈\text{cat}\left|{\rho }_{\text{cat}}\right|\text{cat}〉}$, as shown in Fig. 2(a). The fidelity oscillates in time and reaches the maximum $0.88$ within a very short time. As shown later, the fidelity would not be unit even in absence of any decay effect, which is due to unexpected effects beyond nonlinearity from $H_{\text{eff}}$ under the weak decay condition. So we call below the maximum fidelity we could reach as the optimal fidelity, and the time reaching the optimal fidelity as the optimal time ${\tau }_{op}$.

In order to scrutinize the quantum coherence and interference effects during the cat creating process, we have examined the Wigner function $W\left(\beta \right)=\left(2/\pi \right)\text{Tr}\left[\mathrm{D}\left(-\beta \right)\widehat{\rho }\left(\mathrm{t}\right)\mathrm{D}\left(\beta \right)\mathrm{P}\right]$ with $P={(-1)}^{a^{†}a}$ [37], as plotted in Fig. 2(b), where the interference patterns implies the coherent superposition of the states. Besides, we also employ the probability distribution $P_X{=}_e〈X\left(\theta \right)\left|{\rho }_e\right|X\left(\theta \right){〉}_e$ to portray quantum superposition properties teprl-116-163602,qo, where $|X\left(\theta \right){〉}_e$ is the eigenstate corresponding to the rotated quadrature operator $\widehat{X}\left(\theta \right)=\left(a{e}^{-i\theta }+a^{†}{e}^{i\theta }\right)/\sqrt{2}$. In order to maximize the interference, we choose the rotation angle to be ${\theta }_0=\arg [\alpha ]$. The oscillation in both curves in Fig. 2(c) indicates the quantum interference between the superposition components. In addition, we have also noticed from the values that with the position of the maximal P_X approaching $X\left({\theta }_0\right)=0$, the superposition state is closer to the state $|\text{cat}〉$ and quantum interference becomes stronger.

 

Fig. 3 (a) Optimal fidelity of the cat state as a function of the cat’s size. Inset: the dots denote the oscillating period of the fidelity for $|\alpha {|}^2〉5$ and the line means a constant obtained by fitting these dots. (b) Optimal time $ct{\tau }_{\mathbf{op}}$ in variation with the cat’s size. Inset: the dots represents $\mathbf{\Delta }{\tau }_{\mathbf{op}}$ between the nearest-neighbor ladders and the curve denotes an exponential fitting. (c) Dynamical evolution of the fidelity for three sizes of the cat states around the inflection point $|\alpha {|}^2=5.05$, which is labeled as a red dot in (a). (d) Probability distribution of the rotated quadrature operator P_X at time points $rotect\kappa {\tau }_{\mathbf{op}}=0.02$ ($|\alpha {|}^2=4.8,5$) and $0.025$ ($|\alpha {|}^2=5.1,5.8$). Here the parameters used are the same as in Fig. 2.

Download Full Size | PDF

In general, the fidelity of the prepared optimal cat state is sensitive to its size under the dissipative condition. In contrast to a simple asymptotic decay as usual, the decay in our case shows a subtle behavior (Fig. 3(a)). To understand the physics behind the fidelity, we define ${\tau }_{\text{op}}$ as the time reaching the optimal fidelity of the cat state. We see in Fig. 3(a) that the fidelity decreases rapidly when $|\alpha {|}^2<2$ and ${\tau }_{\text{op}}$ decreases accordingly [see Fig. 3(b)]. Then the fidelity has a small rise, followed by a large descent until $|\alpha {|}^2\approx 5$, during which ${\tau }_{\text{op}}$ remains almost unchanged. For $|\alpha {|}^2〉5$, the fidelity behaves as a periodic oscillation in a linear decay and the corresponding ${\tau }_{\text{op}}$ increases discretely as stairs. Quantitatively, the laws of the changes for the case of $|\alpha {|}^2〉5$ can be described as (some details can be seen in Appendix D),

 

Fig. 4 (a) and (b) Fidelities of the cat state prepared under different flux-qubit dissipation rates and NVE decay rates, where we set $u_0/2\pi =0.1$ MHz, and for convenience of presentation the parameters below are written in units of u₀. We have $\kappa =$1 in (a) and $\gamma =1$ in (b). The solid, dashed and dotted curves in both (a) and (b) correspond to $\tilde{\lambda }=300,200$ and $100$, respectively. Inset: time evolution of the fidelity with $\widetilde{t\lambda }=300$. (c) Probability distribution of the rotated quadrature operator P_X at the time point t₄ of Fig. 2. The dotted, dashed and solid curves correspond to the excitation numbers $n_c=$1, 2 and 5. In (c), $\kappa =1$ and $\gamma =0.4$. Inset: Wigner function $W\left(\beta \right)$ with respect to $n_c=$1 and 5. Other parameters in (a-c) are $pha=2e^{i\pi /4}$, ${ϵ}_p=2$, $\overline{\lambda }=400$ and $\mathbf{\Omega }=50$, respectively.

Download Full Size | PDF

5. The influence from decoherence and dissipation

To check more seriously the influence from the system dissipations, we have numerically investigated the fidelity variation of the prepared cat states with respect to the system decay. As shown in Figs. 4(a) and 4(b), the fidelity of the cat state is damaged by both the flux-qubit dissipation and the NVE decay. As a comparison, one can find that under the same condition the decay of the NVE bring a more serious dissipation than that of flux-qubit due to the fact that the Schrödinger-cat state is prepared in the NVE and thus the NVE decay directly leads to the fidelity decrease. In contrast, the flux-qubit dissipation is associated with an indirect coupling between the flux-qubit and the NVE.

Moreover, considering some imperfect factors in the NVE, we assume the NVE under a finite-temperature reservoir, which modifies the second term in the right-hand side of Eq. (4) as $\kappa \left(n_c+1\right)\mathcal{D}\left[a,\rho \right]/2+\kappa n_c\mathcal{D}\left[a^{er},\rho \right]/2$ with n_c representing the boson number regarding the noise of the reservoir. To explore this imperfection, we employ the Wigner function and the probability P_X in Fig. 4(c). The probability P_X demonstrates less evident fringes with increase of the reservoir noise, although quantum interference of the coherent superposition states is still distinct even in the case of $n_c=5$. Besides, the Wigner function presents a different view for the noise influence on the cat state that with the increase of the noise, the two components $|i\alpha 〉$ and $|-i\alpha 〉$ themselves turn to be more dominant than the interference in between, i.e., the darker ends for the larger n_c [inset of Fig. 4(c)]. This is the reason for the reduced fidelity of the cat state in the more noisy case. Nevertheless, as long as the reservoir noise is not big enough, quantum interference pattern remains existing.

 

Fig. 5 (a) Fidelities of the cat state prepared under different dephasing rates κ_d. (b) A detail dynamical evolution of the fidelities for three different dephasing rates at a coupling $\tilde{\lambda }=300$. (c) The Wigner function $W\left(\beta \right)$ atthe optimal cat state position of in (b) for the case ${\kappa }_d=0.5$. Here $\kappa =1$ and $\gamma =0.4$ and the other parameters are same as in Fig. 4.

Download Full Size | PDF

In addition, inhomogeneous broadening is another key problem in manipulation of the NVEs, which is caused by magnetic dipolar interactions with the nuclear or excess electron spins in diamond [21, 22] and results in the dephasing effects [19]. Previous studies have found that the magnetic field applied along a special direction [100] of each NV-center can largely suppress the inhomogeneous broadening of the NVE [23]. In our case for effectively coupling two flux qubits via the NVE, we should also consider the inhomogeneous broadening resulted from spatial distribution of the NV-centers, which diversifies the coupling strength $g_{1,2}$, is very small and thus could be resorted to dephasing, following the solution in [19]. Under the weak field approximation, the dephasing rate κ_d of the NVE is related to the inhomogeneous width γ_s of the NVE by the relation ${\gamma }_s=2{\kappa }_d$. Therefore, a dephasing term ${\kappa }_d\mathcal{D}\left[a^{†}a,\rho \right]/2$ should be complemented into the master equation (at the right-hand side of Eq. (4)) to quantitatively describe the influence from the inhomogeneous broadening effect. Fig. 5 is the result of numerical simulation of such a master equation.Under the dephasing effect, the interference are still clearly visible in Fig. 5(c). Specially, in the case of ($\tilde{\lambda }=300\kappa$), a small inhomogeneous broadening does not distinctly damage the prepared macroscopic superposition states.

6. Discussion and conclusion

Fast generation and long-time remaining of Schrödinger-cat states are essential to deeper understanding of quantum theory and also to some practical applications. Since they can be measured and characterized experimentally by the quantum spectroscopy [38], the cat states can be employed in the quantum metrology, for example, as a sensitive electrometer [39], for improving the measurement accuracy of phases [40], for engineering a remote amplifier [41] and for other tasks of quantum metrology [42]. In addition, NVEs are practical for scalable quantum information processing, such as used as a magnetic field transducer [43], measuring the oscillating field [44], applied as a magnetometer [45] and nanometre positioning meter [46], detecting the molecular radical reaction [47] and a single nuclear spin [48]. As a result, studying Schrödinger-cat states in the NVEs would help keeping quantum coherence in the NVEs when they work as qubits in quantum computing or as nodes in quantum network.

On the other hand, experimentally detecting the Schrödinger-cat states in spin ensembles is not a trivial task. We can visualize the desired Schrödinger cat state by measuring the Wigner function of the NVE. Earlier experimental works have shown mature techniques to determine the density matrix of the resonators and to directly measure Wigner function in ion traps [49] as well as in microwave superconducting circuits [36, 50]. In general, the Wigner function of Schrödinger-cat states could be determined by the parity of the state [51], by the statistical moments of the field operator [52], or by reconstructing the density matrix of the system using a least square fit to each Husimi Q-function [53]. In our case, based on the self-Kerr nonlinearity, we displace the state of the NVE by a displacement pulse $D\left(\alpha \right)$ and wait for a variable evolution time. Then, by means of state tomography, we may employ the least-square fit to each Husimi Q-function for reconstructing the density matrix of the NVE, and then obtain the Wigner function by straightforward calculation.

Our scheme is feasible using current laboratory techniques [54, 55]. The available dissipation values are $\kappa /2\pi =0.1$ MHz for the NVE [22] and $\gamma /2\pi \in \left[0,100\right]$ kHz for this coupled flux-qubit four-level system [21, 56, 57]. For other parameters, the tunneling energy ${\mathrm{\Omega }}_q/2\pi$ of the flux-qubit is from 1 MHz to 10 GHz [15]. The coupling energy J is related to the effective mutual inductance M and the qubit persistent current, implying that $J/2\pi =$ 2.3 GHz is available if $M=$1.5 pH and $\left|I_i^p\right|=1\mu \mathrm{A}$ [15]. With these values, the level spacing ω, as defined in Fig. 1(b), could match the zero-field splitting D of the NV-center ground state. The coupling strength between the flux-qubit and the NVE is regarding the flux-qubit persistent current strengths and the number of the NV-center, which could reach 116 MHz (about $\overline{\lambda },\tilde{\lambda }/\kappa \sim {10}^3$) for $\left|I_i^p\right|=$ 1 $\mu \mathrm{A}$ and $3.2\times {10}^7$ NV-centers. As such, the giant effective self-Kerr nonlinearity $\chi /\kappa$ could be ${10}^3$ under a weak driving field, implying that we can create a Schrödinger-cat state within a time about $200$ ns.

In conclusion, we have proposed a practical scheme to create macroscopic Schrödinger-cat states in a hybrid system with an NVE interacting with two coupled CJJ-flux qubits. The key point of our work is to try to generate Schrödinger-cat states in a fast and decoherence-suppressed fashion. Numerical results have demonstrated that our scheme works well in a wide parameter range in the presence of real dissipations of the system. The idea can also be used to demonstrate quadrature squeezing and extended to other ensemble systems.

Appendix A establish the model and construct the Hamiltonian of system

1. Constructing a four-level system

The compound Josephson-junction (CJJ) rf-superconducting quantum interference device (rf-SQUID) is employed as a qubit and two such qubits are coupled by a CJJ rf-SQUID coupler (see Fig. 1). The Hamiltonian for this architecture is a quantum Ising spin glass [15, 16] that is written as $H_s={\displaystyle \sum }_{i=1}^2\frac{1}{2}\left[{\omega }_q^i{\sigma }_z^i+{\mathrm{\Omega }}_q^i{\sigma }_x^i\right]+J{\sigma }_z^1{\sigma }_z^2,$ where ${\omega }_q^i/2\pi =2\left|I_i^p\right|{\mathrm{\Phi }}_i^x$ and ${\mathrm{\Omega }}_q^i$ are the bias and tunneling strength of the ith flux-qubit, respectively, and $J=M\left|I_1^p\right|\left|I_2^p\right|$ is the coupling strength. Here $\left|I_i^p\right|$ is the magnitude of the qubit persistent current, ${\mathrm{\Phi }}_i^x$ and M denote the external flux bias and mutual inductance, and they are $insitu$ tunable. ${\sigma }_z^i$ and ${\sigma }_x^i$ are Pauli spin operators. In addition, if the CJJ rf-SQUIDs are identical and the individual tunneling energy can be modulated by local flux tuning at the co-resonance point ${\omega }_q^i=0$ [29], we set ${\mathrm{\Omega }}_q^1={\mathrm{\Omega }}_q^2={\mathrm{\Omega }}_q$ and rewrite the Hamiltonian H_s as $H_s^{\text{'}}=\frac{{\Omega }_q}{2}{\displaystyle \sum }_{i=1}^2{\sigma }_x^i+J{\sigma }_z^1{\sigma }_z^2.$ Thus, the eigenenergies of the Hamiltonian $H_s^{\text{'}}$ can be solved straightforwardly as $\left[E_4,E_3,E_2,E_1]=[\sqrt{J^2+{\mathrm{\Omega }}_q^2},J,-J,-\sqrt{J^2+{\mathrm{\Omega }}_q^2}\right]$, with the corresponding eigenstates denoted by $|4〉$, $|3\overline{a}ngle$, $|2〉$ and $|1〉$. Defining an angle $\cos 2\theta =J/\sqrt{J^2+{\mathrm{\Omega }}_q^2}$, the above four states are presented as

(6)$$\left(\begin{array}{l}|4〉\hfill \\ |3〉\hfill \\ |2〉\hfill \\ |1〉\hfill \end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{llll}\text{cos }\theta \hfill & \text{sin }\theta \hfill & \text{sin }\theta \hfill & \text{cos }\theta \hfill \\ 1\hfill & 0\hfill & 0\hfill & -1\hfill \\ 0\hfill & 1\hfill & -1\hfill & 0\hfill \\ \text{sin }\theta \hfill & -\text{cos }\theta \hfill & -\text{cos }\theta \hfill & \text{sin }\theta \hfill \end{array}\right)\left(\begin{array}{l}|ee〉\hfill \\ |eg〉\hfill \\ |ge〉\hfill \\ |gg〉\hfill \end{array}\right).$$

Defining the transition frequencies as ${\omega }_{i,j}=E_i-E_j$, thus we have ${\omega }_{32}=2J$ and ${\omega }_{42}={\omega }_{31}=\omega =J+\sqrt{J^2+{\mathrm{\Omega }}_q^2}$ with $K\equiv {\omega }_{42}/{\omega }_{31}=1$. Here we mention the exception, i.e., the case of nonzero ${\omega }_q^i$ for which the eigenenergies are solved by moving slightly off the co-resonance point through introducing equal-strength Zeeman terms in $H_s^{\text{'}}$, i.e., $H_s^{\text{'}\text{'}}=\frac{1}{2}{\displaystyle \sum }_{i=1}^2{\omega }_q^i{\sigma }_z^i+H_s^{\text{'}}$ [29]. For this case, K is no longer equal to one and a slight difference will be introduced between ${\omega }_{42}$ and ${\omega }_{31}$. Nevertheless, this case does not change the physical essence of our model. As a result, we will proceed following Eq. (6), that is,

(7)$$\left(\begin{array}{l}|ee〉\hfill \\ |eg〉\hfill \\ |ge〉\hfill \\ |gg〉\hfill \end{array}\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{llll}\text{cos }\theta \hfill & 1\hfill & 0\hfill & \text{sin }\theta \hfill \\ \text{sin }\theta \hfill & 0\hfill & 1\hfill & -\text{cos }\theta \hfill \\ \text{sin }\theta \hfill & 0\hfill & -1\hfill & -\text{cos }\theta \hfill \\ \text{cos }\theta \hfill & -1\hfill & 0\hfill & \text{sin }\theta \hfill \end{array}\right)\left(\begin{array}{l}|4〉\hfill \\ |3〉\hfill \\ |2〉\hfill \\ |1〉\hfill \end{array}\right).$$

Eq. (7) implies,

(8)$$\begin{array}{l}{\sigma }_z^1=\text{cos }\theta \left|4〉〈3\right|+\text{sin }\theta \left|4〉〈2\right|+\text{sin }\theta \left|3〉〈1\right|-\text{cos }\theta \left|2〉〈1\right|+\mathrm{h}.\mathrm{c}.,\hfill \\ {\sigma }_z^2=\text{cos }\theta \left|4〉〈3\right|-\text{sin }\theta \left|4〉〈2\right|+\text{sin }\theta \left|3〉〈1\right|+\text{cos }\theta \left|2〉〈1\right|+\mathrm{h}.\mathrm{c}\dots \hfill \end{array}$$

2. Coupling the four-level system to the NVE

The Hamiltonian for the NVE is given by $H_{\text{NVE}}={\displaystyle \sum }_{i=1}^M[D{S}_{z,i}^2+E\left(S_{x,i}^2-S_{y,i}^2\right)+g_e{\mu }_B{B}_{z,i}{S}_{z,i}]$, where $g_e=2.0028$, ${\mu }_B=14\text{ MHz}\cdot {\text{mT}}^{-1}$, and $D=2.88$ GHz is the zero-field splitting ($E<10$ MHz) [21, 30]. The Zeeman interaction $B_{z,i}{S}_{z,i}$ and strain-induced splitting term $E\left(S_{x,i}^2-S_{y,i}^2\right)$ are negligible. $S_{k,i}$ ($k=x,y,z$) are the spin-1 operators of the NV-center with the eigenstates $|\pm 〉$ and $|0〉$. The magnetic field generated by the circularly persistent current of the jth flux-qubit is assumed to be in parallel with the direction of S_x.

In addition, the coupling between the NVE and the flux-qubits is written as $H_{\text{int}}={\displaystyle \sum }_{i=1}^M{S}_{x,i}\left(g_1{\sigma }_z^1+g_2{\sigma }_z^2\right)$, where the coupling strength g_k can be estimated by $g_k=g_e{\mu }_B{B}_q^k/2$, $B_q^k={\mu }_0{I}_k^p/2R$ with $k=1,2$, ${\mu }_0=4\pi \times {10}^{-7}$ N$\cdot$A${ }^{-2}$ and R is the average distance between the flux-qubits and the individual NV-centers. g_k depends on I_k and R [17], for example, with a typical persistent current value in the flux-qubit $I_p=0.3\mu$A and the distance $R=1.2\mu$m the coupling strength is $g_k/2\pi \sim 4.4$ kHz [21]. Therefore, the total Hamiltonian can be written as $H_1=\frac{{\Omega }_q}{2}{\displaystyle \sum }_{i=1}^2{\sigma }_x^i+J{\sigma }_z^1{\sigma }_z^2+{\displaystyle \sum }_{i=1}^{M}D{S}_{z,i}^2+{\displaystyle \sum }_{i=1}^M{S}_{x,i}\left(g_1{\sigma }_z^1+g_2{\sigma }_z^2\right)$. Considering the relations in Eq. (8), we rewrite the Hamiltonian H₁ as $H_2={\displaystyle \sum }_{i=1}^4{E}_i\left|i〉〈i\right|+{\displaystyle \sum }_{i=1}^{M}D{S}_{z,i}^2+{\displaystyle \sum }_{i=1}^M{S}_{x,i}[\overline{G}\left(\cos \theta \left|4〉+\sin \theta \right|1〉\right)〈3\left|+\tilde{G}\left(\sin \theta \left|4〉-\cos \theta \right|1〉\right)〈2\right|+$h.c.$]$, where $\overline{G}=g_1+g_2$ and $\tilde{G}=g_1-g_2$. Transforming H₂ into a rotating frame via $H_3=e^{-iRt}{H}_2{e}^{iRt}+R$ with $R=D\left|1〉〈1\right|+\left(D-\overline{\omega }\right)\left|2〉〈2\right|-\overline{\omega }\left|4〉〈4\right|+D{\displaystyle \sum }_{i=1}^M|0{〉}_{ii}〈0|$ and $\overline{\omega }=\omega -2J$, we obtain

(9)$$\begin{array}{l}{H}_3=\mathrm{\Delta }|4〉〈4|+\delta |3〉〈3|+{\displaystyle \sum }_{i=1}^M[e^{iDt}{|b〉}_{ii}〈0|+h.c.][\overline{G}(\cos \theta e^{i\overline{\omega }t}|4〉+\\ \sin \theta e^{-iDt}||1〉)〈3|+\tilde{G}(\sin \theta e^{iDt}|4〉-\cos \theta e^{-i\overline{\omega }t}|1〉)〈2|+h.c.],\end{array}$$

where $|b〉=\left(\left|+〉+\right|-〉\right)/\sqrt{2}$, $\mathrm{\Delta }={\omega }_{42}-D$ and $\delta ={\omega }_{31}-D$. Because of ${\omega }_{31}={\omega }_{42}=\omega$ in above discussion, we have $\delta =\mathrm{\Delta }$. Due to the fact of $2D,D+\overline{\omega },D-\overline{\omega }\gg \mathrm{\Delta },\dot{e}lta$, under the rotating-wave approximation [17], Eq. (9) is reduced to the effective Hamiltonian as below,

(10)$$H_4=\mathrm{\Delta }\left|4〉〈4\left|+\delta \right|3〉〈3\right|+\text{sin }\theta {\displaystyle \sum }_{i=1}^M[\left|b{〉}_{ii}〈0\right|\left(\overline{G}\left|1〉〈3\right|+\tilde{G}\left|2〉〈4\right|\right)+h.c.],$$

which is the Hamiltonian H₁ in the main text by replacing $\overline{g}=\sin \theta \overline{G}$ and $\tilde{g}=\sin \theta \tilde{G}$.

We have to emphasize that the effect regarding E is ignored in above treatment, and our treatment below will follow Eq. (10). But if E is large enough, some modifications are necessary. For example, in the case of $\left|\delta \right|=\left|{\omega }_{31}-D-E\right|$ and $\left|\mathrm{\Delta }\right|=\left|{\omega }_{42}-D-E\right|\ll E$, $|b〉$ should be replaced by $|+〉/\sqrt{2}$. In contrast, if $\left|\delta \right|=\left|{\omega }_{31}-D+E\right|$ and $\left|\mathrm{\Delta }\right|=\left|{\omega }_{42}-D+E\right|\ll E$, $|b〉$ should be written as $|-〉/\sqrt{2}$. Meanwhile, for these two cases, D should be replaced by $D+E$ or $D-E$, respectively. With these modifications, following treatments still apply.

3. Deriving the effective self-Kerr nonlinearity Hamiltonian

We investigate the Hamiltonian H₂ (Eq. (3) in the main text) in the dressed states picture. In order to obtain the effective self-Kerr Hamiltonian $H_{eff}$, we will mainly concentrate on the analysis of H₂ excluding the pumping term $H_p={\epsilon }_p\left(a+a^{†}\right)$ as $H_K=\mathrm{\Delta }{\sigma }_{44}+\delta {\sigma }_{33}+\mathrm{\Omega }\left({\sigma }_{32}+{\sigma }_{23}\right)+a(\overline{\lambda }{\sigma }_{31}+\tilde{\lambda }{\sigma }_{42}+$h.c.$)$. Here, we define ${\sigma }_{ij}=|ij|$. Actually, the pumping term H_p as a probe field is weak in our work, thus it is assumed to be virtually negligible.

For convenience, we use $|bosonnumber,flux-qubitstate$ as the notation for the bare states. Then, the general dressed state of the Hamiltonian H_K can be written as $|\mathrm{\Phi }={\displaystyle \sum }_{n=0}^{\infty }{\displaystyle \sum }_{m=1}^4{c}_{nm}|n,m$. Here, $c_{nm}$ is the coefficients of the bare states $|n,m$. On the basis of the level of excitation of the whole system, the Hamiltonian H_K naturally separates into different manifolds, without the driving H_p. Obviously, the ground state of the system H_K is $|g=|0,1$. The bare states $|1,1$, $|0,2$ and $|0,3$, formed the first manifold of the system. Thus, according to Hamiltonian H_K, the dressed states in the first manifold can be written as $|e_0=\left(\overline{\lambda }|0,2-\mathrm{\Omega }|1,1\right)/X$ and $|e_{\pm }=\pm {\mu }_{\pm }|0,3+{\mu }_{\mp }\left(\mathrm{\Omega }|0,2+\overline{\lambda }|1,1\right)/X$. Where, $X=\sqrt{{\mathrm{\Omega }}^2+{\overline{\lambda }}^2}$, ${\mu }_{\pm }=\sqrt{\left(Y\pm \delta \right)/2Y}$, with $Y=\sqrt{{\delta }^2+4X^2}$. The corresponding eigenenergies are ${\epsilon }_0=0$, ${\epsilon }_{\pm }=\left(\delta \pm Y\right)/2$. For the nth manifold ($n\ge 2$), the four bare states of the system are $|n+1,1$, $|n,2$, $|n,3$, and $|n-1,4$.

According to the dressed states $|e_j$ ($j=0,\pm$), we define the polariton creation operators in the first manifold as [58, 59] $P_0^{†}=\left(\overline{\lambda }{\sigma }_{21}-\mathrm{\Omega }{a}^{†}\right)/X$ and $P_{\pm }^{†}=\pm {\mu }_{\pm }{\sigma }_{31}+{\mu }_{\mp }\left(\mathrm{\Omega }{\sigma }_{21}+\overline{\lambda }{a}^{†}\right)/X$, such that $|e_j=P_j^{†}|g$ with $j=0,\pm$. These operators satisfy the commute formula, $\left[P_i,P_j^{†}\right]={\delta }_{ij}$, $\left[P_i,P_j]=[P_i^{†},P_j^{†}\right]=0$ with $i,j=0,\pm$, which means that these polaritons are bosons.

In addition, according to the polariton creation operators $P_j^{†}$, we obtain ${\sigma }_{21}=\mathrm{\Omega }\left({\mu }_{-}{P}_{+}^{†}+{\mu }_{+}{P}_{-}^{†}\right)/X+\overline{\lambda }{P}_0^{†}/X$ and $a^{†}=\overline{\lambda }\left({\mu }_{-}{P}_{+}^{†}+{\mu }_{+}{P}_{-}^{†}\right)/X-\mathrm{\Omega }{P}_0^{†}/X$.

 

Fig. 6 Squeezing spectrum $S\left(\omega \right)$ under different conditions, where (a) the solid (dashed) curve denotes $\mathbf{\Omega }=30$ (36) and (b) $\mathbf{\Omega }=50$ (60) in unit of u₀ with $u_0=$0.1 MHz. Here we choose $\mathbf{\Delta }=\delta =5$ and ${ϵ}_p=2$. Other parameters are the same as in Fig. 1(c) of main text.

Download Full Size | PDF

Besides, for the Hamiltonian H_K, we assume that the detuning of levels $|3$ and $|4$ are sufficiently large ($\mathrm{\Delta }\gg {\gamma }_{42}$ and $\delta \gg \left({\gamma }_{31}+{\gamma }_{32}\right)/2$), so we can adiabatically eliminate levels $|3$ and $|4$. Then, the Hamiltonian H_K can be effectively written as $H_m\simeq \mathrm{\Delta }{\tilde{\lambda }}^2/\left({\gamma }_{42}^2+{\mathrm{\Delta }}^2\right)a^{†}a{\sigma }_{22}+\delta {\overline{\lambda }}^2/\left[\left({\gamma }_{31}^2+{\gamma }_{32}^2\right)+{\delta }^2\right]a^{†}a{\sigma }_{11},$ which was described just by the operators $\left\{a,{\sigma }_{11},{\sigma }_{12},{\sigma }_{22}\right\}$ (and their adjoints). With the expressions ${\sigma }_{21}$ and $a^{†}$ (and their adjoints), the Hamiltonian H_m can be expressed in terms of polariton operators $P_j^{†}$ with $j=0,\pm$. On the other hand, when the condition $\left|\delta \right|$, $\left|\mathrm{\Delta }\right|$, $\left|\overline{\lambda }\right|$, $\left|\tilde{\lambda }\right|\ll \left|{\epsilon }_{\pm }\right|$ is satisfied, the terms containing the operators $P_{\pm }^{†}$ ($P_{\pm }$) will oscillate with much higher frequency than the terms containing only $P_0^{†}$ (P₀). Therefore, under the rotating wave approximation, the resulting effective Kerr Hamiltonian has the simple form as ${\tilde{H}}_{eff}=\tilde{\chi }{(P_0^{†})}^2{P}_0^2$ with $\tilde{\chi }=\left[{\overline{\lambda }}^2/\left({\mathrm{\Omega }}^2+{\overline{\lambda }}^2\right)\right]\left\{\mathrm{\Delta }{\tilde{\lambda }}^2/\left({\gamma }_{42}^2+{\mathrm{\Delta }}^2\right)-\delta {\overline{\lambda }}^2/\left[\left({\gamma }_{31}^2+{\gamma }_{32}^2\right)+{\delta }^2\right]\right\}$. When the weak coupling limit ${\overline{\lambda }}^2/{\mathrm{\Omega }}^2\ll 1$ is satisfied, the effective Hamiltonian becomes $H_{eff}=\chi a^{†2}{a}^2$ with $\chi =\left({\overline{\lambda }}^2/{\mathrm{\Omega }}^2\right)\left\{\mathrm{\Delta }{\tilde{\lambda }}^2/\left({\gamma }_{42}^2+{\mathrm{\Delta }}^2\right)-\delta {\overline{\lambda }}^2/\left[\left({\gamma }_{31}^2+{\gamma }_{32}^2\right)+{\delta }^2\right]\right\}$.

Appendix B the quadrature squeezing

The quadrature squeezing can also be created by the self-Kerr nonlinearity [34, 35]. Define $\epsilon =-2i\chi \text{Tr}\left[{\mathrm{a}}^{†}\mathrm{a}{\rho }_{\text{ss}}\right]$ with ${\rho }_{\text{ss}}$ denoting the steady state solution of Eq. (4). If $|epsilon{|}^2\le {\kappa }^2/3$, the squeezing spectrum has a single peak at $\omega =0$. Otherwise, it splits into two peaks at ${\omega }^2=3|\epsilon {|}^2-{\kappa }^2$ [35]. For the first condition, the squeezing spectrum is written by

(11)$$S\left(\omega \right)=\frac{2{\xi }^2\sqrt{{(1-3{\xi }^2)}^2+16{\xi }^2}-\xi \left[\left({\omega }^2/{\kappa }^2+1-3{\xi }^2\right)\left(1-3{\xi }^2\right)+16{\xi }^2\right]}{\sqrt{{(1-3{\xi }^2)}^2+16{\xi }^2}[\left({\omega }^2/{\kappa }^2+1-3{\xi }^2{)}^2+12{\xi }^2\right]}$$

with the dimensionless parameter $\xi =\left|\epsilon \right|/\kappa$. Shown in Fig. 6(a) is the spectrum with its minimum at $\omega =0$. An optimal squeezing spectrum obtained at $\xi =1/\sqrt{3}$ is $S\left(\omega \right)=-2{\xi }^2/\left({\omega }^4/{\kappa }^4+12{\xi }^2\right)$ with the maximal squeezing $S\left(0\right)=-1/6$. If the second condition is satisfied, the maximum squeezing spectrum is $S\left(\omega \right)=-2{\xi }^2/[\left({\omega }^2/{\kappa }^2+1-3{\xi }^2{)}^2+12{\xi }^2\right)$ [see Fig. 6(b)] and the two peaks appear at ${\omega }_{\pm }=\pm \kappa \sqrt{3{\xi }^2-1}$ with a same squeezing $S\left({\omega }_{\pm }\right)=-1/6$.

 

Fig. 7 (a) Dynamical evolution of the fidelity for different sizes of the cat states. The curves from left to right correspond to the cat states with the size $|\alpha {|}^2$ changing form 2 to 6.5. (b) and (c) are zooming-in plots of Figs. 3(a) and 3(b) of the main text, respectively.

Download Full Size | PDF

Appendix C the ideal Kerr state

The Kerr nonlinearity can be observed in the evolution of a coherent state [36]. An ideal Kerr state arises from the evolution of a coherent state $|\alpha 〉$, which is the eigenstate of the annihilation operator a under the effective Hamiltonian $H_{\text{eff}}=\chi a^{†2}{a}^2$. After evolving for time t, the ideal Kerr state is written as $\left|\psi \left(\alpha ,\theta \right)〉=e^{-|\alpha {|}^2/2}{\displaystyle \sum }_{n=0}^{\infty }\frac{{\alpha }^n}{\sqrt{n!}}{e}^{-in\left(n-1\right)\theta }\right|n〉$ with $\theta =\chi t$. Since $n\left(n-1\right)$ is an even number, the period for $|\psi \left(\alpha ,\theta \right)〉$ is $T=\pi$ that means $\left|\psi \left(\alpha ,\theta +T\right)〉=\right|\psi \left(\alpha ,\theta \right)〉$. Supposing $\theta =\frac{N}{M}T$ with the integer M is indivisible with respect to the integer N, we rewrite the state $|\psi \left(\alpha ,\theta \right)〉$ by a superposition of coherent states as $\left|\psi \left(\alpha ,\theta =\frac{N}{M}T\right)〉={\displaystyle \sum }_{k=0}^{2M-1}{c}_k\right|e^{i{\phi }_k}\alpha 〉.$ Here, ${\phi }_k=k\pi /M$ and $c_k=\frac{1+{(-1)}^{k-N\left(M-1\right)}}{2M}{\displaystyle \sum }_{n=0}^{M-1}\exp \left\{-\frac{i\pi }{M}\left[nk+Nn\left(n-1\right)\right]\right\}$ shows that only M different coherent states exist. For simplicity, we set N = 1 and particularly consider M = 2 that is for the Schrödinger-cat state $|\alpha ,\frac{\pi }{2}〉=\frac{1}{\sqrt{2}}\left(e^{-i\pi /4}\left|i\alpha 〉+e^{i\pi /4}\right|-i\alpha 〉\right)$. The phase difference between two neighboring cat-states is i which could be applied to speculate the effective anharmonicity parameter χ in the evolution based on the master equation.

Appendix D the detail for the dynamical evolution of fidelity in the preparation of different-size cat state

Zooming-in plots around the inflection point $|\alpha {|}^2=5.05$of Fig. 3 are shown in Fig. 7(a) where the change of the optimal fidelity can be explicitly seen in Fig. 7(a). The first oscillating period in Fig. 3 is amplified in Fig. 7(b) and a sharp change in fidelity can be found around $|\alpha {|}^2=5.05$. On the other hand, an amplification for the first stair in Fig. 3 is shown in Fig. 7(c) where a sudden jump at the inflection point $|\alpha {|}^2=5.05$ occurs.

 

Fig. 8 (a) and (c) Optimal fidelity of the cat state as a function of the cat’s size. Inset: the dots denote the oscillating period of the fidelity for $|\alpha {|}^2〉5$ and the line means a linear fitting of these dots. (b) and (d) The time ${\tau }_{\mathbf{op}}$ for obtaining the optimal fidelity of the cat state. Inset: the dots represent the corresponding spacing $\mathbf{\Delta }{\tau }_{\mathbf{op}}$ between the nearest neighbor two ladders and the curve denotes an exponential fitting. The black, red and green data and curves in (a) and (b) correspond to $\gamma =0.2,0.4$ and 1 with $\kappa =1$, and in (c) and (d) they are corresponding to $\kappa =0.5,1$ and 2 with $\gamma =0.4$. The other parameters are same in Fig. 3 of main text.

Download Full Size | PDF

Figure 8 presents some details about the fidelity variation under different flux-qubit dissipation rates γ and NVE decay rates κ. It clear proves that Eq. (5) in the main text hold at different dissipation rate and decoherence rate. Particularly, in Figs. 8(a) and 8(b), the parameters of the oscillation length $\mathrm{\Delta }|\alpha {|}^2$ and optimal time space $\mathrm{\Delta }{\tau }_{\text{op}}$ in Eq. (5) are corresponding to $C=1.91\left(1.92,1.92\right)$, $a_1=3.4\left(3.4,3.4\right)\times {10}^{-3}$, $a_2=0.26,0.26,0.26$ and $a_3=1.96\left(1.96,1.96\right)\times {10}^{-3}$ for $\gamma =0.2,0.4,1$, and they are $C=2.03\left(1.92,1.91\right)$, $a_1=3.4\left(3.39,3.35\right)\times {10}^{-3}$, $a_2=0.26,0.26,0.27$ and $a_3=1.98\left(1.96,1.99\right)\times {10}^{-3}$ for $\kappa =0.5,1,2$ in Figs. 8(c) and 8(d). Therefore, the property denoted by Eq. (5) is determined by the internal dynamic of flux-qubits and independent on the flux-qubits dissipation and NVE decay.

Funding

National Key Research and Development Program of China (Grant No. 2017YFA0304503); National Natural Science Foundation of China (Grants No. 11835011, No. 11804375, No. 11734018, No. 11674360, No. 11574353, and No. 11404377); Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB21010100); China Postdoctoral Science Foundation (Grant No. 2018M642956).

References

1. E. Schrödinger, “Die gegenwärtige situation in der quantenmechanik,” Naturwissenschaften 23, 807–812 (1935). [CrossRef]  

2. C. Monroe, D. M. Meekhof, B. E. King, and D. J. Wineland, “A Schrödinger cat superposition state of an atom,” Science 272, 1131–1136 (1996). [CrossRef]   [PubMed]  

3. A. Ourjoumtsev, H. Jeong, R. Tualle-Brouri, and P. Grangier, “Generation of optical Schrödinger cat from photon number states,” Nature 448, 784–786 (2007). [CrossRef]   [PubMed]  

4. M. Arndt, O. Nairz, J. Vos-Andreae, C. Keller, G. van der Zouw, and A. Zeilinger, “Wave-particle duality of $C_{60}$ molecules,” Nature 401, 680–682 (1999). [CrossRef]  

5. D. Leibfried, E. Knill, S. Seidelin, J. Britton, R. B. Blakestad, J. Chiaverini, D. B. Hume, W. M. Itano, J. D. Jost, C. Langer, R. Ozeri, R. Reichle, and D. J. Wineland, “Creation of a six-atom Schrödinger cat state,” Nature 438, 639–642 (2005). [CrossRef]   [PubMed]  

6. W.-B. Gao, C.-Y. Lu, X.-C. Yao, P. Xu, O. Gühne, A. Goebel, Y.-A. Chen, C.-Z. Peng, Z.-B. Chen, and J.-W. Pan, “Experimental demonstration of a hyper-entangled ten-qubit Schrödinger cat state,” Nat. Phys. 6, 331–335 (2010). [CrossRef]  

7. K. G. Johnson, J. D. Wong-Campos, B. Neyenhuis, J. Mizrahi, and C. Monroe, “Ültrafast creation of large Schrödinger cat states of an atom,” Nat. Commun. 8, 697 (2017). [CrossRef]  

8. M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, “Cavity optomechanics,” Rev. Mod. Phys. 86, 1391–1452 (2014). [CrossRef]  

9. M. Blencowe, “Quantum electromechanical systems,” Phys. Rep. 395, 159–222 (2004). [CrossRef]  

10. O. Romero-Isart, A. C. Pflanzer, F. Blaser, R. Kaltenbaek, N. Kiesel, M. Aspelmeyer, and J. I. Cirac, “Large quantum superpositions and interference of massive nanometer-sized objects,” Phys. Rev. Lett. 107, 020405 (2011). [CrossRef]   [PubMed]  

11. J.-Q. Liao and L. Tian, “Macroscopic quantum superposition in cavity optomechanics,” Phys. Rev. Lett. 116, 163602 (2016). [CrossRef]   [PubMed]  

12. Y. Kubo, C. Grezes, A. Dewes, T. Umeda, J. Isoya, H. Sumiya, N. Morishita, H. Abe, S. Onoda, T. Ohshima, V. Jacques, A. Drau, J.-F. Roch, I. Diniz, A. Auffeves, D. Vion, D. Esteve, and P. Bertet, “Hybrid quantum circuit with a superconducting qubit coupled to a spin ensemble,” Phys. Rev. Lett. 107, 220501 (2011). [CrossRef]   [PubMed]  

13. T. Astner, S. Nevlacsil, N. Peterschofsky, A. Angerer, S. Rotter, S. Putz, J. Schmiedmayer, and J. Majer, “Coherent coupling of remote spin ensembles via a cavity bus,” Phys. Rev. Lett. 118, 140502 (2017). [CrossRef]   [PubMed]  

14. W. L. Yang, Z. Yin, Z. X. Chen, S. Kou, M. Feng, and C. H. Oh, “Quantum simulation of an artificial Abelian gauge field using nitrogen-vacancy-center ensembles coupled to superconducting resonators,” Phys. Rev. A 86, 012307 (2012). [CrossRef]  

15. R. Harris, F. Brito, A. J. Berkley, J. Johansson, M. W. Johnson, T. Lanting, P. Bunyk, E. Ladizinsky, B. Bumble, A. Fung, A. Kaul, A. Kleinsasser, and S. Han, “Synchronization of multiple coupled rf-SQUID flux qubits,” New J. Phys. 11, 123022 (2009). [CrossRef]  

16. R. Harris, T. Lanting, A. J. Berkley, J. Johansson, M. W. Johnson, P. Bunyk, E. Ladizinsky, N. Ladizinsky, T. Oh, and S. Han, “Compound Josephson-junction coupler for flux qubits with minimal crosstalk,” Phys. Rev. B 80, 052506 (2009). [CrossRef]  

17. D. Marcos, M. Wubs, J. M. Taylor, R. Aguado, M. D. Lukin, and A. S. Sørensen, “Coupling nitrogen-vacancy centers in diamond to superconducting flux qubits,” Phys. Rev. Lett. 105, 210501 (2010). [CrossRef]  

18. Z.-L. Xiang, X.-Y. Lü, T.-F. Li, J. Q. You, and F. Nori, “Hybrid quantum circuit consisting of a superconducting flux qubit coupled to a spin ensemble and a transmission-line resonator,” Phys. Rev. B 87, 144516 (2013). [CrossRef]  

19. W. L. Song, W. L. Yang, Z. Q. Yin, C. Y. Chen, and M. Feng, “Controllable quantum dynamics of inhomogeneous nitrogen-vacancy center ensembles coupled to superconducting resonators,” Sci. Rep. 6, 33271 (2016). [CrossRef]   [PubMed]  

20. Y. Maleki and A. M. Zheltikov, “Generating maximally-path-entangled number states in two spin ensembles coupled to a superconducting flux qubit,” Phys. Rev. A 97, 012312 (2018). [CrossRef]  

21. X. Zhu, S. Saito, A. Kemp, K. Kakuyanagi, S. Karimoto, H. Nakano, W. J. Munro, Y. Tokura, M. S. Everitt, K. Nemoto, M. Kasu, N. Mizuochi, and K. Semba, “Coherent coupling of a superconducting flux qubit to an electron spin ensemble in diamond,” Nature 478, 221–224 (2011). [CrossRef]   [PubMed]  

22. X. Zhu, Y. Matsuzaki, R. Amsüss, K. Kakuyanagi, T. S. Oka, N. Mizuochi, K. Nemoto, K. Semba, W. J. Munro, and S. Saito, “Observation of dark states in a superconductor diamond quantum hybrid system,” Nat. Commun. 5, 3424 (2014). [CrossRef]   [PubMed]  

23. S. Saito, X. Zhu, R. Amsüss, Y. Matsuzaki, K. Kakuyanagi, T. S. Oka, N. Mizuochi, K. Nemoto, W. J. Munro, and K. Semba, “Towards realizing a quantum memory for a superconducting qubit: storage and retrieval of quantum states,” Phys. Rev. Lett. 111, 107008 (2013). [CrossRef]  

24. A. Imamoglu, H. Schmidt, G. Woods, and M. Deutsch, “Strongly interacting photons in a nonlinear cavity,” Phys. Rev. Lett. 79, 1467–1470 (1997). [CrossRef]  

25. M. J. Werner and A. Imamoglu, “Photon-photon interactions in cavity electromagnetically induced transparency,” Phys. Rev. A 61, 011801 (1999). [CrossRef]  

26. M. Fleischhauer, A. Imamoglu, and J. P. Marangos, “Electromagnetically induced transparency: optics in coherent media,” Rev. Mod. Phys. 77, 633–673 (2005). [CrossRef]  

27. S Rebić, S M Tan, A S Parkins, and D F Walls, “Large Kerr nonlinearity with a single atom,” J. Opt. B 1, 490–495 (1999). [CrossRef]  

28. Y. Maleki and A. M. Zheltikov, “Witnessing quantum entanglement in ensembles of nitrogen-vacancy centers coupled to a superconducting resonator,” Opt. Express 14, 17849–17858 (2018). [CrossRef]  

29. S. Rebic, J. Twamley, and G. J. Milburn, “Giant Kerr nonlinearities in circuit quantum electrodynamics,” Phys. Rev. Lett. 103, 150503 (2009). [CrossRef]   [PubMed]  

30. F. Jelezko, T. Gaebel, I. Popa, A. Gruber, and J. Wrachtrup, “Observation of coherent oscillations in a single electron spin,” Phys. Rev. Lett. 92, 076401 (2004). [CrossRef]   [PubMed]  

31. G. Chen, H. Zhang, Y. Yang, R. Wang, L. Xiao, and S. Jia, “Qubit-induced high-order nonlinear interaction of the polar molecules in a stripline cavity,” Phys. Rev. A 82, 013601 (2010). [CrossRef]  

32. H. Schmidt and A. Imamoglu, “Giant Kerr nonlinearities obtained by electromagnetically induced transparency,” Opt. Lett. 21, 1936–1938 (1996). [CrossRef]   [PubMed]  

33. P. D. Drummond and D. F. Walls, “Quantum theory of optical bistability. I. nonlinear polarisability model,” J. Phys. A 13, 725–741 (1980). [CrossRef]  

34. M. D. Reid and D. F. Walls, “Squeezing via optical bistability,” Phys. Rev. A 32, 396–401 (1985). [CrossRef]  

35. M. J. Collett and D. F. Walls, “Squeezing spectra for nonlinear optical systems,” Phys. Rev. A 32, 2887–2892 (1985). [CrossRef]  

36. G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg, H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, “Observation of quantum state collapse and revival due to the single-photon Kerr effect,” Nature 495, 205–209 (2013). [CrossRef]   [PubMed]  

37. D. F. Walls and G. J. Milburn, Quantum optics(Springer, 1994). [CrossRef]  

38. M. Kira, S. W. Koch, R. P. Smith, A. E. Hunter, and S. T. Cundiff, “Quantum spectroscopy with Schrödinger-cat states,” Nat. Phys. 7, 799–804 (2011). [CrossRef]  

39. A. Facon, E.-K. Dietsche, D. Grosso, S. Haroche, J.-M. Raimond, M. Brune, and S. Gleyzes, “A sensitive electrometer based on a Rydberg atom in a Schrödinger-cat state,” Nature 535, 262–265 (2016). [CrossRef]   [PubMed]  

40. J. Joo, W. J. Munro, and T. P. Spiller, “Quantum metrology with entangled coherent states,” Phys. Rev. Lett. 107, 083601 (2011). [CrossRef]   [PubMed]  

41. J. S. N. Nielsen, Y. Eto, C. W. Lee, H. Jeong, and M. Sasaki, “Quantum tele-amplification with a continuous-variable superposition state,” Nat. Photon. 7, 439–443 (2013). [CrossRef]  

42. J. Huang, X. Qin, H. Zhong, Y. Ke, and C. Lee, “Quantum metrology with spin cat states under dissipation,” Sci. Rep. 5, 17894 (2015). [CrossRef]   [PubMed]  

43. T. Unden, P. Balasubramanian, D. Louzon, Y. Vinkler, M. B. Plenio, M. Markham, D. Twitchen, A. Stacey, I. Lovchinsky, A. O. Sushkov, M. D. Lukin, A. Retzker, B. Naydenov, L. P. McGuinness, and F. Jelezko, “Quantum metrology enhanced by repetitive quantum error correction,” Phys. Rev. Lett. 116, 230502 (2016). [CrossRef]   [PubMed]  

44. A. Z. Chaudhry, “Ütilizing nitrogen-vacancy centers to measure oscillating magnetic fields,” Phys. Rev. A 90, 042104 (2014). [CrossRef]  

45. T. Wolf, P. Neumann, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, and J. Wrachtrup, “Subpicotesla diamond magnetometry,” Phys. Rev. X 5, 041001 (2015).

46. W. L. Ma, S. S. Li, G. Y. Cao, and R. B. Liu, “Atomic-scale positioning of single spins via multiple nitrogen-vacancy centers,” Phys. Rev. Applied 5, 044016 (2016). [CrossRef]  

47. H. Liu, M. B. Plenio, and J. Cai, “Scheme for detection of single-molecule radical pair reaction using spin in diamond,” Phys. Rev. Lett. 118, 200402 (2017). [CrossRef]   [PubMed]  

48. M. Bruderer, P. F. Acebal, R. Aurich, and M. B. Plenio, “Sensing of single nuclear spins in random thermal motion with proximate nitrogen-vacancy centers,” Phys. Rev. B 93, 115116 (2016). [CrossRef]  

49. D. Leibfried, D. M. Meekhof, B. E. King, C. Monroe, W. M. Itano, and D. J. Wineland, “Experimental determination of the motional quantum state of a trapped atom,” Phys. Rev. Lett. 77, 4281 (1996). [CrossRef]  

50. M. Hofheinz, H. Wang, M. Ansmann, R. C. Bialczak, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, J. Wenner, J. M. Martinis, and A. N. Cleland, “Synthesizing arbitrary quantum states in a superconducting resonator,” Nature 459, 546–549 (2009). [CrossRef]   [PubMed]  

51. B. R. Johnson, M. D. Reed, A. A. Houck, D. I. Schuster, L. S. Bishop, E. Ginossar, J. M. Gambetta, L. DiCarlo, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, “Quantum non-demolition detection of single microwave photons in a circuit,” Nat. Phys. 6, 663 (2010). [CrossRef]  

52. E. Solano, “Selective interactions in trapped ions: state reconstruction and quantum logic,” Phys. Rev. A 71, 013813 (2005). [CrossRef]  

53. S. Haroche and J. M. Raimond, Exploring the quantum: atoms, cavities, and photons(Oxford university, 2006). [CrossRef]  

54. M. Scala, M. S. Kim, G. W. Morley, P. F. Barker, and S. Bos, “Matter-wave interferometry of a levitated thermal nano-oscillator induced and probed by a spin,” Phys. Rev. Lett. 111, 180403 (2013). [CrossRef]   [PubMed]  

55. L. P. Neukirch, E. von Haartman, J. M. Rosenholm, and A. N. Vamivakas, “Multi-dimensional single-spin nano-optomechanics with a levitated nanodiamond,” Nat. Photo. 9, 653–657 (2015). [CrossRef]  

56. J. E. Mooij, T. P. Orlando, L. Levitov, L. Tian, C. H. van der Wal, and S. Lloyd, “Josephson persistent-current qubit,” Science 285, 1036–1039 (1999). [CrossRef]   [PubMed]  

57. R. Harris, M. W. Johnson, S. Han, A. J. Berkley, J. Johansson, P. Bunyk, E. Ladizinsky, S. Govorkov, M. C. Thom, S. Uchaikin, B. Bumble, A. Fung, A. Kaul, A. Kleinsasser, M. H. S. Amin, and D. V. Averin, “Probing noise in flux qubits via macroscopic resonant tunneling,” Phys. Rev. Lett. 101, 117003 (2008). [CrossRef]   [PubMed]  

58. M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, “Strongly interacting polaritons in coupled arrays of cavities,” Nat. Phys. 2, 849–855 (2006). [CrossRef]  

59. K. Cai, Z.-W. Pan, R.-X. Wang, D. Ruan, Z.-Q. Yin, and G.-L. Long, “Single phonon source based on a giant polariton nonlinear effect,” Opt. Lett. 43, 1163–1166 (2018). [CrossRef]   [PubMed]