Collective dynamics and entanglement of two atoms embedded into negative index material

We study the dynamics of two two-level atoms embedded near to the interface of paired meta-material slabs, one of negative permeability and the other of negative permittivity. The interface behaves as a plasmonic waveguide composed of surface-plasmon polariton modes. It is found that significantly different dynamics occur for the resonant and an off-resonant couplings of the plasma field to the atoms. In the case of the resonant coupling, the plasma field does not appear as a dissipative reservoir to the atoms. We adopt the image method and show that the dynamics of the two atoms are completely equivalent to those of a four-atom system. Moreover, two threshold coupling strengths exist, one corresponding to the strength of coupling of the plasma field to the symmetric and the other the antisymmetric mode of the two-atom system. The thresholds distinguish between the non-Markovian (memory preserving) and Markovian (memoryless) regimes of the evolutions that different time scales of the evolution of the memory effects and entanglement can be observed. The Markovian regime is characterized by exponentially decaying whereas the non-Markovian regime by sinusoidally oscillating contributions to the evolution of the probability amplitudes. The solutions predict a large and long living entanglement mediated by the plasma field in both Markovian and non-Markovian regimes of the evolution. We also show that a simultaneous Markovian and non-Markovian regime of the evolution may occur in which the memory effects exist over a finite evolution time. In the case of an off-resonant coupling of the atoms to the plasma field, the atoms interact with each other by exchanging virtual photons which results in the dynamics corresponding to those of two atoms coupled to a common reservoir. In addition, the entanglement is significantly enhanced under the off-resonant coupling.


I. INTRODUCTION
The radiative properties of emitters (e.g. atoms or quantum dots) located inside a dielectric or conducting material can be significantly modified compared to those in vacuum. The modification results from the variation of the density of modes of the EM field which can be adjusted by changing the geometric shape or space period structure of a material [1][2][3]. The radiative properties of emitters can also be modified by locating the atoms close to the surface of a dielectric or conducting material [4][5][6][7][8][9][10][11][12]. In this case, the modification results from the presence of surface EM modes known as plasmon guided (PG) field [13][14][15][16][17][18][19]. A new category of materials has been proposed, so-called meta-materials [20][21][22][23], characterized by specifically designed geometrical structures which drastically modify the density of the EM modes, so the field propagation, also yielding to the PG field [24][25][26][27][28]. Owing to the high local density of modes, emitters may interact strongly with the surface field which affects their radiative properties [29][30][31][32][33][34][35]. For example, when quantum dots are placed at a distance about several tens of nanometers above two dimensional metal surface, strong coupling could be generated between the quantum dots and the collective mode reflected in the presence of the Rabi oscillations [36]. In the structure composed of zero index and left hand materials, maximum quantum interference and a suppression of the atomic decay rate can be achieved between Zeeman levels due to the anisotropy of the EM modes [37]. It * gaox@phy.ccnu.edu.cn has been shown that by changing the strength of the driving field and adjusting position of the quantum dot, the plasma mode on the surface of a metal-nanoparticles material can introduce asymmetrical features into the spectrum [38].
The study of entanglement between distant atoms and controlled the transmission of information between them are vital to the development of quantum information technology [39,40]. A key model to investigate the creation and storage processes of entanglement often follows with a system composed of two-level quantum emitters [41,42]. When atoms coupled to the same PG field, the incoherent spontaneous exchange of photons could occur and results in collective damping [43,44]. It has been revealed when the decay through one of the collective states is deeply depressed, long lived entanglement of the system could be achieved, which only depends on the distance between atoms [45]. By applying nanowire structure, the PG field could be well guided and thus entanglement between quantum dots still exists at several vacuum wave lengths [46]. Xu et al. [47] have shown that the entanglement between two atoms can exist over distances much larger than the resonant wavelength if the space between the atoms is filled with a thin membrane made of single negative (ε < 0 or µ < 0) and left hand materials. However, most of the studies are focus on metal-dielectric structure, and the entanglement only maintains for a small time scale. permeability (MN) and the other of a negative permittivity (EN). We assume that the atoms are located in the MN material. Surface plasmon polaritons (SPP), nonradiative electromagnetic excitations associated with charge density waves propagating along the interface are generated. The SPP propagate in the xand y-directions along the interface between the meta-materials, and rapidly decay in the z-direction. As illustrated in figure 1, two atoms located close to the interface can excite the SPP and the excitation depends on the distance of the atoms from the interface and the polarization of the atomic dipole moments. For a polarization of both atomic dipole moments in the x−z plane, the atoms could excite the SPP which would effectively propagate in the x direction. In other words, the polarization of the atomic dipole moments determines the direction of propagation of the SPP on the interface. Hence, the interface would behave as a directional guiding plasma field mode propagating in the x-direction, formally analogy of a plasmonic waveguide [45,48,49]. We shall demonstrate below that then it would be possible to achieve a strong interaction between two atoms located close to the interface through their coupling to the SPP.
The mathematical approach we adopt here is based on the Green's function method [50]. Our focus is on how the plasma field induced at the interface between the materials changes the dynamics of the atoms, in particular, the population transfer and entanglement. The remarkably simple analytical expressions are derived for the probability amplitudes valid for an arbitrary initial state, arbitrary strengths of the coupling constants of the atoms to the plasma field, and arbitrary distances between the atoms. We find a number of interesting general results. In the first place, we distinguish two different time scales of the evolution of the atomic states, one corresponding to the evolution of the collective symmetric state and the other to the antisymmetric state. Secondly, we find a threshold behavior of the coupling constants which separate the non-Markovian behavior of the system from the Markovian one [51]. The Markovian evolution is usually attributed to a weak coupling of an atom to the field. We show that the collective effects may result in the Markovian evolution even in the limit of a strong coupling of the atoms to the field. Inversely, a non-Markovian evolution can be seen even in the regime of a weak coupling. Thirdly, we find that the plasma field does not appear as a common reservoir to the atoms. In order to explain the behavior of the atoms we adopt the image method and show that the dynamics of the two atoms are completely equivalent to those of a four-atom system. Finally, we consider the case in which the plasma field frequency is offresonant with the atoms and find that in this case the dynamics resemble those of two atoms coupled to a common reservoir. The atoms interact through the exchange of virtual photons resulting in the absence of the images.
The plan of this paper is as follows. In section II we introduce the model and present the explicit analytic expressions for the time dependence of the probability amplitudes of the atoms. Detailed dynamics of the atoms are studied in section III. We assume that a single excitation is present initially in the system and demonstrate how the evolution of the system can be simply understood in terms of the evolution of (Color online) Two atoms located at distances z1 and z2 (z2 > z1) close to the interface between two meta-materials. Each atom excites the SPP propagating along the interface. The excitation size of the SPP depends on the distance of the atom from the interface and polarization of its dipole moment. As illustrated by green and purple ellipses, the excitation area of the surface plasmon polaritons decreases with an decreasing distance zi and their shapes extend more in the x direction than in the y direction if the atomic dipole moments are polarized in the x − z plane.
the atoms and their corresponding images. We then demonstrate in section IV the collective behavior of the atoms in both Markovian and non-Markovian regimes of the evolution. Entangled properties of the atoms are discussed in section V, where we calculate the concurrence for different initial states and different coupling strengths of the atoms to the plasma field. The effect of an off-resonant coupling of the atoms to the plasma field on the collective dynamics and entanglement is discussed in section VI. We summarize our results in section VII. The paper concludes with two Appendices in which we give details of the derivation of the integro-differential equations for the probability amplitudes and the calculations of the integral kernels. Both longitudinal and transverse parts of the Green function are considered in the evaluation of the kernels.

II. ATOMS INTERACTING WITH THE PLASMA FIELD
Perfectly conducting materials are known to generate a strong PG field which is refined in a short regime near the surface [52][53][54]. A plasma field can also be generated at the surface of a meta-material with either negative permittivity (ε < 0) or negative permeability (µ < 0). However, the density of the plasma field near surface only originates from one or several discrete modes, which can be derived by applying the continuous conditions on the boundaries. Recently, Tan et al. [55] have shown that the density of the EM modes can be significantly enhanced at the interface of two meta-materials, one with negative ε and the other with negative µ. Especially, when the materials are perfectly paired, i.e. ε 1 = −ε 2 and µ 1 = −µ 2 , the effective permittivity and permeability, defined as are both zero when the slabs have the same thickness, d 1 = d 2 . In this case, the band gap disappears and the density of modes of the EM field becomes continuous. This means that a large density of the modes exists at the interface between the two perfectly paired negative meta-material slabs, and can be treated as optical topological material [56]. We consider a system composed of two identical atoms located at a distance z 0 from the interface between two different negative index material slabs, µ-negative (MN) and ε-negative (EN) slabs, as shown in figure 2. We assume that the atoms are located in the MN slab and the distance between the atoms, x 21 = x 2 − x 1 , is large compared to the atomic wavelength, x 21 ≫ λ a , so there is no direct interaction between the atoms. Each atom is represented by its ground state |g j , an excited state |e j (j = 1, 2), the atomic transition frequency ω a , and the atomic transition dipole moment p j .
FIG. 2. (Color online) Geometry of the system. The z axis is taken normal to the interface between MN and EN slabs with its origin at the interface. The slabs have thickness d1 and d2, respectively, and are assumed to have infinite extents in the xy plane. Two atoms are embedded in the MN slab at fixed positions (x1, 0, z0) and (x2, 0, z0), where z0 is the distance of the atoms from the interface between the materials. The atomic transition dipole moments p1 and p2 are parallel to each other and oriented in the x − z plane.
The atoms interact with an electromagnetic field via a dipole interaction according to the Hamiltonian [57] H =Ĥ 0 +Ĥ I , is the unperturbed Hamiltonian of the atoms and the field, and is the interaction of the atoms with the field. Here,f † λ (r, ω) andf λ (r, ω) are the creation and annihilation operators which can be viewed as collective excitations of the electromagnetic field, λ = e represents noise polarization of the EN material, λ = m represents noise magnetization of the MN material [12,58],σ † j (σ j ) andσ z are the raising (lowering) and the energy difference operators of atom j. The positive frequency part of the electric field operator at the position r j of the jth atom is given bŷ where ℑ[ε(r, ω)] is the imaginary part of permittivity, ℑ[κ(r, ω)] is the imaginary part of reciprocal of permeability (κ(r, ω) = 1/µ(r, ω)), respectively, and ↔ G (r j , r, ω) is the Green tensor of the field, which characterizes the density of the field modes at the location r j of atom.
For the permittivity and permeability of the slabs, we assume that ε 1 and µ 2 are positive constants, but ε 2 and µ 1 are negative and strongly depend on frequency of electromagnetic field which have following forms [59,60] where ω ep and ω mp are plasma frequencies of the electric and magnetic materials, respectively, ω eo and ω mo are resonance frequencies of the materials, and γ e , γ m are dissipation (losses) parameters of the materials. For clarity of the notation we have omitted the spatial argument. It is clear from Eq. (6) that in the frequency region above the resonance, ω > ω eo (ω > ω mo ) and ω < ω ep (ω < ω mp ), the EN (MN) slab is a single-negative material. Thus, the structure of a meta-material can be designed [61,62].
To study the dynamics of the atoms, we consider the wave function of the system whose the time evolution is governed by the Schrödinger equation If the field was initially at t = 0 in the vacuum state and the atoms shared a single excitation, the wave function of the system at time t > 0, written in the interaction picture is of the form where C 1 (t) is the probability amplitude of the state in which atom 1 is in its excited state |e 1 , atom 2 is in the ground state |g 2 , and the field is in the vacuum state |{0} , C 2 (t) is the probability amplitude of the state in which atom 1 is in its ground state |g 1 , atom 2 is in the excited state |e 2 , and the field is in the vacuum state |{0} , and C λ (r, ω, t) is the probability amplitude of the state in which both atoms are in their ground states, |g 1 , |g 2 , and there is an excitation of the medium-assisted field |1 λ (r, ω) ≡f † λ (r, ω) |{0} . With the interaction (4), the Schrödinger equation (7) transforms into four coupled equations of motion for the probability amplitudes. When the amplitudes C λ (r, ω, t) are eliminated we arrive, as shown in Appendix A, into two coupled integro-differential equations for the probability amplitudes of the atomṡ in which K ii (t, t ′ ) is the integral kernel determined by the imaginary part of the one-point Green tensor, The kernels can be evaluated explicitly, and straightforward but lengthly calculations (for details see Appendix B) lead to the following explicit expressions where δ = ω s − ω a is the detuning of the atomic transition frequency from the plasma field frequency, is the coupling strength of the atoms to the surface plasma field, and determines the strength of the interaction between the atoms resulting from the coupling of the atoms with the same plasma field. Here, Γ A = ω 3 a p 2 a /(3ε 0 π c 3 ) is the spontaneous emission rate of the atoms in free space, assumed the atoms are identical, p a ≡ p 1 = p 2 . The function U (x 21 , z 0 ) depends on the separation between the atoms, x 21 , and also their distance z 0 from the interface. It determines the strength of the coupling between the atoms. In the limit of Thus, for large x 21 , the effects of the coupling between the atoms become negligible and the atoms evolve independently. Notice that at x 21 = 0 the function U (x 21 , z 0 ) is always unity independent of the value of z 0 , the distance of the atoms from the interface. However, the variation of U (x 21 , z 0 ) with x 21 depends strongly on z 0 . This is illustrated in figure 3 which shows U (x 21 , z 0 ) as a function of x 21 for several different values of z 0 . It is seen that for large z 0 , the function U (x 21 , z 0 ) varies slowly with x 21 . In that case, the indirect coupling between the atoms which is provided by the plasma field is effectively quite strong even at large separations. On the other hand, for small z 0 , (z 0 ≪ λ s ), the function U (x 21 , z 0 ) is different from zero only over very small distances x 21 and decays rapidly to zero as x 21 increases.
Thus, a strong coupling of the atoms to the plasma field destroys the collective behavior of the atoms. In the physical terms, the location of the atoms at a small distance z 0 from the interface leads to a strong spatial confinement (localization) of the surface plasmon fields around x 1 and x 2 resulting in a weak overlap of the surface fields produced by the atoms.

III. COLLECTIVE DYNAMICS OF THE ATOMS COUPLED TO SURFACE PLASMA
Having available the explicit expressions for the integral kernels, we can now solve Eqs. (9) and (10) and study the time evolution of the atomic system. If we introduce symmetric and antisymmetric combinations of the probability amplitudes, C s = (C 1 + C 2 )/ √ 2 and C a = (C 1 − C 2 )/ √ 2, corresponding to collective symmetric and antisymmetric states of the two-atom system, we readily find that Eqs. (9) and (10) simplify toĊ where Here, ] are coupling strengths of the symmetric and antisymmetric states to the plasma field, respectively. Clearly, the coupling strengths of the collective states are altered by the atomic interaction with Ω s enhanced and Ω a reduced by U (x 21 , z 0 ). This may have an interesting effect on the dynamics of the atoms that at small distances between the atoms, at which U (x 21 , z 0 ) ≈ 1, the antisymmetric state could be completely decoupled from the interaction with the plasma field leaving only the symmetric state to be strongly coupled to the field.
It is seen from Eqs. (15) and (16) that the equations of motion for the probability amplitudes C s (t) and C a (t) are independent of each other and are similar in form. We therefore need to obtain the solution for C s (t), and then the solution for C a (t) can be obtained simply by replacing Ω s by Ω a . From Eq. (17) it follows that the kernel K s (t, t ′ ) is a function only of the time difference t−t ′ . Therefore, the integro-differential equation (15) can be solved exactly by Laplace transformation. Thus if we get from Eq. (15) or By inverse Laplace transformation we then have the result s is the effective Rabi frequency of the interaction of the atoms with the plasma field, u = γ 2 + 4δ 2 and θ = arctan(2δ/γ). Similarly, for C a (t), we get whereΩ a = 1 16 (ue iθ ) 2 − Ω 2 a . Two important features of the results (21) and (22) should be noted. Firstly, the effective Rabi frequenciesΩ s andΩ a exhibit a threshold effect that depending upon Ω s(a) < u/4 or Ω s(a) > u/4, the Rabi frequencies can be either purely real or purely imaginary. In the other words, the time evolution of the probability amplitudes could be either exponential or sinusoidal. Secondly, we note that below threshold the time evolution of both C s (t) and C a (t) involves two decaying exponentials with a reduced (subradiant) decay constant, 1 4 γ −Ω s(a) , and an enhanced (superradiant) decay constant, 1 4 γ +Ω s(a) . The involvement of the fast and slow decay rates in the evolution of both superradiant and subradiant states seems to contradicts our expectation since, according to Dicke [63] (see also Refs. [64][65][66]), each of the collective amplitudes of a two atom system should decay with a single rate: The symmetric superposition C s (t) should decay with the fast (superradiant) rate while C a (t) should decay with the slow (subradiant) rate. FIG. 4. (Color online) Two atoms located at a distance z0 from the interface between two materials and their images located at a distance z0 behind the interface. The atoms are not directly coupled to each other, but can be coupled by the radiation reflected from the interface. A photon emitted by atom 1 and reflected from the interface towards atom 2 can be viewed as being emitted from the image of the atom 1.
A qualitative understanding of the involvement of both fast and slow decay rates in the evolution of C s (t) and C a (t) may be obtained by considering the interaction of the atoms with the plasma field as the interaction with images of the atoms located at a distance z 0 behind the interface. This is illustrated in Fig. 4, which shows that the two-atom system interacting with the surface plasma field can be seen as a four-qubit system, the two atoms plus two images. The radiation field emitted by either atom 1 or atom 2 and reflected from the interface in the direction normal to the interface can be regarded as the radiation from an image located at a distance z 0 behind the interface. The radiation field emitted by atom 1 and reflected from the interface towards atom 2 can be viewed as being emitted by the image of the atom 1 located at a distance z 0 behind the interface.
To study the evolution of the atoms in terms of the interaction with their images, we write Eqs. (9) and (10) in the following forṁ Following the Fig. 4, the second term on the right-hand side of Eq. (23) may be interpreted as arising from the coupling of the atom 1 to its image, whereas the third term may be interpreted as arising from the coupling of the atom 1 to the image of the atom 2. Thus, we can immediately write Eqs. (23) and (53) asĊ wherẽ are the probability amplitudes of the images of the atom 1 and 2, respectively. The equation of motion for the probability amplitudes of the corresponding images arė We focus on the evolution of the symmetric combinations of the probability amplitudes, which obey the equationṡ It is then straightforward to show that the solution of Eq. (30) for C s (t) is of the same form as Eq. (21). The involvement of the images allows us to write the general solution for C s (t) as a sum of two amplitudes whereD are symmetric and antisymmetric superpositions of the probability amplitudes of the atomic and image states, with A similar treatment can be applied to C a (t), which can be written in the form whereG with .
The reason for the presence of both the superradiant and subradiant terms in Eqs. (21) and (22) is now clear: The superradiant terms are associated with the decay of the symmetric superpositions involving the atomic and image states,D s (t) and G s (t), whereas the subradiant terms are associated with the decay of the antisymmetric superpositionsD a (t) andG a (t).
The slowest decay rate in the system, γ − a = 1 4 γ −Ω a , is the decay rate of the antisymmetric superpositionG a (t) whereas the fastest decay rate, γ + s = 1 4 γ +Ω s , is the decay rate of the symmetric superpositionD s (t).
We may conclude that the interaction of the atoms with the surface plasma field can be viewed as the interaction between the atoms and their corresponding images.

IV. MARKOVIAN AND NON-MARKOVIAN REGIMES OF THE EVOLUTIONS
We have already seen that different locations of the atoms lead to two different Rabi frequenciesΩ s andΩ a determining the evolution of the system and though two different time scales of the evolution. The forms ofΩ s andΩ a show a threshold effect for the Rabi frequencies that depending upon Ω s(a) < u/4 or Ω s(a) > u/4, the time evolution of the probability amplitudes could be either exponential or sinusoidal.
Since Ω s = Ω a , the threshold conditions for the evolution of the symmetric and antisymmetric states do not coincide with each other, and thus we can distinguish between three regions of Ω s and Ω a : (a) Ω s < u/4 and Ω a < u/4, (b) Ω s > u/4 and Ω a < u/4, (c) Ω s > u/4 and Ω a > u/4. Physically, the threshold values ofΩ s andΩ a separate what we can identify as the non-Markovian (memory preserved) regime from the Markovian (memoryless) regime of the evolution [67][68][69][70]. The case (a), in whichΩ s andΩ a are real so the Rabi frequencies contribute to the decay rates, that the exponentially decaying amplitudes of the symmetric and antisymmetric states are a manifestation of a Markovian evolution. In the case (b), the dynamics of the system are partly Markovian and partly non-Markovian. The symmetric state undergoes a non-Markovian whereas the antisymmetric state undergoes a Markovian evolutions. In the case (c), the dynamics of the system are fully non-Markovian that the amplitudes of both symmetric and antisymmetric states undergo an oscillatory evolution, which is a manifestation of a non-Markovian evolution. A non-Markovian evolution is a reversible process characterized by a flow (oscillation) of the information between the atoms and the field, but a Markovian evolution is a irreversible process of a flow (decay) of the information to the field. A possibility to control the non-Markovian dynamics is essential in quantum information technology since it plays a crucial role in preserving quantum memory.

A. Both Ωs and Ωa below threshold
Let us specialize Eqs. (21) and (22) to the case of exact resonance, δ = 0, and first examine the situation when the coupling of the atoms to the plasma field is weak, Ω 0 ≪ γ. In this case, bothΩ s andΩ a could be below threshold. If Ω s(a) < u/4, thenΩ s andΩ a are real and we see from Eqs. (21) and (22) that the Rabi frequencies contribute to the damping rates of the probability amplitudes. Physically, this is the kind of behavior corresponding to a Markovian evolution.
Since Ω s = Ω a , we see that the weak coupling of the atoms to the plasma field may result in the decay of the probability amplitudes C 1 (t) and C 2 (t) with four rates, two enhanced and two reduced rates. Figure 5 shows the time evolution of the populations P 1 (t) = |C 1 (t)| 2 and P 2 (t) = |C 2 (t)| 2 for bothΩ s and Ω a below threshold. At early times, the population P 1 (t) decreases whereas P 2 (t) increases until the populations become equal. At that time, the populations began to decay monotonically. The rate they decay is equal to γ − a , the slowest decay rate of the antisymmetric state. Since the atoms are very strongly coupled to each other, U (x 21 , z 0 ) = 0.95, the decay rate γ − a ≈ 0.002γ. Consequently, the effective decay time of the populations can be very long. The decay of the populations is irreversible so the evolution of the system is Markovian.

B. Ωs above threshold and Ωa below threshold
At large values of U (x 21 , z 0 ) it may happen that Ω a < u/4 and Ω s > u/4 even if the atoms are weakly coupled to the plasma field, i.e. Ω 0 < u/4. For example, when U (x 21 , z 0 ) ≈ 1, we have Ω a ≈ 0 and Ω s ≈ 2Ω 0 . Hence, Ω s can be larger than u/4 even if Ω 0 < u/4 and at the same time Ω a can be smaller than u/4. For Ω s > u/4 the Rabi frequencyΩ s is purely imaginary, and then the time evolution of the probability amplitude C s (t) takes the form whereΩ s = Ω 2 s − 1 16 (ue iθ ) 2 . The temporal evolution of the C s (t) is sinusoidal whereas the temporal evolution of the amplitude C a (t), which is below threshold, is exponential and is given in Eq. (22). In this case, the symmetric mode evolves in the non-Markovian regime whereas the antisymmetric mode evolves in the Markovian regime. It follows that in this case each atom evolves under the simultaneous influence of Markovian and non-Markovian mechanisms. Figure 6 shows the evolution of the populations for Ω s above and Ω a below the threshold of 1 4 γ. At early times, the oscillations of the populations with the Rabi frequency of the symmetric mode are clearly visible. At early times the populations oscillate with the Rabi frequency of the symmetric mode. In other words, the evolution of the populations is reversible but the reversibility occurs in a restricted time range t < 1/γ − a . Beyond t ∼ 1/γ − a the populations decay monotonically that the evolution is irreversible. Thus, we can clearly distinguish between the non-Markovian and Markovian regimes of the evolutions. We see that the upper limit on time of the reversible evolution results from the presence of the interaction between the atoms. Clearly, it is a collective effect. Physically, it is a consequence of the fact that a large part of the population is trapped in the asymmetric state determined by the amplitudeG a (t) thereby lowering the strength of the coupling of the atoms to the plasma field. It is easy to see, since for small distances between the atoms U (x 21 , z 0 ) ≈ 1, we have Ω a ≈ 0, which means that the antisymmetric states decouple from the plasma field. This example also shows that the atoms when behaving collectively can be weakly coupled to the field even if individually they are strongly coupled to the field.

C. Both Ωs and Ωa above threshold
Above the thresholds,Ω s andΩ a are purely imaginary. The time evolution of the probability amplitudes is then given by whereΩ a = Ω 2 a − 1 16 (ue iθ ) 2 . In this case, the time evolution of the probability amplitudes becomes sinusoidal. Such dynamics reflect the reversible property of the system that the evolution is non-Markovian. Figure 7 shows the time evolution of the populations P 1 (t) = |C 1 (t)| 2 and P 2 (t) = |C 2 (t)| 2 when the atoms are strongly coupled to the plasma field with z = 0.05λ s , δ = 0, but are weakly coupled to each other, U (x 21 , z 0 ) = 0.1. Note the presence of two characteristic time scales of the oscillations associated with the presence of two slightly different Rabi frequencies. At short times, t ≪ 1/Ω a , the initially excited atom 1 periodically exchange the excitation with the When the atoms are close to each other the collapses and revivals of the populations are absent. Instead, a periodic localization of the excitation is observed. This is illustrated in figure 8 which shows the evolution of the populations for a small distance between the atoms at which U (x 21 , z 0 ) = 0.8. The manner the populations oscillate is different for P 1 (t) and P 2 (t). We see a periodic localization of the excitation that even at long times the memory effects are still evident. The explanation of this feature follows from the observation that at small distances between the atoms, the time scale of the oscillations of the antisymmetric state is very large, approaches infinity when U (x 21 , z 0 ) → 1. Therefore, the system effectively evolves with a single time scale determined by the Rabi frequency of the symmetric state, t ∼ 1/Ω s .

V. EVOLUTION OF ENTANGLEMENT BETWEEN THE ATOMS
Given the time evolution of the probability amplitudes, we now proceed to evaluate the concurrence, a measure of entanglement between two qubits [65,66]. Following the definition of the concurrence, we find that in terms of the probability amplitudes, C 1 (t) and C 2 (t), the concurrence is given by an expression In terms of the amplitudes of the symmetric and antisymmetric combinations, the concurrence can be written as A positive value of the concurrence, C(t) > 0, indicates entanglement between the atoms, and C(t) = 1 corresponds to maximally entangled atoms. It is clear from Eq. (41) that the atoms are entangled whenever C s (t) = C a (t). Otherwise, the atoms are separable. Thus, to examine the occurrence of entanglement between the atoms we must look at differences in the evolution of the amplitudes C s (t) and C a (t).
If initially, C s (0) = C a (0), then according to the solutions Eqs. (21) and (22), the amplitudes will evolve differently only if U (x 21 , z 0 ) = 0. It then follows that the coupling between the atoms through the plasma field is necessary to create entanglement between the atoms from an initial separable state. The features of the concurrence for the three regions of Ω s and Ω a are illustrated in Figs. 9 -11. Figure 9 shows the effect of increasing interaction strength between the atoms on the concurrence for a weak coupling of the atoms to the plasma field, both Ω s and Ω a below threshold, which corresponds to a Markovian evolution of the system. In figure 9(a) the system starts from the separable state |e 1 |g 2 , whereas in figure 9(b) the initial state of the system is the maximally entangled state |s . We see that even in the weak coupling regime, a large and long living entanglement can be created between the atoms. The entanglement created from the initial separable state increases with an increasing U (x 21 , z 0 ) and attains the maximal value of C = 0.5 for U (x 21 , z 0 ) ≈ 1. The behavior of the concurrence is similar to that noted in the decay of two atoms into a common Markovian reservoir [66].
When the system starts from a maximally entangled state, either |s or |a , the initial entanglement always decays to zero with no entanglement present at long times, as illustrated in figure 9(b). This is readily understood if it is recalled that the symmetric and antisymmetric states evolve independently in time. Thus, if the population of the antisymmetric state was initially zero it will remain zero for all times. In this case, the system of two atoms effectively behave as a single two-level system with the upper state |s and the ground state |g . Then, the initial population of the state |s decays exponentially to the ground state with the rate 2Ω 2 s /γ. Turning now to the case of a strong coupling of the atoms to the plasma field at which Ω s and Ω a are above their thresholds, we show in figure 10 the evolution of the concurrence for weakly (U (x 21 , z 0 ) = 0.1) and strongly (U (x 21 , z 0 ) = 0.8) interacting atoms. The interaction creates a small difference between the frequencies of the oscillation of the symmetric and antisymmetric modes thatΩ s =Ω a . The frequency difference induces beating oscillations of the populations of the atoms, as was seen in figure 7, and one can see from figure 10 that these beating oscillations are rendered visible as beats in the concurrence.
Interesting features of the entanglement also appear when the symmetric mode evolves at Rabi frequency which is above the threshold, Ω s > u/4, and simultaneously the antisymmetric mode evolves at Rabi frequency which is below the threshold, Ω a < u/4. Under this circumstance, the probability amplitude of the symmetric mode is determined by Eq. (38) whereas the amplitude of the antisymmetric mode is given by Eq. (22).  Figure 11 shows the evolution of the concurrence for this special case. We see that the concurrence is zero only at the initial time t = 0. As time progresses the concurrence develops to a nonzero value. The concurrence never becomes zero as time develops, and thus no periodic entanglement quenching occurs. This feature is associated with the fact that with the Rabi frequency Ω a < u/4, the population of the antisymmetric state does not evolve in time leading to a trapping of a part of the atomic populations in their energy states. This is illustrated in figure 12, which shows the time evolution of the population P 1 (t) and P 2 (t) for the same parameters as in figure 11. We see from figure 12(a) that the initial population is periodically transferred between the atoms. However, the transfer is not complete that the populations of the atoms never become zero during the evolution. A part of the population is trapped in the atoms and is not transferred between them. Figure 12(b) shows the evolution of the population P 1 (t) for two initial maximally entangled states, |s and |a . Since in this case P 2 (t) = P 1 (t), we clearly see that the concurrence, if starts from maximally entangled state, it follows the evolution of the population of the atoms.

VI. OFF-RESONANT COUPLING
To this end we have discussed the collective effects induced by the resonant interaction of the atoms with the plasma field. We have established the importance of the images of the atoms in the atomic dynamics. Moreover, we have demonstrated the equivalence of the system with that of four interacting atoms. We now turn to the off-resonant case of the atomic transition frequencies strongly detuned from the plasma frequency, δ ≫ γ, Ω s . Under such condition, the effective Rabi frequencyΩ s can be approximated bỹ Thus, if in Eq. (21) the effective Rabi frequency is replaced by (42), we get, up to terms of order Ω 2 s /δ 2 , A similar expression with s → a gives C a (t). We see that C s (t) is composed of fast and slow oscillating terms varying in time with frequencies δ and Ω 2 s /δ, respectively. Of the two terms it is the one of the small magnitude (Ω 2 s /4δ 2 ) arising from the presence of the images. Thus, the evolution of C s (t) is well determined without much contribution of the images. It is particularly well seen from Eq. (33) that in the limit of δ ≫ Ω s , cos φ ≈ 1 (sin φ ≈ 0) so that the superposition amplitudeD s (t) = 0 andD a (t) is reduced to the atomic amplitudeC s (t). In other words, two atoms significantly detuned from the plasma field are coupled each other by exchanging virtual photons through a short interaction time with the plasma field.
The above considerations are illustrated in figure 13, which shows the time evolution of the atomic populations for δ = 50γ. We see that the atoms exchange the population with frequency 2Ω 2 0 /δ. The fast oscillations seen in the early time of the evolution occur at frequency δ and can be attributed to the involvement of the images in the dynamics of the system. The presence of the fast oscillations only at very early times of the evolution is also consistent with energy-time uncertainty arguments. It is easy to understand. At short times the uncertainty of the energy of the atoms and the plasma field is very large, so one cannot distinguish between ω a and ω s . This results in the presence of the fast oscillation of frequency δ. As time progresses, the frequencies become more distinguishable resulting in the disappearance of the fast oscillations.
It is interesting to contrast the entanglement created at δ = 0 with that created in the resonant case of δ = 0. We have seen in Sec. V that in the resonant case the maximal entanglement which can be created between the atoms from an initial separable state cannot exceed C = 0.5. For offresonant case (δ = 0), however, the entanglement can be significantly enhanced. This is shown in figure 14, which illustrates the time evolution of the concurrence for a large detuning δ. Clearly, at early times of the evolutions the concurrence is larger than 1/2, increases with an increasing U (x 21 , z 0 ) and becoming as large as C = 1. This behavior can be explained in terms of the energy-time uncertainty relation. At short times a large uncertainty in the energy results in a large uncertainty in the localization of the excitation. We see that the increased possibility to distinguish between the frequencies of the atoms and the plasma field results in an enhanced entanglement between the atoms.
As a final remark, we would like to comment about a potential experimental system in which the collective dynamics of the atoms could be observed. It could be done in experiments similar to that of Refs. [23,71,72], where a strong coupling between an artificial atom embodied into a material structure composed of MN and EN meta-materials was observed. The strong coupling was observed as the Rabi oscillations of the temporal evolution of the electric field inside the artificial atom after being excited by a short pulse. The experimental setup could be modified by embodying two artificial atoms into the composed meta-material structure and observe the Rabi oscillation of the electric field of the atoms. The presence of the second atom could lead to the modulation of the Rabi oscillations of the population of the first atom, as seen in figures 5-8, which would be the clear evidence of the collective behavior of the atoms.
To clarify the role of the SPP in the collective behavior of the atoms, we now consider the emission properties near the interface of an atom, represented by its oscillating dipole p i , and calculate the electric field at position r emitted by the atom located at r i . The field is given by [50] Suppose that the dipole moment p i is polarized in the x − z plane, p i = p a (x +z). Then the coupling of atom j, located at an arbitrary position r, to the field produced by the atom i is maximal if the dipole moment p j is parallel to p i . In this case, the magnitude of the field is where Ω 0 and U (x−x i , z i ) are given in Eqs. (13) and (14), respectively. Clearly, the function U (x − x i , z i ) determines the distribution of the field produced by atom i. It is seen from the expression (45) that the field is distributed in the x − z plane and the distribution depends not only on the distance x − x i along the interface but also on the distance z i of the radiating dipole from the interface. As illustrated in figure 3, the variation of the function U (x−x i , z i ) with x−x i depends strongly on the distance z i ≡ z 0 that U (x − x i , z i ), so that the field distribution, decreases with an decreasing distance z i ≡ z 0 . Hence, z 0 should not be too small in order to achieve a strong coupling between distant atoms through their interaction with the SPP propagating along the interface.

VII. CONCLUSIONS
We have studied the dynamics of two two-level atoms located near to the interface of two meta-materials; one of negative permeability and the other of negative permittivity. We have derived analytical expressions for the probability amplitudes of the atomic states valid for an arbitrary initial state, arbitrary strengths of the coupling constants of the atoms to the plasma field, and arbitrary distances between the atoms. We have shown that the effect of the plasma field is to produce several interesting features, such as (1) two different time scales of the evolution of the atomic states, one corresponding to the evolution of the collective symmetric state and the other to the antisymmetric state. The existence of the two evolution time scales results in an entanglement between the atoms even in long times. (2) A threshold behavior of the coupling constants of the atoms to the plasma field which distinguishes between the non-Markovian and Markovian regimes of the evolutions. We have shown that the collective behavior of the atoms may lead to three different regimes of the evolution; fully Markovian, simultaneous Markovian and non-Markovian, and fully non-Markovian evolutions. The three regimes determines three different time scales of the evolution of the memory effects and entanglement. (3) In the case of the resonant coupling of the plasma field to the atoms, the plasma field does not appear as a common reservoir to the atoms. We have adopted the image method and showed that in the resonant case the dynamics of the two atoms are completely equivalent to those of a four-atom system. (4) In the limit of a strong detuning of the plasma field from the atoms the dynamics resemble those of two atoms coupled to a common reservoir. If we introduce the notation we then easily find that Eqs. (50) and (51) simplify to Eqs. (9) and (10).

APPENDIX B
In this Appendix, we evaluate the integral kernels of the integro-differential equations (9) and (10). The kernels involve the imaginary part of the Green tensor ↔ G (r, r i , ω). The tensor when evaluated at an arbitrary space point r, distance R = |r−r i | from the atom located at r i , can be written as [30] where k 1 = √ ε 1 µ 1 ω/c represents the wave number, ε 1 and µ 1 are permittivity and permeability of the MN slab in which the atoms are located, and ↔ I =xx+ȳȳ+zz is the unit dyadic. Using the Weyl's expansion [50] in which β 1 is the z component of the propagation vector, we may express the Green tensor in terms of a two dimensional Fourier transform where ↔ G (k , ω, z, z 0 ) is the Green tensor on a plane (k , z) with constant z coordinate, k = k xx + k yȳ and r = r xx + r yȳ are components, respectively, of the wave vector and the position vector in a plane parallel to the interface, the x − y plane.
The expression (56) allows us to evaluate the the Green tensor in terms of plane waves incident on and reflected from the boundaries between different materials, including the interface between the MN and EN materials and the boundaries between the materials and the exterior regions (vacuum) on either side. Since the atoms are located in the MN material, so that their radiative properties are modified by the field existing inside the material, we will evaluate the Green tensor only at points r inside the MN material. We follow the procedure of Tomaš [50] in evaluating the Green tensor.
The presence of the boundaries results in the field inside the MN material consisting of waves propagating in both the +z and −z directions. Therefore, ↔ G (k , ω, z, z 0 ) can be written in terms of functions U ± q (k , ω, z) defined by imposing boundary conditions in the z direction where the functions U ± q (k , ω, z) describe the electric field in MN slab, with unit strength incident from its upper side (by taking symbol '−') or lower side (by taking symbol '+'), that can be categorized into TM (q = p) and TE (q = s) types of even (ξ p = 1) and odd (ξ s = −1) symmetries in the ±z directions, and Θ(z) is the unit step function. The forms of the functions U ± q (k , ω, z) for the field inside the MN material can be represented in terms of a sum of incident and reflected waves as where e ± p (k ) = (∓β 1k + k z)/k 1 and e ± s (k ) =k ×z are orthonormal polarization vectors of the electric field of p and s polarized waves, respectively;k is the unit vector in the direction of k (k = k k ),x,ȳ andz are unit vectors in the Cartesian coordinates, r q ± are reflection coefficients of the waves propagating in the ±z directions, and results from summing the geometrical series due to the multiple reflections from the boundaries between different materials.
We now proceed to evaluate the components of the Green tensor in Cartesian coordinates, which are given by where n, m = x, y, z, and n, m =x,ȳ,z. Thus, if we apply the explicit forms of the polarization vectors, and use the polar representation for k , k = k (cos φx + sin φȳ), the diagonal components of the Green tensor ↔ G (r, r i , ω) evaluated at a point r near the position r i of the ith atom are then G xx (r, r i , ω) = iµ 1 2(2π) 2 dk k and the off-diagonal components are G xy (r, r i , ω) = G yx (r, r i , ω) = iµ 1 2(2π) 2  cos φ k 2 1 D p r p + e −iβ1(z+z0−2d1) − r p − e iβ1(z+z0) e ik ·(r−ri) , G yz (r, r i , ω) = −G zy (r, r i , ω) = iµ 1 2(2π) 2 dk k 2 × dφ sin φ k 2 1 D p r p + e −iβ1(z+z0−2d1) − r p − e iβ1(z+z0) e ik ·(r−ri) , where R (q) The four terms appearing in Eq. (63) represent waves propagating inside the MN material. Figure 15 illustrates the source of those four terms in the expression (61). The term exp[iβ 1 (z − z 0 )] represents a wave propagated a distance z − z 0 from the source atom located at z 0 . This term is independent of the presence of the interface and the boundaries. Physically, it corresponds to the source field produced by an atom located at z 0 . The term r − exp[iβ 1 (z + z 0 )] represents a wave propagated the distance z + z 0 after the reflection from the bottom interface of the MN material. The term r + exp[−iβ 1 (z + z 0 − 2d 1 )] represents a wave propagated the distance 2d 1 − z − z 0 after the reflection from the upper interface of the MN material. The final term r − r + exp[−iβ 1 (z − z 0 − 2d 1 )] represents a wave propagated the distance 2d 1 + z 0 − z after two reflections, one from the bottom and the other from the upper interfaces.
At this point it should be stressed that the expressions (61) and (62), although evaluated in the presence of the boundaries, they contain terms which are independent of the boundaries. The reason is in the fact that the field is not completely bounded into the area inside the materials. The slabs have finite sizes in the z direction and in the derivation of Eqs. (61) and (62) it has been assumed that there are nonzero transmission coefficients at the boundaries with the exterior vacuum regions. Hence, we may consider the Green tensor as a sum of two terms, a source-field part ↔ GS (r, r i , ω) and a scatteredfield part ↔ GB (r, r i , ω) as [30,50]  field scattered from the interface and boundaries. The source term has the same properties that would apply to the free (unbounded) field.
To evaluate the integrals over φ, which appear in Eqs. (61) and (62) and involving exp[ik · (r − r i )], we assume, for simplicity, that r − r i has only x component so that we can write the dot product in the form k · (r − r i ) = α i cos φ, where α i = k |r − r i |. Hence, we arrive at the following expressions for the diagonal components