Directional dipole radiations and long-range quantum entanglement mediated by hyperbolic metasurfaces

Wei Fang; Wei Fang; Yaping Yang

doi:10.1364/OE.401628

1. Introduction

It is known that in the quantum scheme, the interaction of individual atoms with vacuum fluctuations leads to the spontaneous emissions and Lamb shifts [1]. As pointed out by Dicke, for an atomic sample with smaller size compared to the resonant wavelength, the identical excited atoms could synchronize to emit light coherently, which is often called supperradiance [2]. Since the emerging of this concept, the superradiant property in a variety of systems have been intensely investigated [3–6]. Furthermore, it has been proved that the spontaneous decay of collective state suffers a suppression when supperadiant effect becomes prominent, which could be used to generate long-live entanglement in decoherence-free subspace [7,8]. While initial attention has been focused on the local interactions between close quantum emitters (QEs), entangling arbitrary pairs of distant QEs is still left as a challenge. For example, quantum entanglement with an interemitter separation ranges from subwavelength to several wavelengths were observed in different systems [9–11]. Besides, it has been stressed that the generation of remote quantum entanglement via superradiance can be realized when a pair of QEs strongly coupled to the reservoir [12]. This important point of view indicates the significance in constructing light-matter interactions, which might allow effective controls over the functionality of complex systems. In fact, for a certain class of systems consist of interacting QEs through surfaces, strong light-matter couplings can be achieved [13–15] owing to the excitation of surface plasmon field (SPF) [16,17]. By employing plasmonic nanostructures, not only the efficient energy transfer, but also quantum entanglement between two or more spatially separated QEs can be realized [18–22]. Although great progress has been made, quantum entanglement based on the platforms mentioned above still lacks of selective addressing of individual QE, which is an essential ingredients for implementing quantum computation algorithms. In addition, spatial scalability to macroscopic distances turns out to be unachievable owing to the short propagation length of the SPF and large material losses.

Recently, the development of nanotechnology makes it possible to fabricate artificial materials exhibiting a hyperbolic regime, which often called hyperbolic metamaterials (HMs). Among them, three-dimensional HMs (TDHMs) are highly anisotropic materials that exhibit an inductive response along one (or two) of the optical axes and a capacitive response along the left optical axes, which have attracted a lot of attention since the proposal of the HMs [23–28]. An important feature is that the high-$k$ modes, a kind of electromagnetic field with very large wave vectors and photonic density of states, can propagate inside the bulk of the TDHMs over long distances [24]. Thus the TDHMs are promising for a wide range of applications such as nanoscale resonators [29,30], negative refraction [31,32], hyperlensing and focusing [33,34], emission engineering [35,36] and the enhancement of dipole interactions between QEs [37]. However, the complexity of bulk fabrication, associated with the losses induced by the volumetric effects, hinder the further implementation of TDHMs especially in the quantum domain. To circumvent these issues, two-dimensional version of the HMs – hyperbolic metasurfaces (HMSs), have been theoretically proposed and demonstrated in experiments [38–43]. Particularly, systems with anisotropic electronic and optical properties are promising candidates for the HMSs, such as graphene nanoribbon arrays [38,41], black phosphorus [44,45], van der Waals materials [42,46] and other nanostructures based on lithography and etching technologies [43,47]. In most cases, HMSs are sheets with extreme sub-wavelength thickness that support directional guidings of highly confined plasmon-polaritons [40]. Because the SPF is guided at the surface rather than passing through bulks, HMSs have been predicted to suffer much low energy loss while still exhibiting optical phenomena akin to those in TDHMs. Meanwhile, HMSs have a significant impact on optical designs since their unique structures enable the realization of virtually flat optics via replacing bulky optical components with ultrathin planar elements [41], which in turn provides the possibility in the integration with planarized systems compatible with integrated circuits.

In this paper, we propose a novel scheme to overcome the limitations of traditional three-dimentional plasmonic structures by applying the unique capability of the HMSs that enables spatial scalability to macroscopic distances. In the begining, we investigate the radiation properties of a vertically polarized QE coupled to the HMSs with the Green function method. We focus on the HMSs that support the propagation of the SPF, which are often made of graphene-based periodic structures. The Purcell factor of the QE clearly illustrates different light-matter coupling regions dominated by free-space radiations, plasmon-induced radiations and the quenching effect, respectively. Then by working in the plasmonically coupled region, we show that the radiation energy is channeled as very narrow beams. Meanwhile, the propagating direction of the beams can be controlled by engineering the anisotropy of the HMS, which is experimentally achievable with an electric bias [48–50]. Base on these properties, prominent superradiance between two distant QEs mutually coupled to the SPF can be realized, where an increase in the optical response of the HMS is constructive in supporting long-range superradiance. Furthermore, it is shown that large degree of quantum entanglement, as well as remarkable revivals from sudden death can both be observed at macroscopic distance.

The paper is organized as follows. In Sec. 2, we introduce the model and present the Green function rigorously describing the spatial and spectral distribution of the electromagnetic field in the dissipative reservoir formed by the HMS. Then we study the Purcell factor and radiation pattern of a QE coupled to the HMS. In Sec. 3, we focus on the superradiant behavior of two QEs mutually coupled to the SPF. The master equation for the reduced density operator of the QEs and the general forms of the cross interactions are presented in Sec. 3.1. In order to clarify the influence of the material response on the supported plasmon modes and long-range dipole interactions, both the cross decay and the collective behavior of the QEs system have been investigated, which is the subject of Sec. 3.2. In Sec. 4 we consider different initial states of the system and study the transient entanglement between the QEs with a macroscopic separation. Finally we present our conclusions in Sec. 5 and some essential concepts are clarified in appendixes to better understand the work.

2. Radiation properties of a QE near the HMS

Considering a system consists of a two-level QE coupled to a graphene-based hyperbolic surface. The QE has a distance $l$ above the metasurface, where the transitions can take place between its excited state $\left |e_{1}\right \rangle$ and ground state $\left |g_{1}\right \rangle$. The transition dipole is assumed to be polarized along the $z$ axis and as Fig. 1 indicated, the radiation field could be well guided on the metasurface by working in the surface plasmon region (we will give detail discussions on different light-matter coupling regions later). It should be stressed that there are various ways to realize hyperbolic dispersions supported by surfaces, and the structure mentioned in one of the relevant works [38] has been adopted in our figure to illustrate the metasurface. Owing to the interplay between anisotropic intraband and interband motions, the conductivity tensor of the anisotropic metasurface possesses nonvanished diagonal elements $\overleftrightarrow {\sigma }=\textrm {diag}\left (\sigma _{xx},\sigma _{yy}\right ),\left (\sigma _{xx}\neq \sigma _{yy}\right )$ and it can be tuned on the spot by applying a gate doping [50,51]. In what follows we are interested in the hyperbolic dispersion case, which arises when the surface exhibits an inductive response along one direction and with a capacitive response along the orthogonal one (i.e. $\textrm {Im}[\sigma _{xx}]\textrm {Im}[\sigma _{yy}]<0$, where $\textrm {Im}[\sigma _{jj}](j=x,y)$ denotes the imaginary part of the conductivity).

Fig. 1. A two-level QE near the graphene-based HMS, experiencing the spontaneous emission and mainly emits photons into the directional SPF guided by the surface. The QE has a distance $l$ to the surface with its coordinate denoted by $\mathbf {r}_{1}=\left (0,0,l\right )$, and the transition dipole is assumed to be normal to the surface.

Download Full Size | PDF

It is known that for a distance much larger than the vacuum wave length, the interaction between dipoles through the vacuum modes is negligible. In this scenario, the cross interaction between two QEs could only takes place by mutually coupling to different optical structures [18,19], which is mainly governed by the scattered field. For a point dipole located at $\mathbf {r}_{d}=\left (x_{d},y_{d},z_{d}\right )$, the electric field including the scattered field mediated by the surface can be evaluated via $\mathbf {E}\left (\mathbf {r},\mathbf {r}_{d},\omega _{d}\right )=\omega _{d}^{2}c^{-2}\overset {\leftrightarrow }{\mathbf {G}}\left (\mathbf {r},\mathbf {r}_{d},\omega _{d}\right )\cdot \mathbf {u}_{d}$ [52]. In the above expression, $\mathbf {u}_{d}$ represents the dipole moment, $\omega _{d}$ is the resonant frequency of the dipole transition, and $\mathbf {r}=\left (x,y,z\right )$ denotes the target point where the electric field is observed. Thus it is easy to prove that the field component normal to the surface (denoted by $\mathbf {E}_{z}\left (\mathbf {r},\mathbf {r}_{d},\omega _{d}\right )$) induced by a $z$-oriented dipole is determined by the Green function $\overset {\leftrightarrow }{\mathbf {G}}_{zz}\left (\mathbf {r},\mathbf {r}_{d},\omega _{d}\right )$, which is the $zz$ diagonal element of the relevant Green tensor. After enforcing the boundary conditions and following the methods mentioned in Ref. [53], in the presence of the metasurface the Green function is found to have the form

(1)$$\overset{\leftrightarrow}{\mathbf{G}}_{zz}\left(\mathbf{r},\mathbf{r}_{1},\omega_{a}\right) = \frac{i}{2\pi}\int_{0}^{\infty}\frac{k^{3}e^{i\beta_{0}\left|z-z_{1}\right|}J_{0}\left(k\xi\right)}{2\beta_{0}}dk+\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-ik_{y}\left(y-y_{1}\right)}\left[f_{2}\left(k_{y}\right)+k_{y}^{2}f_{1}\left(k_{y}\right)\right]dk_{y},$$

where for simplicity and without losing generality the metasurface has assumed to be suspended in the vacuum. In practice, owing to the tunnelling effect of the free electrons, the presence of the substrate might bring changes in both the doping level and the scattering rates of transitions [54]. However, these changes only modify the spectral weight as well as the absorption (also emission) intensity and do not really affect the physical mechanism of light-matter couplings. Here $\omega _{a}$ is the transition frequency between the excited state and ground state of the QE, which locates in the infrared region [55,56]. $\xi =\sqrt {\left (x-x_{1}\right )^{2}+\left (y-y_{1}\right )^{2}}$ is the projection of the distance between the QE and the observation point in the $x-y$ plane, $k$ is the radial wave number and $\beta _{0}=\sqrt {k_{0}^{2}-k^{2}}$ is the $z$-component of the wave number in the upper space ($k_{0}=\sqrt {\varepsilon _{0}\mu _{0}}\omega _{a}$, $\varepsilon _{0}$ and $\mu _{0}$ are the vacuum permittivity and permeability, respectively). As implied by Eq. (1), it is obvious that the total Green function at the spatial point $\mathbf {r}$ originates from the superposition of the propagating part directly through the vacuum and the scattered part mediated by the metasurface. For anisotropic case, the inner integrals for the scattered part are

(2)$$f_{1}\left(k_{y}\right) = \frac{i}{2\pi}\int_{-\infty}^{\infty}R_{zz}\left(k_{x},k_{y}\right)\frac{e^{ik_{z}\left(z+z_{1}\right)}}{2k_{z}}e^{-ik_{x}\left(x-x_{1}\right)}dk_{x}$$

and

(3)$$f_{2}\left(k_{y}\right) = \frac{i}{2\pi}\int_{-\infty}^{\infty}R_{zz}\left(k_{x},k_{y}\right)\frac{k_{x}^{2}e^{ik_{z}\left(z+z_{1}\right)}}{2k_{z}}e^{-ik_{x}\left(x-x_{1}\right)}dk_{x},$$

where the reflection coefficient is given by

(4)$$R_{zz}\left(k_{x},k_{y}\right) = \frac{k_{z}^{2}\left(\sigma_{xx}+\sigma_{yy}\right)-k_{z}k_{0}\eta\sigma_{xx}\sigma_{yy}}{2\sigma_{xx}\left(k_{0}^{2}-k_{x}^{2}\right)+2\sigma_{yy}\left(k_{0}^{2}-k_{y}^{2}\right)-k_{0}k_{z}\eta^{-1}\left(4+\eta^{2}\sigma_{xx}\sigma_{yy}\right)}.$$

In the above equation, $k_{x}$ and $k_{y}$ are metasurface-supported wave numbers along $x$ and $y$ directions, where the normal part of the wave number is $k_{z}=\sqrt {k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}$ and the parameter $\eta =\sqrt {\mu _{0}\backslash \varepsilon _{0}}$. In order to investigate the coupling between the QE and the metasurface under different circumstances, it is instructive to inspect the plasmon dispersion relation supported by the two-dimensional surface. As it has been pointed out, the plasmon dispersion relation can be obtained by finding the relevant solutions which would lead to an infinity in the reflection coefficient [57]. For lossless surface where the conductivity is purely imaginary (i.e. $\sigma _{xx}=\sigma _{xx}^{i}i$ and $\sigma _{yy}=\sigma _{yy}^{i}i$) and in the assumption of large wave numbers in the $x$-$y$ plane that $k_{x,y}\gg k_{0}$, we can arrive at the following dispersion relation

(5)$$\frac{k_{x}^{2}}{\sigma_{yy}^{i}}+\frac{k_{y}^{2}}{\sigma_{xx}^{i}} = \frac{i\omega k_{z}}{2}\left(\frac{4\varepsilon_{0}}{\sigma_{xx}^{i}\sigma_{yy}^{i}}-\mu_{0}\right).$$

Obviously, for a certain frequency constant and the conductivity parameters $\sigma _{xx}^{i}\sigma _{yy}^{i}<0$, Eq. (5) actually implies hyperbolic equifrequency surfaces. It should be noted that for anisotropic materials, the directions of the wave front and energy propagation do not coincide with each other. Instead, the transfer direction of the energy is often defined by the group velocity $\left .\nabla _{\mathbf {k}}\omega \left (\mathbf {k}\right )\right |_{\mathbf {k}=\left (k_{x},k_{y},0\right )}$ in the medium [58]. Considering the case of a dipole excites plasmon modes with different wave vectors $\mathbf {k}$, we can expect that the normals to the points on the hyperbola asymptotes would point in the same direction for a given sign of $k_{x}$, with the angle defined by $\theta _{pl}=\arctan \left (\sqrt {-\sigma _{yy}^{i}\backslash \sigma _{xx}^{i}}\right )$ (associated with the result given in Eq. (23)). Thus different from the isotropic case where the SPF equally propagate in all directions along the surface, plasmons supported by the HMS can be guided towards specific directions.

To inspect the mechanisms of the QE-metasurface coupling for an explicit distance, we introduce the Purcell factor ($\textrm {P}_{j}=\Gamma _{j}\backslash \Gamma _{0}$, $j=\textrm {T},\textrm {SP}$) to have a better understanding of the coupling formations between the QE and the electromagnetic environment. Here $\Gamma _{0}$ is the spontaneous decay rate in free vacuum, $\Gamma _{\textrm {T}}$ is the total atomic decay rate in the presence of the metasurface and $\Gamma _{\textrm {SP}}$ is the decay rate corresponds to the radiation into the surface plasmon modes. It has been pointed out that the imaginary part of the classical Green function relates to the mode density of the electromagnetic field [59], and thus can be used to study the spontaneous decay property of the QE. Starting from the definitions the Purcell factors can be expressed in forms of the Green function, as

(6)$$\textrm{P}_{\textrm{T}\left(\textrm{SP}\right)} = 6\pi k_{0}^{-3}\textrm{Im}\left[\overset{\leftrightarrow}{\mathbf{G}}_{zz}^{(sp)}\left(\mathbf{r}_{1},\mathbf{r}_{1},\omega_{a}\right)\right].$$

Here $\textrm {Im}\left [\overset {\leftrightarrow }{\mathbf {G}}_{zz}\left (\mathbf {r}_{1},\mathbf {r}_{1},\omega _{a}\right )\right ]$ and $\textrm {Im}\left [\overset {\leftrightarrow }{\mathbf {G}}_{zz}^{sp}\left (\mathbf {r}_{1},\mathbf {r}_{1},\omega _{a}\right )\right ]$ represent the imaginary parts of the Green functions defined in Eqs. (1) and (8), respectively. In fact the single point Green function reflects the mode density of state at the spatial point of the source, thus determines the spontaneous emission behavior of the QE. As we will see below, by investigating the Purcell factor it is convenient to distinguish different QE-metasurface coupling regions.

It has been demonstrated that with the increase in the carrier concentration of graphene-based heterostructures, the real part of the conductivity suffers decreases due to Pauli blocking, while the imaginary part of the conductivity grows as a result of the enhanced collective motion of free carriers [60]. That is to say, the graphene layer with high carrier concentration in the heterostructures can be served as metasurfaces that exhibit strong optical responses and small material dissipations. Moreover, decreases in the environment temperature can further enhance the material response and reduce the relevant dissipation [61]. In practice, the conductivity of the metasurfaces could possess a sufficiently small real part compared with its imaginary part in the plasmon region, where the latter appears to be highly tunable and ranges from $\mathrm {\mu }\textrm {S}$ to $\textrm {mS}$ [38,41,62]. Firstly consider the QE interacts with the HMS, where the conductivity components are $\sigma _{xx}=10^{-3}+0.5i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$ near the QE’s transition frequency. As shown in Fig. 2, different QE-metasurface coupling regions can be clearly distinguished. It should be noticed that the Purcell factor for the total decay (blue solid line) and the decay through the surface plasmon modes (red dashed-dotted line) are overlapped as the distance $l$ sweeps the interval $\left (10^{-3}\lambda _{a},10^{-2}\lambda _{a}\right )$, which indicates the complete coupling of the QE to the SPF supported by the metasurface. For distance larger than $10^{-2}\lambda _{a}$, the QE gradually decouples from the SPF and decays by radiating photons into the free vacuum modes. Besides, if the QE is close enough to the surface ($l<10^{-3}\lambda _{a}$), material dissipation turns out to be dominant and the corresponding decay rate is largely enhanced by several orders of magnitude owing to the quenching effect [63]. Further study in the inset panel shows that the Purcell factor suffers an increase by enlarging the imaginary part of the conductivity component $\sigma _{xx}$ to $1.5\:\textrm {mS}$, while other parameters are kept unchanged. In this scenario, according to the results given in Eqs. (5) and (20), the increase in the $\sigma _{xx}^{i}$ allows the surface to couple the incident light through the excitation of the surface plasmon modes that possess a broader distribution in the wave numbers. Therefore, the overall local density of states are increased, which in turn strengthens the light-matter coupling in the plasmon region.

Fig. 2. Purcell factor for the QE coupled to the HMS as a function of the distance (normalized to the vacuum wavelength $\lambda _{a}$), where the diagonal elements of the conductivity tensor are $\sigma _{xx}=10^{-3}+0.5i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$. The blue solid line denotes the total Purcell factor of the QE that contains the contributions from the radiation into the free-space, surface plasmon and quenching modes, and the red dashed-dotted line represents the Purcell factor of the QE decays only through the surface plasmon modes. The inset panel shows the Purcell factor for $\sigma _{xx}=10^{-3}+1.5i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$.

Download Full Size | PDF

We then study the spatial distributions of the vertical component of the electric field (with the intensity denoted by $\left |\mathbf {E}_{z}\right |$ and scaled by its vacuum value $\left |\mathbf {E}_{0}\right |$) excited by a $z-$polarized dipole placed at a distance $\lambda _{a}/500$ above from the HMS. It should be mentioned that for the emitter-surface distance we choose, the radiation energy can be well confined near the surface owing to the strong coupling of the QE to the surface plasmon modes (illustrated by the overlap region in Fig. 2). Figure 3(a) demonstrates the scaled radiation field intensity for the case of $\sigma _{xx}=10^{-3}+0.5i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$. Under this circumstance, it is clear that the excited SPF propagates along specific directions on the metasurface, with the angle $\theta _{pl}=45^{\circ }$ as expected. To further inspect how the conductivity of the metasurface determines the directionality of the propagating SPF, the spatial distribution of the scaled field intensity has been investigated with an enhancement in the optical response of the conductivity tensor ($\sigma _{xx}=10^{-3}+1.5i\:\textrm {mS}$). As shown in Fig. 3(b), the plasmon energy is channeled as narrow beams and propagates along the directions with an angle $\theta _{pl}=30^{\circ }$ deviates from the $x$ axis. Moreover, owing to the capability of exciting the surface plasmon modes with larger propagating wave lengths and lower losses (see Fig. 10 for more details), the field intensity can be well maintained during the propagation of the SPF when compared with the former case.

Fig. 3. Spatial distribution of the field intensity $\left |\mathbf {E}_{z}\right |$ (scaled by the vacuum value $\left |\mathbf {E}_{0}\right |$) excited by a normal dipole (green arrow) in the plane $z=l=\lambda _{a}/500$, the HMS is assumed to be placed in the plane $z=0$. The parameters of the surface are: (a) $\sigma _{xx}=10^{-3}+0.5i\:\textrm {mS}$, $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$. (b) $\sigma _{xx}=10^{-3}+1.5i\:\textrm {mS}$, $\sigma _{yy}=10^{-3}-0.5i\:\textrm {mS}$.

Download Full Size | PDF

It is obvious from the above discussions that when a QE strongly coupled to the SPF supported by the HMSs, the propagation of the radiation field exhibits evident directionality along the surface. In the mean time, one can modify the propagating properties of the plasmon beams by engineering the anisotropy of the metasurface. Thus in analogy to the super-Coulombic interaction in three-dimensional hyperbolic material [37], the strong coupling between the QE and the plasmon field, as well as the novel propagating properties of the proposed systems, may probably lead to long-range interactions between QEs in specific directions.

3. Collective behavior of two distant QEs coupled to the HMS

3.1 Master equation for the reduced density operator of the QEs

Consider two QEs near the metasurface, both at the same distance from the surface and decay only through the surface plasmon modes. According to the analysis in the previous section, long-range interaction between the QEs can be expected in specific directions owing to the existence of the guided SPF (the configuration of the system is illustrated by Fig. 4). In the following, we will investigate this interaction. The rotating-wave Hamiltonian of two identical QEs (placed at $\mathbf {r}_{1}$ and $\mathbf {r}_{2}$) interacting with the electromagnetic field in presence of the HMS, in the dipole approximation assumes the form [59]

(7)$$\hat{H} = \hat{H}_{0}+\hat{H}_{I},$$

where

(8)$$\hat{H}_{0} = \hbar\int d\mathbf{r}\int\omega\hat{\mathbf{f}}^{\dagger}\left(\mathbf{r},\omega\right)\hat{\mathbf{f}}\left(\mathbf{r},\omega\right)d\omega+\frac{1}{2}\hbar\omega_{a}\sum_{j=1,2}\hat{\sigma}_{j}^{z}$$

is the unperturbed Hamiltonian of the QEs and the HMS-assisted reservoir, and Hamiltonian

(9)$$\hat{H}_{I} = -\sum_{j=1,2}\left[\hat{\sigma}_{j}^{+}\int\mathbf{p}_{j}\cdot\hat{\mathbf{E}}^{\left(+\right)}\left(\mathbf{r}_{j},\omega\right)d\omega+\textrm{H.c.}\right]$$

denotes the interaction between the QEs and the medium-assisted field. Here the electromagnetic environment has been quantized in order to investigate the dynamics of the QEs subsystem according to the damping theory. The relevant bosonic operators $\hat {\mathbf {f}}^{\dagger }\left (\mathbf {r},\omega \right )$ and $\hat {\mathbf {f}}\left (\mathbf {r},\omega \right )$ represent the creation and annihilation of the elementary excitations in the metasurface-assisted reservoir [59], which include the contributions of both the vacuum and the metasurface environments [64]. Meanwhile, they satisfy the commutation relations $\left [\hat {\mathbf {f}}\left (\mathbf {r}_{n},\omega '\right ),\hat {\mathbf {f}}^{\dagger }\left (\mathbf {r}_{m},\omega \right )\right ]=\delta \left (\mathbf {r}_{n}-\mathbf {r}_{m}\right )\delta \left (\omega '-\omega \right )$ and $\left [\hat {\mathbf {f}}\left (\mathbf {r}_{n},\omega '\right ),\hat {\mathbf {f}}\left (\mathbf {r}_{m},\omega \right )\right ]=0$. The indices $j=1,2$ refer to the QE1 and QE2, with excited states and ground states denoted by $\left |e_{j}\right \rangle$ and $\left |g_{j}\right \rangle$, respectively. In this notation, $\hat {\sigma }_{j}^{+}\left (\hat {\sigma }_{j}^{-}\right )$ represents the raising (lowering) operator with the transition dipole moment $p_{j}$, and $\hat {\sigma }_{j}^{z}$ is the energy difference operator. The positive frequency part of the field can be expressed as

(10)$$\hat{\mathbf{E}}^{\left(+\right)}\left(\mathbf{r}_{j},\omega\right) = i\frac{\omega^{2}}{c^{2}}\sqrt{\frac{\hbar}{\pi\varepsilon_{0}}}\int d\mathbf{r}\sqrt{\textrm{Im}\left[\varepsilon\left(\mathbf{r},\omega\right)\right]}\overset{\leftrightarrow}{\mathbf{G}}\left(\mathbf{r}_{j},\mathbf{r},\omega\right)\cdot\hat{\mathbf{f}}\left(\mathbf{r},\omega\right),$$

where $\textrm {Im}\left [\varepsilon \left (\mathbf {r},\omega \right )\right ]$ denotes the imaginary part of the permittivity refers to the metasurface and $\overset {\leftrightarrow }{\mathbf {G}}\left (\mathbf {r}_{j},\mathbf {r},\omega \right )$ is the Green tensor of the electromagnetic field. Starting from Eqs. (8) and (9), by working in the interaction picture the interaction Hamiltonian can be obtained through the transformation $\hat {V}_{I}\left (t\right )=\exp \left (i\hat {H}_{0}t\backslash \hbar \right )\hat {H}_{I}\exp \left (-i\hat {H}_{0}t\backslash \hbar \right )$. Then applying the second-order perturbation theory [65] by tracing out the field variable through $\dot{\hat{\rho}}_{E} = -\hbar^{-2}\int_{0}^{t}{\textrm{T}{r_{R}}}\left\{ \left[\hat{V}_{I}(t),\left[\hat{V}_{I}(t-\tau), \hat{\rho}(t-\tau)\right]\right]\right\}$, under the Born-Markovian approximation the motion equation for the reduced density operator of the QEs system appeared to be

(11)$$\dot{\hat{\rho}}_{E} = \sum_{m,n=1}^{2}\int_{0}^{\infty}d\omega S\left(\mathbf{r}_{m},\mathbf{r}_{n},\omega\right)\int_{0}^{t}d\tau e^{i\left(\omega-\omega_{a}\right)\tau}\left\{ \left[\hat{\sigma}_{m}^{-}\hat{\rho}_{E}\hat{\sigma}_{n}^{+}-\hat{\rho}_{E}\hat{\sigma}_{m}^{+}\hat{\sigma}_{n}^{-}\right]+\textrm{H.c.}\right\} .$$

Here the photonic reservoir function $S\left (\mathbf {r}_{m},\mathbf {r}_{n},\omega \right )=\omega ^{2}\mathbf {p}_{m}\cdot \textrm {Im}\left [\overset {\leftrightarrow }{\mathbf {G}}\left (\mathbf {r}_{m},\mathbf {r}_{n},\omega \right )\right ]\cdot \mathbf {p}_{n}^{\ast }\backslash \left (\hbar \pi \varepsilon _{0}c^{2}\right )$ relates to the dipole-field ($m=n$) and dipole-dipole ($m\neq n$) couplings, which depends on the mode density of state at the QEs’ positions and the propagating field strength from one QE to the other, respectively. Low environment temperature has been assumed such that the trace over the reservoir degrees can be performed by $\textrm {Tr}_{R}\left [\hat {\mathbf {f}}^{\dagger }\left (\mathbf {r}_{n},\omega \right )\hat {\mathbf {f}}\left (\mathbf {r}_{m},\omega '\right )\right ]=0$ and $\textrm {Tr}_{R}\left [\hat {\mathbf {f}}\left (\mathbf {r}_{n},\omega \right )\hat {\mathbf {f}}^{\dagger }\left (\mathbf {r}_{m},\omega '\right )\right ]=\delta \left (\mathbf {r}_{n}-\mathbf {r}_{m}\right )\delta \left (\omega '-\omega \right )$. In the limit of large time scale, Eq. (11) can be rearranged in the following form

(12)$$\dot{\hat{\rho}}_{E} = \sum_{m,n=1}^{2}\left\{ \frac{\Gamma_{mn}}{2}\left(2\hat{\sigma}_{m}^{-}\hat{\rho}_{E}\hat{\sigma}_{n}^{+}-\hat{\sigma}_{n}^{+}\hat{\sigma}_{m}^{-}\hat{\rho}_{E}-\hat{\rho}_{E}\hat{\sigma}_{n}^{+}\hat{\sigma}_{m}^{-}\right)+i\delta\omega_{mn}\left[\hat{\sigma}_{n}^{+}\hat{\sigma}_{m}^{-},\hat{\rho}_{E}\right]\right\} .$$

In Eq. (12), one-point coefficients $\Gamma _{mm}$ and $\delta \omega _{mm}$ represent the spontaneous decay rates and Lamb shifts of different QEs, where two-point coefficients $\Gamma _{mn}$ and $\delta \omega _{mn}$ denote the cross decay rates and dipole-dipole shifts [66] owing to the mutual coupling with the SPF. They can be generally given in the forms

(13)$$\Gamma_{mn} = 2\pi S\left(\mathbf{r}_{m},\mathbf{r}_{n},\omega_{a}\right)$$

and

(14)$$\delta\omega_{mn} = \mathcal{P}\int_{0}^{\infty}\frac{S\left(\mathbf{r}_{m},\mathbf{r}_{n},\omega\right)}{\omega-\omega_{a}}d\omega,\:\left(m,n=1,2\right)$$

where symbol $\mathcal {P}$ represents the principal part of the relevant integrations. It should be mentioned that both the Lamb and the dipole-dipole shifts are determined by the principal values of the photonic reservoir function $S\left (\mathbf {r}_{m},\mathbf {r}_{n},\omega \right )$, they can be numerically evaluated by transferring the integration to the imaginary axis in the complex frequency plane [67]. However, as we will see below, the shifts mediated by the metasurface is typically at most on the order of the linewidth of the transition if the separation between the QEs is larger than the plasmon wavelength. Meanwhile, the energy shifts suffer a $\pi \backslash 2$ phase shift corresponds to the cross damping and thus can be disregarded when one focuses on the dynamics of two maximally-coupled QEs.

Fig. 4. Schematic diagram of two identical two-level QEs simultaneously coupled to the HMS. By assuming that the radiation energy of both the QEs totally transfers into the surface plasmon modes and channels as the propagating SPF along specific directions, a prominent interaction can be built between the QEs.

Download Full Size | PDF

3.2 Long-range superradiance of the QEs system

Superradiance is a phenomenon that illustrate the cooperative radiation of identical QEs interact with a reservoir, which often attributes to the collective effect of subsystems. It is known that in free space, the strength of superradiance suffers rapid declines with the increasing separation between atoms. Here we introduce the normalized decay factor $R$ to characterize the metasurface-mediated superradiance and subradiance between the QEs, which relates to the ratio of the total decay rate that two QEs interact through the metasurface to the decay rate of two uncoupled QEs in the presence of the metasurface, and has the form

(15)$$R = 1+\frac{\Gamma_{12}+\Gamma_{21}}{\Gamma_{11}+\Gamma_{22}}.$$

In the above expression, $\Gamma _{12}$ ($\Gamma _{21}$) denotes the cross decay rate of the dipole placed at the spatial point $\mathbf {r}_{1}$ ($\mathbf {r}_{2}$) owing to the presence of another dipole placed at $\mathbf {r}_{2}$ ($\mathbf {r}_{1}$), $\Gamma _{11}=\Gamma _{22}$ represents the spontaneous decay rate of the dipole placed at the spatial point $\mathbf {r}_{n}$, their expressions can be acquired through Eq. (13). For our model it is convenient to work in the Dicke state basis, where the QEs system behaves as a four-level system [66] with one excited state $\left |e\right \rangle =\left |e_{1},e_{2}\right \rangle$, one ground state $\left |g\right \rangle =\left |g_{1},g_{2}\right \rangle$ and two intermediate states denoted by the symmetric state $\left |s\right \rangle =\left (\left |g_{1},e_{2}\right \rangle +\left |e_{1},g_{2}\right \rangle \right )\backslash \sqrt {2}$ and antisymmetric state $\left |a\right \rangle =\left (\left |g_{1},e_{2}\right \rangle -\left |e_{1},g_{2}\right \rangle \right )\backslash \sqrt {2}$. Note that the symmetric and antisymmetric states are characterized by the enhanced decay rate $\Gamma _{11}+\Gamma _{12}$ and reduced decay rate $\Gamma _{11}-\Gamma _{12}$, respectively. Since $\Gamma _{12}=\Gamma _{21}$, Eq. (15) simply indicates that different signs of the cross interaction will lead to two distinguished cooperative radiation behaviors denoted by $R$. When $R>1$ ($\Gamma _{12}$ takes positive values), the interaction between the symmetric state and the reservoir is enhanced due to the presence of the metasurface, which leads to the appearance of the superradiance. Correspondingly, when $R<1$ ($\Gamma _{12}$ takes negative values), an inhibition of the dipole-dipole interaction occurs and thus the system is subradiant.

Before discussing the superradiant property of the system, it is beneficial to investigate the cross interaction of the QEs built through their mutual couplings with the SPF. Figure 5(a) displays both the cross decay rate and the dipole-dipole shift (normalized to the spontaneous decay rate, since the QEs have the same distance to the metasurface we have $\Gamma _{11}=\Gamma _{22}$), as a function of the separation between the QEs placed along the plasmon propagating direction. Obviously an oscillatory behavior of the cross decay with the amplitude smaller than the single-QE decay rate can be observed, where the period of the oscillations is affected by the wavelengths of the supported surface plasmon modes. The dipole-dipole shift undergoes similar oscillations for the interemitter separations comparable to the plasmon wavelength, which exhibits smaller amplitude and the relevant oscillating phase has an $\pi \backslash 2$ shift with respect to the cross decay rate. Deviations from the oscillatory behavior occur only when the interemitter separation is well below the vacuum wavelength, indicating that in this regime the QEs are already strongly interacting with each other through the vacuum.

Fig. 5. (a) Normalized cross decay rate $\Gamma _{12}$ and dipole-dipole shift $\delta \omega _{12}$ (scaled by the spontaneous decay rate $\Gamma _{11}$) versus the separation of two QEs ($\mathbf {r}_{12}=\left |\mathbf {r}_{2}-\mathbf {r}_{1}\right |$), the QEs are placed along the propagating direction of the SPF with an angle $\theta =30^{\circ }$ deviates from the $x$ axis and both have a distance $l=\lambda _{a}\backslash 500$ above from the metasurface. The diagonal conductivity elements are $\sigma _{xx}=10^{-4}+3i\:\textrm {mS}$ and $\sigma _{yy}=10^{-4}-1i\:\textrm {mS}$. (b) Normalized cross decay rate versus the separation of two QEs and the conductivity components $\sigma _{xx}=10^{-3}+3\sigma _{o}i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-\sigma _{o}i\:\textrm {mS}$ ($\sigma _{o}\in \left [0.1,3\right ]$), the propagating angle of the plasmon field keeps unchanged when the conductivity of the metasurface varies.

Download Full Size | PDF

Apart from the interemitter separation, we also study the influence of the material response on the cross decay behaviors with a particular set of parameters (see Fig. 5(b)). It should be stressed that the propagating angle of the SPF keeps unchanged ($\theta _{pl}=30^{\circ }$) under these variations. Obviously when both the capacitive and inductive responses of the metasurface are weak (the regime where $\sigma _{o}$ is smaller than $0.5$ in the contourplot), it could barely support long-range interactions between the QEs. The relevant physical scheme can be well understood by combining the results shown in Fig. 10, where only the propagation of the plasmon modes with large wave numbers are permitted. In other words, under this circumstance the platform is poor in supporting the SPF owing to the decrease in plasmon wavelength and large losses in the propagation directions. When the material response is extremely weak, the model is similar with the interaction of two QEs in free space. As illustrated by the figure, prominent cross interactions take place only for interemitter separations much smaller than the vacuum wavelength. Different from the former case where the QEs are in fact weakly coupled to the metasurface, interesting features arise when the optical response of the surface becomes dominant. It is clear that with the increase of the material response, the metasurface gradually plays as an important role in modifying the spontaneous decay behavior of the QEs. As the result, strong interactions between the QEs can be achieved for interemitter separations up to several vacuum wavelength. The reason lies in the existence of the propagating surface plasmon modes with large wavelength, which can be extracted from the results shown in Fig. 10. Meanwhile, low propagation losses of the SPF further enhance the interemitter interaction distance supported by the metasurface.

The discussions above indicate that an adequate tailoring of the metasurface conductivity tensor allows a simple and powerful control of the cross interaction between two QEs, which is mainly determined by the propagating wavelength of the SPF supported by the structure. Base on this property, in Fig. 6 we inspect the superradiant and subradiant behaviors of the QEs mutually coupled to the metasurface. One can see that when the material response is weak (as shown in subplot a), a rapid oscillatory variations between the superradiance and subradiance of the symmetric and antisymmetric transitions can be observed for different interemitter separations. This is due to the small propagating wavelength of the excited SPF illustrated by Figs. 5(b) and 10, where the cross decay rate of the emitters declines to the value around a half of the single-QE decay rate when they are apart from each other with a separation $5\lambda _{a}$. However, when compared with the cases that QEs are interacting through the isotropic surface (dotted line) and the vacuum modes in free space (dashed dotted line), the superradiant phenomenon for QEs mutually coupled to the HMS (solid line) has been enhanced and can be clearly identified at large interemitter separations. On the other hand, when the metasurface exhibits strong optical response (as subplot b demonstrates), a slow oscillatory behavior between the superradiance and subradiance with small amplitude damping rate along the interemitter separation axis can be observed. According to the previous analyses, it is known that the metasurface is a good platform in supporting distant interemitter interactions under this circumstance. Consequently, the QEs can interfere constructively in the symmetric state for the separation up to $10\lambda _{a}$, which is the main reason for the generation of the long-range superradiance.

Fig. 6. Super- and subradiant property of the system. The blue solid curves illustrate the $R$ factor as a function of the separation between two QEs that mutually coupled to the HMS (the propagating angle of the SPF is $\theta _{pl}=30^{\circ }$), where the dotted and dashed-dotted lines represent the cases of two QEs interacting through isotropic surface and free space, respectively. Both the QEs have a distance $l=\lambda _{a}\backslash 500$ above from the surface, with different material parameters: (a) $\sigma _{xx}=10^{-4}+1.5i\:\textrm {mS}$ and $\sigma _{yy}=10^{-4}-0.5i\:\textrm {mS}$ for the HMS and $\sigma =10^{-4}+1.5i\:\textrm {mS}$ for the isotropic surface. (b) $\sigma _{xx}=10^{-4}+9i\:\textrm {mS}$ and $\sigma _{yy}=10^{-4}-3i\:\textrm {mS}$ for the HMS and $\sigma =10^{-4}+9i\:\textrm {mS}$ for the isotropic surface.

Download Full Size | PDF

4. Quantum entanglement at macroscopic distance

Base on the superradiant features of two QEs mediated by the metasurface, in this section we will mainly focus on exploring the quantum entanglement properties of the system. In order to determine the amount of entanglement and the relevant dynamics between the QEs, the widely accepted measure of the entanglement, i.e. the parameter concurrence is adopted here [68]. The values of the concurrence function varies from $0$ to $1$, where the values $0$ and $1$ indicate non-entangled and maximally-entangled features of the system, respectively. We first consider the creation of quantum entanglement between the QEs via their mutual couplings to the SPF supported by the metasurface, with the system initially prepared in the single-excited state $\left |\psi \left (0\right )\right \rangle =\left |e_{1},g_{2}\right \rangle \otimes \left |\left \{ 0\right \} \right \rangle$ (the state $\left |\left \{ 0\right \} \right \rangle$ denotes the state of the vacuum modes with zero excitations). Under this circumstance, the time evolution of the concurrence is given by

(16)$$C\left(t\right) = \sqrt{\left[\rho_{ss}\left(t\right)-\rho_{aa}\left(t\right)\right]^{2}+4\textrm{Im}\left[\rho_{sa}\left(t\right)\right]^{2}}.$$

Here $\rho _{ss\left (aa\right )}\left (t\right )$ denotes the time evolution of the symmetric (antisymmetric) state population, and $\rho _{sa}\left (t\right )$ represents the evolution of one-photon coherence. The time dependence of these density matrix elements can be directly derived through Eq. (12), where their explicit forms are given in Eq. (28). Notice that apart from the population difference between the symmetric and antisymmetric states, the time evolution of the concurrence is also determined by the imaginary part of the coherence between these two states. Since this process is dominated by the off-resonance coupling, it will lead to oscillations in the time-dependent function with frequency difference $2\delta \omega _{12}$ between the states. However, previous analyses reveal that this energy difference vanishes when the QEs are maximally interacted with each other through their mutual couplings to the SPF. Meanwhile, maximums of the concurrence function can be achieved when the QEs system is mostly populated in the symmetric or antisymmetric state, which in turn requires a strong superradiant or subradiant behavior. For these reasons, this frequency shift will be neglected in the following discussions.

The features described above can be easily seen in Fig. 7, where both the concurrence and the evolution of the density matrix elements have been investigated according to Eqs. (16) and (28). As it has been mentioned before, the concurrence actually reflects the entanglement between the QEs. It is clear through Fig. 7(a) that for two distant QEs [55,56] (separated by a distance $10\lambda _{a}$, for HMSs that support the excitations of the SPF range from terahertz to mid-infrared [41,69], this separation can be up to milimeters) initially prepared in the separable state, the concurrence builds up immediately after $\tau =0$. With the increase of time, the concurrence reaches its maximum at some specific time points and then begins to decay. It should be notice that when compared with the isotropic surface (dotted line) and the free space (dashed-dotted line) cases, the concurrence has a prominent enhancement for the QEs interacting with the HMS (solid line). This phenomenon indicates that the excitation and propagation of the long-range SPF can build a strong interaction between two distant QEs via the spontaneous emissions, which finally leads to the creation of the entanglement. The behavior of the time evolution of the concurrence can be well understood by combining the dynamics of the density matrix elements. As illustrated by Fig. 7(b), the symmetric state suffers a rapid depopulation at early times (superradiance), where the population of the antisymmetric state decays much slower (subradiance). Thus the population difference between these two states increases with the time, which is the main reason for the growth in the entanglement. At later times where the symmetric state becomes depopulated, the concurrence decays slowly in time with the subradiant decay rate and overlaps with the population of the antisymmetric state. Owing to the strong interemitter interaction ($\Gamma _{12}=0.95\Gamma _{11}$) supported by the HMS, the lifetime of the entanglement can be effectively prolonged. From the above analyses, it is clear that the maximum concurrence of the spontaneously generated entanglement is always less than the upper bound with $C=0.5$. This is due to the incoherent interaction between the QEs, where the coherence of the system suffers a rapid decrease with the time. As the result, the concurrence is bounded by the population of the dark-like state (the subradiant state). Recently, investigations show that the coherent interaction between QEs can be realized in the non-reciprocal environments, which would lead to the increase in the purity of the systematic states [70]. Consequently, the maximum achievable concurrence could be larger than the reciprocal case and thus it is an efficient method to enhance the entanglement between two QEs.

Fig. 7. Quantum entanglement between distant QEs mediated by the HMS. (a) Time evolution (with the normalized time $\tau =t\Gamma _{11}$) of the concurrence for the QEs initially prepared in the separable state $\left |e_{1},g_{2}\right \rangle$, with a macroscopic separation $r_{12}=\left |\mathbf {r}_{2}-\mathbf {r}_{1}\right |=10\lambda _{a}$ and both have a distance $\lambda _{a}\backslash 500$ above from the surfaces. The concurrence mediated by the HMS (with conductivity elements $\sigma _{xx}=10^{-4}+9i\:\textrm {mS}$ and $\sigma _{yy}=10^{-4}-3i\:\textrm {mS}$) and the isotropic surface (with conductivity $\sigma =10^{-4}+9i\:\textrm {mS}$) are denoted by the blue solid and red dashed-dotted curves, respectively. Also the concurrence for the QEs interacting through the free vacuum modes (red dashed-dotted line) is given for comparison. (b) Time evolution of the populations and the photon coherence for the case that has been illustrated by the blue solid curve in the subplot (a). (c) Maximum entanglement versus the material dissipations and optical responses of the HMS. The interemitter distances where the QEs interact maximally with each other near $r_{12}=10\lambda _{a}$ are considered, and the varied conductivity elements are $\sigma _{xx}=\sigma _{d}+3\sigma _{o}i\:\textrm {mS}$ and $\sigma _{yy}=\sigma _{d}-\sigma _{o}i\:\textrm {mS}$ $\left (\sigma _{d}\in \left [10^{-4},10^{-2}\right ],\sigma _{o}\in \left [0.1,6\right ]\right )$.

Download Full Size | PDF

To develop a brief overview of the quantum entanglement between distant QEs modulated by the metasurface, a contourplot illustrates the dependence of the concurrence on both the material dissipation (real part of the conductivity) and its optical response (imaginary part of the conductivity) has been given in Fig. 7(c). Here we consider the cases that the QEs are maximally interacted with each other for any set of the material parameters, where the specific separation can be determined by finding the relevant peak in the oscillating curve of the factor $R$ near $r_{12}=10\lambda _{a}$. It is clear that the concurrence is negligible when the material response is weak (i.e. both the imaginary parts of $\sigma _{xx}$ and $\sigma _{yy}$ are small), which indicates an absence of the entanglement between distant QEs. On the other hand, for strong material responses, prominent increases in the concurrence can be observed and the generated quantum entanglement approaches the degree of $45\%$ with the corresponding pure state of the system. The physics associated with these phenomena are easily understood by combining the results in Figs. 6 and 10, where the surface exhibits strong optical response behaves as an excellent platform in supporting long-range interemitter interactions through their mutual couplings to the SPF. Also we can see that for HMSs with strong optical responses, the entanglement becomes more robust to the material dissipations, which mainly caused by the decline in the propagation losses of the SPF.

In real experimental systems, quantum entanglement degrades with time due to the inevitable coupling to the reservoir, where various kinds of entanglement dynamics can be observed. In some scenarios, the concurrence vanishes at a finite time rather than decays asymptotically, which is known as the sudden death of entanglement. Commonly, in the presence of strong interactions between two QEs which are closely spaced, quantum entanglement can be revived from death after a finite time interval. In what follows we focus on these special dynamics of the entanglement at macroscopic distances. Here consider the system initially prepared in an entangled state $\left |\varphi \left (0\right )\right \rangle =\left (\sqrt {\eta }\left |e\right \rangle +\sqrt {1-\eta }\left |g\right \rangle \right )\otimes \left |\left \{ 0\right \} \right \rangle$, where $\eta$ is positive and satisfies $0\leq \eta \leq 1$. The state corresponds to an excitation of the system into a linear superposition of the states in which both or neither of the atoms is excited, and it reduces to the maximally entangled Bell state in the special case of $\eta =1\backslash 2$. Under this circumstance, the concurrence is determined by $C\left (t\right )=\textrm {max}\left \{ 0,C_{coh}\left (t\right ),C_{dis}\left (t\right )\right \}$, where the coherence and the population distribution dominated criterions are

(17)$$C_{coh}\left(t\right) = 2\left|\rho_{ge}\left(t\right)\right|-\left[\rho_{ss}\left(t\right)+\rho_{aa}\left(t\right)\right]$$

and

(18)$$C_{dis}\left(t\right) = \left|\rho_{ss}\left(t\right)-\rho_{aa}\left(t\right)\right|-2\sqrt{\rho_{ee}\left(t\right)\rho_{gg}\left(t\right)}.$$

It should be pointed out that at $t=0$, only the populations $\rho _{ee}\left (0\right )=\eta$, $\rho _{gg}\left (0\right )=1-\eta$ and two-photon coherence $\rho _{ge}\left (0\right )=\sqrt {\eta \left (1-\eta \right )}$ are nonvanished. Thus the coherence contributes to the concurrence since $C_{coh}\left (0\right )>0$ and $C_{dis}\left (0\right )<0$. Then according to Eq. (10), populations in the symmetric state ($\rho _{ss}$, which is superradiant for large $\Gamma _{12}$) and antisymmetric state ($\rho _{aa}$, which is subradiant for large $\Gamma _{12}$) will be spontaneously accumulated during the evolution, which may lead to nonvanished concurrence governed by $C_{dis}\left (t\right )$. This kind of the dynamical evolution is often related to spontaneously generated entanglement that originates from the spontaneous decay of the system.

Figure 8 shows the dynamical evolutions of the concurrence for two QEs separated by macroscopic distance and initially prepared in the state $\left |\varphi \left (0\right )\right \rangle$. Depending on the cross decay rate (which determines the superradiant and subradiant properties of the system) and the coefficient $\eta$, different dynamics can be clearly identified. As illustrated by the subplot (a), the condition for the appearance of the sudden death of the entanglement is $\eta >1/2$ when the QEs are independent from each other. It is easily verified that in this case the symmetric and the antisymmetric states are equally populated ($\rho _{ss}=\rho _{aa}$ for all times) owing to the vanished cross interaction, thus the criterion $C_{dis}\left (t\right )$ keeps negative and only the coherence dominated term contributes to the concurrence. The relevant death time for the entanglement can be derived via the condition $C_{coh}\left (t\right )=0$ and turns out to be $T_{d}=\Gamma _{11}^{-1}\ln \left [\left (\eta +\sqrt {\eta \left (1-\eta \right )}\right )\backslash \left (2\eta -1\right )\right ]$, which has solutions for $\eta >1\backslash 2$ and $\eta \neq 1$. The situation changes when the QEs become interactive, that is, with nonvanished cross decay rate. As shown in Fig. 8(b), the revival of entanglement appears for moderate cross interaction strength. Moreover, two revivals of the entanglement after finite death times occur when the QEs are strongly interacting with each other and the initial state of the system is largely inverted (region $\textrm {II}$ in Fig. 8(c), marked by white dashed lines with $0.86<\eta <0.9$). The underlying physics is much more complicated than that in the noninteracting case and we will discuss it later. Finally, to gain an insight into the dynamics of the entanglement, the time evolution of the concurrence varies with the interemitter separation has been plotted for an initially inverted system. As expected, different dynamics can be distinguished. It is clear that for separations where moderate interactions between the QEs take place, a finite interval of complete disentanglement is followed by an entanglement revival at large times (the revival feature is highlighted by yellow dashed line). However, when the interaction between the QEs becomes prominent ($\Gamma _{12}\approx \Gamma _{11}$), a new interval exhibits the feature of the entanglement revival emerges (highlighted by red dashed line) and thus there are totally two revivals in the entanglement for such cases.

Fig. 8. Time evolution of the concurrence for two QEs separated by macroscopic distances and initially prepared in the state $\left |\varphi \left (0\right )\right \rangle$, other parameters are $l=\lambda _{a}/500$, $\sigma _{xx}=10^{-4}+9i\:\textrm {mS}$ and $\sigma _{yy}=10^{-4}-3i\:\textrm {mS}$. In subplots (a), (b) and (c) the dependence of the concurrence on $\eta$ has also been investigated, where different cross interaction strengths are considered: (a) Noninteracting QEs system with $r_{12}=9.75\lambda _{a}$ and $\Gamma _{12}\approx 0$. (b) Moderate interacting QEs system with $r_{12}=9.83\lambda _{a}$ and $\Gamma _{12}\approx 0.5\Gamma _{11}$. (c) Strong interacting QEs system with $r_{12}=10\lambda _{a}$ and $\Gamma _{12}\approx \Gamma _{11}$. The gray regions illustrate an absence of the entanglement between the QEs, where the white dashed lines are introduced to distinguish regions with different concurrence dynamics. Regions denoted by $\textrm {I}$ and $\textrm {II}$ demonstrate the dynamics that the entanglement revives once and twice, respectively. Subplot (d) displays the time evolution of the concurrence as a function of the interemitter separation for $\eta =0.89$, where regions correspond to type $\textrm {I}$ (yellow dashed line) and the first revival in type $\textrm {II}$ (red dashed line) are highlighted.

Download Full Size | PDF

The origin of the sudden death and revivals of the entanglement can be well understood by studying the evolutions of the collective state populations and photon coherences. First focus on the moderate coupling case with the interemitter separation $r_{12}=9.83\lambda _{a}$ and the cross decay rate $\Gamma _{12}\approx 0.5\Gamma _{11}$, which relates to the results that have been discussed in Fig. 8(b). As it can be seen from Fig. 9(a), the reason for the revival of the entanglement lies in the unequal decays of the populations $\rho _{aa}$ and $\rho _{ee}$. At early times the two-photon coherence $\rho _{ge}$ is larger than other density matrix elements such that the criterion $C_{coh}\left (t\right )>0$, which leads to the presence of the entanglement. Then the entanglement suffers a sharp decrease with the growth in the collective state populations and vanishes at the time where $\rho _{ss}+\rho _{aa}>2\left |\rho _{ge}\right |$ (cross the intersection point of the blue and red solid lines). Apparently, sudden death of the entanglement is caused by the significant population accumulation in both the symmetric and antisymmetric states. At long times, most of the density matrix elements that include excitations decay to zero except the population of the antisymmetric state, which remains large due to its subradiant feature. Therefore when the criterion $C_{dis}\left (t\right )$ becomes positive (with the critical time point $T_{dis}=2\backslash \left (3\Gamma _{12}\right )\ln \left [4\Gamma _{11}\backslash \left (\sqrt {\eta }\left (\Gamma _{11}-\Gamma _{12}\right )\right )\right ]$ marked by the yellow dot) the entanglement revives from death.

Fig. 9. The transient behavior of the coherence- and population-dependent terms in both the criterions $C_{coh}\left (t\right )$ and $C_{dis}\left (t\right )$, where different cases are considered: (a) $r_{12}=9.83\lambda _{a}$, $\Gamma _{12}\approx 0.5\Gamma _{11}$ and $\eta =0.89$. (b) $r_{12}=10\lambda _{a}$, $\Gamma _{12}\approx \Gamma _{11}$ and $\eta =0.89$. For each panel, the inset provides the evolution of the concurrence under the same conditions.

Download Full Size | PDF

As it has been shown in Fig. 8(c), the entanglement could revive twice when the interaction between the QEs is strong. Under this circumstance, by considering a largely inverted initial state of the strongly interacting system (with $\eta =0.89$, $r_{12}=10\lambda _{a}$ and $\Gamma _{12}\approx \Gamma _{11}$), the dynamical behavior of the density matrix elements is displayed in Fig. 9(b). It is shown that at early times the death of the entanglement is mainly due to the population accumulation in the symmetric state, where the first revival of the entanglement is caused by the fast decay of the same state. Within this death-revival period, the coherence $\rho _{ge}$ plays a key role in the production of the entanglement. Thus one can treat the first revival as a re-built of the initial entanglement, which has been unmasked by destroying the symmetric state population. The origin of the appearance of the second revival, which leads to the long time entanglement, is equivalent to the case that has been discussed in Fig. 9(a). On such a time scale, the system is mostly populated in the antisymmetric state and the entanglement decays asymptotically with the subradiant decay rate.

5. Conclusions

In conclusion, we have investigated the dipole radiation properties of one QE, also the collective behavior and quantum entanglement dynamics of two distant QEs coupled to the HMS. It has been shown that by adjusting the distance of the QE to the HMS, the radiation energy of the QE can well transfer into the surface plasmon modes, which indicates a strong coupling between the QE and the HMS. In this case, the excited SPF channels as narrow rays propagating along specific directions, where the control over the propagating SPF can be achieved by engineering the anisotropy of the HMS. Further studies demonstrate that our structure allows strong cross interactions between distant QEs mutually coupled to the HMS. As the result, prominent superradiance and subradiance occur even for macroscopic interemitter separations. Finally we have shown that remarkable quantum entanglement between QEs with macroscopic separations can be realized owing to the existence of spatially scalable interaction channels supported by the HMS. For QEs systems initially prepared in the superpositions of the ground state and the excited state, strong revivals from sudden death can also be observed.

As we have pointed out, HMSs often possess high tunability in their optical properties [48,62]. Meanwhile, the supported plasmon modes have been predicted to suffer much low energy loss than that in TDHMs, which originates from the unique surface-guided property offered by the HMSs [42,60]. It is known that for quantum computing and quantum information processing that aim at on-chip architecture [71], the efficient channeling of the signal from the source qubit toward the desired target qubit is especially crucial for practical applications. However, even in the plasmonic waveguides remarkable entanglement can only be observed for the interemitter separations within several wavelengths [18,19]. Different from normal plasmonic structures, our results show that notable quantum entanglement between the QEs could be generated for a macroscopic separation up to $10\lambda _{a}$ (reaches the milimeter magnitude for terahertz transitions) or even larger by coupling the QEs to the HMSs. This intriguing property allows the preparation of the QEs with macroscopic separations, which dramatically reduces the difficulty in fabricating the hybrid systems consist of two dimensional materials and semiconductor quantum dots [72]. Thus compared with traditional methods, our proposal provides a more promising route in achieving on-chip platforms for two- or many-body quantum entanglement, which may have potential applications in quantum computing and quantum information processing.

Appendix A: Surface plasmons supported by the HMSs

In order to understand the behavior of the surface modes, it is instructive to inspect the relevant dispersion relation, which can be acquired by investigating the zero(s) of the denominator of the reflection coefficient given in Eq. (4)

(19)$$2\sigma_{xx}\left(k_{0}^{2}-k_{x}^{2}\right)+2\sigma_{yy}\left(k_{0}^{2}-k_{y}^{2}\right)-k_{0}k_{z}\eta^{-1}\left(4+\eta^{2}\sigma_{xx}\sigma_{yy}\right) = 0.$$

Then starting from Eq. (19) and after arrangements one can arrive at the result shown in Eq. (5) by assuming a low loss surface with purely imaginary conductivity. Rigorously, the root of the Eq. (19) for a given wave number $k_{y}$ reveals to be

(20)$$k_{x}^{sp} = \pm\sqrt{\alpha\pm2i\beta\left(\beta^{2}+\alpha+1\right)^{\frac{1}{2}}},$$

where the parameters are

(21)$$\alpha = k_{0}^{2}-2\beta^{2}-\left(k_{y}^{2}-k_{0}^{2}\right)\sigma_{xx}^{-1}\sigma_{yy},$$

(22)$$\beta = k_{0}\eta^{-1}\left(\sigma_{xx}^{-1}+\eta^{2}\sigma_{yy}\backslash 4\right).$$

It should be noted that the signum outside the square root in Eq. (20) indicates the forward or backward propagation of the surface plasmon waves along the $x-$axis, whereas the signum inside the square root denotes the possible plasmonic modes supported by the HMSs. The proper choice of the plasmonic modes would lead to a decaying surface wave travelling away from the dipole, which is determined by the imaginary part of the wave number $k_{x}^{sp}$. On the other hand, Eq. (20) interestingly predicts the propagation angle of the surface waves for large wave numbers ($k_{y}\gg k_{0}$), where the wave numbers of the plasmon modes satisfy the relation

(23)$$\frac{k_{x}^{sp}}{k_{y}} = \pm\sqrt{-\frac{\sigma_{yy}}{\sigma_{xx}}}.$$

In Fig. 10, both the variations in the real and imaginary parts of the supported plasmonic wave number $k_{x}^{sp}$ are studied for specific sets of parameters, where the propagation angle of the SPF is $\theta _{pl}=30^{\circ }$. It is clear that when the optical response of the HMS is sufficiently weak (for $\sigma _{o}<0.2$), only the SPF with large wave numbers can be excited. Meanwhile, the SPF suffers considerable propagation losses owing to the large values of $\textrm {Im}\left [k_{x}^{sp}\right ]$. However, the circumstance is different when the material response becomes prominent. In this case the HMSs support the propagation of the SPF with small wave numbers and low propagation losses, which can be clearly identified for $\sigma _{o}>1.5$. Thus in general, one can draw conclusions that the HMSs with strong optical response are advantaged in supporting the long-range propagation of the SPF.

Fig. 10. Normalized surface plasmon wave number $k_{x}^{sp}$ versus the in-plane wave number $k_{y}$ and the conductivity components of the HMS $\sigma _{xx}=10^{-3}+3\sigma _{o}i\:\textrm {mS}$ and $\sigma _{yy}=10^{-3}-\sigma _{o}i\:\textrm {mS}$, where the absolute values of different parts are plotted: (a) The real part of the wave number $k_{x}^{sp}$, which relates to the propagating wavelength of the SPF. (b) The imaginary part of the $k_{x}^{sp}$ that relate to the propagation loss of the SPF.

Download Full Size | PDF

Then continue with Eq. (20), by combining the results given in Eqs. (1)$-$(4), in the plasmon region (the dipole decays only through the surface plasmon modes supported by the HMSs) the inner integrals appeared in the Green function can be approximately evaluated as

(24)$$f_{1}^{sp}\left(k_{y}\right) = R_{zz}^{sp}\left(k_{x}^{sp},k_{y}\right)\left(2k_{z}^{sp}\right)^{-1}e^{ik_{z}^{sp}\left(z+z_{1}\right)}e^{-ik_{x}^{sp}\left(x-x_{1}\right)},$$

(25)$$f_{2}^{sp}\left(k_{y}\right) = R_{zz}^{sp}\left(k_{x}^{sp},k_{y}\right)\left(2k_{z}^{sp}\right)^{-1}\left(k_{x}^{sp}\right)^{2}e^{ik_{z}^{sp}\left(z+z_{1}\right)}e^{-ik_{x}^{sp}\left(x-x_{1}\right)}.$$

Here $k_{z}^{sp}=\sqrt {k_{0}^{2}-\left (k_{x}^{sp}\right )^{2}-k_{y}^{2}}$ represents the $z-$component wave number, notice that owing to the dominant contribution of the SPF to the spontaneous decay the branch cut terms of the integrals have been ignored in the above equations. In this scenario, the reflection coefficient denoted by $R_{zz}^{sp}\left (k_{x}^{sp},k_{y}\right )$ can be redefined just by replacing the denominator of the reflection coefficient in Eq. (4) with its partial derivative over $k_{x}$ (and then replace $k_{x}$ by $k_{x}^{sp}$). Finally with the combination of the Eqs. (1), (24) and (25), one can obtain the $zz$ component of the Green tensor

(26)$$\overset{\leftrightarrow}{\mathbf{G}}_{zz}^{sp}\left(\mathbf{r},\mathbf{r}_{1},\omega_{a}\right) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-ik_{y}\left(y-y_{1}\right)}\left[f_{2}^{sp}\left(k_{y}\right)+k_{y}^{2}f_{1}^{sp}\left(k_{y}\right)\right]dk_{y}.$$

Appendix B: Time evolution of the reduced density matrix elements for the QEs

From the master equation given in Eq. (12), in Dicke basis the motion equations for the populations of the collective states (diagonal elements) and the photon coherences (off-diagonal elements) are found to be in the following forms

(27)$$\begin{aligned} &\dot{\rho}_{ee}\left(t\right) = -2\Gamma_{11}\rho_{ee}\left(t\right),\\ &\dot{\rho}_{ss}\left(t\right) = -\left(\Gamma_{11}+\Gamma_{12}\right)\left[\rho_{ss}\left(t\right)-\rho_{ee}\left(t\right)\right],\\ &\dot{\rho}_{aa}\left(t\right) = -\left(\Gamma_{11}-\Gamma_{12}\right)\left[\rho_{aa}\left(t\right)-\rho_{ee}\left(t\right)\right],\\ &\dot{\rho}_{gg}\left(t\right) = \left(\Gamma_{11}+\Gamma_{12}\right)\rho_{ss}\left(t\right)+\left(\Gamma_{11}-\Gamma_{12}\right)\rho_{aa}\left(t\right),\\ &\dot{\rho}_{as}\left(t\right) = -\left(\Gamma_{11}+2i\delta\omega_{12}\right)\rho_{as}\left(t\right),\\ &\dot{\rho}_{ge}\left(t\right) = -\Gamma_{11}\rho_{ge}\left(t\right).\end{aligned}$$

Obviously, the transitions including the symmetric state $\left |e\right \rangle \rightarrow \left |s\right \rangle \rightarrow \left |g\right \rangle$ take place with an enhanced decay rate $\Gamma _{11}+\Gamma _{12}$ (the relevant transition channels are superradiant), where the transitions relate to the antisymmetric state $\left |e\right \rangle \rightarrow \left |a\right \rangle \rightarrow \left |g\right \rangle$ occur with a reduced decay rate $\Gamma _{11}-\Gamma _{12}$ (the relevant transition channels are subradiant). Solving the above equations by substituting the initial conditions of the state populations and photon coherences, we can acquire the time-dependent solutions of the density matrix elements for arbitrary initial conditions

(28)$$\begin{aligned} &\rho_{ee}\left(t\right) = \rho_{ee}\left(0\right)e^{-2\Gamma_{11}t},\\ &\rho_{ss}\left(t\right) = \rho_{ss}\left(0\right)e^{-\left(\Gamma_{11}+\Gamma_{12}\right)t}+\rho_{ee}\left(0\right)\frac{\Gamma_{11}+\Gamma_{12}}{\Gamma_{11}-\Gamma_{12}}\left[e^{-\left(\Gamma_{11}+\Gamma_{12}\right)t}-e^{-2\Gamma_{11}t}\right],\\ &\rho_{aa}\left(t\right) = \rho_{aa}\left(0\right)e^{-\left(\Gamma_{11}-\Gamma_{12}\right)t}+\rho_{ee}\left(0\right)\frac{\Gamma_{11}-\Gamma_{12}}{\Gamma_{11}+\Gamma_{12}}\left[e^{-\left(\Gamma_{11}-\Gamma_{12}\right)t}-e^{-2\Gamma_{11}t}\right],\\ &\rho_{sa}\left(t\right) = \rho_{sa}\left(0\right)e^{-\left(\Gamma_{11}+2i\delta\omega_{12}\right)t},\\ &\rho_{ge}\left(t\right) = \rho_{ge}\left(0\right)e^{-\Gamma_{11}t}.\end{aligned}$$

It is clear that the dynamical evolution of the QEs system strongly depends on the initial values of each quantum states.

Funding

National Natural Science Foundation of China (11874287, 11947044); Huaqiao University (605-50Y19046); Science and Technology Commission of Shanghai Municipality (18JC1410900).

Disclosures

The authors declare no conflicts of interest.

References

1. W. E. Lamb and R. C. Retherford, “Fine structure of the hydrogen atom by a microwave method,” Phys. Rev. 72(3), 241–243 (1947). [CrossRef]

2. R. H. Dicke, “Coherence in spontaneous emission process,” Phys. Rev. 93(1), 99–110 (1954). [CrossRef]

3. A. F. van Loo, A. Fedorov, K. Lalumière, B. C. Sanders, A. Blais, and A. Wallraff, “Photon-mediated interactions between distant artificial atoms,” Science 342(6165), 1494–1496 (2013). [CrossRef]

4. A. Goban, C. L. Hung, J. D. Hood, S. P. Yu, J. A. Muniz, O. Painter, and H. J. Kimble, “Superradiance for atoms trapped along a photonic crystal waveguide,” Phys. Rev. Lett. 115(6), 063601 (2015). [CrossRef]

5. Z. Meir, O. Schwartz, E. Shahmoon, D. Dron, and R. Ozeri, “Cooperative Lamb shift in a mesoscopic atomic array,” Phys. Rev. Lett. 113(19), 193002 (2014). [CrossRef]

6. M. O. Araújo, I. Krešić, R. Kaiser, and W. Guerin, “Superradiance in a large and dilute cloud of cold atoms in the linear-optics regime,” Phys. Rev. Lett. 117(7), 073002 (2016). [CrossRef]

7. P. G. Brooke, “Decoherence-free quantum-information processing using dipole-coupled qubits,” Phys. Rev. A 75(2), 022320 (2007). [CrossRef]

8. A. González-Tudela, V. Paulisch, D. E. Chang, H. J. Kimble, and J. I. Cirac, “Deterministic generation of arbitrary photonic states assisted by dissipation,” Phys. Rev. Lett. 115(16), 163603 (2015). [CrossRef]

9. M. F. Riedel, P. Böhi, Y. Li, T. W. Hänsch, A. Sinatra, and P. Treutlein, “Atom-chip-based generation of entanglement for quantum metrology,” Nature 464(7292), 1170–1173 (2010). [CrossRef]

10. T. Wilk, A. Gaëtan, C. Evellin, J. Wolters, Y. Miroshnychenko, P. Grangier, and A. Browaeys, “Entanglement of two individual neutral atoms using Rydberg blockade,” Phys. Rev. Lett. 104(1), 010502 (2010). [CrossRef]

11. M. Saffman, T. G. Walker, and K. Mølmer, “Quantum information with Rydberg atoms,” Rev. Mod. Phys. 82(3), 2313–2363 (2010). [CrossRef]

12. G. Y. Chen, Y. N. Chen, and D. S. Chun, “Spontaneous emission of quantum dot excitons into surface plasmons in a nanowire,” Opt. Lett. 33(19), 2212–2214 (2008). [CrossRef]

13. A. Delga, J. Feist, J. Bravo-Abad, and F. J. Garcia-Vidal, “Quantum emitters near a metal nanoparticle: strong coupling and quenching,” Phys. Rev. Lett. 112(25), 253601 (2014). [CrossRef]

14. F. H. L. Koppens, D. E. Chang, and F. J. G. García de Abajo, “Graphene plasmonics: a platform for strong light-matter interactions,” Nano Lett. 11(8), 3370–3377 (2011). [CrossRef]

15. D. E. Gómez, K. C. Vernon, P. Mulvaney, and T. J. Davis, “Surface plasmon mediated strong exciton-photon coupling in semiconductor nanocrystals,” Nano Lett. 10(1), 274–278 (2010). [CrossRef]

16. W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon subwavelength optics,” Nature 424(6950), 824–830 (2003). [CrossRef]

17. Y. Luo, A. I. Fernandez-Dominguez, A. Wiener, S. A. Maier, and J. B. Pendry, “Surface plasmons and nonlocality: a simple model,” Phys. Rev. Lett. 111(9), 093901 (2013). [CrossRef]

18. A. Gonzalez-Tudela, D. Martin-Cano, E. Moreno, L. Martin-Moreno, C. Tejedor, and F. J. Garcia-Vidal, “Entanglement of two qubits mediated by one-dimensional plasmonic waveguides,” Phys. Rev. Lett. 106(2), 020501 (2011). [CrossRef]

19. M. Hensen, T. Heilpern, S. K. Gray, and W. Pfeiffer, “Strong coupling and entanglement of quantum emitters embedded in a nanoantenna-enhanced plasmonic cavity,” ACS Photonics 5(1), 240–248 (2018). [CrossRef]

20. R. C. Ge and S. Hughes, “Quantum dynamics of two quantum dots coupled through localized plasmons: an intuitive and accurate quantum optics approach using quasinormal modes,” Phys. Rev. B 92(20), 205420 (2015). [CrossRef]

21. G. Y. Chen, N. Lambert, C. H. Chou, Y. N. Chen, and F. Nori, “Surface plasmons in a metal nanowire coupled to colloidal quantum dots: scattering properties and quantum entanglement,” Phys. Rev. B 84(4), 045310 (2011). [CrossRef]

22. M. Otten, R. A. Shah, N. F. Scherer, M. S. Min, M. Pelton, and S. K. Gray, “Entanglement of two, three, or four plasmonically coupled quantum dots,” Phys. Rev. B 92(12), 125432 (2015). [CrossRef]

23. H. N. S. Krishnamoorthy, Z. Jacob, E. Narimanov, I. Kretzschmar, and V. M. Menon, “Topological transitions in metamaterials,” Science 336(6078), 205–209 (2012). [CrossRef]

24. A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, “Hyperbolic metamaterials,” Nat. Photonics 7(12), 948–957 (2013). [CrossRef]

25. M. A. K. Othman, C. Guclu, and F. Capolino, “Graphene-based tunable hyperbolic metamaterials and enhanced near-field absorption,” Opt. Express 21(6), 7614–7632 (2013). [CrossRef]

26. L. Ferrari, C. H. Wu, D. Lepage, X. Zhang, and Z. W. Liu, “Hyperbolic metamaterials and their applications,” Prog. Quantum Electron. 40, 1–40 (2015). [CrossRef]

27. S. M. Prokes, O. J. Glembocki, J. E. Livenere, T. U. Tumkur, J. K. Kitur, G. Zhu, B. Wells, V. A. Podolskiy, and M. A. Noginov, “Hyperbolic and plasmonic properties of Silicon/Ag aligned nanowire arrays,” Opt. Express 21(12), 14962–14974 (2013). [CrossRef]

28. X. D. Yang, J. Yao, J. Rho, X. B. Yin, and X. Zhang, “Experimental realization of three-dimensional indefinite cavities at the nanoscale with anomalous scaling laws,” Nat. Photonics 6(7), 450–454 (2012). [CrossRef]

29. J. D. Caldwell, A. V. Kretinin, Y. G. Chen, V. Giannini, M. M. Fogler, Y. Francescato, C. T. Ellis, J. G. Tischler, C. R. Woods, A. J. Giles, M. Hong, K. Watanabe, T. Taniguchi, S. A. Maier, and K. S. Novoselov, “Sub-diffractional volume-confined polaritons in the natural hyperbolic material hexagonal boron nitride,” Nat. Commun. 5(1), 5221 (2014). [CrossRef]

30. S. V. Zhukovsky, A. A. Orlov, V. E. Babicheva, A. V. Lavrinenko, and J. E. Sipe, “Photonic-band-gap engineering for volume plasmon polaritons in multiscale multilayer hyperbolic metamaterials,” Phys. Rev. A 90(1), 013801 (2014). [CrossRef]

31. V. P. Drachev, V. A. Podolskiy, and A. V. Kildshev, “Hyperbolic metamaterials: new physics behind a classical problem,” Opt. Express 21(12), 15048–15064 (2013). [CrossRef]

32. K. V. Sreekanth, A. D. Luca, and G. Strangi, “Negative refraction in graphene-based hyperbolic metamaterials,” Appl. Phys. Lett. 103(2), 023107 (2013). [CrossRef]

33. P. N. Li, M. Lewin, A. V. Kretinin, J. D. Caldwell, K. S. Novoselov, T. Taniguchi, K. Watanabe, F. Gaussmann, and T. Taubner, “Hyperbolic phonon-polaritons in boron nitride for near-field optical imaging and focusing,” Nat. Commun. 6(1), 7507 (2015). [CrossRef]

34. A. Fang, T. Koschny, and C. M. Soukoulis, “Optical anisotropic metamaterials: Negative refraction and focusing,” Phys. Rev. B 79(24), 245127 (2009). [CrossRef]

35. M. A. Noginov, H. Li, Y. A. Barnakov, D. Dryden, G. Nataraj, G. Zhu, C. E. Bonner, M. Mayy, Z. Jacob, and E. E. Narimanov, “Controlling spontaneous emission with metamaterials,” Opt. Lett. 35(11), 1863–1865 (2010). [CrossRef]

36. A. N. Poddubny, P. A. Belov, and Y. S. Kivshar, “Spontaneous radiation of a finite-size dipole emitter in hyperbolic media,” Phys. Rev. A 84(2), 023807 (2011). [CrossRef]

37. C. L. Cortes and Z. Jacob, “Super-Coulombic atom-atom interactions in hyperbolic media,” Nat. Commun. 8(1), 14144 (2017). [CrossRef]

38. I. Trushkov and I. Iorsh, “Two-dimensional hyperbolic medium for electrons and photons based on the array of tunnel-coupled graphene nanoribbons,” Phys. Rev. B 92(4), 045305 (2015). [CrossRef]

39. D. Correas-Serrano, J. S. Gomez-Diaz, M. Tymchenko, and A. Alù, “Nonlocal response of hyperbolic metasurfaces,” Opt. Express 23(23), 29434–29448 (2015). [CrossRef]

40. O. Y. Yermakov, A. I. Ovcharenko, M. Song, A. A. Bogdanov, I. V. Iorsh, and Y. S. Kivshar, “Hybrid waves localized at hyperbolic metasurfaces,” Phys. Rev. B 91(23), 235423 (2015). [CrossRef]

41. J. S. Gomez-Diaz and A. Alù, “Flatland optics with hyperbolic metasurfaces,” ACS Photonics 3(12), 2211–2224 (2016). [CrossRef]

42. W. L. Ma, P. Alonso-González, S. J. Li, A. Y. Nikitin, J. Yuan, J. Martín-Sánchez, J. Taboada-Gutiérrez, I. Amenabar, P. N. Li, S. Vélez, C. Tollan, Z. G. Dai, Y. P. Zhang, S. Sriram, K. Kalantar-Zadeh, S. T. Lee, R. Hillenbrand, and Q. L. Bao, “In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal,” Nature 562(7728), 557–562 (2018). [CrossRef]

43. O. Y. Yermakov, A. A. Hurshkainen, D. A. Dobrykh, P. V. Kapitanova, I. V. Iorsh, S. B. Glybovski, and A. A. Bogdanov, “Experimental observation of hybrid TE-TM polarized surface waves supported by a hyperbolic metasurface,” Phys. Rev. B 98(19), 195404 (2018). [CrossRef]

44. T. Low, R. Roldán, H. Wang, F. N. Xia, P. Avouris, L. M. Moreno, and F. Guinea, “Plasmons and screening in monolayer and multilayer black phosphorus,” Phys. Rev. Lett. 113(10), 106802 (2014). [CrossRef]

45. X. Huang, Y. Q. Cai, X. W. Feng, W. C. Tan, D. M. N. Hasan, L. Chen, N. Chen, L. Wang, L. Huang, and T. J. Duffin, “Black phosphorus carbide as a tunable anisotropic plasmonic metasurface,” ACS Photonics 5(8), 3116–3123 (2018). [CrossRef]

46. P. N. Li, I. Dolado, F. J. Alfaro-Mozaz, F. Casanova, L. E. Hueso, S. Liu, J. H. Edgar, A. Y. Nikitin, S. Vélez, and R. Hillenbrand, “Infrared hyperbolic metasurface based on nanostructured van der Waals materials,” Science 359(6378), 892–896 (2018). [CrossRef]

47. A. A. High, R. C. Devlin, A. Dibos, M. Polking, D. S. Wild, J. Perczel, N. P. de Leon, M. D. Lukin, and H. Park, “Visible-frequency hyperbolic metasurface,” Nature 522(7555), 192–196 (2015). [CrossRef]

48. A. Fallahi and J. Perruisseau-Carrier, “Design of tunable biperiodic graphene metasurfaces,” Phys. Rev. B 86(19), 195408 (2012). [CrossRef]

49. P. Y. Chen, J. Soric, Y. R. Padooru, H. M. Bernety, A. B. Yakovlev, and A. Alù, “Nanostructured graphene metasurface for tunable terahertz cloaking,” New J. Phys. 15(12), 123029 (2013). [CrossRef]

50. D. K. Efetov and P. Kim, “Controlling electron-phonon interactions in graphene at ultrahigh carrier densities,” Phys. Rev. Lett. 105(25), 256805 (2010). [CrossRef]

51. L. Ren, Q. Zhang, J. Yao, Z. Z. Sun, R. Kaneko, Z. Yan, S. Nanot, Z. Jin, I. Kawayama, and M. Tonouchi, “Terahertz and infrared spectroscopy of gated large-area graphene,” Nano Lett. 12(7), 3711–3715 (2012). [CrossRef]

52. M. S. Tomǎs, “Green function for multilayers: light scattering in planar cavities,” Phys. Rev. A 51(3), 2545–2559 (1995). [CrossRef]

53. J. S. Bagby and D. P. Nyquist, “Dyadic Green-functions for integrated electronic and optical circuits,” IEEE Trans. Microwave Theory Tech. 35(2), 207–210 (1987). [CrossRef]

54. M. Buscema, G. A. Steele, H. S. J. van der Zant, and A. Castellanos-Gomez, “The effect of the substrate on the Raman and photoluminescence emission of single-layer MoS2,” Nano Res. 7(4), 561–571 (2014). [CrossRef]

55. R. Xie, D. Battaglia, and X. Peng, “Colloidal InP nanocrystals as efficient emitters covering blue to near-infrared,” J. Am. Chem. Soc. 129(50), 15432–15433 (2007). [CrossRef]

56. Y. Zhang, K. Shibata, N. Nagai, C. Ndebeka-Bandou, G. Bastard, and K. Hirakawa, “Terahertz intersublevel transitions in single self-assembled InAs quantum dots with variable electron numbers,” Nano Lett. 15(2), 1166–1170 (2015). [CrossRef]

57. M. C. Chew, Waves and fields in inhomogeneous media (IEEE Press, 1995).

58. A. Eroglu, Wave propagation and radiation in gyrotropic and anisotropic media (Springer Press, 2010).

59. H. T. Dung, S. Y. Buhmann, L. Knöll, D. G. Welsch, S. Scheel, and J. Kästel, “Electromagnetic-field quantization and spontaneous decay in left-handed media,” Phys. Rev. A 68(4), 043816 (2003). [CrossRef]

60. A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. L. Koppens, “Highly confined low-loss plasmons in graphene-boron nitride heterostructures,” Nat. Mater. 14(4), 421–425 (2015). [CrossRef]

61. G. X. Ni, A. S. McLeod, Z. Sun, L. Wang, L. Xiong, K. W. Post, S. S. Sunku, B. Y. Jiang, J. Hone, C. R. Dean, M. M. Fogler, and D. N. Basov, “Fundamental limits to graphene plasmonics,” Nature 557(7706), 530–533 (2018). [CrossRef]

62. F. Wang, Y. B. Zhang, C. S. Tian, C. Girit, A. Zettl, M. Crommie, and Y. R. Shen, “Gate-variable optical transitions in graphene,” Science 320(5873), 206–209 (2008). [CrossRef]

63. S. Kühn, U. Håkanson, L. Rogobete, and V. Sandoghdar, “Enhancement of single-molecule fluorescence using a gold nanoparticle as an optical nanoantenna,” Phys. Rev. Lett. 97(1), 017402 (2006). [CrossRef]

64. B. Huttner and S. M. Barnett, “Quantization of the electromagnetic field in dielectrics,” Phys. Rev. A 46(7), 4306–4322 (1992). [CrossRef]

65. M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, 1997).

66. R. H. Lehmberg, “Radiation from an N-atom system. II. Spontaneous emission from a pair of atoms,” Phys. Rev. A 2(3), 889–896 (1970). [CrossRef]

67. S. Y. Buhmann, L. Knöll, D. G. Welsch, and H. T. Dung, “Casimir-Polder forces:A nonperturbative approach,” Phys. Rev. A 70(5), 052117 (2004). [CrossRef]

68. Z. Ficek and R. Tanaś., “Entanglement induced by spontaneous emission in spatially extended two-atom systems,” J. Mod. Opt. 50(18), 2765–2779 (2003). [CrossRef]

69. S. Dai, Q. Ma, M. K. Liu, T. Andersen, Z. Fei, M. D. Goldflam, M. Wagner, K. Watanabe, T. Taniguchi, M. Thiemens, F. Keilmann, G. C. A. M. Janssen, S.-E. Zhu, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, “Graphene on hexagonal boron nitride as a tunable hyperbolic metamaterial,” Nat. Nanotechnol. 10(8), 682–686 (2015). [CrossRef]

70. B. Zhang, S. J. You, and M. Lu, “Enhancement of spontaneous entanglement generation via coherent quantum feedback,” Phys. Rev. A 101(3), 032335 (2020). [CrossRef]

71. T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L. O’Brien, “Quantum computers,” Science 464(7285), 45–53 (2010). [CrossRef]

72. S. Rodt, S. Reitzenstein, and T. Heindel, “Deterministically fabricated solid-state quantum-light sources,” J. Phys.: Condens. Matter 32(15), 153003 (2020). [CrossRef]

Directional dipole radiations and long-range quantum entanglement mediated by hyperbolic metasurfaces

Abstract

1. Introduction

2. Radiation properties of a QE near the HMS

3. Collective behavior of two distant QEs coupled to the HMS

3.1 Master equation for the reduced density operator of the QEs

3.2 Long-range superradiance of the QEs system

4. Quantum entanglement at macroscopic distance

5. Conclusions

Appendix A: Surface plasmons supported by the HMSs

Appendix B: Time evolution of the reduced density matrix elements for the QEs

Funding

Disclosures

References

Cited By

Figures (10)

Equations (28)

Optics Express