Enhanced phonon blockade in a weakly-coupled hybrid system via mechanical parametric amplification

We consider a hybrid spin-mechanical setup, as schematically illustrated in Fig. 1, where a single NV center is embedded near the freely vibrating end of a singly-clamped diamond cantilever of dimensions ($\ell,w,t$). Two cylindrical nanomagnets are symmetrically arranged on the two sides of the cantilever to produce strong magnetic field gradient at short distances Ma et al. (2016); Cai et al. (2017); Sánchez Muñoz et al. (2018); Yin et al. (2019). The electronic spin of the NV center can be coupled to the cantilever resonator based on the fact that the bending motion of the cantilever modulates the local magnetic field sensed by the NV spin Arcizet et al. (2011); Kepesidis et al. (2016). The cantilever is electrically pumped by a periodic drive (under the cantilever) which modulates the mechanical spring constant in time Rugar and Grütter (1991); Lemonde et al. (2016); Li et al. (2020a), generating the process of MPA. In addition, a weak mechanical drive is applied to the cantilever for inducing the PB effect.

Following the specific quantization process as in Refs. Lemonde et al. (2016); Li et al. (2020a), the Hamiltonian of the cantilever oscillator takes the form

where $\hat{a}$ ($\hat{a}^{\dagger}$) is the annihilation (creation) operator for the fundamental vibrational mode of frequency $\omega_{m}$, $\omega_{p}$ and $\Omega_{p}$ are the pump frequency and amplitude, $\omega_{L}$ and $\varepsilon_{L}$ are the frequency and strength of the linear mechanical drive.

The electronic ground state of the single NV center we studied is a $S=1$ spin triplet with basis states $|m_{s}\rangle$, where $m_{s}=0,\pm 1$. The ground-state energy-level structure of this NV center is shown in the inset of Fig. 1. The zero-field splitting between the degenerate sublevels $|m_{s}=\pm 1\rangle$ and $|m_{s}=0\rangle$ is $D=2\pi\times 2.87$ GHz. The strong magnetic field B generated by the nanomagnets at the position of the NV center (see below) causes a Zeeman splitting much larger than the zero-field splitting $D$, which can be canceled by introducing a static magnetic field $\textbf{B}_{\text{static}}$ (in the opposite direction to B) Rabl et al. (2009). The magnetic fields of these two parts form a superposed field of magnitude $B_{s}$ such that $\mu_{B}B_{s}\ll D$, thus giving a Zeeman splitting $\delta_{z}=2g_{e}\mu_{B}B_{s}$, where $g_{e}\simeq 2$ is the Land$\acute{\mathrm{e}}$ g-factor of NV centers and $\mu_{B}=14\ \mathrm{GHz}/\mathrm{T}$ the Bohr magneton. On this basis, two microwave (MW) fields polarized in the $x$ direction, $B_{x}^{\pm}(t)=B_{0}^{\pm}\cos\omega_{\pm}t$, are introduced to induce the transitions between states $|\pm 1\rangle$ and $|0\rangle$, whose Rabi frequencies $\Omega_{\pm}\equiv g_{e}\mu_{B}B_{0}^{\pm}$ and (red) detunings $\Delta_{\pm}\equiv D\pm\delta_{z}/2-\omega_{\pm}$ are supposed to be identical in what follows for simplicity, i.e., $\Omega_{\pm}=\Omega$ and $\Delta_{\pm}=\Delta$. In a rotating frame with respect to the MW frequencies $\omega_{\pm}$ and under the rotating-wave approximation (RWA), the Hamiltonian of the single NV center is written as Rabl et al. (2009)

Next, we turn our attention to the spin-phonon interaction. In principle, the single NV spin surrounded by a magnetic field $\textbf{B}(\textbf{r})$ couples to the mechanical oscillator via the position-dependent Zeeman shift. The corresponding magnetic interaction Hamiltonian is generally written as $\hat{H}_{\mathrm{int}}=\mu_{B}g_{e}\hat{\textbf{S}}\cdot\textbf{B}(\textbf{r}_{0})$, where $\hat{\textbf{S}}$ denotes the spin operator and $\textbf{r}_{0}$ the position of the NV center. Assume that the cantilever resonator (also the NV spin) oscillates only along the $z$ direction, which allows us to expand the Hamiltonian $\hat{H}_{\mathrm{int}}$ up to second order in terms of the cantilever displacement $z$, resulting in $\hat{H}_{\mathrm{int}}\simeq\mu_{B}g_{e}\hat{\textbf{S}}\cdot\left[\partial\mathbf{B}/\partial z(0)\hat{z}+\frac{1}{2}\partial^{2}\mathbf{B}/\partial z^{2}(0)\hat{z}^{2}\right]$. Due to the specific arrangement of nanomagnets, the generated magnetic field yields an extremum at the position of the NV center, thus nullifying the first derivative $\partial\mathbf{B}/\partial z(0)$, that is, the first-order magnetic gradient is equal to zero. This can be verified in Fig. 2, where we show in Fig. 2(a) the magnetic field distribution in the $xz$ plane generated by two symmetrically-placed nanomagnets, and the magnetic strength along the $z$ axis for $x=0$ is given in Fig. 2(b). Here, the two cylindrical magnets are of the same size, with a diameter of 30 nm and height of 40 nm, and the distance between them is adjustable, ranging from tens to hundreds of nanometers; The magnetic substance is selected as Dy due to a high saturation magnetization up to $\mu_{0}M_{\mathrm{s}}=3.7$ T Sánchez Muñoz et al. (2018). For the diamond cantilever, its dimension is set as ($\ell,w,t$)=(4,0.1,0.02) $\mu$m in our simulations. Note that the components of different dimensions as simulated here could be fabricated individually with modern nanofabrication techniques; See Appendix A for a brief discussion regarding experimental realizations of device fabrication and arrangement. As seen in Fig. 2(b), the magnetic strength displays an extremum (precisely, a local minimum) at $z=0$ where the NV center is located. Also, it shows that the magnetic strength between the nanomagnets decreases with the increase of the magnet separation $D_{z}$.

The resulting interaction Hamiltonian $\hat{H}_{\mathrm{int}}$, after eliminating the first-order magnetic gradient, is given by $\hat{H}_{\mathrm{int}}\simeq\frac{1}{2}\mu_{B}g_{e}G\hat{z}^{2}\hat{S}_{z}$, where $G=\partial^{2}B_{z}/\partial z^{2}(0)$ represents the second-order magnetic gradient. Note here that the mechanical oscillator only couples to the $z$-component of the NV spin due to null second derivatives of the magnetic field along the $x$ and $y$ directions. Expressing the position operator $\hat{z}$ with $\hat{z}=z_{\mathrm{zpf}}\left(\hat{a}^{\dagger}+\hat{a}\right)$, we end up with

where $g=\frac{1}{2}\mu_{B}g_{e}z_{\mathrm{zpf}}^{2}G$ is the two-phonon coupling rate. Since $g$ is proportional to the square of the zero-point motion of the oscillator $z_{\mathrm{zpf}}$, which ranges from tens to hundreds of femtometers for various systems Ghadimi et al. (2018); Norte et al. (2016); Weber et al. (2016); Singh et al. (2014); Arcizet et al. (2011), the quadratic coupling strength is intrinsically much smaller than the single-phonon coupling strength ($g_{0}$) that is proportional to the zero-point fluctuation ($g_{0}\propto z_{\mathrm{zpf}}$) Rabl et al. (2009, 2010); Li et al. (2020a); Zhou et al. (2018); Li et al. (2015, 2016). Also, the magnitude of the second-order magnetic gradient determines linearly the quadratic coupling strength, which is controllable over a large range by adjusting the magnets separation (see below).

To sum up, we obtain the total Hamiltonian for the studied hybrid system

To diagonalize the spin-only part (i.e., $\hat{H}_{\mathrm{NV}}$), we then transform $\hat{H}_{\mathrm{Total}}$ to a dressed-state basis spanned by $\{|\mathcal{D}\rangle=(|+1\rangle-|-1\rangle)/\sqrt{2},|\mathcal{G}\rangle=\cos\theta|0\rangle-\sin\theta|\mathcal{B}\rangle,|\mathcal{E}\rangle=\cos\theta|\mathcal{B}\rangle+\sin\theta|0\rangle\}$, with $|\mathcal{B}\rangle=(|+1\rangle+|-1\rangle)/\sqrt{2}$ and $\tan(2\theta)=\sqrt{2}\Omega/\Delta$. The resulting system Hamiltonian in a rotating frame with respect to $\omega_{p}$, under the RWA, is simplified as (see Appendix C for a detailed derivation)

where $\delta_{m(L)}=\omega_{m(L)}-\omega_{p}$, $\delta_{ed}=\omega_{ed}-2\omega_{p}$, and $\hat{\sigma}_{+}=\hat{\sigma}_{-}^{\dagger}\equiv|\mathcal{E}\rangle\langle\mathcal{D}|$. Notably, Eq. (5) has the form of a two-phonon Jaynes-Cummings (JC) Hamiltonian in the presence of both linear and nonlinear (two-phonon) mechanical drive, which suggests degenerate two-phonon exchange between the mechanical mode and an effective TLS.

So far, we have derived the driven, two-phonon JC Hamiltonian for the hybrid spin-mechanical system. In this subsection, we proceed to estimate the experimentally achievable two-phonon coupling rate as well as other system parameters. Based on the numerical simulations in Figs. 2(a) and 2(b), it is straightforward to acquire the values of the second-order magnetic gradient $G=\partial^{2}B_{z}/\partial z^{2}(0)$ for different magnet separations. Then, the two-phonon coupling $g$ depends directly on the value of the zero-point fluctuation of the oscillator $z_{\mathrm{zpf}}=\sqrt{\hbar/2M\omega_{m}}$, according to $g=\frac{1}{2}\mu_{B}g_{e}z_{\mathrm{zpf}}^{2}G$. The fundamental frequency $\omega_{m}$ and effective mass $M$ of the oscillator are estimated as $\omega_{m}\approx 3.516\times\left(t/\ell^{2}\right)\sqrt{E/12\rho}\approx 2\pi\times 3.8$ MHz and $M=\rho\ell wt/4\approx 7\times 10^{-18}$ kg, with the Young’s modulus and mass density of diamond being $E\approx 1.22\times 10^{12}\mathrm{~{}Pa}$ and $\rho\approx 3.52\times 10^{3}\mathrm{~{}kg/m^{3}}$, respectively, which gives $z_{\mathrm{zpf}}\approx 563$ fm. In Fig. 2(c) we plot the calculated two-phonon coupling rate $g$ versus the magnet separation $D_{z}$ ranging from 40 to 100 nm. It shows that a small magnet separation is desirable for achieving a large coupling rate. Despite this, we will take hereafter a relatively large value of $D_{z}$ from the perspective that a larger magnet separation would be more beneficial to actual operation in experiments. As an example, we take $D_{z}=80$ nm, giving $G\approx 7.9\times 10^{14}\mathrm{~{}T/m^{2}}$ and therefore $g\approx 2\pi\times 1.1$ Hz. A small misalignment in the field orientation due to imperfectly positioned magnets could affect slightly the second-order magnetic gradient and therefore the two-phonon coupling rate, leading to a relative error of less than $2\%$ in the presence of a $10^{\circ}$ misalignment (see Appendix B for details).

The time-dependent pump on the mechanical oscillator introduces nonlinear drive to the mechanical mode, producing the effect of degenerate parametric amplification. Specifically, by performing a unitary squeezing transformation $\hat{U}_{s}=\exp\left[r_{p}\left(\hat{a}^{2}-\hat{a}^{\dagger 2}\right)/2\right]$ Leroux et al. (2018); Wang et al. (2020b, c); Chen et al. (2021), where the squeezing parameter $r_{p}$ is defined via $r_{p}=(1/2)\operatorname{arctanh}\left(\Omega_{p}/\delta_{m}\right)$, the system Hamiltonian (5) can be transformed to the squeezed frame, yielding (see Appendix D)

where $\delta_{s}=\delta_{m}/\cosh(2r_{p})$ is the squeezed-oscillator frequency, $g_{\mathrm{eff}}=g\cosh^{2}(r_{p})$ and $\varepsilon_{L}^{\prime}=\varepsilon_{L}\cosh(r_{p})$ are the effective two-phonon spin-mechanical coupling strength and linear mechanical driving strength, respectively. The exponential coupling enhancement with respect to the original coupling strength, $g_{\mathrm{eff}}/g=\cosh^{2}(r_{p})$, is visualized in Fig. 3, where one sees that $g_{\mathrm{eff}}$ can be two orders of magnitude larger than $g$ for a modest squeezing parameter $r_{p}=3$. It should be mentioned that the scheme proposed here is distinct from the exponential enhancement associated with the single-phonon (photon) spin-resonator interaction using mechanical (optical) parametric amplification Li et al. (2020a); Lü et al. (2015); Qin et al. (2018); Leroux et al. (2018); Wang et al. (2019a, b); Chen et al. (2019); Qin et al. (2019); Wang et al. (2020b, c, d); Chen et al. (2021).

Utilizing MPA to enhance the spin-mechanical coupling also amplifies the mechanical noise inevitably, which could break down any nonclassical behavior in the system due to the amplified dissipation. An effective strategy against this adverse effect is to use the squeezed-vacuum-reservoir technique, which has been studied extensively in recent reports Li et al. (2020a); Lü et al. (2015); Qin et al. (2018); Leroux et al. (2018); Wang et al. (2019a, b); Chen et al. (2019); Qin et al. (2019); Wang et al. (2020b, c); Chen et al. (2021). To be specific, by integrating an auxiliary, broadband microwave resonator to the original hybrid system, which acts as an engineered squeezed-vacuum reservoir, the squeezing-enhanced mechanical noise could be effectively suppressed, ensuring that the squeezed mechanical mode equivalently interacts with a thermal vacuum reservoir (see Appendix E for details). The dynamics of the system is thus governed by the standard Lindblad master equation which in the squeezed frame takes the form Li et al. (2020a)

where $\mathcal{D}[\hat{o}]\hat{\rho}=\hat{o}\hat{\rho}\hat{o}^{\dagger}-\hat{o}^{\dagger}\hat{o}\hat{\rho}/2-\hat{\rho}\hat{o}^{\dagger}\hat{o}/2$, $\gamma_{m_{\mathrm{eff}}}$ is the engineered effective mechanical dissipation rate as a consequence of the coupling between the mechanical mode and the auxiliary reservoir, and $\gamma_{z}$ is the dephasing rate of the spin qubit. The quality factor of the oscillator $Q$ is associated with $\gamma_{m_{\mathrm{eff}}}$ via the relation $Q=n_{\mathrm{th}}\omega_{m}/\gamma_{m_{\mathrm{eff}}}$, with $n_{\mathrm{th}}=\left[\exp\left(\hbar\omega_{m}/k_{B}T\right)-1\right]^{-1}$ being the average number of thermal phonons in the oscillator at the temperature $T$. In the present work the environmental temperature is assumed to be about tens of millikelvin Tsaturyan et al. (2017), resulting in tens of mean thermal phonon number for the mechanical frequency of tens of MHz. When it comes to the NV center, its dephasing time, given by $T_{2}=1/\gamma_{z}$, is assumed to be on the order of milliseconds in our calculations. While the spin properties of shallow NV centers are not as good as those for NV centers in bulk diamond Bar-Gill et al. (2013); De Lange et al. (2010); Kolkowitz et al. (2012), many experimental efforts have been devoted to improve the coherence properties of shallow NV centers Mamin et al. (2013); Ishikawa et al. (2012); Ohashi et al. (2013); Ohno et al. (2012); Myers et al. (2014, 2014), which make the extension of coherence time possible with the development and exploration of future experimental techniques. For clarity, we also evaluate the effect of spin dephasing on the PB performance by simulating both steady-state second-order correlation function and mean phonon number with varying spin dephasing rate over a large range (see below).

The master equation (7) gives an effective cooperativity $C^{\prime}=g_{\mathrm{eff}}^{2}/(\gamma_{m_{\mathrm{eff}}}\gamma_{z})$, which can be strengthened significantly with increasing the squeezing parameter $r_{p}$, as shown in Fig. 3. Here, the cooperativity enhancement is approximated as $C^{\prime}/C\approx\cosh^{4}(r_{p})$. This result is justified for the case where the effective mechanical dissipation rate is almost unchanged after introducing MPA. The enhancement in the effective cooperativity provides great potentials for various applications in quantum information technologies. In the following, we explore the possibilities of using the designed device to induce strong PB and investigate the statistical characteristics of phonons. As we will show, both analytically and numerically, the enhanced cooperativity contributes essentially to the observation of strong PB in the hybrid system that works in the weak-coupling regime.

We start by introducing two criteria for PB. The first criterion is based on a comparison of the phonon-number distribution with standard Poissonian distribution. Specifically, when $n$-PB occurs, the phonon-number distribution $P(m)$ (with normalization $\sum_{m=0}^{\infty}P(m)=1$) satisfies Hamsen et al. (2017); Huang et al. (2018)

where $\mathcal{P}(m)=\langle\hat{m}\rangle^{m}e^{-\langle\hat{m}\rangle}/m!$ is the Poissonian distribution with the same average phonon number $\langle\hat{m}\rangle$ as the mechanical oscillator. This indicates a sub-Poissonian phonon-number statistics for $(n+1)$ phonons with the simultaneous super-Poissonian statistics of the first $n$ phonons. Also, the relative deviation of a given phonon-number distribution from the corresponding Poissonian distribution is given by $[P(m)-\mathcal{P}(m)]/\mathcal{P}(m)$.

The second criterion is based on calculating the phonon correlation functions that characterize the statistical properties of phonons. The normalized equal-time $\mu$th-order correlation function is defined as

which is related to the probability of simultaneously measuring $\mu$ phonons. In the steady state, the equal-time $\mu$th-order correlation function becomes $g^{(\mu)}(0)_{ss}=\text{lim}_{t\rightarrow\infty}g^{(\mu)}(t,0)$. Note that the larger value of $g^{(\mu)}(0)_{ss}>1$, the higher probability of $\mu$-phonon bunching; The smaller value of $g^{(\mu)}(0)_{ss}<1$, the higher probability of $\mu$-phonon antibunching. Particularly, the steady-state, equal-time second-order correlation function is given by

In the case of a finite-time delay between two phonon detection, one can apply the steady-state, delayed-time second-order correlation function

The criterion for $n$-phonon blockade with respect to the phonon-number distribution, given in Eq. (8), can be translated into the following conditions Hamsen et al. (2017); Huang et al. (2018)

For single PB with $n=1$, we obtain

For the system with a weak linear driving of the mechanical resonator, the mean phonon number is very small, i.e., $\langle\hat{m}\rangle\ll 1$, which approximatively gives $f\rightarrow 1$ and $f^{(1)}\rightarrow 1$. Then the condition simplifies to the usual criterion of single PB, $g^{(2)}(0)_{ss}<1$. The two criteria, given respectively in Eqs. (8) and (13)-(14), will facilitate us to identify the occurrence of single PB in the studied system.

To analyze the steady state of our system and obtain the analytical expression of the second-order correlation function, we introduce an effective non-Hermitian Hamiltonian involving the system dissipation

where

is a reformulation of Eq. (6) after assuming $\delta_{ed}=2\delta_{s}$ and $\delta_{s}-\delta_{L}=\delta$ with $\delta$ being a small laser-drive detuning. In the weak-driving regime $(\varepsilon_{L}^{\prime}\ll g_{\mathrm{eff}})$, by truncating the infinite-dimensional Hilbert space to the two-phonon excitation subspace, the wave function of the system could be expressed with the ansatz

with the coefficients being probability amplitudes of the basis states. Substituting the state $|\psi\rangle_{\mathrm{s}}$ and the Hamiltonian $\hat{H}_{\mathrm{nH}}^{\mathrm{S}}$ to the Schr$\ddot{\mathrm{o}}$dinger equation $i|\dot{\psi}\rangle_{\mathrm{s}}=\hat{H}_{\mathrm{nH}}^{\mathrm{S}}|\psi\rangle_{\mathrm{s}}$, we obtain the following equations of motion for the probability amplitudes

By solving the above equations, the steady-state solutions of the subsystem are readily achieved (as given detailedly in Appendix F). Then, we obtain the analytical expression of the equal-time second-order correlation function

with $\zeta=(\gamma_{\mathrm{z}}/2)^{2}+\gamma_{m_{\mathrm{{eff}}}}^{2}-4g_{\mathrm{eff}}^{2}$. As seen in Eq. (19), the value of $g^{(2)}(0)_{ss}$ depends directly on $\delta$: for $\delta=0$, $g^{(2)}(0)_{ss}$ reaches its minimum

This implies that strong single PB (or, two-phonon antibunching effect) could occur when the mechanical mode is resonantly driven. Further, for $\varepsilon_{L}^{\prime}\ll\gamma_{m_{\mathrm{{eff}}}}$, by expressing $g_{\mathrm{eff}}^{2}/(\gamma_{m_{\mathrm{{eff}}}}\gamma_{\mathrm{z}})$ as $C^{\prime}$ in Eq. (20), we have

This explicitly shows that the remarkably enhanced cooperativity in our scheme could make the quantum effect of PB as strong as possible.

To confirm the analytical results, we now numerically study the full quantum dynamics of the system. The numerical simulations were performed by solving the master equation (7) using the Quantum Toolbox in Python (QuTiP) Johansson et al. (2013, 2012). We consider the system with weak bare coupling: $\gamma_{m_{\mathrm{{eff}}}}=2g$ Li et al. (2016). Fig. 5(a) depicts the first- and second-order correlation functions versus the driving detuning $\delta/g_{\mathrm{eff}}$. For $\delta/g_{\mathrm{eff}}$ fluctuating within a small range, the value of $g^{(2)}(0)_{ss}$ keeps below 1, corresponding to the region where PB occurs. In particular, it is shown that $g^{(2)}(0)_{ss}$ has a dip at $\delta=0$ (i.e., a resonant drive), which is consistent with the analytical description above. We also see that, at $\delta=0$, $g^{(2)}(0)_{ss}$ is smaller than $f$ defined in the criterion given in Eq. (13), while $g^{(1)}(0)_{ss}$ is greater than $f^{(1)}$ defined in the criterion given in Eq. (14), as shown in Figs. 5(b) and 5(c), respectively. Moreover, by comparing the phonon-number distribution $P(m)$, calculated numerically via $P(m)=\left\langle m\left|\hat{\rho}_{ss}\right|m\right\rangle$ with $\hat{\rho}_{ss}$ the steady-state solutions of the master equation, with the Poisson distribution $\mathcal{P}(m)$, we find that the single-phonon probability is enhanced as $P(1)>\mathcal{P}(1)$, while the two-phonon probability is suppressed as $P(2)<\mathcal{P}(2)$, as depicted in Figs. 5(d) and 5(e), respectively. Fig 5(f) shows the deviations of the phonon distribution to the standard Poisson distribution with the same mean phonon number when $\delta=0$, which characterizes the second-order sub-Poissonian statistics (i.e., two-phonon antibunching). The above results provide sufficient and unambiguous evidence for single PB.

For a weak mechanical drive $\varepsilon_{L}^{\prime}$, the system is well truncated into a two-phonon excitation subspace spanned by four basis states given in Eq. (17). However, this approximation becomes invalid when higher levels are excited in the case of a relatively strong drive. Thus, in order to estimate the effect of state truncation on the quality of PB, we introduce a fidelity describing the sum of the steady-state probabilities of the four basis states Miranowicz et al. (2013); Yin et al. (2019); Wang et al. (2016), given by

In Fig. 5(g) we plot the fidelity $F\left(\rho_{ss}\right)$ versus the driving strength $\varepsilon_{L}^{\prime}/\gamma_{m_{\mathrm{{eff}}}}$. For $\varepsilon_{L}^{\prime}/\gamma_{m_{\mathrm{{eff}}}}<0.5$, the value of $F\left(\rho_{ss}\right)$ keeps close to 1, whereas further increasing the driving strength leads to a rapid decrease of $F\left(\rho_{ss}\right)$. In this case, the two-phonon (or even multi-phonon) state will be populated slightly, as seen in Fig. 5(h), and correspondingly, the effect of PB will be spoiled, resulting in an increased $g^{(2)}(0)_{ss}$ (see the dashed curve in Fig. 5(g)). Thus, we conclude that the occurrence of strong PB requires the mechanical driving strength to be as small as possible to suppress the two-phonon excitation. Nonetheless, a large single-phonon probability determining the efficiency of single-phonon emission is also crucial for practical applications such as the single-phonon source.

We now proceed to show how to enhance the PB emerging in the weakly coupled system by modulating MPA. Specifically, we plot in Fig. 6(a) the logarithmic second-order correlation function $g^{(2)}(0)_{ss}$ versus the squeezing parameter $r_{p}$ as well as the $Q$ factor of the mechanical oscillator. As seen, for a fixed value of $Q$, the phonon antibunching effect can be gradually strengthened with increasing $r_{p}$. In particular, when $r_{p}\rightarrow 3$, $g^{(2)}(0)_{ss}<10^{-2}$ is achievable even for $Q$ as low as $10^{7}$ (corresponding to $\gamma_{m_{\mathrm{{eff}}}}/g\approx 20$). Note that here, in the region where strong PB occurs, the single-phonon probability can not be very high due to a low ratio $\varepsilon_{L}^{\prime}/\gamma_{m_{\mathrm{{eff}}}}$. This is illustrated in Fig. 6(b), where we see that the single-phonon probability $P(1)$ is lower than 0.1 for $r_{p}=3$ and $Q=10^{7}$, which can be increased to exceed 0.1 with increasing $Q$. In addition, one sees in Fig. 6(a) that for larger $Q$, e.g., $Q=10^{8}$, the reduction of $g^{(2)}(0)_{ss}$ with increasing $r_{p}\to 3$ slows down gradually, in contrast to the case with smaller $Q$. This is because for $Q=10^{8}$, the effective driving strength $\varepsilon_{L}^{\prime}$ enhanced by the strong MPA exceeds the mechanical decay rate $\gamma_{m_{\mathrm{{eff}}}}$, thus weakening the effect of PB in accordance with our previous analysis. In this case (with a relatively high-$Q$ mechanical oscillator), it is desirable to apply a weaker bare drive $\varepsilon_{L}$ to ensure strong PB being triggered. To this end, we plot in Fig. 6(c) the logarithmic correlation $g^{(2)}(0)_{ss}$ as functions of the mechanical $Q$ factor and the driving strength $\varepsilon_{L}$ for a fixed squeezing parameter $r_{p}=3$. One clearly sees that strong-PB region with $g^{(2)}(0)_{ss}<10^{-3}$ emerges around $Q=10^{8}$ when applying a weak drive with strength $\varepsilon_{L}<0.1g$. Correspondingly, the single-phonon probability $P(1)>0.1$ can be achieved in the region of $g^{(2)}(0)_{ss}<10^{-3}$, as we exemplify with point A in Figs. 6(c) and 6(d). Also, point B in both plots exemplifies a typical case where both strong PB and a large single-phonon probability could coexist in our weak-coupling system with a low-$Q$ oscillator. Note that the performance of PB in the general hybrid system with single-phonon spin-mechanical coupling is relatively weak due to the limitation on the strength of MPA (see Appendix G for a detailed discussion). Thus, our two-phonon-based scheme may provide an attractive platform to access PB with better performance and large-scale tunability.

It is also interesting to examine the effect of the adjustable MPA on the phonon statistics when the system is originally in the strong-coupling regime. As shown in Fig. 6(a), we observe a small region of phonon bunching (or, super-Poissonian phonon statistics) characterized by $g^{(2)}(0)_{ss}>1$, which appears for $\gamma_{m_{\mathrm{{eff}}}}<g$ and a small $r_{p}$. By gradually enhancing the mechanical amplification, phonon antibunching with $g^{(2)}(0)_{ss}<1$ emerges, and finally strong PB with $g^{(2)}(0)_{ss}<0.1$ is achieved. This reveals that our approach provides a flexible tunability of the phonon statistics, and thus could enable a number of applications in the quantum control of phonon as well as phonon-based quantum networks.

In the above analysis we have taken a fixed spin dephasing rate $\gamma_{z}=2\pi\times 10$ Hz, corresponding to a dephasing time $T_{2}\approx 16$ ms. This actually gives a future demonstration based on an assumption that the coherence properties of the shallow NV center we are considering can be greatly improved. To further evaluate the effect of varying $\gamma_{z}$ on the PB and also show the PB performance within the extent of current technology, we plot in Fig. 7(a) the steady-state second-order correlation function $g^{(2)}(0)_{ss}$ and mean phonon number $\langle\hat{a}^{\dagger}\hat{a}\rangle_{ss}$ versus $\gamma_{z}/(2\pi)\in\left[1,1500\right]$ Hz. The hollow markers locate $\gamma_{z}=2\pi\times 10$ Hz ($T_{2}\approx 16$ ms), which reproduce the results obtained by point A in Figs. 6(c) and 6(d). With increasing $\gamma_{z}$ the value of $g^{(2)}(0)_{ss}$ increases sharply at first and then slows down, while the value of $\langle\hat{a}^{\dagger}\hat{a}\rangle_{ss}$ is nearly independent of varying $\gamma_{z}$. For $\gamma_{z}/(2\pi)=1$ kHz or $T_{2}^{\ast}\approx 160$ $\mu$s, which are parameters accessible to current technology Ishikawa et al. (2012), PB with $g^{(2)}(0)_{ss}<10^{-1}$ and $\langle\hat{a}^{\dagger}\hat{a}\rangle_{ss}>0.1$ can be observed (see solid markers in Fig. 7(a)).

Fig. 7(b) shows the effects of thermal phonons on the performance of the enhanced PB. It is clear that increasing the mean thermal phonon number $n_{\mathrm{th}}$ does not destroy the PB effect significantly, as $g^{(2)}(0)_{ss}$ can still stay far below $10^{-2}$ when $n_{\mathrm{th}}$ increases to 200. By contrast, a large $n_{\mathrm{th}}$ has a greater effect on the mean phonon number in the mechanical oscillator, reducing $\langle\hat{a}^{\dagger}\hat{a}\rangle_{ss}$ to be about $0.013$ for $n_{\mathrm{th}}=200$. This suggests that the present protocol could yield better performance by operating at cryogenic temperatures.