Single-Photon Pulse Induced Transient Entanglement Force

We show that a single photon pulse (SPP) incident on two interacting two-level atoms induces a transient entanglement force between them. After absorption of a multi-mode Fock state pulse, the time-dependent atomic interaction mediated by the vacuum fluctuations changes from the van der Waals interaction to the resonant dipole-dipole interaction (RDDI). We explicitly show that the RDDI force induced by the SPP fundamentally arises from the two-body transient entanglement between the atoms. This SPP induced entanglement force can be continuously tuned from being repulsive to attractive by varying the polarization of the pulse. We further demonstrate that the entanglement force can be enhanced by more than three orders of magnitude if the atomic interactions are mediated by graphene plasmons. These results demonstrate the potential of shaped SPPs as a powerful tool to manipulate this entanglement force and also provides a new approach to witness transient atom-atom entanglement.

Single photon induced forces and torques correspond to the fundamental limit of optical linear momentum and angular momentum exchange with atoms [1]. Their direct detection is an open challenge since state-of-the-art quantum detectors are only sensitive to energy and arrival time of single photons [2]. Recent advances in temporal shaping of single photon scattering from atoms has shed light on the role of the temporal waveform of Fock states [3]. In light of these developments, it is an open question how single photon waveforms influence dipoledipole interactions (DDI) between atoms. Of particular interest is the exploration whether single photon shaped waveforms incident on interacting atoms can lead to experimentally observable transient effects.
During the last two decades, many techniques have been utilized to enhance the strength of the DDI and the corresponding force [4], such as utilizing micro-cavity [5][6][7], surface plasmons [8][9][10], and hyperbolic materials [11]. Especially, the strong DDI induced large energy shift in highly excited atoms (e.g. Rydberg) has been proposed as the mechanism for "Rydberg blockade", which provides a novel approach for quantum information processing [12,13] and simulation of quantum phase transition [14,15]. However, single-photon pulse as a tool to manipulate the transient dipole-dipole force has not been explored.
In this letter, we show the existence of a unique transient entanglement force between two neutral atoms induced by a single photon pulse. With the help of our defined force operator, we explicitly show that the resonant dipole-dipole interaction (RDDI) force fundamentally arises from two-body entanglement, which is significantly different from the van der Waals force. Our theoretical framework combines quantum theories of singlephoton pulse (SPP) scattering [16][17][18][19] and the macroscopic QED approach of dipole-dipole interaction [20][21][22][23]. We thus show that the quantum statistics of the incident (Fock-state vs coherent-state) pulses lead to significant differences in the induced RDDI entanglement forces. After absorption of a SPP, the inter-atomic force changes from the extremely weak van der Waals force [4,24,25] to the RDDI force [26,27] with the amplitude enhanced by ∼ 10 orders of magnitude.
We propose an experiment to detect this SPP induced force with two levitated neutral atoms (see Fig. 1), which are separated with distance r ∼ 1µm by optical tweezers operating at the magic wavelength [28][29][30]. Even with this enhancement, detection of such a weak transient RDDI force is still a difficult challenge. Therefore, we we demonstrate that the SPP induced RDDI force can be significantly enhanced by placing the atoms near a graphene layer with the assistance of graphenebased surface-plasmon polaritons. By investigating the full quantum dynamics of single-photon absorption, we predict optimum entanglement generation mechanisms conducive to experimental inquiry. Finally, we argue that the proposed effect can be differentiated from previously known dipolar interactions since the SPP induced entanglement force can be tuned from repulsive to attractive by tuning the polarization of the incident pulse.
Dipole-Dipole Interaction Force Operator-With the help of the Hellmann-Feynman theorem [32], we define a quantum operator to characterize the force generated by the coherent part of the DDI [33], where F mn (r) ≡ −∂U mn (r)/∂r is determined by the atom-atom interactionÛ (r) = mn U mn (r) |m n| induced by electromagnetic vacuum fluctuations [21,34] and |m ∈ {|gg , |eg , |ge , |ee } for a two-level-atom pair. The DDI force is always along the axis joining the two atoms. Our defined force operator allows us directly to classify the DDI force into two categories: (1) van der Waals force between two atoms in a directproduct state, such as the force for two ground-state Two atoms on top of a graphene layer (z0 is the height). These two atoms (the yellow spheres) are levitated by two separated optical tweezers. The relative displacement between the two atoms is r = x2 − x1 = r ex, which is along x-axis. The linearly polarized SPP propagates along y-direction, with polarization being parallel ( with θ = 0) or perpendicular (⊥ with θ = π/2) to r. For two ground-state atoms, the van der Waals force mediated by the vacuum fluctuations is extremely small (∼ 10 −35 N for r ≈ 1 µm, far beyond the state-of-art of the force sensitivity [31]). After absorption of a SPP, the atom-atom interaction changes to the resonant dipole-dipole interaction (RDDI) as shown in (c). The corresponding force is enhanced more than 10 orders to ∼ 10 −21 N. We emphasize that this RDDI force for atoms on states |Ψ ± = (|eg ± |eg )/ √ 2 is an entanglement force, which is fundamentally different from the van der Waals force. atomsF vdW = F gg,gg |gg gg|; (2) RDDI force for entangled atoms, e.g., We will show how to control this force with a SPP later. We emphasize that the latter RDDI force fundamentally arises from two-body entanglement [35]. The eigenvectors of the force operatorF RDDI (r) are the two Bell states with eigenvalues ±F eg,ge (r). This immediately reveals that, to maximize the RDDI force, one needs to prepare the atom pair in one of these two entangled states. We also note that, the RDDI force presents a readout of two-body entanglement. This entanglement force between transition dipoles is fundamentally different from van der Waals force [25] and the force generated by the permanent dipole-dipole interaction [20]. We emphasize that the maximum possible RDDI force (the eigenvalue of the force operator) is determined by the atom-atom distance r. However, the exact time-dependent envelope of the RDDI force in a specific dynamical process is determined by the atomic two-body entanglement. Dynamical Entanglement Force-The master equation method has been broadly applied to study the DDI and entanglement between neutral atoms [21,34,36,37]. We now incorporate the SPP-absorption dynamics with the traditional master equation to show the time-dependent entanglement force induced by a SPP [33], We have introduced an extra qubit degree of freedom ρ PN (t) to characterize the photon number degree (see more details in [18,33]). The initial value ofρ(t) is given byρ(0) =Î PN ⊗ ρ(0), whereÎ PN is the two-dimensional identity matrix and ρ(0) = |gg gg| denotes the initial state of the atom pair. The quantum pumping from a SPP is characterized by, where γ jj = γ 0 is the spontaneous decay rate of the atoms in vacuum. The coefficient η j characterizes the pumping efficiency, which is determined by the effective scattering cross section of the jth atom. The wave-packet amplitude of a Gaussian SPP is given by with center frequency ω 0 and pulse length τ f [19]. The time that the center of the pulse arrives at the jth atom is given by t j = k 0 · x j /ω 0 (|k 0 | = ω 0 /c). The absorption of the pulse is characterized by the Pauli matrixτ − of the extra qubit degree. The interatomic RDDI are included in the regular time-independent Lindblad superoperator [21,34] where ω 0 is the energy splitting of the two-level atoms, and the energy shifts δ ij = U eg,ge (r)/ and decay rates γ ij are given in the supplementary material [33]. Single-photon pulse induced transient entanglement force between two Rb atoms (D1 transition from 5 2 S 1/2 → 5 2 P 1/2 ). The force reaches its maximum when the photon absorption probability is largest. The magnitude of the RDDI force oscillates with atom-atom distance around r ∼ 1 µm. Here, the time is in the unit of 1/γ0 (γ0 is spontaneous decay rate of the atom in free space). Perpendicularly polarized pulse (⊥) is selected and its pulse length is set as γ0τ f = 0.63. The pumping efficiency is set as η1 = η2 = 1/ √ 2. The exact data of the Rb atom is given in the supplemental material.
The time-dependent RDDI entanglement force, F RDDI (r, t) = T r[ρ(t)F RDDI (r)], induced by a SPP for different atom distance is displayed in Fig. 2. For a fixed inter-atomic distance, the RDDI force increases after the pulse excites the atoms and decreases with time when atoms re-emit the photon. We can also see the amplitude of the RDDI force oscillates with atom distance r, due to the oscillation in the matrix elements F ge,eg (r) of the RDDI force operator. The van der Waals force has been neglected here as it is negligibly small [33]. The impulse force from the incident pulse is estimated to be F imp ≈ ω 0 /cτ f ∼ 10 −20 N with center frequency ω 0 ≈ 2π × 3.77 × 10 14 Hz and pulse length τ f ∼ 30 ns. But this force is along y-axis, which is perpendicular to the inter-atomic force in x-direction and can be relieved by the trapping force in y-axis. Thus, the only relevant force along the axis joining the two atoms is the RDDI entanglement force.
Quantum entanglement fundamentally determines the time-dependent RDDI force induced by a SPP. Here, we use the concurrence to quantitatively characterize the two-qubit entanglement [38]. As shown in Fig. 3 (a), for fixed atom-atom distance r = 1.2 µm, the concurrence C(t)(the dashed-pink line) and the RDDI force F RDDI (t) (the solid-blue line), as well as the excitation probability of the first atom P 1e (t) (the dotted-red line), reach their maxima simultaneously for homogeneous pumping case (η 1 = η 2 ). But for the local pumping of the first atom case with η 1 = 1 and η 2 = 0 [see The excitation probability P1e(t) first reaches its maximum and then the force and the concurrence reach their maximum in the local pumping case with η1 = 1, η2 = 0, and τ f γ0/2π = 0.75. Thus, the two-body entanglement is generated via the dipole-dipole interaction. Here, the atom-atom distance is fixed as r = 1.2 µm. In the double-y-axis figure, FRDDI(t) is associated with the left yaxis and both C(t) and P1e(t) are associated with the right y-axis. Fig. 3 (b)], C(t) and F RDDI (t) reach their peaks at the time, which is later than the time when P 1e (t) reaches its maximum. Thus, it is the entanglement instead of the total excitation probability that maximizes the RDDI force. We also see that there are two ways to generate the quantum entanglement between the atoms: (1) homogeneous pumping to the symmetric state |Ψ + directly by the SPP; (2) local pumping of single atom to state |eg and then the RDDI evolves the atoms to entangled states. Here, we show that the first one is more efficient for entanglement generation. The total photon absorption probability P e,tot (t) for both homogeneous [P e,tot (t) = 2P 1e (t) in Fig. 3(a)] and local pumping cases [P e,tot (t) = P 1e (t) in Fig. 3(b)] are almost the same. But the entanglement and the RDDI entanglement force under homogeneous pumping are much larger than that of local pumping case. This is because the projection of the atomic state ρ(t) on the entangled state |Ψ + under homogeneous pumping is much larger.
The existing theory [4,20,21,34] can not describe the quantum pulse induced DDI force. Now, we show that the force induced by a Fock-state pulse is significantly different from the one induced by a coherent-state pulse. As explained in Ref. [19], the absorption probability of Fock-state SPP by a two-level atom is much higher than that of coherent-state pulse. Thus, the corresponding force is larger as shown by the blue lines in Fig. 4. However, there exists an optimal pulse length τ f,opt to reach the largest excitation probability of the atoms for Fock-state pulses [18]. For fixed pulse length τ f γ 0 = 0.3, the maximum entanglement force decreases with photon number in Fig. 4 (a), as the total excitation probability decreases [18]. But the force induced by coherent pulse always increases with the mean photon number [see Fig. 4 (b)]. In an experiment, larger entanglement force can be obtained by optimizing the pulse length to increase the atomic excitation probability for given atomic transition frequency and DDI strength [33]. The entanglement force can be enhanced significantly by engineering the nanophotonic environment near the atoms. As a practical illustration, we demonstrate this enhancement by placing the atoms near a graphene layer as depicted in Fig. 1 (b). The surface plasmon polaritons of graphene have been previously shown to allow conventionally forbidden atomic transitions [39] in addition to enhancing other well-known physical effects such as decay rate of emitters [40] and Förster energy transfer rate [41]. Here, we find that the time dependent entanglement force can be enhanced significantly by placing the atoms near a graphene layer. Figure 5 (a) demonstrates the distance dependence of the engtanglement force. For atomic transition frequency close to graphene surface plasmon polaritons (exact data provided in [33]), the enhancement factor is larger than 1000 at atom-surface distance z 0 = 10 nm (red curve). While the graphene-based surface plasmon polaritons occur in the terahertz to near infrared band [40,42], similar enhancement at optical frequencies are feasible with other plasmonic materials such as gold and silver [43].
Precise Control of The Entanglement Force-Now, we show SPP as a novel tool to precisely control the atomic force: (1) a more than ten orders of DDI force amplitude change can be induced by a SPP; (2) the induced entanglement force can be continuously tuned from being repulsive to attractive by varying the polarization of the pulse. For relevant inter-atomic separations (r ∼ 1 µm), Force (N) Time * + ( C + = 10nm * + 0 1 = 0.001 C + = 20nm * + 0 1 = 0.0019 C + = 50nm * + 0 1 = 0.0075 Figure 5. (a) Single-photon pulse induced entanglement force between two atoms placed near a graphene-layer interface. Here, the forces have been normalized by the eigenvalue of the force operatorFRDDI(r) in free space. (b) Eigen value of the force operatorFRDDI(r) for two Rb atoms in free space as a function of atomic distance r. The induced RDDI forces FRDDI are different for parallel ( ) and perpendicular (⊥) polarizations of single-photon pulses, as shown by the solid-green (F ) and dotted-blue (F ⊥ ) curves. In the subplot, we show the force FRDDI with r = 0.8 µm (marked by the vertical dashed line) for different polarization angle (θ with respect to x-axis) of the pulse in xz-plane. This clearly shows the change in sign of the force from repulsive to attractive.
the van der Waals force is around ∼ 5 × 10 −35 N (see Fig. 1 in [33]), which is far beyond the state-of-art force sensitivity (∼ 10 −24 N/ √ Hz) [31]. As the van der Waals force arises from higher-order process, thus it is much weaker than the RDDI force. After absorption of a SPP, the RDDI force dominates with a greatly enhanced amplitude ∼ 10 −22 N. This force can be further enhanced upto 10 −19 N with surface plasmons-plaritons making its experimental detection possible.
For atomic transition between states connected by linearly polarized light, the direction of the corresponding transition dipole is determined by the polarization of the incident pulse. As shown in Fig. 5 (b), both the forces induced by parallelly ( ) and perpendicularly (⊥) polarized pulses oscillate with the atomic distance around r ∼ 1 µm. But these two forces have a phase shift and usually have opposite signs (especially in the near region r < 0.5 µm). Thus, we can control the force to be repulsive or attractive by changing only the polarization of the pulse. More importantly, we can continuously tune the value of the RDDI force via tuning the pulse polarization angle θ in xz-plane with fixed atom-atom distance (r) (see the subplot).
Looking ahead, our work provides a natural platform to investigate photoassociation in chemical reactions and bioprocesses [30]. By generalizing the force operator to multi-atom case, we can also study the role of the manybody entanglement in the collective force of neutral atom ensemble [19,44].
Multiply both side with m|, we have In most case, due to the non-adiabatic transition terms in the square brackets, there does not exist a well defined force operator for a microscopic system, such as the exchanging interaction in a condensed-matter lattice. As the distance between the two atoms is much larger than the size the the atoms, thus the atomic wave function is not dependent on the relative distance r and the second term at the right-hand-side disappears (i.e., l|∂n/∂r = 0). In the atomic Hamiltonian, only the dipole-dipole interaction part depends on the inter-atomic distance r. As the corresponding force is always along the co-axis line, we can define a scalar operator for this force as,F We note that this force operator only works for weak atom-field coupling case. If the two atoms strongly coupled to a resonant cavity field, one can not eliminate the degree of the cavity mode to obtain an effective interaction Hamiltonian as shown in Eq. (3). In this case, the inter-atomic force is not only dependent on atom-atom separation, but also the position of each atom [2]. More important, the magnitude of the forces experienced by the two atoms can be different, which violates Newton's third law for a macroscopic body. We do not consider this case in this paper.
Different elements in the operatorF (r) correspond to different virtual processes generated forces. We emphasize that only the anti-diagonal elements of the two-body interaction in (3) can be mediated by second-order processes [3] and all the other terms result mainly from fourth order processes. Thus, the corresponding forces are weak. In this paper, we only focus on two forces. The first one is the van der Waals (vdW) force F vdW ∝ F gg,gg (r) between two ground-state atoms, which mainly arises from fourth-order process [4,5] and usually is extremely small. An incident single-photon pulse (SPP) can pump the atom pair to an entangled state. In this case, the interaction changes to the resonant dipole-dipole interaction (RDDI), which plays the key role in energy transfer between different molecules in chemical and biological processes. As the RDDI is mediated by second-order processes, the corresponding force F RDDI ∼ F ge,eg (r) between the two atoms will be greatly enhanced. In the following, we present the approach to calculate the elements ofÛ (r) andF (r).

II. MODEL HAMILTONIAN FOR ATOM-VACUUM FIELD INTERACTION
The Hamiltonian of the total system is given bŷ where the Hamiltonian of the field modes in an arbitrary linear (non-magnetic) media is given by [6,7] and the ladder operators of the eigen modes satisfy the commutation relations and The Hamiltonian of the two atoms is where ω a,j is the energy splitting of the jth atom andσ + j = (σ − j ) † = |e j g j | is the Puali matrix. There are two forms of Hamiltonian to describe the interaction between the atoms and the electromagnetic field. One is the minimum coupling and the other one is the multipolar coupling [4]. The difference and relation between these two forms of interaction can be found in [4,8]. Here, we use the multiploar interaction Hamiltonian.
where µ j,eg is the electric dipole transition element of the jth atom. The electric field operator can be expanded with the eigen modes of the field aŝ whereÊ The function ← → G (x, x , ω) is the classical Green tensor obeying the equation Here, we assume that the media is a non-magnetic material with constant permeability µ 0 = 1 and the frequency dependent complex dielectric constant ε(x, ω). The Green tensor has the properties We will show that both the van der Waals interaction and the resonant dipole-dipole interaction can be easily obtained with the Green tensor.

III. VAN DER WAALS INTERACTION
The van der Waals interaction between two atoms has between well studied. A detailed calculation of the coherent van der Waals interaction in free space is presented in Ref. [4]. Here, we only present the more general form of van der Waals interaction between two identical atoms obtained by Safari and his collaborators [5], The incoherent part of van der Waals interaction has been neglected, as it is usually negligible small compared to the spontaneous decay rate of the atoms.

IV. HEISENBERG-EQUATION METHOD TO CALCULATE THE RESONANT DIPOLE-DIPOLE INTERACTION
To calculate the resonant dipole-dipole interaction (RDDI), we take one of the atoms as the source and the other as the test dipole. The interaction between the two atoms is given by the energy of the second atom in the field generated by the source (the first one) [9]. By taking only the interaction between the field and the first dipole into account, the motion equation of the noise operatorf (x, ω) in Heisenberg picture is given by with formal solution Then, we substitute this noise operator to the electric field and split the electric field in to two partsÊ(x, is the free field and the field due to the source is given bŷ where we have used the relation in Eq. (18) and replaced t − τ with s. The extra energy of the system by adding the second dipole is given bŷ As shown in Ref. [9], the effective coupling between the two atoms is already of second-order of atom-field coupling. Thus, we can replace the atomic operator in the integral witĥ where we have assumed ω a,1 = ω a,2 = ω 0 . Finally, we can write the effective interaction Hamiltonian between the two atoms in the Heisenberg picture aŝ containing both coherent and incoherent interaction parts V mn (r) = U mn (r) + iγ mn (r)/2. Here, Next, we will take the Markov approximation by extending the upper limit in the integral over s to ∞ and using Then real and imaginary (principle integral) parts of the integral will contribute to the cooperative decay rates and energy shift respectively. Using µ 0 = 1/ε 0 c 2 , we obtain the cooperative decay rates Only γ eg,ge and γ ge,eg will be involed in the following. The energy shift are given by where the second line we have used the Onsager reciprocity condition x 2 , ω) and Im ← → G (x 1 , x 2 , ω) = −Im ← → G * (x 1 , x 2 , ω), the last step we have used the Kramers-Kronig relations (we have assume that µ j are real). Similarly, we have The coherent parts of the DDI U eg,ge and U ge,eg contribute most to the entanglement force. The eigen function of the corresponding force operatorF RDDI = F eg,ge |eg ge| + F ge,eg |ge ge|, are the two Bell states |Ψ ± = (|eg + |ge )/ √ 2. We also note that similar entanglement force also exist for the other two Bell states (|gg ± |ee ) / √ 2. However, the corresponding force element F gg,ee (r) |gg ee| oscillates (with frequency 2ω 0 ) very fast in time and the time-averaged force is zero. Thus, the dipole-dipole force induced by transitions |gg ↔ |ee has been neglected in the following. The mean value of the RDDI force is also zero for two atoms in a single-excitation product state (i.e., |eg or |ge ). For atoms in the state |eg , only the van der Waals force F eg,eg |eg eg|, which is similar to the force between two ground-state atoms F gg,gg |gg gg|, has non-zero mean value.

V. MASTER-EQUATION METHOD TO CALCULATE THE RESONANT DIPOLE-DIPOLE INTERACTION
In this section, we calculate the RDDI by re-deriving the Lindblad form master equation for a two-level-atom pair. This master equation can also be found in Refs. [10][11][12]. Here, we present some non-trivial details. In the following, for simplicity, we consider two identical atom case µ j,eg = µ j,eg = µ = d 0 e d . For atomic states connected by linearly polarized light, the direction of the transition dipoles e d are determined by the polarization of the incident pulse. This makes it possible to precisely control the RDDI force by tuning the polarization of the pulse as shown in the main text.
To derive the master equation fo the atomic paire, we first split the total Hamiltonian intoĤ =Ĥ 0 +V , wherê H 0 =Ĥ F + j=1,2Ĥ A,j andV = jĤ AF,j . In the interaction picture with respect toĤ 0 , the atom-field coupling becomes time-dependent The density matrix of the whole system satisfies the Liouville-von Neumann equation Its formal solution is given by Substituting this formal solution back to the right hand side (r.h.s.) of Eq. (46), we have After tracing off the degrees of freedom of the field, we obtain the master equation for the atomic density matrix ρ(t) ≡ Tr F [ρ tot (t)]. Without loss of generality, the initial state of the whole system is assumed to be in a product state ρ tot (0) = ρ(0) ⊗ ρ F (0). Here, we only consider the vacuum fluctuations, thus the initial state of the field is ρ F (0) = |0 0|. It can be easily verified that the first term at the r.h.s. of Eq. (48) does not contribute to the motion equation of ρ(t). Then, we have where we have taken the Born approximation [13]. In the following, we will handle the energy conserving terms (such asσ + iσ − j ) and energy non-conserving terms [such asσ + iσ + j exp(2iω 0 t)], respectively.

A. Energy conserving terms
We first consider the energy conserving terms. After tracing off the degrees of the field, we obtain where which exactly return to the results obtained with the Heisenberg equation in Eqs. (36) and (43). We emphasize that the non-rotating wave terms in the Hamiltonian (10) are very important to obtain the full energy shift. Otherwise, only the first term in Eq. (52) exists under the rotating wave approximation. This can also be found in the seventh chapter of the book [4].
When i = j,L +− characterizes the spontaneous decay of the atoms induced by the local field fluctuations with rates γ 11 = γ 22 . The Lamb shift δ ii (usually diverges without renormalization) will be neglected in the following.
When i = j,L +− characterizes the coherent (δ ij -terms) and incoherent (γ ij -terms) coupling between the two atoms mediated by the non-local field fluctuations as explained in the main text. Especially, the relative-position dependent energy shift δ 12 (r) can generate an entanglement force between the two atoms, which can be measured in experiment with the state-of-art force detection sensitivity.

B. The energy non-conserving terms
Now we consider the energy non-conserving terms. Similarly, after tracing off the field degree, we obtain Here, we have neglected the terms with i = j to obtain the traditional Lindblad master equation. Furthermore, these terms oscillate with frequency 2ω 0 in the interaction picture, thus its time-averaged effect is negligible small. 85 Rb Transition Frequency ω0 Transition Dipole element d0 Spontaneous Decay γ0 Life time 1/γ0 D1 (5 2 S 1/2 → 5 2 P 1/2 ) 2π × 3.77 × 10 14 Hz 2.54 × 10 −29 C · m 2π × 5.75 × 10 6 Hz 27.68 × 10 −9 s Table I. The data of the 85 Rb atom used in this paper coming from Ref. [14]. We note that the spontaneous decay rate can be obtained directily from Eq. (63) with ω0 and d0.
Finally, we obtain the the full master equation for two atoms in the Schrödinger picture We can see that, different from the energy conserving terms, the energy non-conserving terms cannot be expressed in the Lindblad form [10]. Thus, the coherent coupling part (δ ± ij ) cannot be simply explained as the energy shifts of the atomic level. But these terms oscillate with frequency 2ω 0 in the interaction picture. In most case, it is safe to neglect the energy non-conserving terms (the second line).

VI. DIPOLE-DIPOLE FORCE IN FREE SPACE
In this section, we present the details to calculate the strength of the van der Waals interaction and the RDDI for two atom in vacuum. We first recover the well known results for the dipole-dipole interaction obtain by mode expansion. It is easy to find that if we let r = x 2 − x 1 = (r, 0, 0), only the diagonal elements of the free space Green tensor are nonzero [7,15], It is straightforward to verify that, for the single point Green's function, the real part diverges, but the imaginary part does not, Then, we can obtain the well known spontaneous decay rate of an atoms in free space, We will take γ 0 = 1 as the unit of frequency and 1/γ 0 as the unit of time in this paper. As shown in the next section, both the coherent and incoherent dipole-dipole interaction can be greatly enhanced by engineering the electromagnetic environment to change the Green tensor.

A. van der Waals interaction
As the ground-state atoms can be excited by arbitrarily polarized virtual photons. Thus, to calculate the van der Waals interaction, we need average out the polarization angle by taking the spherically symmetric polarizability tensor [see Eq. (49) in Ref. [5]]. Finally, the van der Waals interaction between two ground-state atom is given by Using the method presented in [4] (see Chaps. 7.5 and 7.6), we can verify that: U gg,gg (r) ∼ 1/r 6 , ur 1 1/r 7 , ur 1 .
Thus, the corresponding force F vdW (r) will be of scale∼ 1/r 7 in the near region and ∼ 1/r 8 in the far region. As shown by the pink curve in Fig. 1, the van der Waals force F vdW (r) deviate from the line 1/r 7 (the thin black line) slightly in the far region. decay rates γ 12, for parallelly polarized atoms and γ 12,⊥ for perpendicularly polarized atoms behave differently, but both of them converges to the spontaneous decay rate γ 0 in the near region [see the subplot in Fig. 2(a)]. The coherent part of the RDDI diverges in the near region. More importantly, δ 12, and δ 12,⊥ usually have opposite signs, especially in the near region. This lays the foundation to tune the RDDI force by tuning the polarization of the pulse as explained in the main text. The matrix element of the force operatorF RDDI are given by and The numerical results are displayed in Fig. 1. In the near region, the RDDI force decreases with 1/r 4 . In the far region, F RDDI, decreases with 1/r 2 (green solid line) and F RDDI,⊥ vanishes with scaling 1/r (blue dotted line).

VII. DIPOLE-DIPOLE FORCE NEAR PLANAR INTERFACE
where α(ω) = 2πσ(ω)/ε 0 c is the dimensionless in-plane conductivity of the graphene. The optical conductivity of a graphene layer can be split into intra-band and inter-band contributions σ(ω) = σ intra (ω) + σ inter (ω) with [19,20] and where τ D is the relaxation time in the Drude model, E F the graphene's Fermi energy, and the function .
The RDDI strength for two atoms on top of a graphene layer is given by Then, the eigen value of the RDDI force operator on the state |Ψ + is obtained as F (r) = −∂U eg,ge (r)/∂r. In Fig. 3, to show the enhancement in the RDDI force due to the graphene layer, we re-scale F (r) with the eigen value F vacuum (r 0 ) of the corresponding RDDI force operator in vacuum at r 0 = 1.05µm (denoted by the vertical black line). Comparing with the subplot, we see that more than three order enhancement in the force can be obtained if the atoms are very close to the graphene layer (z 0 = 10 nm). We also see that this enhancement decreases fast with the hight of the atoms z 0 and vanishes for z 0 > 500 nm.
In the main text, the corresponding time-dependent entanglement force induced by a SPP has been displayed. The inter-atomic distance is set as r = 1.05µm as marked by the dark vertical line in Fig. 3 and the atom-interface distance is set as z 0 = 10, 20, 50 nm. The pulse length τ f has been optimized to get the maximum entanglement force as both the local spontaneous decay rate γ ii and the cooperative decay rates γ ij defined in Eq. 51 have also been greatly enhanced by the graphene layer.

VIII. DYNAMICAL FORCE AND TWO-BODY ENTANGLEMENT
In this section, we study the dynamics of a two-level-atom pair. Different from the previous literatures, we prepare the atom pair in the ground state |gg instead of a single-excited state (e.g., |eg ). In 2012, Ben et al. derived a powerful time-dependent master equation for n-photon broadband pulse interacting with an arbitrary quantum system. Here, we generalize this method to calculate the dynamical RDDI force. More important, we explicitly show that the RDDI force fundamentally arises from two-body entanglement.
The total master equation including the single-photon pumping process is given by, whereρ(t) = ρ PN (t) ⊗ ρ(t) is an effective density matrix and we have introduced an extra qubit degree of freedom ρ PN (t) to characterize the photon number degree (see more details in Ref. [21]). The initial value ofρ(t) is given by  Figure 3. The eigen value F ⊥ (r) of the resonant dipole-dipole interaction (RDDI) force operator on state |Ψ + for two atoms on top of a graphene layer. Different curves denote different atom-interface distance z0. In the subplot, we display the details of the curve for free-space case and the curves with z0 = 200 nm and z0 = 500 nm. Here, the electric dipole moments (along z-direction) of the atoms are perpendicular to the relative displacement r and F ⊥ (r) has been re-scaled by the eigen value Fvacuum(r0) of the corresponding RDDI force operator in vacuum at r0 = 1.05µm (denoted by the vertical black line). The Fermi energy of the graphene is set as EF = 1.0 eV and the relaxation time is taken as τD = 10 −13 s. To obtained a large enhancement in the RDDI force, the energy splitting of the two-level atoms is set as ω0 = 0.7 eV different from the optical transition in Rb atoms as shown in previous section. The graphene layer is considered to lie on an ε(ω0) = 2.5 substrate. ρ(0) =Î PN ⊗ ρ(0), whereÎ PN is the two-dimensional identity matrix and ρ(0) = |gg gg| is the initial state of the atom pair. The the first term at right hand side (r.h.s.) of Eq. (93) characterizes the free evolution of the atom pair without the pumpinĝ The second term characterizes the pumping of the SPP, with Pauli matricesτ − characterizing the absorption of the SPP. The parameter η j characterizes the pumping efficiency of the jth atom determined by its effective scattering cross section, t j = (x j · e p )/c is the time of the center of the pulse arriving the jth atom, and is the Fourier transform of the pulse spectrum function. For a Gaussian SPP, its wave packet amplitude in the time-space domain is given by, In the main text, we assume the pulse propagates along the x-axis and arrives at the two atoms at the same time t 1 = t 2 = 0. The pumping efficiency η j in practice shoule be much smaller than 1 [22,23], but its can be enhanced by adding a mode converter [24]. In our simulation, we take η 1 = η 2 = 1/ √ 2 for the homogeneous pumping case and η 1 = 1, η 2 = 0 for the locally pumping case.
This effective master equation method can be straightforwardly generalized to n-photon Fock-State pulse case by replacing the Pauli matrixτ ± in Eqs.
and replacing the 2 × 2 identity matrixÎ PN with the (n + 1) × (n + 1) identity matrix. Actually,ρ(t) is not a real density matrix of a physical system, as Trρ(0) = n for n-photon Fock-state pulse. Thus, only its projection on the specific subspace has physical meaning. The expected value of any atomic operatorÔ is given by whereP is the projection operator of the extra qubit degree with the only non-zero element P 11 = 1. We also note that, to handle the coherent-state pulse case, we only need to replace all the photon related operators (i.e.,τ ± ,Î PN , andP ) with the constant 1. This powerful time-dependent master equation can be used to uniformly study the quantum photon pulse scattering process. The dynamical RDDI force has been shown in the main text. Here, we show the time-dependent concurrence C(r, t) of the atom pair for different atom-atom distance in Fig. 4 (a). The concurrence is calculated with the method presented in the seminal work [25] by extracting the density matrix ρ(t) of the atom pair fromρ(t) via whereÎ a is the identity matrix of the atomic degrees of freedom. The concurrence C(r, t) has similar dynamical behavior as the the RDDI force shown in the main text. For the homogeneous pumping case η 1 = η 2 = 1/ √ 2, the two-body entanglement is mainly generated by the single-photon pulse. Thus, the concurrence always reaches its maximum at the peak of the absorption probability. The oscillation with the atom-atom distance r results from the oscillation in the incoherent coupling strength γ 12,⊥ (r), as the maximum absorption probability is dependent on the decay rate of the atomic system [21][22][23]. To verify this, we plot the excitation probability P 1e (r, t) of the first atom as a function of time and atom-atom distance with fixed photon pulse in Fig. 4 (b). The envelopes of the concurrence and P 1e (r, t) have the same shape. We can also enhance the RDDI force by changing the pulse length τ f to optimize the two-body entanglement (see Fig. 5). Here, we see that, for homogeneous pumping case with η 1 = η 2 = 1/ √ 2, the optimal pulse length maximizes the local excitation probability of the first atom P 1e , the inter-atomic force F RDDI , and the concurrence C simultaneously [see Fig. 5(a)]. But, for local pumping case with η 1 = 1 and η 2 = 0, only the pulse length optimizing C maximizes the RDDI force [see Fig. 5(b)]. A shorter pulse optimizes the photon absorption probability P 1e , but the entanglement and the force are suppressed due to the low entanglement generation rate via the weak RDDI coupling and the fast spontaneous decay rates of the atoms. Thus, the homogeneous pumping is a more efficient way to generate the entanglement force.