Floquet Many-body Engineering: Topological and Many-body Physics in Phase Space Lattices

Hamiltonians which are inaccessible in static systems can be engineered in periodically driven many-body systems, i.e., Floquet many-body systems. We propose to use interacting particles in a one-dimensional (1D) harmonic potential with periodic kicking to investigate two-dimensional (2D) topological and many-body physics. Depending on the driving parameters, the Floquet Hamiltonian of single kicked harmonic oscillator has various lattice structures in phase space. The noncommutative geometry of phase space gives rise to the topology of the system. We investigate the effective interactions of particles in phase space and find that the point-like contact interaction in quasi-1D real space becomes a long-rang Coulomb-like interaction in phase space, while the hardcore interaction in pure-1D real space becomes a confinement quark-like potential in phase space. We also find that the Floquet exchange interaction does not disappear even in the classical limit, and can be viewed as an effective long-range spin-spin interaction induced by collision. Our proposal may provide platforms to explore new physics and exotic phases by \textit{Floquet many-body engineering}.


I. INTRODUCTION
Since the concept of topological order was first introduced into condensed matter physics in 1973 [1], topological phenomena have been intensively investigated in the past decades. Today, topology lies at the heart of many research fields, e.g., quantum Hall physics [2], topological insulators/superconductors [3,4], and many more. The origin of topology in physics comes from the geometric phase factor of a quantum state when it moves along an enclosed path. In quantum Hall physics, the geometric phase is induced by the applied magnetic field and the resulting energy spectrum, also called the Hofstadter's butterfly [5], is a fractal; while the band can be characterized by its topological invariant (Chern number or TKNN invariant), which relates to the quantized Hall conductance directly [6]. In topological insulators/superconductors, the spin-orbit coupling takes the role of an effective magnetic field [7,8] resulting in the geometric quantum phase factor. For the ultracold atoms in optical lattice [9,10], the geometric quantum phase (Berry phase [11]) is generated by shaking the lattice, which creates an artificial gauge field [12][13][14][15][16][17].
An alternative way to study topological physics is employing the noncommutativity of phase space in quantum mechanics. In a noncommutative space, the concept of point is meaningless due to the commutative relationship * lzguo@tfp.uni-karlsruhe.de [X,P] = iλ. Instead, we should define a coherent state |α which is the eigenstate of the lowering operater, i.e., a|α = α|α withâ ≡ (X + iP)/ √ 2λ. As shown in Fig. 1, we observe that a coherent state moving along a closed path in phase space naturally acquires an additional quantum phase factor e iS /λ , where S is the enclosed area [18]. This observation reveals the origin of topology in the study of many dynamical systems, e.g., the kicked harmonic oscillator (KHO) [19][20][21][22] and the kicked Harper model (KHM) [23,24]. The energy spectra of these dynamical systems exhibit butterfly structure and band topology similar to quantum Hall systems [25,26]. In the strong kicking strength regime, the dynamical systems become chaotic and exhibits many novel behaviors such as dynamical localization, which has an intimate relation with the topology of bands [27][28][29].
In many-body physics of equilibrium systems, many exotic phases of matter emerge when interaction makes the system strongly correlated. It is the interplay between topology and interaction that gives rise to the fractional quantum Hall effect [30][31][32], and many other fascinating phenomena [33][34][35][36], like fractional charge and anyons [37][38][39][40][41][42]. Alternatively, it is also possible to engineer novel phases in periodically driven systems, i.e., the Floquet systems. The Hamiltonian of a Floquet system is a periodic function in time, i.e., H(t) = H(t + T ). The Floquet theory [43,44] allows us to describe stroboscopic time-evolution for every period by a time-independent Hamiltonian which is called the Floquet Hamiltonian H F arXiv:1710.09716v2 [quant-ph] 3 Nov 2017 and is defined by e − i T H F ≡ T e − i T 0 H(t)dt , or equivalently Here, T is the chosen stroboscopic time step and T is the time-ordering operator. Exotic Floquet Hamiltonians [45][46][47][48][49][50][51] which are inaccessible in static systems can be engineered from Eq. (1) and a range of novel physical phenomena, such as phase space crystals [52,53], Anderson localization (or many-body localization) in time domain [54][55][56][57] and spontaneous breaking of discrete timetranslation symmetry (Floquet time crystals) [58][59][60][61][62][63][64][65][66], can be created by Floquet engineering [67][68][69]. While most work focus on the single-particle physics of (dissipative) Floquet systems, the possible new physics by Floquet many-body engineering has become an active research direction in recent years. Unlike the static manybody systems, the generic nonintegrable Floquet manybody systems are expected to be heated up, by the driving field, to a trivial stationary state with infinite temperature [70][71][72]. However, before reaching the long-time featureless infinite-temperature state, there is a prethermal state with exponentially long lifetime for high driving frequencies, and therefore existing a prethermal dynamics which can be described by the time-independent Floquet Hamiltonian (1) [73][74][75][76][77][78][79][80][81][82]. By introducing disorder as in many-body localized systems [83] or coupling the Floquet many-body system to a cold bath [84], it is also possible to protect the metastable prethermal state.
In this paper, we investigate cold atoms trapped in 1D harmonic potential with a stroboscopically applied optical lattice. The equation of motion of a single atom corresponds to the kicked harmonic oscillator (KHO) and we find that that intriguing 2D topological and many-body physics emerges in phase space. The Floquet Hamiltonian of a single KHO, in the rotating wave approximation (RWA), forms various lattice structures in phase space depending on the driving parameters. The full dissipative quantum dynamics shows that the stationary state forms a lattice structure in phase space but with a finite size limited by the dissipation rate. Furthermore, we consider the interaction between cold atoms and find that the point-like contact interaction of cold atoms in real space becomes a long-range Coulomb-like interaction in phase space. More interestingly, the hard-core interaction of cold atoms in real space becomes a long-range potential which increases linearly with the distance in phase space, i.e., a quark-like confinement potential. We also find the Floquet exchange interaction has Coulomblike long-range behavior, which does not disappear in the classical limit and becomes an effective spin-spin interaction.

II. MODEL AND HAMILTONIAN
We start from interacting cold atoms trapped by an elongated three-dimensional harmonic potential, with the ) with two complex numbers ξ 1 , ξ 2 determining the moving path. The geometric phase factor is given by e i 1 λ S with S = 1 2 Im[ξ 2 ξ * 1 ] being the area of the enclosed path (blue area). radial motion cooled down to the ground state. In this way, the spatial motion of the atoms is restricted to the remaining axial direction. In general, the onedimensional system is described by where V(x i −x j ) is the two-body interaction, which is typically contact or hard-core interactions in the context of cold atoms [10,[85][86][87][88][89][90]. The H s (x i ,p i , t) is the single-particle Hamiltonian which can be explicitly timedependent. Here, the single-particle Hamiltonian is the quantum kicked harmonic oscillator, which is described by H s = 1 where ω 0 is the axial harmonic frequency and m is the atom's mass. The periodic term is implemented by a stroboscopic optical lattice, which can be created by two counter-propagating laser beams with off-resonant frequency far away from internal electronic transitions [9,10]. Parameters k and K 0 are the wave vector of the laser beams and the kicking strength, respectively. Parameter T d is the time period between adjacent kicking pluses. We scale the coordinate and momentum by the units of √ /(λmω 0 ) and √ m ω 0 /λ with the parameter λ ≡ k 2 /mω 0 , respectively. Finally, we have the dimensionless single-particle Hamiltonian scaled by ω 0 /λ where K = λK 0 / ω 0 is the dimensionless kicking strength, τ = ω 0 T d is the dimensionless kicking period and the time t has also been scaled by ω −1 0 . The commutation relationship between the coordinate and the momentum is now [x i ,p j ] = iλδ i j , where the dimensionless parameter λ plays the role of an effective Planck constant. Thus, the semiclassical regime corresponds to the limit λ → 0. Accordingly, the two-body interaction will be given by the new scaled dimensionless observables as V(x i −x j ). Our remaining paper is organized as follows. In Sec. III, we discuss the single-particle physics neglecting interaction of paticles. We first introduce the topological band theory of phase space lattices in Sec. III A. Then, in Sec III B, we investigate the dissipative quantum dynamics of a KHO in a realistic environment and show how a lattice structure is formed in phase space. In Sec. IV, we consider the interactions and investigate the many-body dynamics. We first develop a general theory of transforming a given real space interaction potential to a phase space interaction potential in Sec IV A. Then, in Sec IV B, we apply our theory of phase space interaction to the special cases of contact interaction and hard-core interaction of cold atoms, and give the analytical expressions of corresponding phase space interactions. In Sec IV C and IV D, we investigate the many-body dynamics in the classical limit and discuss the concept of dynamical crystals. Finally, we summarize our results in Sec. V.

III. PHASE SPACE LATTICES
In this section, we investigate the single-particle Hamiltonian of the quantum KHO, i.e., Eq. (3), in the resonant condition that the kicking period satisfies τ = 2π/q 0 with q 0 an integer. When the kicking strength is weak |K| 1, the single-particle dynamics is still dominated by the fast harmonic oscillation. Then we transform into an appropriately chosen rotating frame generated by the free time- 2λ. We transform the coordinates and momenta of particles by Here, the operatorsX i andP i describe the dynamics of the i-th atom's phase and amplitude. For the harmonic oscillator,X i andP i are fixed and correspond to the initial state ofx i (t) andp i (t). In our case, however, the phase and amplitude of KHO are slightly changed every harmonic time period due to the weak kicking. The time-evolution ofX i andP i is slow compared to the fast global harmonic oscillation and can be obtained stroboscopically from the time-evolution ofx i (t) andp i (t) every time period of 2π.
The canonical transformation of the single-particle Hamiltonian is given byÔ(t)H sÔ In RWA, we drop the fast oscillating terms and arrive at the timeindependent Hamiltonian (see detailed derivation in Appendix A) Here, we have dropped the index of the operators since we are considering single-particle physics. Another way of deriving H RWA (X,P) is based on the series expansion of the Floquet Hamiltonian (1) in order of the kicking strength K. By replacing the Planck constant by a dimensionless one λ and choosing the stroboscopic time step T = q 0 τ, we have the time-evolution operator in one stroboscopic time step T e − i λ q 0 τ 0 H s (t)dt =F q 0 with the Floquet operator for one periodF ≡ e −i(x 2 +p 2 )τ/2λ e iK cosx/λ . In the Appendix A, we show that H RWA (X,P) is indeed the first order expansion of the Floquet operatorF q 0 .
To display the symmetries of H RWA (X,P) in phase space, we calculate the averaged H RWA (X,P) in the coherent state representation (see the details in Appendix B), i.e., α|H RWA (X,P)|α = e −λ/4 H RWA (X, P).
Here, the coherent state |α is defined by the eigenstate of the annihilation operatorâ ≡ (X + iP)/ √ 2λ, i.e., a|α = α|α . The averaged position and momentum are . The quantity H RWA (X, P) has the same expression as Eq. (5) but replacing the operatorsX andP by the averaged values X and P respectively. In Fig. 2, we plot H RWA (X, P) in phase space for different q 0 . We see that the H RWA (X, P) has a square lattice structure for q 0 = 4, hexagonal lattice structure for q 0 = 3 or q 0 = 6, and even quasicrystal structure for q 0 = 5 or q 0 ≥ 7. The translational symmetry of the Hamiltonian (5) in phase space gives rise to the band structure of its spectrum.

A. Band Structure and Topology
We will deal with the case of square lattice (q 0 = 4) in detail but the results can be readily generalized to the case of hexagonal lattice (q 0 = 3 or 6). For q 0 = 4, the effective Hamiltonian (5) is further simplified as This Hamiltonian is closely related to the established Harper's equation, which is a tight binding model governing the motion of noninteracting electrons in the presence of a two-dimensional periodic potential and a uniformly threading magnetic field [5,91]. The H sq (X,P) is invariant under discrete translation in phase space by two op-eratorsT 1 ≡ e i 2π λP andT 2 ≡ e i 2π λX , i.e., The translation operatorsT 1 andT 2 generate an invariance group G of H sq [92], which is a nonabelian group due to the identity [T r 1 ,T s 2 ] =T 1T2 (1 − e −i4π 2 rs/λ ) with integer powers r, s ∈ Z. However, the group G has abelian subgroups G a generated byT r 1 andT s 2 if 2πrs/λ ∈ Z, which means the value of the parameter λ/2π needs to be a rational number, i.e., λ/2π = p/q, where p and q are coprime integers. Here, we choose the abelian subgroup G a generated by the following two generators (r = 1, s = p) Therefore, we can find the common eigenstates of commutative operatorsT X andT P with eigenvalues given by e i2πk X and e i2πpk P respectively. The boundaries of the two dimensional Brillouin zone are defined by 0 ≤ k X ≤ 1 and 0 ≤ k P ≤ 1/p, where k X and k P are quasimomentum and quasicoordinate, respectively [93]. The corresponding eigenvalues of the Hamiltonian H sq are also called quasienergies.
The discrete translational symmetry in phase space allows us to determine the quasienergy spectrum numerically in Zak's kq-representation (see the instruction in Appendix C or Ref. [93]). Given the parameters of λ/2π = p/q, k X and k P , the eigenvalues E of H sq are determined by the following polynomial equation (see the derivation in Appendix D or Ref. [94]) The left hand side of Eq. (9) takes values in the range [−2, 2] when the quasimomentum k X and quasicoordinator k P run over the whole Brillouin zone. The right hand side of Eq. (9) is a periodic function of λ with period 2π. Therefore, the quasienergy spectrum is also a periodic function of λ with period 2π. In Fig. 3(a), we plot the quaienergy spectrum for λ/2π ∈ [0, 1], showing a Hofstadter butterfly structure identical to that in quantum Hall systems. In Fig. 3(b), (d) and (e), we plot the quasienergy band structures in the two-dimensional Brillouin zone (k X , k p ) for λ/2π = 1/2, 1/3 and 2/3, respectively. For the given parameter λ/2π = 1/2, we can obtain the analytical solutions from Eq. (9), i.e., The two bands touch each other at the central point of the Brillouin zone, i.e., (k X = 1 2 , k P = 1 2 ), where the dispersion relationship becomes linear near the touching point, i.e., E ≈ ± πK √ 2 |k| with |k| ≡ (k X − 1 2 ) 2 + (k P − 1 2 ) 2 , as shown in Fig. 3(c). In general, the two innermost bands always touch each other for even integer q. We also see that the quasienergy band structure is two-fold degenerate for λ/2π = 2/3 while there is no degeneracy for λ/2π = 1/2 and λ/2π = 1/3. In fact, for each rational λ/2π = p/q (remembering p, q are coprime integers), the spectrum contains q bands and each band has a p-fold degeneracy due to the fact that the invariance group G can be expressed as the coset sum p r=1T r 1 G a [92]. We denote the quasienergy states by |ψ b,k with k ≡ (k X , k P ) and b the band index counting from the bottom. To visualize the quasienergy states, we define the Husimi Q-function of a given eigenstate in phase space [95] where |α is the coherent state introduced at the beginning in this section. In Figs. 4(a) and (b), we plot the Q-functions of eigenstates |ψ 1,(0,0) and |ψ 2,(0,0) for λ 2π = 1 2 , which are the ground-like states of the lower band and upper band, respectively. Comparing the Q-functions of the two states to the phase space lattices shown in Fig. 2(a), we see that the Q-function of eigenstate |ψ 1,(0,0) mostly occupies the negative phase space lattice H sq (X, P) < 0 while the Q-function of eigenstate |ψ 2,(0,0) is shifted by π along both X and P directions in phase space, mostly occupying the positive phase space lattice H sq (X, P) > 0. In fact, our system has a chiral symmetry defined by the chiral operatorT c ≡ e i π λX e i π λP , i.e.,T c H sqT † c = −H sq . Thus, for a given eigenstate |ψ b,k , there must be another eigen-stateT † c |ψ b,k with opposite quasienergy. In Figs. 4(c) and (d), we plot the Q-functions of eigenstates |ψ 1,( 1 4 , 1 2 ) and |ψ 1,( 3 4 , 1 2 ) for λ 2π = 2 3 , which are the degenerate states of the lower band shown in Fig. 3(e). We see that the period along X-direction is 2π while the period along Pdirection is 4π. In fact, this degeneracy depends on the discrete translation operators we choose in Eq. (8). For λ 2π = p q , the period of any Q-function of eigenstate is 2πp-period in P-direction and 2π-period in X-direction. Theses p-degenerate states are given by with the same quasicoordinator but different quasimomenta. In the case of λ 2π = 2 3 , the two degenerate states of |ψ b,k andT 2 |ψ b,k in the X-representation are Therefore, the Q-functions of the two degenerate states are the same in the X-dimension. But due to the relationship ψ b,k |P|ψ b,k = T 2 ψ b,k |P|T 2 ψ b,k + 2π, the Q-functions are shifted by 2π in the P-dimension.
The underlying topology of a quasienergy band is de- fined by the Chern number [96] where the contour C is integrated over the boundary of the Brillouin zone. The Chern number associated with a gap is subtle here. For the equilibrium systems, the Chern number of a gap is defined by the sum of the Chern numbers of the energy bands below the gap. However, in our present work, we are dealing with a Floquet system far from equilibrium. The general statistic law of the Floquet states for the long-time stationary state is an ongoing research topic [97][98][99]. Actually, the positive and negative sublattices shown in Fig. 2(a) make no difference in the frame of Floquet theory. We assume that the statistic mechanics near the ground state of each sublattice can be described by an effective Floquet-Gibbs states [99]. Therefore, we define the Chern number of a gap below (above) the zero energy line as the sum of Chern numbers of all the quasienergy bands below (above) the gap. As shown in Fig. 3(a), the Chern number of some gaps are calculated and labelled symmetrically with respect to zero energy line.

B. Full Dissipative Quantum Dynamics
In the above section, our analysis is based on the rotating wave approximation where the kicking strength needs to be weak |K| 1. In this section, we will investigate the full quantum dynamics of KHO based on the original full Hamiltonian (3) and confirm the validity of the rotating wave approximation, which is used to derive the effective Hamiltonian (6). From a practical point of view, the oscillators are inevitably in contact with the environment, which is conventionally modeled by a harmonic bath model. The coupling with the environment results in dissipation or decoherence of the quantum system. Here, we describe the dissipative dynamics of the quantum KHO by the following master equation, where κ characterizes the dissipation rate and n 0 is the Bose-Einstein distribution of the thermal bath. The dissipative dynamics is described by the Lindblad superoperator defined by D[Ô]ρ ≡ÔρÔ † − 1 2 (Ô †Ô ρ+ρÔ †Ô ), wherê O is an arbitrary operator. The two Lindblad terms in Eq. (12) represents relaxation and heating processes respectively. We notice that some authors also choose the non-Lindblad Caldeira-Leggett master equation to describe the dissipative dynamics [96]. Here, we choose the Lindblad master equation (12) since it can give the correct thermal equilibrium state of harmonic oscillator without kicking force while the non-Lindblad Caldeira-Leggett master equation cannot [100].
As the kicks act as delta-functions, we can separate the dissipative dynamics from the kicking dynamics. In order to solve the dissipative dynamics, we define the characteristic function of the Wigner distribution by [95] w(s, k) ≡ dx e ixk/λ x + s 2 |ρ|x − s 2 , Then the master equation (12) without kicking can be transformed into the following Fokker-Planck equation [96] ∂ t w+( The dissipative dynamics between two successive kicks is solvable from the above Fokker-Planck equation. Given the initial state at the moment right after n − 1 kicks w(s, k; τ + n−1 ), where τ + n−1 = (n − 1)τ + ∆ with a positive infinitesimal increment ∆, the final state at the moment right before n kicks w(s, k; τ − n ) with τ − n = nτ − ∆, is given by the following map [96] w(s, k; with s r = e − κτ 2 (k sin τ+s cos τ) and k r = e − κτ 2 (k cos τ−s sin τ). The kicking dynamics is an instantaneous unitary transformation ρ → ρ =Û K ρÛ † K withÛ K ≡ e −iKτ/λ cosX . In Appendix E, we prove that the corresponding map of the The initial state is a coherent state with its center locating on (0, 0) and (−π, π) as marked by two dashed circles in figures (a) and (b) respectively. The Q-functions are plotted after 3000 kicks. Other parameters: characteristic function of the Wigner distribution at the time τ n = nτ is given by where the J j are the jth-order cylindrical Bessel function. Hence, the full dynamics of the quantum KHO in contact with a thermal bath is realized by applying the two maps (14) and (15) sequentially. From the characteristic function w(s, k), it is direct to obtain the corresponding Husimi-Q function (see Appendix F). In Fig. 5(a), we evolve the dynamics of the system starting from the ground state of the harmonic oscillator. We then plot the Husimi Q-functions of the states after 1000 and 3000 kicks in Fig. 5(b) and Fig. 5(c) respectively. We see clearly that a final state with square lattice structure in phase space forms gradually revealing the underlying square structure of Hamiltonian (6). Interestingly, we find that the transient state shown in Fig. 5(b) has no reflection symmetries with respect to X and P although the RWA Hamiltonian (6) has. There is a chiral feature as marked by the dashed lines along the backbone of the quasiprobability distribution. This chirality is a reflection of the topological property of our system and the noncommutative geometry [101,102] of the phase space. As approaching the stationary state, the chirality disappears in the end.
Without dissipation, the quantum KHO will experience unbounded diffusion for resonant condition, where the average energy of the harmonic oscillator increases infinitely due to the energy pump from kicking [25]. When dissipation is present, the diffusion process approaches a nonequilibrium stationary state with a finite size in phase space depending on the driving strength and dissipation rate. In Fig. 6, we plot the average energy of the KHO λ â †â = X 2 +P 2 /2 as a function of the kicking number for different dissipation rates. We see that the smaller the dissipation rate is, the lager the phase space lattice is in the long-time limit. If the dissipation rate is so strong that the system can relax to its ground state during the successive kicks, the lattice state cannot be formed in phase space. Therefore, in order to create a phase space lattice with enough large size, the dissipation rate has to be much weaker than the kicking strength, i.e., κ |K|.
In Figs. 2(a) and (b), we also notice that there are actually two identical square lattices with a relative shift in phase space, which support eigenstates with positive and negative quasienergy respectively. In Figs. 7(a) and (b), we plot the two Husimi Q-functions evolving from two coherent states with different initial positions in phase space, i.e., (X, P) = (0, 0) and (X, P) = (−π, π) respectively. We see that a state initially prepared on one sublattice stays on that lattice during the evolution and has negligible occupation on the other sublattice. This is different from the static potential, where the minimum points correspond to stable state while maximum points correspond to unstable state. Since we are working on a dynamical system far from equilibrium, both minimum and maximum points of the Hamiltonian in phase space are stable; only the saddle points are unstable. This is the reason why we define the Chern number of the gaps symmetrically with respect to the zero line for the Hofstadter's spectrum in Fig. 3(a).

IV. MANY-BODY DYNAMICS
In the above discussion, we have neglected the interaction terms in the original Hamiltonian (2). From this section, we will consider the interactions between particles. Using the free time-evolution operatorÔ(t) = exp(iλ iâ † iâ i t) defined at the beginning in Sec. III, the total Hamiltonian in the rotating frame is given by the canonical transformation, i.e., In the RWA, we drop the fast oscillating terms and arrive at the time-independent Hamiltonian Here, H RWA (X i ,P i ) is the single-particle RWA Hamiltonian given by Eq. (5). The RWA interaction potential In general, U(X i ,P i ;X j ,P j ) is defined in phase space and depends on both coordinates and momenta of two particles. Thus, we call U(X i ,P i ;X j ,P j ) the phase space interaction potential. We aim to determine the explicit form of U(X i ,P i ;X j ,P j ) in this section.

A. Phase Space Interaction Potential
For two arbitrary particles, we introduce the operatorŝ X c ≡ (X 1 +X 2 )/2,P c ≡ (P 1 +P 2 )/2 representing the coordinator and momentum of two particles' center of mass, and the operators ∆X ≡X 1 −X 2 , ∆P ≡P 1 −P 2 representing their relative displacement in phase space. We further define the operator of phase space distance bŷ It is important to notice that the background of the phase space interaction potential is a noncommutative space. From the commutation relationship [X i ,P j ] = iλδ i j , we have [∆X, ∆P] = i(2λ), and [R 2 ,X c ] = 0 which means the motion of two particles' center of mass and their relative motion are independent. Thus, we write the common eigenstate of commutative operatorsR 2 andX c as a product state Ψ (X 1 , where the wave function f (X c ) is the state of two particles' center of mass and the wave function Φ(X 1 − X 2 ) describes their relative motion. Reminiscent of the Hamiltonian operator of a harmonic oscillator, the eigenvalues of operator R 2 are given by 4λ(N +1/2) with N = 0, 1, 2, ···. Therefore, the eigenvalues of the operatorR are given by For each N, the corresponding eigenstate is given by where H N (•) is the Hermite polynomial of degree N. We choose functions δ(X c − C), i.e., the eigenstate of oper-atorX c , as the basis of two particles' center of mass. Therefore, we use the Dirac notation |N, C to represent the total eigenstate, which is determined by two good quantum numbers C and R N , i.e.,X c |N, C = C|N, C and R|N, C = R N |N, C . In the coordinate representation, the total eigenstate has the explicit form X 1 , There is a fundamental difference between the commutative real space and the noncommutative phase space. The concept of point is meaningless in noncommutative space. Instead, we are only allowed to define the coherent state |α as the point in noncommutative geometry. Similarly, the concept of distance also needs to be reexamined. The distance of two particles in real space is a continuous variable from zero to infinity. However, the distance in phase space is a quantized variable and has a lower limit √ 2λ as seen from Eq. (18). Here, we actually provide a description for the quantization of the noncommutative background.
We now start to determine the phase space interaction potential U(X 1 ,P 1 ;X 2 ,P 2 ). From the transformation (4), the relative displacement of two particles in the rotating frame isÔ(t)(x 1 −x 2 )Ô † (t) = ∆P sin t + ∆X cos t. Therefore, for a given real space interaction potential, we haveÔV(x 1 −x 2 )Ô † = +∞ −∞ dqV qQ with the Fourier coefficients V q = 1 2π +∞ −∞ dxV(x)e −iqx and the operator Q ≡ exp iq(∆P sin t + ∆X cos t) . The matrix element of the operatorQ in theR-representation is given by the Laguerre polynomials [53,103] where |N and |M are the eigenstates of the operatorR given by Eq. (19). In the RWA, we only keep the timeindependent diagonal elements of the matrix (20), i.e., N|Q|N with N = 0, 1, 2, ···. Thus, given an arbitrary real space interaction potential V(x 1 − x 2 ), we find a compact expression for the phase space interaction potential In the eigenbasis |N, C , we have Here, the interaction potential U(R) takes the value U(R N ) with R N = 2 √ λ(N + 1/2). If the two paricles have spins, their spatial state is either antisymmetric or symmetric depending on the symmetry of total spin state, i.e., We sperate the average phase space interaction potential by U ± = ψ ± (X 1 , X 2 )|U(R)|ψ ± (X 1 , X 2 ) ≡ U c ± U e . Here, we have defined the direct interaction U c and the the exchange interaction U e , respectively, The direct interaction part U c = 1 2 ( U + + U − ) corresponds to the classical interaction while the exchange interaction part U e = 1 2 ( U + − U − ) is a pure quantum effect without classical counterpart, which we call U e the Floquet exchange interaction for our system. In theR-representation, they have been calculated in the Appendix G, i.e., with the overlap integral In the Appendix H, we have given the overlap integral I N (R) for the two displaced coherent states ϕ(X) = 1 √ πλ 1 2 e − 1 2λ (X−R/2) 2 and φ(X) = 1 √ πλ 1 2 e − 1 2λ (X+R/2) 2 , where R is the distance between the centers of two coherent states in phase space. Below, we will calculate the analytical expressions of U c (R) and U e (R) for contact and hardcore interactions of ultracold atoms.

B. Applications
In this section, we apply our general theory of phase space interaction to the special cases of contact interaction and hardcore interaction for ultracold atoms. We show that, in quasi-1D, the point-like contact interaction in real space becomes a long-range Coulomb-like interaction in phase space. In pure 1D, the hardcore interaction in real space produces a quark-like confinement interaction potential in phase space, which increases linearly with the phase space distance of two atoms.

Contact Interaction
In the experiments, the ultracold atoms are confined in one dimension if the transverse trapping frequency ω ⊥ is much larger than the longitudinal trapping frequency ω z . If the characteristic length of transverse trapping l ⊥ ≡ √ /(mω ⊥ ) is much larger than the cold atom's size, i.e., in the quasi-1D, the effective interaction between cold atoms is described by the contact interaction V(x 1 − x 2 ) = εδ(x 1 − x 2 ), where ε is the interaction strength [10,85,86].
and Eq. (21), the phase space interaction potential can be calculated .
Here, N = 0, 1, 2, 3 · ·· and Γ(•) is the gamma function. We see that U(R N ) is zero for odd integer N and finite for even integer N. The wave function of two atoms' relative motion Φ N (∆X) is antisymmetric for odd N, which means the probability amplitude is zero when the two atoms contact each other. The result is that the total average interaction of Φ N (∆X) is zero for odd N. The direct phase space interaction U c (R) and Floquet exchange interaction U e (R) of two cold atoms, which are described by two coherent states, can be calculated from Eq.s (23), (25) and (H8) Here, R is the distance in phase space between the centers of two coherent states and I 0 (•) is the zeroth order modified Bessel function of the first kind. In the large distance limit, we use the asymptotic behavior of the special function I 0 (z) ∼ e z / √ 2πz for z 1 and have In Fig. 8(a), we plot the U c (R) as function of R and its long-range asymptotic behavior. We see that a pointlike contact interaction indeed becomes a long-range Coulomb-like interaction in the long-distance limit, which is consistent with the pure classical analysis [51]. As shown in Eq. (26), we also find that the Floquet exchange interaction is equal to the direct phase space interaction, i.e., U e (R) = U c (R), and does not disappear even in the classical limit λ → 0, which cannot happen in a static system. Usually, the effective spin-spin interaction in Heisenberg model comes from the quantum exchange interaction between the nearest-neighbouring electrons and cannot be explained by classical dynamics. One should always keep in mind that we are investigating the effective stroboscopic dynamics and the two atoms indeed collide with each other during every stroboscopic time step. The phase space interaction is actually the time-averaged real space interaction in one harmonic period. The spin-spin interaction in Heisenberg model is a short-range interaction due to the exponentially small wave function overlap of two next-nearest-neighbouring electrons. However, here in our system, the Floquet exchange interaction has long-range behavior following Coulomb's law. In the classical limit λ → 0, the long-range Floquet exchange interaction can be viewed as an effective long-range spin-spin interaction induced by collision of two atoms. The equality U e (R) = U c (R) comes from the δ function modelling the contact interaction and the fact that the spatial antisymmetric state of two atoms has zero probability to touch each other. If the interaction potential between cold atoms is different from the δ-function model, it is possible to tune the phase space interaction U c (R) and collision-induced spinspin interaction U e (R) independently in the experiments.

Hardcore Interaction
If the characteristic length of transverse trapping l ⊥ is much smaller than the cold atom's size, which is called pure-1D, the contact interaction is no longer valid for the description of interaction between cold atoms [87][88][89][90]. In this situation, the atom can be viewed as a hardcore particle with a radius a, which means the interaction potential between the two atoms is infinite when their distance is smaller than 2a and zero when the distance is larger than 2a. Our theory of phase space interaction can be applied to the small hardcore limit a √ λ. Since the two atoms can not contact each other due to the hardcore interaction, the engenstates of phase space distance operatorR have to be zero at zero distance, which means that only the odd eigenstates Φ 2m+1 (∆X) with m ∈ N satisfy this condition. The even eigenstate should be reconstructed as Φ 2m (∆X) = sgn(∆X)Φ 2m+1 (∆X) with sgn(•) the sign function. The eigenstates Φ 2m (∆X) and Φ 2m+1 (∆X) are degenerate with the same eigenvalue R 2m = R 2m+1 = 2 √ λ(2m + 3/2). In the Appendix I, we calculate the phase space interaction potential of the hardcore interaction potential for odd integers N = 2m + 1 Here, [ N 2 ] means the closest integer number less than N 2 . For even integers N = 2m, we have U(R 2m ) = U(R 2m+1 ).
Here, we find that U(R N ) can be approximated very well by the linear relationship U(R N ) ≈ 2aπ −1 R N .
The direct phase space interaction U c (R) and the Floquet exchange interaction U e (R) of two coherent states can be calculated from Eq.s (23), (28) Here, R is phase space distance between the centers of two coherent states and the overlap integral I 2m+1 is given by Eq. (H8). The zero Floquet exchange interaction comes from the degeneracy of the symmetric and antisymmetric states. Thus, there is no collision-induced spin-spin interaction for the hardcore interaction. Using the linear approximation U(R N ) ≈ 2aπ −1 R N and Eq. (29), we have In the long-distance limit, we have the asymptotic expression of Eq. (30), i.e., This is consistent again with the classical analysis [51]. In Fig. 8(b), we plot the direct phase space interaction potential U c (R) as a functions of phase space distance R. We see that the linear relationship (31) (blue dashed line) is a very good approximation of Eq. (29) (black solid curve) and Eq. (30) (red dotted-dashed curve). It is interesting to find that the linear phase space interaction potential (31) mimics the interaction potential between quarks in QCD [104,105]. Actually, this surprising behavior of hardcore atoms can be understood in a simple picture. Since the two atoms have a tiny hardcore radius a, they prefer to oscillate in a synchronized way, i.e., in phase. If the atoms are out of phase due to the finite phase space distance R, they are more likely to collide with each other during the oscillation. The collision effect becomes stronger as the phase space distance R is larger, resulting in a confinement potential in the end.

C. Classical Many-body Dynamics
Although it is very difficult to numerically simulate the quantum many-body dynamics from the original many-body Hamiltonian (2), we can simulate the classical many-body dynamics and verify our theory of phase space interaction. From now on, we consider the classical dynamics of spinless atoms and replace all the operators by their corresponding classical quantities. The time evolution of the original coordinates, x i (t) and p i (t), of a single atom are given by the canonical equations of motion (EOM) from Eq. (2) As seen from Eq. (4), the values of X i (t) and P i (t) can be obtained from the time evolution of x i (t) and p i (t) stroboscopically every time period of 2π. In this sense, the X i (t) and P i (t) define the time evolution of the amplitude and phase of an oscillating atom in the discrete time domain t = 2πm with m = 0, 1, 2, · · ·. This method is called Poincaré map [106,107].
In the rotating frame, we write the RWA many-body Hamiltonian explicitly for q 0 = 4 where R i j = (X i − X j ) 2 + (P i − P j ) 2 is the classical phase space distance of two arbitrary atoms. Depending on the original interaction potential, the phase space interaction potential U c (R i j ) takes the form of either Eq. (27) or Eq. (31). The EOM of X i (t) and P i (t) is described by dX i /dt = ∂H T RWA /∂P i , dP i /dt = −∂H T RWA /∂X i . Using Eq. (33), we have the explicit form of EOM Using the above two methods, we can calculate the dynamics of many interacting atoms, and compare them to verify our phase space interaction theory. In Fig. 9, we investigate the dynamics of three interacting particles. For convenience, we introduce the complex position of the jth atom in phase space Z j (t) ≡ X j (t) + iP j (t). As shown in Fig. 9(a), we set the three atoms initially at the local equilibrium points of single particle Hamiltonian, i.e., Z 1 (0) = 2πi, Z 2 (0) = −2π − 2πi and Z 3 (0) = 2π − 2πi. If the displacement of each atom in phase space is small, i.e., |Z i (t) − Z i (0)| 1, we can linearize the EOM (34), and have the following solution where e ji ≡ [Z j (0) − Z i (0)]/|Z j (0) − Z i (0)| is the unit vector directing from ith atom's initial position to jth atom's initial position in phase space as shown by the black arrows in Fig. 9(a). The above linear solution indicates that each atom oscillates harmonically around a shifted equilibrium point Z i (0) − 2K −1 j i (dU c /dR i j )e ji and with the amplitude ∆Z i = 2K −1 | j i (dU c /dR i j )e ji |.
In Figs. 9(b)-(d), we show the three-body dynamics with the real space contact interaction potential The corresponding phase space interaction potential is given by U c (R i j ) = ε π 1 R i j . We plot the phase space trajectories of the first atom for different interaction strengths ε = 0.194 and ε = 0.775 in Fig. 9(b), and the corresponding time evolutions of positions and momenta in Figs. 9(c) and (d). We see that the dynamics given by the three methods agree with each other very well for weak interaction ε = 0.194 as shown in Fig. 9(c). However, for a larger interaction ε = 0.775, the linear solution (35) breaks down while the RWA EOM (34) is still a very good approximation as shown in Fig. 9(d).
In Figs. 9(e) and (f), we show the three-body dynamics with the hardcore interaction, which is mod- eled approximately by an inverse power-law potential n with a high power n = 20 in our numerical simulations. The corresponding phase space interaction potential is given by U c (R i j ) = 2aπ −1 R i j . As shown in Fig. 9(e), for a small hardcore radius a = 0.05, the phase space interaction is already strong enough to make the three particles overcome the potential barrier of phase space lattice and exhibit global motions. In Fig. 9(f), we compare the results from Poincaré map (black dots) and the RWA EOM (34) (red lines), which agree with each other very well.

D. Dynamical Crystals
In Fig. 10, we show the dynamics of seven interacting atoms for contact interaction with ε = 0.194. The seven atoms are located initially at the equilibrium points with zero momenta as shown by the seven red dots. It can be seen that the two atoms at the ends oscillate with the largest amplitude. If the interaction is weak enough, the atoms only oscillate locally around their equilibrium points. If the interaction is strong enough, the two edge atoms can overcome the potential barrier of the phase space lattice, and destroy the crystal state. The existence of the crystal state is guaranteed by the condition that the oscillating amplitude of the edge atom ∆Z edge is smaller than the radius of the dashed circle indicated in Fig. 10, i.e., ∆Z edge < ( We can estimate the critical condition from the linear solution (35). For the contact interaction, the oscillating amplitude of the edge atom converges for infinite atoms, i.e., ∆Z edge = ε 12πK . Therefore, the critical interaction strength for the existence of 1D crystal state in phase space is given by For hardcore interaction with U c (R i j ) = 2aπ −1 R i j , the oscillating amplitudes of the two edge atoms can be estimated by ∆Z edge = 4a πK (N − 1), where N is the number of atoms. The oscillating amplitude ∆Z edge increases linearly with the number of atoms, which means it is impossible to create an infinitely long 1D crystal state with hardcore interaction. For a given kicking strength K and hardcore radius a, the critical atom number for the existence of 1D crystal state is We call the stable crystal state in phase space formed by many atoms, the dynamical crystals.
One should distinguish the concepts of phase space lattice discussed in Sec. III and the dynamical crystal introduced here. Phase space lattice refers to the periodic structure in phase space of the single-particle Hamiltonian (5) without consideration of atomic interaction, while the dynamical crystal refers to the many-body state formed by interacting atoms. In the experiments, the dynamical crystal can be realized by two basic steps: first, prepare the initial state of atoms via applying a very strong static optical lattice; then, suddenly turn off the strong static optical lattice and add a weak optical lattice stroboscopically.
If the atoms have spins and tightly bound at the their fixed points by the phase space lattice, the direct phase space interaction U c (R) does not play a role in the dynamics. However, as discussed in Sec. IV B 1, the contact interaction can induce a Coulomb-like Floquet exchange interaction U e (R) = ε π 1 R , and thus the system shown in Fig. 10 can be modelled by a 1D spin chains with isotropic spin-spin interaction. The famous Mermin-Wagner theorem claims that, at any nonzero temperature, a one-or two-dimensional isotropic Heisenberg model with finiterange exchange interaction can be neither ferromagnetic nor antiferromagnetic [108]. This theorem clearly excludes a variety of types of long-range ordering in low dimensions, and is crucial to the search for low-dimensional magnetic materials in the recent years [109][110][111]. Here in our model, the collision-induced spin-spin interaction has a Coulomb-like long-range behavior, which is beyond the definition of finite-range interaction in Mermin-Wagner theorem [112]. Hence, the dynamical crystals actually provide a possible platform to test the Mermin-Wagner theorem and search for other new phenomena with longrange interactions such as causality and quantum criticality [113][114][115][116][117][118][119], nonlocal order [120][121][122], etc.

V. SUMMARY
In summary, we have studied the possibility to create new physics by Floquet many-body engineering in the dynamical system of kicked interacting particles in 1D harmonic potential. Our system exhibits surprisingly rich topological and many-body physics in 2D phase space. In the weak kicking strength regime K 1, the singleparticle RWA Hamiltonian has various lattice structures in phase space depending on the kicking period. The topological physics comes from the noncommutative geometry of phase space, which naturally provides a geometric quantum phase. We analyzed the topological quasienergy band structure of the square phase space lattice. We investigated the full dissipative quantum dynamics of a single kicked harmonic oscillator using master equation and the Fokker-Planck equation. The time evolution of the Husimi Q-functions confirms that the nonequilibrium stationary state is indeed a lattice state in phase space, but has a finite size due to the dissipation.
For the many-body dynamics, we made several findings and predictions based on the theory of phase space interaction potential. We found that the original contact interaction becomes a long-range Coulomb-like interaction in phase space, and the hardcore interaction becomes a quark-like confinement interaction in phase space. For the contact interaction, we predicted that the long-range Floquet exchange interaction does not disappear even in the classical limit, and can be viewed as collision-induced spin-spin interaction. We investigated the classical many-body dynamics and proposed the concept of dynamical crystals. We found that the contact interaction can create an infinitely long 1D dynamical crystal but the hardcore interaction cannot.
Finally, we point out that our method can increase the speed of numerical simulation significantly. For example, in simulating the dynamics of seven interacting atoms in Sec. IV D, the method of Poincaré map based on the original Hamiltonian costs more than ten hours using Wolfram Mathematica while the method based on the phase space interaction only needs one second. The reason is that our method only needs to calculate the dynamics on the stroboscopic time points by averaging the dynamics between stroboscopic steps using the phase space interaction potential.
Recently, we learned of a related study [123] in which the authors also discussed that the contact interactions between atoms can result in exotic long-range interactions in the effective description of the resonantly driven many-body system.
where we have neglected the particle index and simplified the notation of the summation in the kicking part. As discussed in the main text, we transform Eq. (A1) into the rotating frame by employing the unitary transforma-tionÔ = exp i(â †â + 1/2)t and using the relationship (4) where we defineM(X,P) = exp i X cos t +P sin t . The harmonic term in Eq. (A2) disappears due to the resonant condition.
The element ofM(X,P) in the basis of Fock states is evaluated to be [53,103] l|M(X,P) Here the L k−l l are the generalized Laguerre polynomials. Inserting Eq. (A3) intoĤ RF we have The sum of the Dirac δ-functions in Eq. (A4) obeys the following identity [19] n δ( t τ − n) = n cos 2πnt τ . (A5) Making use of this relation and dropping all terms relevant to t inĤ RF (rotating wave approximation) we get The sum over n can be formulated as the sum over k = l + nq 0 with the help of the formula that is, Using Eq. (A3), we have the final effective Hamiltonian e iX cos(2π j/q 0 )+iP sin(2π j/q 0 ) + h.c.
Another way to derive the effective Hamiltnoian Eq. (5) is to start from the Floquet operator in one harmonic oscillation Following the same procedure in [25] , we can reformulate Eq. (A9) aŝ ExpandingF q 0 into a power series of the kicking strength K and keeping the terms in the first order, we again get the effective HamiltonianĤ RWA .

Appendix C: Zak's kq Representation
The kq representation introduced in Ref. [93] is the complete orthonormal basis constructed by the common eigenstates of the translation operators in both X and P directions. In general, for the translation operation in X direction e iPA/λ , the "shortest" commutative translation operator in P direction is e iXB/λ with B = 2πλ/A. Given the dimensionless Planck constant λ = p/q, where p and q are coprime integers, we choose A = 2π/q and B = 2πp for constructing our Zak's kq representation. In the coordinate representation, the basis for given quantum numbers k X and k P is [93] with l ∈ Z. State function (C1) is composed of a series of Dirac's delta function with shifted phases. Note that the quantum numbers k X and k P take value in the region 0 ≤ k X < q, 0 ≤ k P < 1/p, which is q-times larger than the Brillioun zone. Since function (C1) is the eigenstate of translation operator e i 2π qλP , it is also the eignestate of the translation operatorT X = e i 2π qλP q , defined by Eq. (8) in the main text, with the the same quasimomentum k X but different Brillioun range 0 ≤ k X < 1. States (C1) with q different quasimomenta, i.e., φ (k X ,k P ) (X), φ (k X +1,k P ) (X), · · · and φ (k X +q−1,k P ) (X), can be treated as q degenerate states of operatorT X , which guarantee the completeness of the Zak's basis.
Appendix D: Derivation of Eq. (9) We outline how to derive Eq. (9) in the main text. The quasienergy state |ψ b,k can be spanned as |ψ b,k = q−1 m=0 u m |φ m (X) , where m ≡ (k X + m, k P ) with 0 ≤ k X < 1 and φ m (X) is the basis (C1) of the kq-representation. In this subspace, the eigenequation is simply [94] 2 cos λ(k X + m)u m + e −iλk P u m−1 + e iλk P u m+1 = 4E K u m (D1) together with the boundary condition u 0 = u q . To eliminate the dependence on k P , we make the substitution u m = e imλk P u m , which leads to and the boundary conditionū 0 = e −i2πpk Pū q . This equation can be formulated as where U m = (ū m ,ū m+1 ) and the matrix T m is From this recursive relation it is easy to find that U 1 = T(k X )U q with T(k X ) = q m=1 T m . As U 1 = e −i2πpk P U q , we have the secular equation det(T(k X ) − e −i2πpk P ) = 0, which can be expanded as Due to the cyclic permutation invariance of matrix trace, we have Tr T(k X + 1) = Tr T(k X ). Thus the expansion of Tr T(k X ) in terms of power series of e iλk X is simplely The coefficient C q can be directly evaluated by extracting the term with the highest power, resulting in C q = −1.
where we have used the Jacobi-Anger relation in the second line.
Appendix F: Obtaining Q(X, P) from w(s, k) The definition of the characteristic function w(s, k) we used in the main text can be formulated into an elegant form [95] w(η, η * ) = Tr ρe ηa † −η * a (F1) with η = (s + ik)/ √ 2λ, whereas the characteristic function of the Husimi distribution is given by, Clearly, the two characteristic functions are related through C Q (η, η * ) = w(η, η * )e −|η| 2 /2 . Once C Q (η, η * ) is obtained, the Husimi distribution Q(X, P) can be retrieved by the Fourier transform Appendix G: Expressions for U c and U e First, we calculate the matrix element of U(R) for two given functions f (X 1 , X 2 ) and g(X 1 , X 2 ) as following Here, we have used the property ofX c |C = C|C and the resulting f (X c )|C = f (X c )δ(X c − C)dX c = f (C) for any f (X c ) in the representation of coordinate of center of mass. Then, we apply the result (G1) to calculate U c and U e defined in the main text, i.e., the direct integral and the exchange integral Here, the overlap integral I N is given by Appendix H: I N for coherent states Now, we assume the two states ϕ(X) and φ(X) are two displaced squeezed coherent states described by The product state ϕ(X 1 )φ(X 2 ) = ϕ(C + 1 2 ∆X)φ(C − 1 2 ∆X) can be calculated from Eq. (H1) ϕ(X 1 )φ(X 2 ) = β √ π e −β 2 C 2 e −β 2 ( 1 2 ∆X−r m ) 2 .
Appendix I: U(R N ) for hardcore interaction We derive U(R N ) for hard-core interaction in Eq. (28). In the rest frame, assuming the interaction potential between two atoms is V(x 1 −x 2 ), the eigen equation of energy is given by Herr, E T is the total energy. We introduce the coordinate of central of mass x c ≡ (x 1 + x 2 )/2 and the relative coordinate x ≡ x 1 − x 2 . By separating the eigenstate into a product state Ψ = φ(x c )ψ(x), we have the eigen equation for the motion of center of mass − λ 2 2M with M = 2 the total mass and E c the energy of center of mass motion. The eigen equation for the relative motion is − λ 2 2µ with µ = 1/2 the reduced mass, E = E T − E c is the energy of relative motion. The solutions of Eq. (I2) are just the harmonic motions. We now try to find the solutions of Eq. (I3). Without interaction V(x 1 −x 2 ), the eigen problem is determined by H r φ n (x) = E n φ n (x) with H r ≡ − λ 2 2µ The eigenstates are given by where the parameter ζ = √ 1/(2λ). With consideration of hard-core interaction, i.e., V(x 1 −x 2 ) = +∞ for |x 1 −x 2 | < 2a and V(x 1 −x 2 ) = 0 for |x 1 −x 2 | > 2a, the boundary condition requires that wavefunction must be zero at x ∈ [−a, a]. For odd integer n, we assume the approximate eigenstates are just repulsed outside the hard-core region, i.e., n ∈ N (I5) For even integer n, the wave functions φ n (x), however, do not satisfy the hard-core boundary condition and the continuity condition. Therefore, we construct the symmetric eingenstates from antisymmetric states n ∈ N (I6) The energy levels to the first order correction are φ 2n |H r |φ 2n = φ 2n+1 |H r |φ 2n+1 = λ(2n + 1) + 1 2 a 2 + a ζ 1 (2n + 1)!2 2n+1 √ π [n+ 1 2 ],[n+ 1 2 ] k,l=0 (−1) k+l [(2n + 1)!] 2 (2n + 1 − k − l)! k!l!(2n + 1 − 2k)!(2n + 1 − 2l)! 2 2(2n+1−k−l) , where we have used µ = 1/2 in the last step. In fact, one can prove that U(R N ) is the first order correction, from the weak interaction, to the N-th energy level of the harmonic trapping potential. Relabelling N ≡ 2n + 1, we have from Eq. (I7) which is exactly Eq. (28) in the Sec. IV B 2.