Robust entanglement with 3D nonreciprocal photonic topological insulators

TL;DR

This paper demonstrates that entanglement between two-level systems near a 3D topological insulator interface is highly robust due to nonreciprocal surface plasmons, with potential applications in quantum technologies.

Contribution

It develops a master equation approach using the Green function to analyze entanglement in nonreciprocal, inhomogeneous environments, highlighting topological protection effects.

Findings

Entanglement remains robust despite material defects.

Nonreciprocal environment can prevent entanglement at close distances.

Entanglement is insensitive to variations in topological material parameters.

Abstract

We investigate spontaneous and pumped entanglement of two level systems in the vicinity of a photonic topological insulator interface, which supports a nonreciprocal (unidirectional), scattering-immune and topologically-protected surface plasmon polariton in the bandgap of the bulk material. To this end, we derive a master equation for qubit interactions in a general three-dimensional, nonreciprocal, inhomogeneous and lossy environment. The environment is represented exactly, via the photonic Green function. The resulting entanglement is shown to be extremely robust to defects occurring in the material system, such that strong entanglement is maintained even if the interface exhibits electrically-large and geometrically sharp discontinuities. Alternatively, depending on the initial excitation state, using a non-reciprocal environment allows two qubits to remain unentangled even for very close spacing. The topological nature of the material is manifest in the insensitivity of the entanglement to variations in the material parameters that preserve the gap Chern number. Our formulation and results should be useful for both fundamental investigations of quantum dynamics in nonreciprocal environments, and technological applications related to entanglement in two-level systems.

Equations

$$\left[\nabla \times {\mu }^{-1}(\mathbf{r},\omega )\nabla \times -\frac{{\omega }^2}{c^2}\epsilon (\mathbf{r},\omega )\right]\mathrm{E}(\mathbf{r},\omega )=i\omega {\mu }_0{\mathrm{j}}_s(\mathbf{r},\omega ),$$

$$\left[\nabla \times \nabla \times -\frac{{\omega }^2}{c^2}\epsilon (\mathbf{r},\omega )\right]\mathrm{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\mathrm{I}\delta (\mathbf{r},{\mathbf{r}}^{\prime })$$

$${\hat{\mathrm{P}}}_s(\mathbf{r},\omega )=-i\sqrt{\frac{\hslash {\epsilon }_0}{\pi }}\mathrm{T}(\mathbf{r},\omega )\cdot \hat{\mathrm{f}}(\mathbf{r},\omega ),$$

$$\mathrm{T}(\mathbf{r},\omega )\cdot {\mathrm{T}}^{†}(\mathbf{r},\omega )=\frac{1}{2i}\left[\mathbf{\epsilon }(\mathbf{r},\omega )-{\mathbf{\epsilon }}^{†}(\mathbf{r},\omega )\right],$$

$$\hat{\mathrm{E}}(\mathbf{r},\omega )=i\sqrt{\frac{\hslash }{\pi {\epsilon }_0}}\frac{{\omega }^2}{c^2}\int d{\mathbf{r}}^{\prime }\mathrm{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )\cdot \mathrm{T}(\mathbf{r},\omega )\cdot \hat{\mathrm{f}}({\mathbf{r}}^{\prime },\omega ),$$

$$2i\frac{{\omega }^2}{c^2}\int d^3{\mathrm{r}}^{\prime \prime }\mathrm{G}(\mathbf{r},{\mathbf{r}}^{\prime \prime },\omega )\cdot \mathrm{T}({\mathbf{r}}^{\prime \prime },\omega )\cdot {\mathrm{T}}^{†}({\mathbf{r}}^{\prime \prime },\omega ){\mathrm{G}}^{†}({\mathbf{r}}^{\prime },{\mathbf{r}}^{\prime \prime },\omega )$$

$$=\mathrm{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )-{\mathrm{G}}^{†}({\mathbf{r}}^{\prime },\mathbf{r},\omega ).$$

$$\mathrm{H}$$

$$-\sum _i{\int }_0^{\infty }d\omega ({\hat{\mathrm{d}}}_i\cdot \mathrm{E}({\mathbf{r}}_i,\omega )+H.c.),$$

$${\rho }_{T,I}={\rho }_I(0){\mathrm{R}}_0+\frac{-i}{\hslash }{\int }_0^{t}d{t}^{\prime }[{\mathrm{H}}_I(t^{\prime }),{\rho }_{T,I}(t^{\prime })]$$

$${\mathrm{H}}_I(t)=-\sum _i\left({\int }_0^{\infty }d\omega {\sigma }_i^{†}(t){\mathrm{d}}_i\cdot \mathrm{E}({\mathbf{r}}_i,\omega )e^{-i(\omega -{\omega }_l)t}+H.c.\right)$$

$${\partial }_t{\rho }_I$$

$$-\frac{1}{{\hslash }^2}{\int }_0^{t}d{t}^{\prime }{\mathrm{Tr}}_R\{\left[{\mathrm{H}}_I(t),\left[{\mathrm{H}}_I(t^{\prime }),{\rho }_{T,I}(t^{\prime })\right]\right]\}.$$

$${\partial }_t{\rho }_I=-\frac{1}{{\hslash }^2}{\int }_0^{\infty }d{t}^{\prime }{\mathrm{Tr}}_R\left\{\left[{\mathrm{H}}_I(t),\left[{\mathrm{H}}_I(t-t^{\prime }),{\rho }_I(t){\mathrm{R}}_0\right]\right]\right\}.$$

$${\mathrm{Tr}}_R\left\{{\mathrm{H}}_I(t){\mathrm{H}}_I(t-t^{\prime }){\rho }_I(t){\mathrm{R}}_0\right\}=\sum _{i,j}{\mathrm{d}}_{\alpha i}{\mathrm{d}}_{\beta j}$$

$${\int }_0^{\infty }d\omega e^{i({\omega }_0-\omega )t^{\prime }}{\sigma }_i^{†}{\sigma }_j{\rho }_{T,I}{\mathrm{Tr}}_R\left({\hat{\mathrm{E}}}_{\alpha }({\mathbf{r}}_i,\omega ){\hat{\mathrm{E}}}_{\beta }^{†}({\mathbf{r}}_j,\omega ){\mathrm{R}}_0\right)$$

$${\sigma }_i=|g_i⟩⟨e_i|,{{\sigma }_i}^{†}=|e_i⟩⟨g_i|$$

$${\mathrm{Tr}}_R\left({\mathrm{E}}_{\alpha }({\mathbf{r}}_i,\omega ){\mathrm{E}}_{\beta }^{†}({\mathbf{r}}_j,\omega ){\mathrm{R}}_0\right)=\frac{\hslash }{\pi {\epsilon }_0}\frac{{\omega }^4}{c^4}$$

$$\int d^3\mathbf{r}{\mathrm{G}}_{\alpha \gamma }({\mathbf{r}}_i,\mathbf{r},\omega )\left[\frac{{\epsilon }_{\gamma {\gamma }^{\prime }}(\mathbf{r},\omega )-{\epsilon }_{\gamma {\gamma }^{\prime }}^{†}(\mathbf{r},\omega )}{2i}\right]{\mathrm{G}}_{{\gamma }^{\prime }\beta }^{\ast }({\mathbf{r}}_j,\mathbf{r},\omega )$$

$$=\frac{\hslash }{2i\pi {\epsilon }_0}\frac{{\omega }^2}{c^2}\left({\mathrm{G}}_{\alpha ,\beta }({\mathbf{r}}_i,{\mathbf{r}}_j,\omega )-{\mathrm{G}}_{\beta ,\alpha }^{\ast }({\mathbf{r}}_j,{\mathbf{r}}_i,\omega )\right).$$

$${\mathrm{Tr}}_R\left\{{\mathrm{H}}_I(t){\mathrm{H}}_I(t-t^{\prime }){\rho }_I(t){\mathrm{R}}_0\right\}=\frac{\hslash }{2i\pi {\epsilon }_0{c}^2}\sum _{i,j}{\sigma }_i^{†}{\sigma }_j{\rho }_I(t)$$

$${\int }_0^{\infty }\left({\mathrm{d}}_{\alpha i}{\mathrm{G}}_{\alpha \beta }({\mathbf{r}}_i,{\mathbf{r}}_j,\omega ){\mathrm{d}}_{\beta j}-{\mathrm{d}}_{\beta j}{\mathrm{G}}_{\beta \alpha }^{\ast }({\mathbf{r}}_j,{\mathbf{r}}_i,\omega ){\mathrm{d}}_{\alpha i}\right)$$

$${\omega }^{2}d\omega e^{i({\omega }_0-\omega )t^{\prime }}.$$

$${\mathrm{Tr}}_R\left\{{\mathrm{H}}_I(t-t^{\prime }){\rho }_I(t){\mathrm{R}}_0{\mathrm{H}}_I(t)\right\}=\frac{\hslash }{2i\pi {\epsilon }_0{c}^2}\sum _{i,j}{\sigma }_j{\rho }_I(t){\sigma }_i^{†}$$

$${\int }_0^{\infty }\left({\mathrm{d}}_{\beta i}{\mathrm{G}}_{\beta \alpha }({\mathbf{r}}_i,{\mathbf{r}}_j,\omega ){\mathrm{d}}_{\alpha j}-{\mathrm{d}}_{\alpha j}{\mathrm{G}}_{\alpha \beta }^{\ast }({\mathbf{r}}_j,{\mathbf{r}}_i,\omega ){\mathrm{d}}_{\beta i}\right)$$

$${\omega }^{2}d\omega e^{i({\omega }_0-\omega )t^{\prime }}.$$

$$\mathcal{P}{\int }_{-\infty }^{\infty }\frac{\mathrm{Re}{\mathrm{G}}_{\alpha \beta }}{\omega -{\omega }_0}d\omega =-\pi \mathrm{Im}{\mathrm{G}}_{\alpha \beta }$$

$$\mathcal{P}{\int }_{-\infty }^{\infty }\frac{\mathrm{Im}{\mathrm{G}}_{\alpha \beta }}{\omega -{\omega }_0}d\omega =\pi \mathrm{Re}{\mathrm{G}}_{\alpha \beta }.$$

$${\partial }_t{\rho }_s(t)=-\frac{i}{\hslash }\left[{\mathrm{H}}_s+{\mathrm{V}}^{AF},{\rho }_s(t)\right]+\mathcal{L}\rho (t),$$

$$\mathcal{L}{\rho }_s(t)$$

$$\sum _i\frac{{\Gamma }_{ii}({\omega }_0)}{2}\left(2{\sigma }_i{\rho }_s(t){\sigma }_i^{†}-{\sigma }_i^{†}{\sigma }_i{\rho }_s(t)-{\rho }_s(t){\sigma }_i^{†}{\sigma }_i\right)$$

$$+\sum _{i,j}^{i\ne j}\frac{{\Gamma }_{ij}({\omega }_0)}{2}\left(\left[{\sigma }_j{\rho }_s(t),{\sigma }_i^{†}\right]+\left[{\sigma }_i,{\rho }_s(t){\sigma }_j^{†}\right]\right)$$

$$+\sum _{i,j}^{i\ne j}{g}_{ij}({\omega }_0)\left(\left[{\sigma }_j{\rho }_s(t),-i{\sigma }_i^{†}\right]+\left[i{\sigma }_i,{\rho }_s(t){\sigma }_j^{†}\right]\right).$$

$${\Gamma }_{ij}({\omega }_0)=\frac{2{\omega }_0^2}{{\epsilon }_0\hslash c^2}\sum _{\alpha ,\beta =x,y,z}{\mathrm{d}}_{\alpha i}\mathrm{Im}\left({\mathrm{G}}_{\alpha \beta }({\mathbf{r}}_i,{\mathbf{r}}_j,{\omega }_0)\right){\mathrm{d}}_{\beta j},$$

$$g_{ij}({\omega }_0)=\frac{{\omega }_0^2}{{\epsilon }_0\hslash c^2}\sum _{\alpha ,\beta =x,y,z}{\mathrm{d}}_{\alpha i}\mathrm{Re}({\mathrm{G}}_{\alpha \beta }({\mathbf{r}}_i,{\mathbf{r}}_j,\omega )){\mathrm{d}}_{\beta j}.$$

$${\mathrm{H}}_s=\sum _i\hslash \Delta {\omega }_i{\sigma }_i^{†}{\sigma }_i,$$

$${\mathrm{V}}^{AF}$$

$$-\hslash \left({\Omega }_2{e}^{-i{\Delta }_{l}t}{\sigma }_2^{†}+{\Omega }_2^{\ast }{e}^{i{\Delta }_{l}t}{\sigma }_2\right)$$

Full text

I Introduction

Entanglement as a quantum resource is important for a range of emerging applications, including quantum computing 5 and quantum cryptography 6 . A main obstacle to the development of entanglement-based systems is decoherence associated with the unavoidable coupling between a quantum system and the degrees of freedom of the surrounding environment 10 . However, reservoir engineering methods have changed the idea of trying to minimize coupling to the environment to one of modifying the properties of the environment in order to achieve a desired state. These methods include using dissipative dynamics Plenio02 ; Hartmann ; Verstraete ; Krauter ; Lucas ; Sarlette , recently extended to systems out of thermal equilibrium Bellomo12 ; Bellomo13 ; Bellomo13EPL ; Bellomo13NJP ; Bellomo15 , as well as, e.g., exploiting the effect of measurements and feedback to achieve a desired final state Mancini ; Stevenson .

Another emerging resource for reservoir engineering is the use of nonreciprocal environments Zoller3 . In particular, there has been considerable investigation of quantum spin networks in chiral waveguides Garcia2 ; Zoller0 ; Zoller1 ; Zoller2 ; MLS ; SJ ; MS . The previous work on spin dynamics in quantum chiral environments has focused on one-dimensional (1D) waveguide models. Here, we investigate two-level (spin) quibit interactions mediated by uni-directional surface-plasmon-polaritions (SPPs) at the interface of a photonic topological insulator (PTI) and a topological-trivial material.

PTIs represent a broad class of materials that are attracting wide interest for both fundamental and applied reasons Haldane ; Rechtsman ; Lu ; Kh . Perhaps their most celebrated aspect is their ability to support SPPs that are unidirectional, propagate in the bulk bandgap, and are topologically protected from backscattering at discontinuities Joannopoulos ; Veronis ; Joannopoulos2 ; Rechtsman1 ; Arthur ; Mario2 ; Hassani2 ; Hassani3 . PTIs can be broadly divided into two classes, (i) those with broken time reversal (TR) symmetry, which are photonic analogs of quantum Hall insulators (photonic quantum Hall effect (PQHE)), and (ii) those that are time-reversal-invariant but have broken inversion symmetry, which are photonic analogs of electronic topological insulators/quantum spin Hall insulators (photonic quantum spin Hall Effect (PQSHE)). Although as a specific example we consider PTIs of the PQHE type, the formulation presented here is general.

In this work, we develop a master equation (ME) for three-dimensional (3D), nonreciprocal, inhomogeneous and lossy environments, based on the macroscopic canonical quantization scheme described in Welsch0 ; Welsch1 ; Trung , extended to nonreciprocal media Buh . In Section II.1 we derive the master equation, and in Section II.2 we present the equations for concurrence as a measure of entanglement. In Section III we consider the topological aspect of concurrence for a PQHE-type PTI system consisting of a plasma continuum. Then, qubit entanglement dynamics are examined for several waveguiding systems. We focus of the aspects unique to the topological and nonreciprocal environment, such as the preservation of entanglement in the presence of large defects. Three appendices present a discussion of various approximations used in the development of the ME, a comparison with previous 1D chiral MEs, and a derivation of the unidirectional concurrence.

II Theoretical model

In this section, we first derive a general ME valid for both reciprocal and nonreciprocal, inhomogeneous and lossy environments. This form is valid for 3D, 2D and 1D systems since it is expressed in terms of the electromagnetic Green function. Then, we present concurrence expressions for the unidirectional case. The physical system we will consider is that of two qubits at the interface of a PTI and another (eventually topologically trivial) medium, as depicted in Fig. 1, although the development is completely general.

II.1 Master equation for general 3D nonreciprocal environments

We consider qubits with transition frequency $\omega_{0}$ interacting through a general nonreciprocal environment. For a derivation in the reciprocal case, see SteveH .

The classical electric field satisfies

[TABLE]

where $c$ is the vacuum speed of light, $\mu(\mathbf{r},\omega)$ and $\boldsymbol{\varepsilon}(\mathbf{r},\omega)$ are the material permeability and permittivity, and $\boldsymbol{\mathrm{j}}_{s}(\mathbf{r},\omega)$ is the noise current. In this work, we suppose that the medium is non-magnetic, $\mathbf{\mu}(\mathbf{r},\omega)=\mathbf{I}$, where $\boldsymbol{\mathrm{I}}$ is the unit dyad, but that the permittivity is a tensorial quantity. By defining the noise current in terms of polarization as $\boldsymbol{\mathrm{j}}_{s}=-i\omega\boldsymbol{\mathrm{P}}_{s}$, which is associated with material absorption by the fluctuation-dissipation theorem, the electric field Green tensor is the solution of

[TABLE]

and the electric field is $\boldsymbol{\mathrm{E}}(\mathbf{r},\omega)=(\omega^{2}/c^{2}\varepsilon_{0})\int_{V}d\mathbf{r}^{\prime}\boldsymbol{\mathrm{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)\cdot\boldsymbol{\mathrm{P}}_{s}(\mathbf{r}^{\prime},\omega)$. Following the standard macroscopic canonical quantization Welsch0 ; Welsch1 ; Trung , the noise polarization can be expressed in term of the bosonic field annihilation operator as Buh

[TABLE]

where

[TABLE]

and, for the special case of a symmetric permittivity tensor (e.g., a reciprocal medium), $\boldsymbol{\mathrm{T}}(\mathbf{r},\omega)=\sqrt{\mathrm{Im}\boldsymbol{\varepsilon}(\mathbf{r},\omega)}$. The bosonic field operators $\hat{\boldsymbol{\mathrm{f}}}(\mathbf{r},\omega)$ obey the commutation relations $[\hat{\mathrm{f}}_{j}(\mathbf{r},\omega),\hat{\mathrm{f}}^{\dagger}_{j^{\prime}}(\mathbf{r}^{\prime},\omega^{\prime})]=\delta_{jj^{\prime}}\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(\omega-\omega^{\prime})$ and $[\hat{\mathrm{f}}_{j}(\mathbf{r},\omega),\hat{\mathrm{f}}_{j^{\prime}}(\mathbf{r}^{\prime},\omega^{\prime})]=0$. The noise polarization operator generates the electric field operator

[TABLE]

where $\boldsymbol{\mathrm{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)$ is the classical electric field Green tensor. The nonreciprocal Green tensor has the following useful property Buh

[TABLE]

Under the dipole approximation, the governing Hamiltonian of a system of qubits (two level atoms) interacting with the surrounding environment can be written as

[TABLE]

where the right side can be decomposed into the reservoir Hamiltonian $\mathrm{H}_{r}$ (first term), the qubit Hamiltonian $\mathrm{H}_{s}$ (second term), and the interaction Hamiltonian $\mathrm{H}_{sr}$ (third term). We can modify the total Hamiltonian to include the coherent drive (external laser pump) Hamiltonian $\mathrm{V}^{AF}$, given later ((II.1)). We transform to a frame rotating with the laser frequency $\omega_{l}$ ($\mathrm{H}\rightarrow\hat{\mathrm{U}}^{\dagger}(t)\mathrm{H}\hat{\mathrm{U}}(t),~{}\hat{\mathrm{U}}(t)=e^{-i\omega_{l}\sum_{i}{\sigma}^{\dagger}_{i}{\sigma}_{i}t}$) and write the total density matrix of the qubit system and reservoir according to the Schr$\ddot{\mathrm{o}}$dinger equation $\partial_{t}\rho_{T}=-i[\mathrm{H},\rho_{T}]/\hbar$, then we transform to the interaction picture ($\hat{\mathrm{O}}_{I}=\hat{\mathrm{U}}^{\dagger}(t)\mathrm{H}\hat{\mathrm{U}}(t),~{}\hat{\mathrm{U}}(t)=e^{-i(\mathrm{H}_{s}+\mathrm{H}_{r})t/\hbar}$) where $\partial_{t}\rho_{T,I}=-i[\mathrm{H}_{I},\rho_{T,I}]$ with $\mathrm{H}_{I}=\mathrm{H}_{sr,I}$. We integrate to find

[TABLE]

where $\mathrm{R}_{0}$ is the initial reservoir density matrix. In the interaction picture, by considering $\Gamma_{ii}\ll\omega$ for optical frequencies we make the rotating wave approximation (RWA) in $\mathrm{H}_{I}$ and drop the rapidly varying counter-rotating terms proportional to ${\sigma}^{\dagger}(t^{\prime})\boldsymbol{\mathrm{f}}^{\dagger}(\mathbf{r}^{\prime},\omega)e^{i(\omega_{l}+\omega)t^{\prime}}$ and its Hermitian conjugate. The interaction Hamiltonian in the interaction picture reduces to

[TABLE]

To find the system density matrix we insert (8) into the interaction picture Schr$\ddot{\mathrm{o}}$dinger equation and trace over the reservoir,

[TABLE]

Aside from the rotating wave approximation, we apply a number of other approximations to the density matrix to simplify this further (see Appendix I). We first take the mean initial system reservoir coupling to be zero such that $\mathrm{Tr}_{R}\{\frac{-i}{\hbar}\left[\mathrm{H}_{I},\rho_{I}(0)\mathrm{R}_{0,I}\right]\}=0$. Then we apply the Born approximation, which states that the reservoir will be largely unaffected by its interaction by the system. Next, we assume that the evolution of the density matrix only depends on its current state (Born-Markov approximation). The Born-Markov approximation comes from the assumption that the reservoir relaxation time is much faster than the relaxation time of the system, and so the memory effect of the reservoir can be ignored. Lastly, we make a second Markov approximation, extending the upper limit of the time integral to infinity to produce a fully Markovian equation. With these simplifications we have

[TABLE]

We suppose that the atomic transition frequency of the qubits is $\omega_{0}$. Then, for the first term in (11) we have

[TABLE]

where

[TABLE]

are the atomic lowering/raising operators describing energy level transitions for each qubit, and where it is supposed that one of the qubits is polarized along $\alpha$ and the other one is polarized along $\beta$. Considering (5) for the nonreciprocal Green tensor and $\mathrm{Tr}_{R}\{\hat{\boldsymbol{\mathrm{f}}}(r,\omega)\hat{\boldsymbol{\mathrm{f}}}^{\dagger}(r^{\prime},\omega^{\prime})\mathrm{R}_{0}\}=(\bar{n}(\omega)+1)\delta(r-r^{\prime})\delta(\omega-\omega^{\prime})$ with zero thermal photon occupation $\bar{n}(\omega)=0$, it can be easily shown that

[TABLE]

Thus, we have

[TABLE]

Following the same procedure for the second term in (11),

[TABLE]

Replacing (II.1) and (II.1) in (11) and performing the time integral over $t^{\prime}$ gives the evolution of the density matrix in the interaction picture, where we have used the Kramers-Kronig relation

[TABLE]

Transforming back to the Schr$\ddot{\mathrm{o}}$dinger picture, we obtain the master equation for the two-level system dynamics

[TABLE]

where

[TABLE]

Equation (II.1) is one of the primary results of this work, and is applicable to both reciprocal and nonreciprocal environments and an arbitrary number of qubits. In (II.1), $\mathcal{L}$ is the Lindblad superoperator for the general nonreciprocal medium, involving the dissipative decay rate, $\Gamma_{ij}(\omega_{0})$, and the coherent coupling terms, $g_{ij}(\omega_{0})$, in terms of the electromagnetic Green dyadic,

[TABLE]

The Hamiltonian of the decoupled qubits is

[TABLE]

where $\Delta\omega_{i}=\omega_{0}-\omega_{l}-\delta_{i}$, with $\delta_{i}=g_{ii}$ being the Lamb shift and $\omega_{l}$ is the laser frequency of an external source. The Lamb shift for optical emitters is in general on the order of a few GHz, therefore the effect of the Lamb-shift for optical frequencies is small $(\omega_{i}\sim 10^{15}~{}\mathrm{Hz},~{}\delta_{i}\sim 10^{9}~{}\mathrm{Hz})$, and can be ignored, or assumed to be accounted for in the definition of the transition frequency $\omega_{0}$. In (18), the term

[TABLE]

represents the external coherent drive applied to each qubit at laser frequency $\omega_{l}$. Due to its large amplitude we treat the drive field as a $c$-number where $\Omega_{i}=\mathbf{d}_{i}\cdot\mathbf{E}_{0}^{i}/\hbar$ is a Rabi frequency and $\Delta_{l}=\omega_{0}-\omega_{l}$ is the detuning parameter.

For the reciprocal case where $\Gamma_{ij}=\Gamma_{ji}$ and $g_{ij}=g_{ji}$ it can be shown that (II.1) is the well-known reciprocal (bidirectional) master equation ME_O ; Garcia . In the reciprocal case, some terms associated with $g_{ij}=g_{ji}$ cancel each other out and are eliminated from the dissipative term. For example, ${\sigma}_{i}{\rho}_{s}(t){\sigma}_{i}^{\dagger},i\neq j$, appears in the nonreciprocal case but is absent in the reciprocal case.

For a system of two qubits, (II.1) can be written in the simple form

[TABLE]

A comparison with previous 1D chiral ME formulations is provided in Appendix II.

II.2 Transient entanglement: Unidirectional SPP-assisted qubit communication

In this work, all numerical results are computed using the master equation (18) with the general 3D Lindblad superoperator (II.1), where the Green tensor for complicated environments is obtained numerically. However, as shown in Appendix III, if the system of qubits are communicating through a strongly nonreciprocal environment e.g., ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{1},{\boldsymbol{\mathrm{r}}}_{2})=0$ ($\Gamma_{12}=g_{12}=0$) and ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{2},{\boldsymbol{\mathrm{r}}}_{1})\neq 0$, then the concurrence (as a measure of entanglement Wootters ) is

[TABLE]

where it has been assumed that the qubits are both polarized along the $y$-axis. This is the general unidirectional result. Concurrence reaches its maximum value at $t=1/\Gamma_{11}$, such that $\mathcal{C}_{\mathrm{max}}=2\sqrt{\tilde{\Gamma}^{2}_{12}/4+\tilde{g}^{2}_{21}}/e$, where $\tilde{\Gamma}_{21}$ and $\tilde{g}_{21}$ are rates normalized by $\Gamma_{11}$.

Although the Hamiltonian in nonreciprocal systems is non-Hermitian, it can be seen (Appendix III) that the density matrix is Hermitian, probability conservation holds ($\mathrm{Tr}(\rho)=1$), and that diagonal elements of the density operator can be interpreted as population densities, as for Hermitian Hamiltonians.

For two identical qubits interacting through a reciprocal medium,

[TABLE]

One of the main differences between the concurrence in the reciprocal case (II.2), and in the unidirectional case (II.2), is the presence of the sinusoidal term in (II.2). When $g_{12}$ is strong enough, this sinusoidal term causes oscillations in the transient concurrence related to photons being recycled between the two qubits, with a period that corresponds to the round trip time of the coupled qubits through the reciprocal medium (Rabi oscillations). For the unidirectional case (II.2) Rabi oscillations can not occur.

It was shown in Garcia that for qubits coupled to an infinite reciprocal waveguide system, the positions of $\Gamma_{ij}$ maxima/minima correspond to positions of $g_{ij}$ minima/maxima (for finite waveguides, see H1 ). Thus, in general, coherent and dissipative regimes become dominant at different separations between emitters. It was further shown in Garcia that for an infinite reciprocal plasmonic waveguide the best entanglement was obtained when $\Gamma_{ij}$ was large and $g_{ij}$ was small (forming the dissipative regime), which forces a restriction on the positioning of the qubits in the reciprocal case. However, in the unidirectional case the qubit positioning is unimportant, as detailed in Garcia2 , and the qubits can be anywhere in the coherent or dissipative regimes, which is a practical advantage of these unidirectional systems.

In order to demonstrate this difference between reciprocal and unidirectional systems, we consider two cases of pure dissipative and pure coherent qubit communication for a model system where we simply assign the Green function values. Fig. 2 shows that the pure dissipative regime is dominant for the reciprocal case while the dissipative or coherent nature of the qubit communication is unimportant for the unidirectional case.

III Numerical results

A unidirectional SPP can be provided by the interface between a PTI and a topologically-trivial material. When operated in a common bandgap of the two materials (or if the trivial medium is opaque), the SPP is unidirectional, topologically protected from back-scattering, and diffraction-immune, providing an ideal implementation of a strongly nonreciprocal system for qubit interactions. Although here we implement a PTI as a PQHE using a continuum plasma Mario2 -Hassani3 , many other implementations of PTIs are possible, of both PQHE and PQSHE types, and qualitatively would behave in a similar manner.

III.1 Continuum photonic topological insulator realization of a nonreciprocal surface plasmon polariton environment

We assume a magnetized plasma having the permittivity tensor

[TABLE]

where

[TABLE]

and where $\omega_{c}=\left(q_{e}/m_{e}\right)\mathrm{B}_{z}$ is the cyclotron frequency ($\mathrm{B}_{z}$ is the applied bias field), $\omega_{p}^{2}=N_{e}q_{e}^{2}/\varepsilon_{0}m_{e}$ is the squared plasma frequency ($N_{e}$ is the free electron density and $q_{e}$ and $m_{e}$ are the electron charge and mass, respectively), and $\nu$ is the collision frequency. Initially, we set $\nu=0$ to focus on the effect of unidirectionality, but later the effect of loss is considered. The magnetized plasma is able to support a bulk TE mode with dispersion $k^{2}_{\mathrm{TE}}=\varepsilon_{33}(\omega/c)^{2}$ and a bulk TM mode with dispersion $k^{2}_{\mathrm{TM}}=\varepsilon_{\mathrm{eff}}(\omega/c)^{2}$, where $\varepsilon_{\mathrm{eff}}=(\varepsilon^{2}_{11}-\varepsilon^{2}_{12})/\varepsilon_{11}$ Arthur . Both bulk modes are reciprocal. The Chern number of the bulk TE mode is trivial, and so TE modes are not considered further in this work. The Chern numbers of the bulk TM modes are nonzero, and at the interface of the magnetized plasma and a topologically-trivial (simple) medium the gap Chern number is $\mathrm{C_{gap}=1}$ Mario2 ; Hassani2 , indicating the presence of one nonreciprocal, backscattering-immune TM-SPP that crosses the bandgap (bandstructure is shown later, in Fig. 4).

III.2 Entanglement evaluation in different environments

We first consider the behavior of the concurrence for qubits in several different environments, and establish that the best entanglement occurs for a PTI/opaque medium interface. Fig. 3a, shows a comparison of concurrence between four different cases of two qubits interacting through: 1) vacuum, 2) at the interface of a gold half-space and vacuum, 3) at the interface of a magnetized plasma and vacuum, and 4) at the interface of a magnetized plasma and an opaque medium. Here and in the following, the Green function is calculated numerically COMSOL .

The system of qubits were initially prepared in state $\left|4\right>=\left|e_{1}\right>\otimes\left|g_{2}\right>$, such that the left qubit is initially in the excited state while the right qubit is in the ground state. It can be seen that the interface between the magnetized plasma and the opaque medium has higher concurrence then the other cases, due to the existence of a strong SPP and the fact that there can be no radiation into either bulk half-space. Thus, in the following, we focus on the magnetized plasma/opaque medium geometry.

Figure 3b shows the excited unidirectional SPP at the interface of the magnetized plasma and the opaque medium, demonstrating the unidirectional nature of the SPP, and Fig. 3c shows the case of pumped concurrence, where the qubit depopulation is compensated by applying an external laser source in resonance with the atomic transition frequency. The pump intensity must be chosen carefully, as illustrated in Fig. 3d, which shows the steady state concurrence for a wide range of laser intensities (a laser pump can be applied to the qubits via, e.g., a fiber penetrating into the material). It can be seen that the laser intensity can not be too large, otherwise the qubits will interact mostly with the laser. Ideally, the pump should be strong enough to keep the system interacting, but weak enough for the qubit interaction to dominate the dynamics. It is clear from Fig. 3d that unequal pumping leads to larger steady state concurrence.

III.3 Topological aspect of entanglement

In this section we briefly show the topological aspect of entanglement in a PTI system. Figure 4 shows the reciprocal bulk bands (solid blue) for the biased plasma, and the unidirectional gap-crossing SPP (dashed red) dispersion for a biased-plasma/opaque medium interface, for different values of bias. For $\omega_{c}>0$ the gap Chern number is -1 Mario2 ; Hassani2 , and there is a positive-traveling SPP ($v_{g}=d\omega/dk>0$), topologically-protected against backscattering. At $\omega_{c}=0$ the gap closes, the material becomes topologically-trival (gap Chern number is 0), and there exists a reciprocal SPP. For $\omega_{c}<0$ the gaps reopens, the gap Chern number is 1, and there is a negative-traveling SPP ($v_{g}<0$), topologically-protected against backscattering.

Figure 5 shows the concurrence when the left dot has an initial excitation (state $\left|4\right>=\left|e_{1},g_{2}\right>$). The concurrence is rather insensitive to the bias as long as the topology does not change, however, when the gap closes and reopens the concurrence vanishes.

III.4 Preserving entanglement in the presence of large defects

Perhaps the most important aspect of using PTIs for entanglement is the possibility of robust SPPs, topologically-immune to backscattering (and immune to diffraction if operated in the bulk bandgap) in the presence of any arbitrary large obstacle or defect. To examine this, we compare two cases: 1) the interface between an opaque medium and a biased plasma, and 2) the interface between the same opaque medium and an unbiased plasma.

In the nonreciprocal case, this unidirectional and scattering-immune SPP provides the ability to preserve the entangled state of two qubits in plasmonic systems even in the presence of very non-ideal interfaces. Figure 6 shows the transient concurrence for the cases of biased/unbiased plasmas with flat and defected interfaces. Although for the flat interface the biased plasma provides better concurrence then the reciprocal (unbiased) case, this could be perhaps altered by adjustment of the two material half-space properties. However, the point is that in the presence of a defect, as shown in the right panel, the reciprocal SPP suffers from a strong reflection at the defect, as expected, whereas the nonreciprocal SPP (biased plasma) detours around the defect, leading to the same concurrence as without the defect.

III.5 Finite-width waveguide

The previous results were for an infinitely-wide interface. In this section we examine the effect of lateral confinement of the SPP Hassani2 . Figure 7a shows the finite-width waveguide geometry. In order to efficiently confine the SPP along the propagation axis, the plasma is extended past the interface to form partially-extended sidewalls. Only partial side walls are needed to prevent radiation in space, since the SPP is confined to the vicinity of the interface.

Lateral confinement of the unidirectional SPP improves both the transient and steady state (pumped) concurrence. Fig. 7b shows the transient and steady state concurrence of two qubits initially prepared in state $\left|4\right>$. In comparison to Fig. 3a, it can be seen that lateral confinement increases both the maximum transient concurrence and the steady state concurrence. Figure 7c shows the dynamics of the qubits under external pumping, where $\rho_{11},~{}\rho_{22},~{}\rho_{33},~{}\mathrm{and}~{}\rho_{44}$ are the probabilities of finding both qubits to be in ground state, both qubits in the excited state, the first qubit in the ground state and the second qubit in the excited state, and vice versa, respectively. Figure 7d shows the steady state concurrence for a wide range of pump values. The behavior is similar to the case of the infinite interface, Fig. 3d, except that the range of pump values that result in large steady state concurrence is extended, and the maximum achievable steady state concurrence is larger in the case of the finite-width waveguide.

In Fig. 8, qubit concurrence is shown for a finite-width waveguide having a defect which spans the entire waveguide width. It can be seen that the concurrence is minimally affected by the defect. Although not shown, as with Fig. 6, in the reciprocal (unbiased) case the defect eliminates the concurrence.

III.6 Effect of different initial state preparations

An interesting behavior of the concurrence arising from having a unidirectional SPP is that, e.g., if the medium supports only a right going SPP, then the initially excited qubit should be the left qubit, otherwise the qubits remain unentangled, as shown in Fig. 9a for the unpumped case $\Omega_{1}=\Omega_{2}=0$. Figure 9b shows the dynamics of the qubits for this unpumped case. It can be seen that $\rho_{33}$, which is the probability of finding the right qubit in the excited state and the left qubit in the ground state, starts from 1 and then drops rapidly. However, $\rho_{44}$, which is the probability of finding the excitation being in the left qubit with the right qubit in the ground state, is always zero, meaning that the excitation lost from the right qubit never gets captured by the left qubit. This behavior is particular to a unidirectional environment, and allows for keeping two qubits disentangled at any qubit separation, even if one of them carries an excitation.

However, by applying an external pump we can achieve non-zero concurrence, as also depicted in Fig. 9a

The pump is turned on at $\mathrm{t}=0$, and instead of immediately becoming non-zero, the concurrence remains zero for a period of time, then starts raising as a sudden birth in concurrence and reaches a non-zero steady state value. This delayed sudden-birth is quite different from the pumped reciprocal case.

It is also possible to consider different initial states which can give other possible unidirectional SPP assisted dynamical evolutions. Figure 10 shows the case of the initial state being the maximally entangled Bell state $\left|\Psi_{\mathrm{Bell}}\right>=(\left|1\right>+\left|2\right>)/\sqrt{2}$.

We consider that the qubits are interacting through the finite-width waveguide depicted in Fig. 7a. Figure 10a shows the time evolution of the concurrence for both pumped and non-pumped cases. In contrast to the previous cases, the concurrence starts from one due to the maximum degree of entanglement of the initial Bell state. For the non-pumped case the concurrence diminishes in time as the system becomes disentangled, resulting in a sudden death of entanglement. It remains zero for a period of time, then the entanglement experience a rebirth before decaying exponentially at long times. For the externally pumped case, the concurrence exponentially decays but the qubits do not become completely disentangled. Fig. 10b shows the dynamics of the qubits for the pumped case. The population probabilities $\rho_{11}$ and $\rho_{22}$ start from 0.5 due to the Bell state preparation. An interesting behavior in the qubit dynamics is the unequal steady state values $\rho_{33}$ and $\rho_{44}$ values under pumping with equal intensities $|\Omega_{1}|=|\Omega_{2}|$ (in the reciprocal case, $\rho_{33}=\rho_{44}$).

III.7 Lossy biased plasma

In a lossy medium the SPP loses power as it propagates along the interface, resulting in weaker qubit entanglement. In order to study the effect of loss, we suppose the qubits are interacting through an infinitely-wide interface as considered in Fig. 3, but for three different collision frequencies; $\nu=0$, $\nu/2\pi=270$ MHz and $\nu/2\pi=500$ MHz. Qubits are initially prepared in the state $\left|4\right>=\left|e_{1}\right>\otimes\left|g_{2}\right>$. Figure 11, left panel, shows the transient concurrence. Increasing the collision frequency reduces the concurrence, and for collision frequencies greater than 500 MHz loss dominates the system and an entangled state is not achievable for this relatively wide qubit separation of $1.6\lambda_{0}$.

The right panel of Fig. 11 shows the steady state concurrence of the pumped system, versus pumping intensity. In comparison to the lossless case (Fig. 3d), the range of pump intensities that give non-zero steady state concurrence has decreased, and the maximum achievable concurrence value is diminished compared to the lossless case.

IV Conclusions

We have derived a master equation for qubit dynamics in a general three-dimensional, nonreciprocal, inhomogeneous and lossy environment. Spontaneous and pumped entanglement were investigated for two qubits in the vicinity of a photonic topological insulator interface, which supports a nonreciprocal (unidirectional), scattering-immune surface plasmon polariton in the bandgap of the bulk material. We have illustrated the topological nature of the entanglement, and it was shown that large defects in the interface do not impact entanglement for the PTI case, whereas a defect has considerable effect for the reciprocal case. Several initial qubit states were considered, as well as the influence of pump intensity and material loss. Particularities arising from the unidirectional nature of the qubit communication were highlighted.

Acknowledgments

The authors would like to thank Stefan Buhmann for helpful discussions.

APPENDIX

Appendix I: Master equation approximations

Here we briefly discuss the approximations used in the derivation of the master equation (18)-(II.1).

The first approximation made in the derivation is the rotating wave approximation (RWA) where in the interaction picture we drop the rapidly varying counter-rotating terms in $\mathrm{H}_{I}$. This approximation is valid for $\Gamma_{ii}\ll\omega$. The qubit transition frequency is $\omega/2\pi=200$ THz, and we assume a dipole moment $\mathrm{d}=60$ D. For an interface made of lossless biased plasma, $\Gamma_{ii}/2\pi\sim 450$ MHz, and for the lossy biased plasma with $\nu/2\pi=500$ MHz, $\Gamma_{ii}/2\pi\sim 2$ GHz. For the non biased plasma-opaque medium interface (interface supporting reciprocal SPP) $\Gamma_{ii}/2\pi\sim 75$ MHz. In all cases the condition for the validity of the RWA is strongly met.

We also applied the Born-Markov approximation, which comes from the assumption that the reservoir relaxation time, $\tau_{R}$, is much faster than the relaxation time of the qubit system $\tau_{S}=1/\Gamma_{ii}$. This allows for the expansion of the exact equation of motion for the density matrix up to second order, and makes the quantum master equation local in time. For a nonreciprocal medium the fluctuation-dissipation theorem Buh is

[TABLE]

where $\mathrm{N}(\omega,\mathrm{T})=2/\left(\mathrm{exp}(\hbar\omega/k_{B}\mathrm{T})-1\right)$ for negative frequencies and $\mathrm{N}(\omega,\mathrm{T})=1+2/\left(\mathrm{exp}(\hbar\omega/k_{B}\mathrm{T})-1\right)$ for positive frequencies, where $k_{B}$ is Boltzmann’s constant. Regarding $\boldsymbol{\mathrm{E}}(\mathbf{r},\omega)=(\omega^{2}/c^{2}\varepsilon_{0})\int_{V}d\mathbf{r}^{\prime}\boldsymbol{\mathrm{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)\cdot\boldsymbol{\mathrm{P}}(\mathbf{r}^{\prime},\omega)$, it can be shown that

[TABLE]

which reduces in the reciprocal case to

[TABLE]

The bath relaxation time can be estimated by looking at the decay time of the correlation

[TABLE]

The Green function consists of homogeneous (vacuum) and scattered terms, and $\tau_{R}$ will be dominated by the slower scattered field contribution (for the vacuum term, $\tau_{R}(T)=\hbar/\pi k_{B}T$ BP , so that $\tau_{R}(300K)\sim 10$ fs). For the scattered part of the Green function for an interface made of non-biased plasma-opaque medium (interface supporting reciprocal SPP), using the Green function in DN , $\tau_{R}\sim 10^{-11}$ s for $\nu=500$ MHz and $\nu=270$ MHz, whereas $\tau_{S}=1/\Gamma_{ii}\sim 10^{-8}$ s, so that we can ignore the reservoir relaxation time.

Appendix II: Comparison with previous 1D chiral theory

Here we discuss the relation between the general ME we derived in terms of the exact electromagnetic Green function, resulting, for two qubits, in the Lindblad (II.1), and the 1D phenomenological ME for two level systems coupled to a 1D chiral reservoir presented in Garcia2 ; Zoller2 (see also Zoller3 ; Zoller1 ). The 1D chiral theory is based on the notion of right and left, defining couplings $\gamma_{\textrm{R,L}}$, whereas the theory presented here is based on qubit interactions $\Gamma_{ij}$; note that $\Gamma_{ij}$ plays the role of a $\Gamma_{\textrm{right}}$ if $x_{i}>x_{j}$, but plays the role of $\Gamma_{\textrm{left}}$ if $x_{i}<x_{j}$. To facilitate the comparison with the 1D chiral theory we will assume two qubits with positions $x_{1}$ and $x_{2}$, with $x_{2}>x_{1}$. In Garcia2 ; Zoller2 phenomenological quantities $\gamma_{iR},\gamma_{1L}$ for $i=1,2$ are utalized, and setting $\gamma_{1R}=\gamma_{2R}=\gamma_{R}$ and $\gamma_{1L}=\gamma_{2L}=\gamma_{L}$, the 1D chiral Lindblad superoperator is

[TABLE]

where $k_{L,R}=\omega_{0}/v_{gL,R}$, with $v_{g}$ being the group velocity of the guided photons.

If we assume now a plasmonic environment, the total emission of the source can be divided into several decay channels: $\Gamma_{11}=\Gamma_{\mathrm{r}}+\Gamma_{\mathrm{nr}}+\Gamma_{\mathrm{SPP}}$, where $\Gamma_{\mathrm{r}}$ represents free-space radiation, $\Gamma_{\mathrm{nr}}$ represents losses in the material (quenching), and $\Gamma_{\mathrm{SPP}}$ represents excitation of SPPs. Material absorption and radiation do not contribute to strong qubit-qubit interactions, and therefore we are interested in systems with strong decay through the plasmon channel, $\Gamma_{\mathrm{SPP}}$, where the fraction of all emissions that are coupled to plasmons is expressed by $\beta_{ij}=\Gamma_{ij,\mathrm{SPP}}/\Gamma_{11}$, with $i\neq j$.

Assuming a plasmonic environment with a preferred propagation axis, here taken as $x$, in order to connect our formulation with previous 1D chiral formulations Garcia2 ; Zoller2 we introduce a particular 1D plasmonic version of (II.1),

[TABLE]

where ${\mathrm{k}}_{\mathrm{ij}}={\mathrm{k}}_{\text{{spp},}\mathrm{ij}}=\mathrm{k}_{\mathrm{spp,ij}}^{\prime}+i\mathrm{k}_{\mathrm{spp,ij}}^{\prime\prime}$ are the SPP wavenumbers. In the systems considered here the bulk modes are reciprocal, whereas the interface SPP is strongly nonreciprocal (unidirectional). Thus, to compare with the 1D chiral ME it is sensible to consider the SPP (nonreciprocal) contribution.

As defined in (Appendix II: Comparison with previous 1D chiral theory), $\Gamma_{ij,\mathrm{SPP}}$ is discontinuous at $x_{i}=x_{j}$ in the nonreciprocal case, i.e., $\Gamma_{ij,\mathrm{SPP}}=2\beta_{ij}\Gamma_{11}$ as $|x_{i}\rightarrow x_{j}|$, whereas at $x_{i}=x_{j}$, $\Gamma_{ij,\mathrm{SPP}}=(\beta_{12}+\beta_{21})\Gamma_{11}$. As we show below, the SPP contribution in the considered PTI system is indeed discontinuous at $x_{i}=x_{j}$. However, the exact $\Gamma_{ij}$, which contains both the SPP and radiation continuum, is continuous at the source point even in the nonreciprocal case. As another example of this, a 3D analytical Green function for a nonreciprocal bulk medium is provided in AL (see their Eq. (117)), where $\Gamma_{ij}$ is also seen to be continuous.

Equating (II.1) in the 1D case (i.e., using (Appendix II: Comparison with previous 1D chiral theory)) and (Appendix II: Comparison with previous 1D chiral theory) term by term, the two Lindblad superoperators will be equal if

[TABLE]

If we now make the assignments

[TABLE]

then (35)-(37) are satisfied and (II.1) becomes strictly equal to (Appendix II: Comparison with previous 1D chiral theory). It is worth stressing that physically the two formulations still differ, since (39) is not exact (phase velocity and group velocity are different quantities). Nonetheless it is interesting to try to connect the phenomenological parameters in the model (Appendix II: Comparison with previous 1D chiral theory) to the corresponding ones in (II.1), which are obtained in terms of the Green function, and hence can be computed for arbitrary environments.

Using the rates defined in (Appendix II: Comparison with previous 1D chiral theory), (II.2) reduces to

[TABLE]

which is distance-independent in the lossless case, as noted in Garcia2 (using (38), (40) is the same as Eq. (6) in Garcia2 ).

Discontinuity of the SPP Here we show that for the strongly nonreciprocal (unidirectional) case, and for a general nonreciprocal case, near the source point the SPP contribution to the Green function is discontinuous. We also show that for nonreciprocal systems, $\Gamma_{21}>\Gamma_{11}$ can occur.

To avoid analytical complications of the general 3D case, we first assume a simple 2D model of a $z$-directed and $z$-invariant magnetic current source located at $x=0,y=d$ inside a biased plasma half-space, adjacent to an opaque half-space occupying $y<0$, as depicted in Fig. 12a. The resulting magnetic field in the plasma is bulkedge ; SR

[TABLE]

where $\mathrm{A}_{0}=i\omega\varepsilon_{0}\varepsilon_{\text{eff}}\mathrm{I}_{m}$, with $\mathrm{I}_{m}$ the magnetic current (set to unity) and $\mathrm{R_{0}}$ accounts for the interface,

[TABLE]

where $\gamma_{p}=\sqrt{\mathrm{k}_{x}^{2}-\varepsilon_{\text{eff}}\mathrm{k}_{0}^{2}}$, $\gamma_{m}=\sqrt{\mathrm{k}_{x}^{2}-\varepsilon_{m}\mathrm{k}_{0}^{2}}$ and $\varepsilon_{m}$ is the permittivity of the metal (opaque medium). The field in the absence of the interface is

[TABLE]

where $\mathrm{H}_{0}^{(1)}$ is the Hankel function of the first kind and order zero and $\rho=\sqrt{x^{2}+(y-d)^{2}}$. The source-point singularity is contained in $\text{Im}(\mathrm{H}_{0}^{(1)})$, and $\Gamma\sim\text{Im}(\textrm{G}_{yy})\sim\text{Re}(\textrm{E}_{y})\sim\text{Re}(\textrm{H}_{z})$.

The interface reflection coefficient $\mathrm{R}_{0}$, leading to the scattered field, contains pole singularities at the SPP wavenumbers (e.g., the denominator of $\textrm{R}_{0}$ is the SPP dispersion equation). For $|\varepsilon_{m}|\rightarrow\infty$ (perfect conductor), there is one pole at $k_{\mathrm{spp},x}=\pm k_{0}\sqrt{\varepsilon_{11}}$ for $\omega_{c}\gtrless 0$. For $|\varepsilon_{m}|$ finite the dispersion equation must be solved numerically, and the plasma may be strongly nonreciprocal, supporting a unidirectional SPP (operating in the bulk bandgap), nonreciprocal, supporting SPPs traveling in opposite directions with unequal wavenumbers (operating above the bulk bandgap), or, in the unbiased (no bandgap) case, reciprocal.

Complex-plane analysis of the magnetic field leads to its evaluation as the sum of a branch cut integral (continuous spectrum) and a discrete residue (SPP) contribution, the latter being

[TABLE]

where Res*(±)* is the residue of $\mathrm{R}_{0}$ evaluated at $\mathrm{k}_{x}=\mathrm{k}_{x,\text{SPP}}^{\left(\pm\right)}$, and $\gamma_{p}^{\left(\pm\right)}=\sqrt{\left(\mathrm{k}_{x,\text{SPP}}^{\left(\pm\right)}\right)^{2}-\epsilon_{eff}\mathrm{k}_{0}^{2}}$, where $\mathrm{k}_{x,\text{SPP}}^{\left(\pm\right)}$ is the SPP pole for $\mathrm{k}_{x}\gtrless 0$ (forward propagating or backward propagating), and where $\theta(x)$ is the Heaviside step function. In the strongly nonreciprocal (unidirectional) case, only one pole is present, leading to only one term in (44).

Figure 12b shows the magnetic field in the bulk bandgap for $\omega_{c}>0$ obtained by numerical evaluation of the Sommerfeld integral (41), and by assuming only the residue component (44) (since we operate in the bulk bandgap and the gap Chern number is $-1$, then there is one unidirectional SPP). The opaque medium is topologically-trivial, and is an unbiased plasma having $\varepsilon=-2$. As shown in the close-up Fig. 12c, the residue accurately approximates the field except very close to the source, where the real-part of the residue ($\propto\Gamma_{\textrm{SPP}}$) has an unphysical discontinuity, indicated by the two black dots. In this case, the radiation continuum compensates for the discontinuity of the residue, such that the real-part of the full Sommerfeld integral ($\propto\Gamma$), is continuous, and the SPP peak is pushed away from the source point.

As a result of the importance of the radiation continuum near the source, at some points $\mathrm{H}_{z}(x=0)<\mathrm{H}_{z}(x>0)$, so that $\Gamma_{21}$ exceeds $\Gamma_{11}$. Figure 12d shows the unbiased (reciprocal) case for the full Sommerfeld integral, where the field peak occurs at $x=0$ and $\Gamma_{21}<\Gamma_{11}$ at all points. In general, there is a quadrature relationship between the dissipative and coherent rates.

Figure 13 shows the magnetic field at a frequency outside the bandgap, where we have two SPPs propagating in opposite directions with unequal wavenumbers. As with the unidirectional case, the residue shows a discontinuity at the source point.

Considering now the 3D case of an electric dipole source at the interface, Fig. 14 shows the dissipative decay and coherent rates (II.1) along the interface, computed using the finite element method (COMSOL, COMSOL ). In this case, it is impossible to separate the discrete and continuum contributions to the field. It can be seen that, as predicted by the previous analytical 2D model, it occurs that $\Gamma$ is nearly discontinuous at the source point (the discontinuity of the discrete spectrum is softened by the radiation continuum), and that $\Gamma_{21}>\Gamma_{11}$ at some points. The coherent rate becomes unbounded at the source due to the well-known divergence of the real part of the Green function.

Appendix III: Concurrence in the unidirectional case

In this section we derive the concurrence for a unidirectional system.

Suppose that the system of qubits are communicating through a strongly nonreciprocal environment, so that the communication is strictly unidirectional, such as occurs for SPPs at PTI interfaces. Assuming that ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{1},{\boldsymbol{\mathrm{r}}}_{2})$ and ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{2},{\boldsymbol{\mathrm{r}}}_{1})$ are the dyadic Green function propagators along two opposite directions, the unidirectionality assumption leads to, e.g., ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{1},{\boldsymbol{\mathrm{r}}}_{2})=0$ ($\Gamma_{12}=g_{12}=0$) and ${\boldsymbol{\mathrm{G}}}({\boldsymbol{\mathrm{r}}}_{2},{\boldsymbol{\mathrm{r}}}_{1})\neq 0$.

Under this unidirectionality assumption, the 3D Lindblad superoperator (II.1) reduces to

[TABLE]

where it has been assumed that $\Gamma_{11}=\Gamma_{22}$.

Defining the basis

[TABLE]

and considering the system of qubits to be initially prepared in the state $\left|4\right>=\left|e_{1}\right>\otimes\left|g_{2}\right>$, it can be shown that for the non-pumped case the non-zero components of the density matrix in (Appendix III: Concurrence in the unidirectional case) are ($\rho$=$\rho_{s}$)

[TABLE]

where $\gamma=\Gamma_{21}/2+ig_{21}$. For all times the density matrix is block diagonal. Concurrence for arbitrary materials can be calculated as Wootters

[TABLE]

where $u_{i}$ are arranged in descending order of the eigenvalues of the matrix ${\rho}(t){\rho}^{y}(t)$, where ${\rho}^{y}(t)=\sigma_{y}\otimes\sigma_{y}{\rho}^{\star}(t)\sigma_{y}\otimes\sigma_{y}$ is the spin-flip density matrix with $\sigma_{y}$ being the Pauli matrix. We have

[TABLE]

such that $x=|\rho_{34}|^{2}+\rho_{33}\rho_{44}$, $y=2\rho_{34}\rho_{33}$ and $z=2\rho_{43}\rho_{44}$ and

[TABLE]

which leads to (II.2).