Kapitza-Dirac blockade: A universal tool for the deterministic preparation of non-Gaussian oscillator states

Harmonic oscillators count among the most fundamental quantum systems with important applications in molecular physics, nanoparticle trapping, and quantum information processing. Their equidistant energy level spacing is often a desired feature, but at the same time a challenge if the goal is to deterministically populate specific eigenstates. Here, we show how interference in the transition amplitudes in a bichromatic laser field can suppress the sequential climbing of harmonic oscillator states (Kapitza-Dirac blockade) and achieve selective excitation of energy eigenstates, Schr\"{o}dinger cats and other non-Gaussian states. This technique can transform the harmonic oscillator into a coherent two-level system or be used to build a large-momentum-transfer beam splitter for matter-waves. To illustrate the universality of the concept, we discuss feasible experiments that cover many orders of magnitude in mass, from single electrons over large molecules to dielectric nanoparticles.

The harmonic oscillator is a paradigmatic text book example of fundamental quantum physics and it has remained at the heart of modern research. Quantum harmonic oscillators have been realized with single electrons [1,2], single ions [3], ultra-cold quantum gases [4], and dielectric nanoparticles [5]. For all these systems, cooling to the oscillator ground state has been successfully demonstrated. Our present proposal is motivated by the challenge to prepare highly non-classical states, mesoscopic Schrödinger cat states and large-momentum transfer beam splitters, independent of detailed oscillator properties.
Throughout the last two decades, macroscopic quantum superposition was realized in widely different systems [6]. Neutrons were delocalized over 10 cm [7]. Atoms were put in superpositions on the half-meter scale [8] or in momentum states separated by more than 1000hk [9], and molecules in excess of 25,000 Da were delocalized over hundred times their size [10]. Lately, it has been proposed to prepare dielectric nanoparticles in distinct position states [6,11] to test the nature of quantum collapse [12], quantum decoherence [13][14][15], or even the quantum nature of gravity [16,17].
Here we propose to exploit Kapitza-Dirac blockade as a universal tool for preparing Schrödinger cat or non-Gaussian states in general, with single electrons, molecules and nanoparticles, differing in mass by more than nine orders of magnitude.
In the original proposal by Kapitza and Dirac [18], electrons were assumed to be coherently scattered by the ponderomotive force of a standing light wave. It took several advances in electron beam and laser technology to realize this idea 70 years later, both in the Raman-Nath and in the Bragg regime [19,20]. The idea of optical phase gratings has also been extended to atoms [21][22][23] and molecules [24,25], where the interaction is mediated by the optical dipole force between the laser electric field and the particle's polarizability. Recently, the Kapitza-Dirac effect has been used in electron microscopy to facilitate all-optical phase masks [26,27].
The inelastic Kapitza-Dirac (KD) effect was first discussed in [28]. It differs from its elastic counterpart not only by the use of laser fields with different frequencies, but also by the presence of a harmonic trap. The latter modifies the conditions for energy and momentum conservation. The inelastic KD-effect is similar to stimulated Raman scattering [29][30][31][32] but operates without internal states, so that absorption, spontaneous emission and decoherence in objects with broad resonance lines can be avoided. It can thus be applied to particles that do not exhibit any internal states at all.
In this letter, we introduce Kapitza-Dirac blockade as  Kapitza-Dirac blockade and quantum state control in a 1D harmonic oscillator. A collimated particle populates the ground state of a 1D harmonic trap. A pair of counter-propagating bichromatic KD-laser fields E1,2 interacts with the particle and changes its energy and momentum while it is in the trap. a new feature emerging from the quantum mechanical treatment of the inelastic KD-effect. By using judiciously chosen laser frequencies, the Kapitza-Dirac blockade can drive a controlled parametric resonance while blocking undesired transitions almost entirely.
We start for simplicity with a particle traveling through a 1D harmonic trap (see FIG. 1). The trap can be realized via the ponderomotive potential on electrons or via the optical dipole potential on polarizable particles. The scheme preferably starts from the harmonic oscillator ground state, which can be populated with high probability by filling the trap with a tightly collimated particle beam.
The two KD-laser pulses with frequencies ω 1,2 (ω 1 > ω 2 ) and a 1/e pulse duration τ KD are assumed to propagate counter to each other along the x-axis. Their polarization should be chosen to avoid wave-mixing with the trapping lasers.
The transition amplitude g n (η) has zero-crossings (see FIG. 3(a)), which implies that the transition |n → |n+2 can be suppressed at certain values of the momentum detuning δ p because of destructive interference between transition amplitudes starting from different k-values (see FIG. 2). This Kapitza-Dirac blockade is a powerful tool, as it can stop the sequential excitation in the energy ladder and allows us to prepare non-Gaussian harmonic oscillator states. The energy-momentum conservation Eq. (2) implies that the Lamb-Dicke parameter η = (N m + δ p )/2 is independent of any oscillator details. In consequence, the Kapitza-Dirac blockade is independent of the specific oscillator realization, and only the central KD-frequency ω KD = ck 0 η is oscillator dependent.
As a first example, we propose to prepare a single energy eigenstate |n = 2 , starting from the ground state |n = 0 . Setting δ p = 0.83 suppresses the |n = 2 → |n = 4 transition down to < 0.2 % of its maximum value (see FIG. 3(c)). As a result, transitions starting in |n = 0 will end deterministically in |n = 2 without populating |n = 4 or other excited states in the energy ladder (see FIG. 4(a)). For a Gaussian-shaped time envelope, complete population inversion occurs when where Ω R is the system specific Rabi frequency [33]. Note, that τ KD should be sufficiently long to keep the pulse bandwidth below the trap frequency Ω 0 in order to avoid off-resonant excitation. When the pulse duration is doubled, the Rabi cycle is completed to the ground state (see FIG. 4(b)). The Kapitza-Dirac blockade has thus transformed the harmonic oscillator into an effective two-level system -with promising applications in quantum information processing. Kapitza-Dirac blockade can also prepare a Schrödinger cat state. We show this with a heuristic example of δ p = −0.60 which suppresses the transition |n = 12 → |n = 14 to < 0.1 % of its maximum value (see FIG. 3(b)). The oscillator then undergoes sequential excitation from |n = 0 up to |n = 12 in steps of ∆n = 2, but |n = 14 is not excited. The pulse intensity I KD and duration τ KD are adjusted together to maximize the population distribution at n max = 8 while avoiding off-resonant excitation. The width of the final population distribution is sub-Poissonian (see FIG. 4(c)), due to the Kapitza-Dirac blockade. This leads to an amplitude-squeezed Schrödinger cat state which is identified by inspecting the Wigner function in FIG. 4(d). The maximum spatial and momentum separation of the cat state, ∆x cat and ∆p cat , are Substituting Eq. (2) to Eq. (5) with N m = 2, the number of photon recoils is which is again independent of the oscillator properties. If we take n max = 650 and δ p = −1.8, the maxi-mum momentum separation is ∆p cat ≈ 1000hk 532 , where k 532 ≡ 2π/532 nm. Kapitza-Dirac blockade is therefore a promising tool for realizing an all-optical large-momentum-transfer (LMT) beam splitter [8,9,36,37]. Moreover, the amplitude-squeezed Schrödinger cat state leaves the divided beam rather well-collimated. A Kapitza-Dirac-LMT beam splitter used in conjunction with an optical Bragg grating [20,25,38] could facilitate large-area matter-wave interferometry, in the future.
Also multi-component cat states [39][40][41] can be prepared by properly choosing N m and δ p [33]. An example is given in FIG. 5 where a 3-component cat state is prepared. The KD-laser frequencies ω 1,2 in this case are determined by taking N m = 3 and δ p = −1.57, which suppresses the transition |n = 18 → |n = 21 . As in multi-slit diffraction, interference among multiple components of the cat state produces sharper fringes compared to those between each pair (see FIG. 5(b)). Thus, a multi-component cat state can be more sensitive in probing quantum decoherence or collapse effects.
Similarly, Gaussian states such as vacuum squeezed state (N m = 2) or coherent state (N m = 1) can be prepared with appropriate I KD and δ p .
Kapitza-Dirac blockade holds universally. Here we start the experimental discussion with the example of an electron beam 1D-trapped by the ponderomotive potential of a standing light wave [42]. Close to the potential minimum, the potential can be approximated as a harmonic trap U p (x) ≈ q 2 e I S /2 0 c 3 m e x 2 , where m e and q e are the electron mass and charge, and I S is the standing wave intensity. The ponderomotive trap frequency is The intensity of each trapping laser that makes the standing wave is I TL = I S /4. The trap ground state can be populated by a well-collimated electron beam with transverse kinetic energy mv 2 x /2 hΩ 0 /2. If we consider the LMT beam splitter example above, using parameters as in Table I, and if we tune the dynamics such that the electron reaches maximum momentum separation when it leaves the trap (see FIG. 1), it will leave in two distinct wave packets separated by 171 µm in real space, already 1 mm behind the trap [43][44][45].
All aspects discussed above can equally be realized with neutral massive particles: in Table I we discuss as the second example a porphyrin derivative (TPPF84) with high mass and high vapor pressure [46]. The third example is a silicon dioxide (SiO 2 ) nanoparticle with low absorption of infrared trapping light.
The experimental schemes for all three systems can be similar to that of FIG. 1. For molecules and nanoparticles a harmonic trap can be realized by the dipole potential of a standing wave. With λ TL as the trapping laser wavelength, the dipole trap frequency is where α and m are the particle's polarizability and mass.
There are three criteria for choosing experimental parameters. First, the laser conditions should satisfy τ TL , where τ is the 1/e pulse duration, W z is the 1/e beam diameter along the z-axis, and v z is the particle speed in the z-direction. Second, the number of Rayleigh scattered photons should be small, n sca < 1, to avoid decoherence and dephasing. Third, the central KD-frequency ω KD should be in the visible or the infrared regime because ultraviolet light would be absorbed in most materials. This implies that lower trap frequencies are preferred for more massive particles. The proposed experimental parameters for large amplitude-squeezed Schrödinger cats of electrons, molecules, and nanoparticles are listed in Table I, designed according to the empirical formulas where |n bk is the blockade state and the transition |n bk → |n bk + 2 is suppressed. The system specific coefficient D is m e ω KD /q 2 e for electrons or 1/αω KD for molecules and nanoparticles. The momentum transfer ∆p cat is well comparable with the state-of-the-art in matter-wave beam splitting [8,9,25,36,37], but going beyond it as the concept can be applied universally to any 1D trapped particle that scatters light coherently. Additionally, the use of the Kapitza-Dirac blockade yields narrow momentum distributions and avoids overlap between the supposedly distinct wave packets.
Wy-KD (µm) 100 10 100 The simulations were performed on a supercomputer using time-dependent Schrödinger equations [33,47]. The results are shown in FIG. 6. The momentum transfer ∆p cat is 1000hk 532 for electrons, 650hk 532 for molecules, and 2900hk 532 for nanoparticles. The highest eigenstate available for excitation is determined by the onset of trap anharmonicity [33]. In the nanoparticle simulation, the maximal spatial separation in the trap is ∆x cat ≈ 1.2 µm. The required beam velocities for molecules and nanoparticles are two orders of magnitude lower than the stateof-the-art of free beams. However, cooling of molecules and nanoparticles is a rapidly advancing field and the proposed parameters are within reach [48][49][50][51][52].
In conclusion, we propose to use the Kapitza-Dirac blockade for manipulating the motional quantum states of single electrons, large molecules, and dielectric nanoparticles. The state preparation scheme is universal, applicable for different particles and independent of trap details and dimensions.
Our simulations demonstrate the experimental feasibility of preparing various non-Gaussian states and large amplitude-squeezed Schrödinger cat states. The latter has applications for all-optical LMT beam splitters in matter-wave interferometry. All these results explicitly rely on the coherent but destructive interference in the harmonic oscillator transition amplitudes in the presence of bichromatic light fields: the Kapitza-Dirac blockade.  The proposed control scheme can also be employed for trapped ions and neutral atoms without invoking their internal states. In 2D or 3D harmonic traps one also finds entanglement among different motional degrees of freedom [33]. This can complement existing methods in quantum computing and quantum simulation.
The authors thank Aephraim M. Steinberg, Peter W. Milonni, Christopher Monroe, and Uroš Delić for advice and discussions. W. C. Huang wishes to give a special thanks to Yanshuo Li for supports and helpful discussions. This work utilized high-performance computing resources from the Holland Computing Center of the University of Nebraska. Funding for this work comes from NSF EPS-1430519 and NSF PHY-1602755.
This supplementary material is aimed for providing a detailed account for the quantum mechanical treatment of the inelastic Kapitza-Dirac (KD) effect. Additionally, we extend the analysis to 2D and 3D harmonic oscillators and show that not only the Kapitza-Dirac blockade persists to work with added dimensionality, but the possibility of entanglement also arises.
where p = −ih∇ and A = A TL + A 1 + A 2 . We use q e = −e and m e as electron charge and mass. The vector potential A TL of the standing wave is where k TL = ω TL /c and I S is the intensity of the standing wave at the center of the waist. Note that the individual trapping laser beam in the standing wave has a pulse intensity I TL = I S /4 and the 1/e pulse duration τ TL . We denote x as the particle position operator hereafter. The pair of counter-propagating running waves are called KD-lasers, and their vector potential A 1,2 are where A 1,2 = 2I KD / 0 cω 2 1,2 , k 1,2 = ω 1,2 /c, I KD is the intensity of individual KD-lasers at the center of the waist, and τ KD is the 1/e pulse duration. Here we approximate all laser fields to be homogenous over the area of interaction in the y-z plane. This approximation is justified when (1) the laser beam sizes in the y-direction are much larger than the electron beam size along the same direction, and (2) the laser beam sizes in the z-direction satisfy the relation where w z is the 1/e beam diameter along the z-axis, and v z is the electron's speed along the z-axis. The subscripts -TL and -KD denote the trapping and the KD-lasers respectively. The first term in Eq. (1) is the electron kinetic energy. The second term in Eq. (1) is simplified to be p · A = (−ih∇) · A + A · (−ih∇) = A · p in the Coulomb gauge ∇ · A = 0. Thus it can be combined with the third term A · p. At the first-order, the A · p term corresponds to single-photon scatterings which do not conserve energy and momentum for electrons. However, at the second order the A · p term can give rise to wave mixing and result in resonant scattering. From a classical analysis [1], it can be shown that the combined action of the A · p term and the A 2 term in Eq. (1) is equivalent to keeping only the A 2 term in Eq. (1). Since the polarizations of the standing wave and the running waves are orthogonal A TL ·A 1,2 = 0, the A 2 term can be expanded as A 2 = A 2 TL +A 2 1 +A 2 2 +2A 1 ·A 2 and the Hamiltonian in Eq. (1) can be rewritten aŝ Assuming that the standing wave is much stronger than the running waves, |A TL | |A 1,2 |, the Hamiltonian can be separated into an unperturbed partĤ and an interaction partĤ Using Eq. (2), the unperturbed Hamiltonian can be written asĤ 0 = p 2 /2m e + U (x, t), where The first term in Eq. (8) is time-dependent, while the second term with cos (2ω TL t) can be discarded because it only causes non-resonant driving and the effect averages out. Thus, we obtain the ponderomotive potential for the well-known elastic KD-effect, where the Gaussian-shaped time envelope is approximated as exp (−2t 2 /τ 2 TL ) ≈ 1 because it is slow-varying during the interaction time τ KD . Depending on the relation between the potential strength U 0 = q 2 e I S /2 0 cω 2 TL m e , the recoil energy E R =h 2 (2k TL ) 2 /2m e , and the transit time ∆t = W z-TL /v z , there are three parameter regimes of interaction: [5]. We assume that the experimental parameters are chosen for interaction in the channeling regime, so that the electron is transversely confined by the ponderomotive potential as it travels through the standing wave. Close to the potential minimum x m = λ TL /4, where λ TL = 2π/k TL , the ponderomotive potential can be approximated as a harmonic trap where x is repurposed to denote the displacement from the potential minimum x m , and the trap frequency is Therefore, the unperturbed Hamiltonian in Eq. (6) becomeŝ where (p x , p y , p z ) ≡ −ih(∂/∂ x , ∂/∂ y , ∂/∂ z ). Thus, for parameters in the channeling regime the x-motion of the electron is confined in an 1D harmonic trap while the y-and z-motions freely propagate. The general solution to the electron's wavefunction is where φ(y, z, t) is the wavefunction of a free wavepacket, and Ω n ≡ Ω 0 (n + 1/2) characterizes the excitation of the harmonic oscillator in the x-direction. We are interested in solving the probability amplitude C n (t) for each eigenstate |n of the harmonic oscillator. Using Eq. (3), the interaction Hamiltonian in Eq. (7) can be expanded aŝ and the corresponding Schrödinger equation for C n (t) is where Ω mn ≡ Ω m − Ω n . To analyze the excitation dynamics, we group Eq. (15) into even and odd transitions, The terms in the first summation represent even transitions, while those in the second summation represent odd transitions. Assuming that ω 1,2 are many orders of magnitude higher than Ω 0 , only the ω 1 − ω 2 term in Eq. (14) can resonantly excite the oscillator, which leads to a simplification of Eq. (14), where cos ((ω 1 − ω 2 )t) and sin ((ω 1 − ω 2 )t) are expanded. In Eq. (17), the first term with cos ((k 1 + k 2 )x) gives rise to even transitions |n → |n ± 2k , where k is a positive integer. This is because the position operator of a harmonic oscillator can be expressed in terms of the rising and lowering operators, x = x 0 (b † +b), where x 0 = h/2m e Ω 0 , b † |n = √ n + 1|n + 1 , andb|n = √ n|n − 1 . The expansion of cos ((k 1 + k 2 )x) contains only even powers ofb † andb.
The resonant excitation frequency for even transitions |n → |n±2k is ω 1 −ω 2 = 2kΩ 0 because e ±i(ω1−ω2)t = e ±i2kΩ0t in Eq. (17) cancel with e ∓i2nΩ0t in Eq. (16) for n = k. Similarly, the sin ((k 1 + k 2 )x) term gives rise to the odd transitions, and the corresponding excitation frequency is ω 1 − ω 2 = (2k + 1)Ω 0 . The above observation can also be derived from the basis of parametric resonance. Let us switch to the interaction picture and replace the position operator x in Eq. (17) by x(t), To obtain the frequency components of cos ((k 1 + k 2 )x) and sin ((k 1 + k 2 )x) in Eq. (17), we first find the frequency components in x m (t) using the binomial expansion, The frequency components, |m − 2n|Ω 0 , are either even or odd depending on the power m. Since there are only even powers of x(t) in cos ((k 1 + k 2 )x), the frequency components in cos ((k 1 + k 2 )x) are all even. Similarly, the frequency components in sin ((k 1 + k 2 )x) are all odd. Given the excitation frequency ω 1 − ω 2 = 2kΩ 0 , there will be some frequency components in cos ((k 1 + k 2 )x) that cancel with e ±i(ω1−ω2)t and lead to parametric excitation. On the other hand, because cos ((k 1 + k 2 )x) contains only even powers ofb † andb, the transitions involved in parametric excitation at ω 1 − ω 2 = 2kΩ 0 are all even. A rigorous and general proof can be obtained by examining the transition matrix element of the interaction Hamiltonian in Eq. (17), where We can find m| cos ((k 1 + k 2 )x)|n and m| sin ((k 1 + k 2 )x)|n by first evaluating the matrix element where φ m (x) is an eigen-wavefunction of the harmonic oscillator in real space. Using the convolution theorem where ψ 1,2 (k) = (1/2π) ∞ −∞ dxφ 1,2 (x)e −ikx , Eq. (22) can be transformed to Here we remark that Eq. (24) depicts the essential element in both elastic and inelastic KD-effect -a periodic potential gives rise to momentum kicks. While Eq. (24) is general and valid for both free particles and harmonic oscillators, in the case of harmonic oscillators, there is one additional property due to the parity of the oscillator eigenstates |m and |n , The physical interpretation of Eq. (25) is that for harmonic oscillators the transition amplitudes for stimulated emission (the integral with ψ n (k + (k 1 + k 2 ))) and absorption (the integral with ψ n (k − (k 1 + k 2 ))) are equal in strength but differ by a phase factor of the parity change (−1) m+n . Using Eqs. (24) and (25), the transition matrix element in Eq. (20) can be evaluated as For even transitions m = n + 2k, where k is a positive integer, Eq. (26) becomes For odd transition m = n + 2k + 1, Eq. (26) becomes If ω 1 − ω 2 = 2kΩ 0 , it can be shown from Eq. (16) that only even transitions are resonantly driven. Their transition matrix elements can be evaluated by Eq. (27) which involves only the cos ((k 1 + k 2 )x) term in Eq. (17). Similarly, in the case of ω 1 − ω 2 = (2k + 1)Ω 0 only the sin ((k 1 + k 2 )x) term in Eq. (17) is involved in the resonant excitation. With Eqs. (16) and (26), we can summarize the conditions for resonant excitation in terms of energy-momentum conservation, where k 0 = m e Ω 0 /2h, N m is a positive integer, and δ p ≥ −N m is the momentum detuning discussed in the main text.
The energy-momentum relation in Eq. (29) and the transition amplitudes in Eqs. (27) and (28) together give a clear physical picture of the inelastic KD-effect as an oscillator receiving a momentum kickh(k 1 + k 2 ) and making a transition to a higher energy level. This picture allows us to give a physical interpretation for the Kapitza-Dirac blockade mechanism as described in FIG. 2 of the main text. However, the integral form of Eqs. (27) and (28) are impractical to implement in numerical simulations, especially when the number of the involved eigenstates is large. Additionally, to identify the exact value of the momentum detuning δ p for Kapitza-Dirac blockade, it is desirable to derive the analytic forms for Eqs. (27) and (28). Cahill and Glauber [2] have shown that for a harmonic oscillator, where η ≡ (k 1 + k 2 )x 0 is the Lamb-Dicke parameter, and L (m−n) n (y) is the generalized Laguerre polynomial. This equation gives Eqs. (27) and (28) the desired analytic forms, So far we have been discussing the general theory of the inelastic KD-effect. Let us apply these results to the specific cases discussed in the main text. In the case of ω 1 − ω 2 = 2Ω 0 , which is used in the main text to populate a single eigenstate |n = 2 and to prepare a Schrödinger cat state, the interaction Hamiltonian in Eq. (17) is reduced tô where f (t) is given in Eq. (21). Note that this interaction Hamiltonian is the effective KD-potential V KD (x, t) in the main text. Using Eq. (31), we define the dimensionless transition amplitude Therefore, according to Eq. (16) the corresponding Schrödinger equation is where hermiticity is used to give n + 2|Ĥ int |n = n|Ĥ int |n + 2 . We note that even in the Lamb-Dicke regime η 1, the inelastic KD-effect as described in Eq. (34) can still be driven. For ω 1 − ω 2 = 3Ω 0 , which is used in the main text to prepare a 3-component Schrödinger cat state, the interaction Hamiltonian in Eq. (17) is reduced tô The dimensionless transition amplitude ξ n (η) ≡ n + 3| sin ((k 1 + k 2 )x)|n = − n!/(n + 3)!η 3 L n (η 2 )e −η 2 /2 is used to identify the proper momentum detuing δ p = −1.57 for the Kapitza-Dirac blockade at |n = 18 . The corresponding Schrödinger equation is In FIG. 2, we give an illustration for the transition processes in Eqs. (34) and (36).
(b) A polarizable particle in a 1D dipole trap Next, we consider the interaction between two pairs of counter-propagating laser fields and a polarizable neutral particle (i.e. atoms, molecules, and dielectric nanoparticles). The proposed experimental setup resembles that in FIG. 1 with the only difference that the lasers are continuous-wave instead of pulsed. This is because the time scale of dynamics as determined by the trap frequency is much longer for massive particles compared to electrons. Polarizable neutral particles interact with laser fields through the dipole force, so the quantum Hamiltonian iŝ where α and M are the polarizability and mass of the particle, and E = E TL + E 1 + E 2 . For most large molecules and dielectric nanoparticles, the optical polarizability is approximately equal to the static polarizability. For atoms, the near-resonant enhancement in the optical polarizability can be used as a tuning parameter. Since the polarizations of the standing wave and the running waves are orthogonal E TL · E 1,2 = 0, the E 2 term can be expanded as E 2 = E 2 TL + E 2 1 + E 2 2 + 2E 1 · E 2 and the Hamiltonian in Eq. (37) can be rewritten aŝ The electric field of the trapping laser is The electric fields of the KD-lasers are where E 1,2 = 2I KD / 0 c. The KD-laser pulse duration τ KD controls the interaction time, and it can be set by an acoustic optical modulator. As in Eqs. (2) and (3), we approximate both the trapping and the KD-lasers to be homogeneous in the y-z plane, assuming that the length scale of intensity variation is much larger than the particle beam size along the same direction. We also assume that the KD-laser fields satisfy the criteria given in Eq. (4).
Assuming that the trapping laser is much stronger than the KD-lasers, |E TL | |E 1,2 |, the Hamiltonian can be separated into an unperturbed partĤ and an interaction partĤ Using Eq. (39), the unperturbed Hamiltonian can be written asĤ The first term in Eq. (43) is time-independent. Assuming ω TL to be far from any electronic resonance of the particle, the second term with cos (2ω TL t) can be discarded because it only causes non-resonant driving and the effect averages out. As a result, we obtain the dipole potential from Eq. (43), Close to the potential minimum x m = 0, the dipole potential can be approximated as a harmonic trap where a constant in the expansion is dropped out, and the trap frequency is Thus, the unperturbed Hamiltonian in Eq. (41) becomeŝ Since Eq. (47) is the same as Eq. (12), the general solution to the particles' wavefunction is also described by Eq. (13). We can examine the interaction between the KD-lasers and the oscillator by evaluating Eq. (42), which has the same form as Eq. (14). Comparing Eqs. (47) and (48) to Eqs. (12) and (14), we see that the only difference is a direct replacement of coefficients, Therefore, all the analysis and results from Eq. (15) to Eq. (36) apply for polarizable neutral particles.
(c) Extension to 2D and 3D harmonic traps We now extend our discussion to particles in 2D and 3D harmonic traps, which are common for atoms, ions, molecules, and dielectric nanoparticles. The analysis presented here focuses on 3D harmonic oscillators but the results can be adapted for 2D harmonic oscillators. Let us consider a non-degenerate 3D harmonic oscillator with trap frequencies (Ω x , Ω y , Ω z ) along the x-, y-, and z-axes. We extend the unperturbed Hamiltonian in Eq. (47) to 3D, The oscillator is assumed to be initially in the ground state |0 x |0 y |0 z , where |n x,y,z are the eigenstates of the 3D harmonic oscillator along the x-, y-, and z-axes. We define these three degrees of freedoms as the x-, y-, and z-oscillators hereafter. If the KD-lasers propagate along the x-axis, only the x-oscillator can be resonantly excited.
The general solution to the oscillator's wavefunction is where Ω n ≡ Ω x (n + 1/2). Since the dynamics of C n (t) pertains to only the x-oscillator, the analysis from the previous subsection applies here. When the KD-lasers travel in an oblique angle with respect to the potential axes, different degrees of freedom of the 3D harmonic oscillator can interact with the KD-laser fields at the same time and this leads to entanglement. As an example, we consider KD-laser fields traveling in the xy-plane with wavevectors k 1,2 = ± (cos (θ), sin (θ), 0) ω 1,2 /c, where π/2 ≥ θ ≥ 0 is the angle between the wavevector k 1,2 and the x-axis. A schematic is given in FIG. 3. The electric fields of the KD-lasers are modified from Eq. (40) to be where r = (x, y, z) is the 3D position operator of the particle. Because the interaction involve both the x-and y-oscillators, the general solution to the particle's wavefunction is modified from Eq. (51) to be where |m, n ≡ |m x |n y , Ω (m) x ≡ Ω x (m + 1/2), and Ω (n) y ≡ Ω y (n + 1/2). The Schrödinger equation for C m,n (t) is Using Eq. (52), the interaction Hamiltonian in Eq. (48) is modified to bê Assuming ω 1,2 to be many orders of magnitude higher than the trap frequencies Ω x,y,z , only the ω 1 − ω 2 term in Eq. (55) can resonantly excite the harmonic oscillator. Therefore, Eq. (55) can be simplified aŝ where cos ((ω 1 − ω 2 )t) and sin ((ω 1 − ω 2 )t) are expanded. The transition matrix element in Eq. (54) is thus m , n |Ĥ int |m, n = f α (t) m , n | cos ((k 1 + k 2 ) · r)|m, n + h α (t) m , n | sin ((k 1 + k 2 ) · r)|m, n , where Using Eq. (24), we can evaluate the matrix element where ψ n (k) and φ n (k) are the eigen-wavefunctions of the x-and y-oscillators in the momentum space, and (k 1 + k 2 ) x,y ≡ (k 1 + k 2 ) · x,y are the projection of k 1 + k 2 along the x-and y-axes. Using Eq. (59), we transform Eq. (57) into which is a generalization of Eq. (26). With Eqs. (54) and (60), we can summarize the conditions for resonant excitation in terms of energy-momentum conservation, where n x,y are positive integers, k x0 = M Ω x /2h, k y0 = M Ω y /2h, and δ px,py are the momentum detuning for the x-and the y-oscillator respectively. The central KD-laser frequency ω KD and the laser propagation angle θ can be obtained from Eq. (61), so the KD-laser frequencies are ω 1,2 = ω KD ± (n x Ω x + n y Ω y )/2. We can see from Eqs. for m > m and n > n. Note that the Lamb-Dicke parameters here are η x = (n x + δ px )/2 and η y = (n y + δ py )/2. Assuming m = m + n x and n = n + n y , the matrix element in Eq. (57) have the analytic forms m + n x , n + n y | cos ((k 1 + k 2 ) · r)|m, n = 1 where F (nx,ny) m,n (η x , η y ) ≡ m! (m + n x )! n! (n + n y )! (η x ) nx (η y ) ny e −(η 2 For ω 1 − ω 2 = n x Ω x + n y Ω y , the Schrödinger equation in Eq. (54) can be simplified to include only the resonant terms, ih dC m,n (t) dt = m, n|Ĥ int |m − n x , n − n y C m−nx,n−ny (t)e i(nxΩx+nyΩy)t + m, n|Ĥ int |m + n x , n + n y C m+nx,n+ny (t)e −i(nxΩx+nyΩy)t .
II. Wigner function and wave function in the discrete limit Given a wavefunction ψ(x) in the configuration space, the Wigner function is [4] W (q, p) = 1 2πh where q and p are position and momentum respectively. Defining k ≡ p/h and f q (x) ≡ ψ * (q − x/2) ψ (q + x/2), we can rewrite Eq. (75) as a Fourier transform of f q (x), As (x, k) is a conjugated pair in Fourier transform, the integral can be discretized using the relation where ∆n = 1, and the discreteness of x implies a finite size of the one-dimensional k-space L k = 2π/∆x with ∆x = x n+1 − x n . Using Eq. (77), we can transform Eq. (76) into W (q, k) = 1 hL k n f q (x n )e ikxn .
The Wigner functions presented in FIG. 4 and FIG. 5 of the main text are computed through Eq. (79). Additionally, we note that the probability distributions |ψ(x)| 2 and |ϕ(p)| 2 in the position and momentum space can be obtained as marginals of the Wigner function, where L p =hL k . In our numerical simulation of inelastic KD-effect, the time-dependent Schrödinger equations in Eq. (34) and Eq. (36) are solved for a N -state harmonic oscillator. The real and imaginary part of the Schrödinger equation are separated, so there are 2N real-number differential equations to be solved simultaneously. We employ the adaptive 5th order Cash-Karp Runge-Kutta method as the equation solver in our FORTRAN code [6]. The oscillator is assumed to be in the ground state initially. The integration time is 5 times longer than the 1/e pulse duration of the KD-lasers. The integration stepsize is 1/100 of the natural period 2π/Ω 0 . The most computationally intense simulation is for a nanoparticle oscillator with N = 6000 eigenstates (see FIG. 6(c) in the main text). The simulation took 400 GB memory and 90 hours of computation time on a single core of a supercomputer (Crane, Holland Computing Center). Once we obtain the probability amplitude C n (t) for each eigenstate, we compute the oscillator's wavefunction using MATLAB, and the wavefunction is subsequently used for computing the Wigner function in Eq. (79). We note that for states n > 170 the value of the normalization factor 1/2 n n! in the analytical formula of the oscillator eigen-wavefunction is too large for the computer to compute. To circumvent this problem, we solve the eigen-wavefunction numerically using the time-independent Schrödinger equation in the range of 0 ≤ x ≤ x t , where x t = (n + 1/2)2h/M Ω 0 is the turning point. In the tunneling regime x ≥ x t , we patch an ansatz function to the eigen-wavefunction where φ (x) ≡ dφ(x)/dx and x 0 = h/2M Ω 0 . The eigen-wavefunction in the range of x < 0 is obtained through symmetry φ(x < 0) = (−1) n φ(−x), where n is the quantum number of the eigenstate. The numerically solved eigenwavefunctions are in very good agreement with the analytic formula. The data size of the nanoparticle oscillator eigen-wavefunctions is 24 GB.