Engineering special nonclassical quantum states are of paramount importance in quantum information science, metrology, and exploring fundamental physics [1–5]. In particular, $n$-photon states play an essential role in a wide range of quantum technologies, including high-NOON states [6], quantum communication [7], lithography [8], spectroscopy [9,10], and biological sensing [11,12]. Methods for generating $n$-photon states were proposed theoretically in cavity quantum electrodynamics (QED) [13–15], Rydberg atomic ensembles [16,17], and atom-coupled photonic waveguides [18–20]. Among them, $n$-photon bundle states generated through either Mollow physics [13,14] or deterministic parametric downconversion [15] are of particular interest since they possess special statistic properties. However, due to the intrinsic weak scattering interactions between photons, the experimental realization of $n$-photon states still remains a challenge.

On the other hand, the ability to manipulate individual phonon allows the experimental creation of $n$-phonon states in both circuit quantum acoustodynamics [21,22] and macroscopic mechanical resonators [23–25]. These low-energy and long-lived novel phonon states could facilitate the study of the decoherence mechanisms [26] and the building of the quantum memories and transducers [27–29]. More interestingly, it was also proposed that $n$-phonon bundle states can be generated in acoustic cavity QED by employing the Stokes processes [30] and hybrid system of nitrogen-vacancy centers and nanomechanical resonators via sideband engineering [31].

We note that existing schemes [13–15,30] for creating $n$-quanta bundle states operate in the far dispersive regime where the frequency of the photonic/phononic mode is far detuned from the transition frequency of the two-level system. Hence the resulting steady-state photon/phonon numbers are typically very small. Moreover, in order to resolve states with distinct phonon numbers, many schemes [21–25,30] also require that the two-level system is strongly coupled to the phononic modes, which demands high-finesse acoustic cavities or mechanical resonators with long coherence time and poses a challenge to the current experiments.

In this work, we propose to generate $n$-phonon bundle states by utilizing the motional degrees of freedom of a trapped alkaline-earth atom. A position-dependent clock laser is introduced to couple the atom’s center-of-mass motion to its electronic ground and the long-lived excited states, which leads to a generalized quantum Rabi model (QRM) with unprecedented tunability. We then investigate the $n$-phonon bundle states prepared in the resonant regime where the frequency of the motional mode is in resonance with the two-level atom. We show that due to the strong anharmonicity of the energy spectrum, distinct motional $n$-phonon bundle states can be well resolved in the atom–phonon resonance regime. Compared to the existing schemes for generating $n$-quanta bundle states, the system proposed here has the following advantages. (i) Since our scheme operates in the resonant regime, the typical average steady-state phonon number is much larger than that of the schemes operating in the far dispersive regime. (ii) A strong coupling regime can be readily achieved here as both the motional state and the atomic internal states possess long lifetimes. (iii) In our configuration, the effective pump field is provided by the high tunable clock detuning, and our system also facilitates the study of Mollow physics without suffering heating and decoherence. Therefore, the proposed system can be used as a high-quality source for multiphonon states.

Without loss of generality, we consider a single ${}^{87}\mathrm{Sr}$ atom trapped in a one-dimensional (1D) harmonic potential along the $x$ direction. Figure 1(a) illustrates the level structure and laser configuration of the system. An ultranarrow clock laser with wavelength ${\lambda }_C=698\mathrm{nm}$ drives the single-photon transition between the ground state ${{}^1\mathrm{S}}_0$ ($|g⟩$) and the long-lived excited state ${{}^3\mathrm{P}}_0$ ($|e⟩$) with the atomic clock transition frequency ${\omega }_a$. Since the lifetime of the excited state is roughly 160 s [32], the spontaneous emission and decoherence of this state can be safely ignored. We assume that the Rabi frequency of the clock laser takes a position-dependent form, ${\mathrm{\Omega }}_l(x)=({\mathrm{\Omega }}_{0}x-i\mathrm{\Omega }/2)e^{-i\kappa x}$, where ${\mathrm{\Omega }}_0$ and $\mathrm{\Omega }$ are coupling strengths and $\kappa =k_C\cos \varphi$ is the effective laser wave vector that is tunable by varying the tilting angle $\varphi$ of the clock laser. As shall be shown, the position-dependent Rabi coupling of ${\mathrm{\Omega }}_l(x)$ is crucial to the success of our scheme and it can be experimentally generated by tailoring the clock laser by using a spatial light modulator [33–37]. We further assume that the 1D harmonic potential is state independent which can be generated by a trap laser at the “magic” wavelength ${\lambda }_L=813\mathrm{nm}$. After performing a gauge transformation $|g⟩\to e^{i\kappa x/2}|g⟩$ and $|e⟩\to e^{-i\kappa x/2}|e⟩$ [38], the resulting Hamiltonian of the system is $$\mathcal{H}/\hslash =\omega {\widehat{a}}^{†}\widehat{a}+\frac{\mathrm{\Omega }}{2}{\sigma }_y+\delta {\sigma }_z+\frac{g_x}{2}({\widehat{a}}^{†}+\widehat{a}){\sigma }_x+\frac{i{g}_k}{2}({\widehat{a}}^{†}-\widehat{a}){\sigma }_z,$$(1)where $\omega$ is the frequency of the harmonic trap which, unlike the detuning of an optical cavity mode, is always positive; $\widehat{a}$ is the annihilation operator of the bosonic phonon for the external motional mode; ${\sigma }_{x,y,z}$ are the Pauli matrices; $\delta$ is the single-photon detuning from the bare clock transition; and $g_x=2{\mathrm{\Omega }}_0{x}_0$ and $g_k=\kappa \omega x_0$ with $x_0$ being the zero-point fluctuation amplitude of the harmonic oscillator. Here, $g_x$ and $g_k$ can be understood as the strengths of the spin–orbit coupling (SOC) in the real ($\sim x{\sigma }_x$) and momentum ($\sim p_x{\sigma }_z$) spaces [39], respectively.

To transfer Eq. (1) into a more familiar form, we introduce a spin rotation, $e^{i\pi {\sigma }_x/4}$, which rotates spin operators according to ${\sigma }_x\to {\sigma }_x$, ${\sigma }_y\to {\sigma }_z$, and ${\sigma }_z\to -{\sigma }_y$. As a result, Hamiltonian of Eq. (1) becomes a generalized QRM, $${\mathcal{H}}^{\prime }/\hslash =\omega {\widehat{a}}^{†}\widehat{a}+\frac{\mathrm{\Omega }}{2}{\sigma }_z-\delta {\sigma }_y+\frac{g_x}{2}({\widehat{a}}^{†}+\widehat{a}){\sigma }_x-\frac{i{g}_k}{2}({\widehat{a}}^{†}-\widehat{a}){\sigma }_y,$$(2)which is an essential building block in quantum information and quantum optics [40–43]. Apparently, ${\mathcal{H}}^{\prime }$ reduces to the standard QRM by setting $\delta =0$ and letting either $g_x$ or $g_k$ be zero. Furthermore, Hamiltonian of Eq. (2) has an extremely high tunability. For instance, it turns into the Jaynes–Cummings model (JCM) when $g_k=g_x$ and the anti-JCM when $g_k=-g_x$. Comparing with the experiment using the standard sideband transitions [43], the emerged Hamiltonian of the JCM or the anti-JCM is exact without rotating-wave approximation, and correspondingly the high-frequency term is completely eliminated in Eq. (2). In particular, the $\delta {\sigma }_y$ term in Eq. (2) which represents an effective classical pumping field of the atom is physically realized through the clock shift $\delta$. Consequently, the strong pumping regime ($\delta \gg \omega$) can be readily achieved without suffering severe laser-induced heating, in contrast to the Raman-induced SOC in atomic gases [44].

For a complete description of the system, we should also take into account the dissipation of the phonons. As a result, the dynamics of the atom–phonon system is now described by the master equation $$\frac{\mathrm{d}\rho }{\mathrm{d}t}=-i[{\mathcal{H}}^{\prime },\rho ]+\frac{{\kappa }_e}{2}\mathcal{D}[{\sigma }_{-}]\rho +\frac{{\kappa }_d}{2}\mathcal{D}[\widehat{a}]\rho +{\gamma }_d\mathcal{D}[{\widehat{a}}^{†}\widehat{a}]\rho ,$$(3)where $\rho$ is the density matrix of the atom–phonon system; ${\kappa }_e$ is the atomic spontaneous emission rate of the clock state with ${\sigma }_{-}=({\sigma }_x-i{\sigma }_y)/2$; ${\kappa }_d$ and ${\gamma }_d$ are, respectively, the decay rate and the dephasing factor of the phonons; and $\mathcal{D}[\widehat{o}]\rho =2\widehat{o}\rho {\widehat{o}}^{†}-{\widehat{o}}^{†}\widehat{o}\rho -\rho {\widehat{o}}^{†}\widehat{o}$ denotes the Lindblad type of dissipation. For a given set of parameters, Eq. (3) can be numerically evolved in the basis $\{|n,\sigma ⟩\}$ until a steady-state density matrix is obtained, where $|n⟩$ is the Fock state of the phonon and $|\sigma ⟩=|e⟩$ or $|g⟩$ represents the atomic state. We should emphasize that the realized effective Hamiltonian of Eq. (2) of the generalized QRM is originating from the rotating-wave approximation since $|g_{x,k}/{\omega }_a{|<10}^{-10}$. Thus, the standard Lindblad Eq. (3) is valid even at the coupling strength $g_{x,k}$ comparable to the effective atom (phonon) detuning $\mathrm{\Omega }$ ($\omega$) [45]. Therefore, the dissipative terms of the two-level atom and phonon emerged by external environment in our mode can be treated as independent with safely neglecting atom–phonon coupling [46,47].

To characterize the statistic properties of the phonons, we introduce the generalized $k$th-order correlation function, $$g_n^{(k)}({\tau }_1,\dots ,{\tau }_n)=\frac{⟨\prod _{i=1}^k{[{\widehat{a}}^{†}({\tau }_i)]}^n\prod _{i=1}^k{[\widehat{a}({\tau }_i)]}^n⟩}{\prod _{i=1}^k⟨{[{\widehat{a}}^{†}({\tau }_i)]}^n{[\widehat{a}({\tau }_i)]}^n⟩},$$(4)first introduced in Ref. [13], where ${\tau }_1\le \dots \le {\tau }_n$. As can be seen, this definition generalizes the standard $k$th-order correlation function [48] for isolated phonons to bundles of $n$ phonons. We note that the equal time correlation function $g_n^{(k)}(0)$ for ${\tau }_1=\dots ={\tau }_k$ is straightforwardly calculated using the steady-state density matrix solution of Eq. (3); while the multitime correlation function $g_n^{(k)}({\tau }_1,\dots ,{\tau }_k)$ with ${\tau }_1<\dots <{\tau }_k$ can be obtained by utilizing the quantum regression theorem [49]. Now, an $n$-phonon bundle state ($n>1$) should satisfy two conditions: $g_1^{(2)}(0)>g_1^{(2)}(\tau )$ to ensure bunching with respect to single phonon and $g_n^{(2)}(0)<g_n^{(2)}(\tau )$ to secure antibunching between separated bundles of phonons [13,15]. We point out that the correlation functions of $g_1^{(2)}$ and $g_n^{(2)}$ can be directly extracted from the steady-state phonon-number distribution $p(q)=\mathrm{tr}(|q⟩⟨q|\rho )$ which, for the motional phonons, can be measured via spin state-resolved projective measurement [43,50]. While for acoustic phonons, although $g_1^{(2)}$ can be measured via a Hanbury Brown and Twiss interferometer [51], experimental measuring of the multiphonon $g_n^{(2)}$ is still a challenging task.

Finally, we specify the parameters used in numerical simulations. For the ${}^{87}\mathrm{Sr}$ atom, the single-photon recoil energy $E_C/h=4.68\mathrm{kHz}$ sets up the typical energy scale for the system, where $h$ is the Planck constant. In current experiments, the phonon frequency $\omega$ is tunable and can be up to about $(2\pi )100\mathrm{kHz}$ [32]. As a result, we realize a tunable $g_k/h$ bound between $-21.6\mathrm{kHz}$ and 21.6 kHz. For convenience, we fix the value of the SOC strength in real space at $g_x/h=9\mathrm{kHz}$ such that it is of the same order as $g_k/h$, which yields the atomic decay for the clock state ${\kappa }_e/g_x=1.1\times {10}^{-8}$ [32]. We further fix the values of the decay rate and the dephasing factor at ${\kappa }_d/g_x=0.005$ and ${\gamma }_d/g_x=0.0005$, respectively. Now, the free parameters of the system reduce to phonon frequency $\omega$, Rabi frequency $\mathrm{\Omega }$, laser detuning $\delta$, and SOC strength $g_k$. In below, we study the statistic properties of the phonon state by varying these parameters.

Before presenting our results on phonon statistics, it is instructive to explore the energy spectrum of the system. In Fig. 1(b), we demonstrate the familiar level structure of a JCM ($g_k=g_x$) with $\mathrm{\Omega }=\omega$ and $\delta =0$. The eigenenergies of the $n$th pair of dressed states $|n,\pm ⟩$ are $E_{n,\pm }=\hslash (n\omega \pm g_x\sqrt{n})$, where $n=\mathrm{1,}\mathrm{2,}\dots$ and “$+$” (“$-$”) denotes the upper (lower) branch. Particularly, $n$-phonon resonance occurs when the lower dressed state $|n,-⟩$ is tuned on resonance with the vacuum state of the system, $\omega ={\omega }_n$, where ${\omega }_n=g_x/\sqrt{n}$ is the $n$-phonon resonance frequency. The energy spectrum for the more general case is presented in Appendix B, where ${\omega }_n$ should also depend on other parameters.

Let us first consider the single phonon states by fixing the phonon frequency at $\omega =g_x$. Figures 2(a) and 2(b) display, respectively, the equal time second-order correlation function $g_1^{(2)}(0)$ and steady-state phonon number $n_s=\mathrm{tr}({\widehat{a}}^{†}\widehat{a}\rho )$ in the $\mathrm{\Omega }{g}_k$ parameter plane with $\delta /g_x=0.005$. As can be seen, when $g_k$ and $\mathrm{\Omega }$ are changed, the values of $g_1^{(2)}(0)$ and $n_s$ vary over a wide range. Particularly, around $(\mathrm{\Omega },g_k)=(\pm 1,\pm 1)g_x$, the system reaches the strong sub-Poissonian statistics region with $g_1^{(2)}(0)\sim {10}^{-4}$. In these regions, a considerably large number of phonons ($n_s>0.2$) are observed as well. In Fig. 2(c), we further plot the interval dependence of the second-order correlation function $g_1^{(2)}(\tau )$ at $(\mathrm{\Omega },g_k)=(\pm 1,\pm 1)g_x$, which shows phonon antibunching since $g_1^{(2)}(0)<g_1^{(2)}(\tau )$. The evidence clearly shows that strong phonon blockade (SPB) is achieved in these regions.

The SPB around $(\mathrm{\Omega },g_k)=(\mathrm{1,1})g_x$ can be understood by noting that the energy spectrum of the resulting JCM at $g_k=g_x$ is highly anharmonic in the strong coupling regime $g_x/\kappa \gg 1$ [Fig. 1(b)]. Therefore, the condition for one phonon excitation at single-phonon resonance ($\mathrm{\Omega }=g_x$) will block the excitation of a second phonon. At first sight, the SPB around $(\mathrm{\Omega },g_k)=(-1,-1)g_x$ may seem strange, as the anti-JCM realized at $g_k=-g_x$ breaks the conservation of the number of the total excitations due to the counter-rotating terms. To explain this, we note that an anti-JCM is equivalent to a JCM under the unitary transformation $e^{-i\pi {\sigma }_x/2}$, i.e., ${\sigma }_{\pm }\to {\sigma }_{\mp }$ and ${\sigma }_z\to -{\sigma }_z$. Consequently, an SPB should occur at $-\mathrm{\Omega }=g_x$. We point out that, because of the small phonon decay in the clock transition, the lifetime of the phonon blockade ($\tau \sim {\kappa }_d^{-1}$) in our system can be very long.

To reveal more details, we also plot, in Fig. 2(d), the distributions of $g_1^{(2)}(0)$ and $n_s$ along the line $g_k=\mathrm{\Omega }$ on the $\mathrm{\Omega }{g}_k$ plane. As can be seen, both $g_1^{(2)}(0)$ and $n_s$ possess two local minima close to $g_k=\pm g_x$. In addition, $g_1^{(2)}(0)$ grows rapidly from these minima when $g_k$ is slightly tuned away from $\pm g_x$. $g_1^{(2)}(0)$ approaches unity for a coherent phonon state at $g_k=0$ where Hamiltonian of Eq. (2) reduces to a standard QRM.

We now turn to study the properties of phonon states in the atom–phonon resonance regime, i.e., $\mathrm{\Omega }=\omega$. In addition, due to the equivalence of two SPB regions in Fig. 2(a), we shall focus, without loss of generality, on the SPB region with $g_k=g_x$. The atom–phonon resonance condition can be satisfied by tuning external atom trap potential and/or power of the classical laser field in the experiment. In particular, it should be noted that our results remain qualitatively unchanged even deviating from the resonance condition, corresponding to the slight shift in the position of $n$-phonon resonance ${\omega }_n$.

Figure 3(a) shows the phonon number $n_s$ as a function of the phonon frequency $\omega$ and clock shift $\delta$. In particular, the dashed lines plot the $\delta$ dependence of the $n$-phonon resonance frequency ${\omega }_n$ obtained by numerically diagonalizing ${\mathcal{H}}^{\prime }$. An immediate observation is that, for a given $\delta$, $n_s(\omega )$ exhibits multiple peaks at exactly the phonon resonance frequencies ${\omega }_n$, signaling the existence of possible multiphonon states. Remarkably, because our scheme operates in the resonant regime, the average phonon number of the phonon states generated here is much larger than that of the phonon states produced in the dispersive regime. In addition, there exists a large area on the $\omega \delta$ plane in which the average phonon number is higher than 0.2.

To explore the statistic properties of the phonon emissions, we first map out, in Fig. 3(b), $g_1^{(2)}(0)$ on the $\omega \delta$ parameter plane around $\omega ={\omega }_1$. As can be seen, phonon blockade is realized for the whole parameter space covered by Fig. 3(b). There even exists a large area on the $\omega \delta$ plane such that the strong phonon blockade condition, say $g_1^{(2)}{(0)<10}^{-2}$, is satisfied. Then combined with requirement of large phonon emission number, say $n_s>0.2$, our system can be used as a high-quality single motional phonon source operating in a parameter regime that is easily accessible to current experiments. To further quantify the quality of the single-phonon states around ${\omega }_1$, we introduce $\tilde{p}(q)\equiv qp(q)/n_s$, which measures the fraction of $q$-phonon states among the total emitted phonons. In Fig. 3(c), we plot the distribution of $\tilde{p}(1)$ on the $\omega \delta$ plane. As can be seen, for the parameter region of our interest, nearly 100% of phonon emission is of the single-phonon nature. It should also be noted that $\tilde{p}(1)$ alone is not sufficient to judge the quality of the single-phonon source, because, by comparing Figs. 3(b) and 3(c), the main feature of $\tilde{p}(1)$ is inconsistent with that of $g_1^{(2)}(0)$ when $n_s$ is small.

Finally, we study the properties of the emitted phonon in the multiphonon resonance regime. As shown in Fig. 3(a), only for sufficiently large pumping, i.e., $\delta /g\ge 0.015$ and 0.094 for two- and three-phonon resonances, respectively, the phonon numbers $n_s$ around ${\omega }_2$ and ${\omega }_3$ become significant ($>0.2$). We shall then explore the statistic properties of the multiphonon resonances in these regions. Figure 4 summarizes the main properties of the motional multiphonon states. Due to the similarity between the two- and three-phonon states, we shall only discuss the two-phonon emissions. In Fig. 4(a), we map out the third-order correlation function in the two-phonon resonance regime of the $\omega \delta$ plane. As can be seen, there exists a large parameter regime with $g_1^{(3)}(0)<1$ in which the $(n+1)$-phonon emission is blockaded. Meanwhile, the minimum value of $g_1^{(3)}(0)$ is achieved at $\omega ={\omega }_2$, where the highest phonon numbers for the multiphonon resonance is reached. This indicates that at the strongest three-phonon blockade, we have the highest phonon number.

To further confirm the bundle-emission nature of the phonon states, we plot, in Fig. 4(c), the typical interval dependence of the correlation functions $g_1^{(2)}(\tau )$ and $g_2^{(2)}(\tau )$. As can be seen, the criteria for motional two-phonon bundle states, $g_1^{(2)}(0)>g_1^{(2)}(\tau )$ and $g_2^{(2)}(0)<g_2^{(2)}(\tau )$, are indeed satisfied. Another observation from Fig. 4(c) is that the decay times for both $g_1^{(2)}(\tau )$ and $g_2^{(2)}(\tau )$ are proportional to ${\kappa }_d^{-1}$, which indicates that the decay of the bunching for single phonon and the decay of the antibunching for the separated bundles of phonons are of the same time scale for the two-phonon bundle states. The two-phonon nature of the emission is also demonstrated by the phonon-number distribution $\tilde{p}(q)$ shown in Fig. 4(e). Indeed, for the two-phonon emission case, $\tilde{p}(q)$ becomes negligibly small for $q>2$. We point out that, as shown in the right panels of Fig. 4, the three-phonon emission possesses similar statistic properties as those of two-phonon emission. The observed phonon probability for steady-state $n$-phonon bundle states exhibits a monotonical decreasing distribution with increasing $q$ for $1\le q\le n$, arising in dynamical processes of bundle emissions [13]. In general, the $n$-phonon bundle states contain the various dynamical processes of emissions in Fock state $|q⟩$ with a distinguishable short temporal window $\sim 1/(q{\kappa }_d)$. Thus, the multiphonon state ($|q⟩$) can be directly extracted with a high phonon probability by choosing a very short temporal window [13,30]. Moreover, we should note that the generated motional $n$-phonon bundle states are essentially different from the experimentally observed $n$-quanta blockade, i.e., two-photon blockade [52], where the quantum statistic for the latter only characterizes the single photon but not for separated bundles of photons with satisfying $n$-photon bunching $g_1^{(n)}(0)>1$ and $(n+1)$-photon antibunching $g_1^{(n+1)}(0)<1$ as well.

We remark that the underlying reason that $n$-phonon bundle states are well-resolved in the parameters space is the strong anharmonicity of the energy spectrum such that the condition $|{\omega }_n-{\omega }_m|\gg {\kappa }_d$ is satisfied for any $n\ne m$. This is in striking contrast to the proposed $n$-photon bundle emission utilizing Mollow physics [53] in which the $(n+1)$th-order process of quantum states is used [13,14]. In addition, the whole process must be operated in the far dispersive regime and under a strong pump field, which, in terms of our model, requires $|\omega -\mathrm{\Omega }|\gg g_x\sqrt{n+1}$ and $|\delta |\gg g_x\sqrt{n+1}$, simultaneously. Finally, we check that our numerical results are not affected when we include the weak atomic decay of the clock state (${\kappa }_e/{\kappa }_d\ll 1$) and phonon dephasing (${\gamma }_d/{\kappa }_d\ll 1$), albeit the strong damping of ${\gamma }_d$ and ${\gamma }_e$ could induce a significant decoherence for realization of nonclassical quantum states [54–56].

Based upon the currently available techniques in experiments, we have proposed to generate motional $n$-phonon bundle states using a trapped alkaline-earth atom driven by a clock laser. Since our system works in the resonant regime, the steady-state phonon number contained in the bundle states is 3 orders of magnitude larger than those obtained in the earlier theoretical schemes operating in the far dispersive regime. Moreover, the quality of the $n$-phonon bundle states is also demonstrated by the strong antibunching for the separated bundles of phonons and bunching for single phonons. Finally, we emphasize that the nonclassical nature of the long-lived motional $n$-phonon bundle states can be quantum-state transferred to the photonic mode by applying a readout cavity field with phonon–photon beam-splitter interaction [57,58]. Moreover, as the scheme mitigates the laser-induced heating, it could inspire an interesting opportunity of exploration of Dicke phase transition for the external motional modes [59], superradiances from the clock transition [60], and novel quantum states of matters hindered by heating [61–63]. Our proposal for trapped single atoms could be equivalently applied to the hybrid spin–mechanical systems [64,65]. In particular, the parameter in our model is outside the Lamb–Dicke regime with $|\kappa x_0|\sim 0.2$, which could provide a versatile platform for exploring long-lived mesoscopic entanglement for trapped atoms [66]. Furthermore, we could expect that the proposed system provides versatile applications in quantum metrology limited by decoherence [50,67] and in fundamental tests of quantum physics [68].

Acknowledgment. We are grateful to Yue Chang for insightful discussions.