Non-Markovian effects in waveguide-mediated entanglement

During the last few years, many efforts have been undertaken to control the interaction between quantum emitters and the electromagnetic (EM) field. One of the principal motivations that drive the interest in this research area lies in quantum information. The development of efficient quantum information devices requires the control of light–matter interaction at a quantum level and, in particular, the precise handling of the entanglement resource. Many promising applications of entangled light–matter states have been proposed, such as quantum teleportation, quantum cryptography or simple quantum logic gates [1–3]. The fundamental system all these applications rely on is the two-qubit ensemble. Specifically, much interest is focused on generating entanglement by means of coupling the qubits to a common EM environment. In order to pursue this objective, the properties of the surrounding EM modes are crucial. As the free-space interaction between emitters is not strong enough for entanglement purposes, additional structures have to be included. These structures modify the density of EM states through the Purcell effect, giving rise to collective phenomena. Two main ways have been proposed to achieve these modifications: optical cavities [4, 5], in which one or more discrete modes are responsible for the dynamics, and guided modes [6–8], in which qubits are coupled to an EM continuum. Hybrid systems combining both cavities and waveguides have also been proposed [9–11]. Much theoretical and experimental work has been carried out to improve light–matter coupling in various systems, such as photonic crystal cavities [12, 13], photonic crystal waveguides [14] and dielectric slot [8] and plasmonic waveguides [15–17]. Despite the fact that cavity systems have been studied in more detail, waveguides are also good candidates for a large-scale implementation of quantum devices. The feasibility of experimental realization and the possibility of generating long-lived entanglement [17–19] strongly support this idea. Moreover, waveguides offer an easy way to couple not only qubits but also photons themselves [20, 21], thus providing an excellent workspace for quantum computation applications.

As commented above, an intense light–matter coupling is extensively sought by researchers. However, the assumption that the entanglement properties will be enhanced for greater qubit–photon coupling is not always correct, as we demonstrate in this work. Previous investigations have reported the possibility of entanglement generation between qubits coupled to plasmonic waveguides [17, 18, 22]. These works made use of the master equation formalism, in which the EM degrees of freedom are traced out. After a so-called Markov approximation [23, 24], the system under study is described by the density matrix of the two-qubit subspace. The dynamics of this density matrix are determined by a master equation through a coherent and a dissipative term. Both contributions contain the effect of the traced EM modes which mediate the interaction. Although very successful in three-dimensional quantum optics calculations, this method has some limitations when the dimensionality of the system is lower. Specifically, if photons are confined in one-dimensional (1D) waveguides they undergo successive reflections between the qubits. This, as we will show, can play a key role in the system dynamics, so that averaging the effect of the EM field in this manner is not always a valid approximation. In this paper, instead of the master equation, we solve the Hamiltonian completely. The real-space formalism we use has been developed by other authors [25], and is widely used in photon scattering problems [22, 26, 27]. The set of eigenstates found in these works has successfully described the system for a large variety of initial states. Nevertheless, we will show that for certain initial conditions some extra eigenstates must be added in order to describe the system dynamics. We present a detailed analysis of the qubit populations and the degree of entanglement for different regimes. Our results show that the dynamics of the system is Markovian for weak waveguide–qubit coupling. In this context, we call the dynamics Markovian when the system is well described by a master equation formalism, being the EM modes traced out. On the other hand, when the coupling is increased a new non-Markovian regime appears, in which the entanglement generation is worsened. The different behavior of the populations in this strongly coupled regime could be used for studying the strong coupling phenomenon in waveguide QED systems.

This paper is organized as follows. In section 2 we solve Schrödinger's equation for the lossless case. This section shows in an easy way how new localized eigenstates appear. In section 3 we diagonalize the Hamiltonian in the presence of losses. Once we have solved Schrödinger's equation, we study the evolution of the qubit populations in section 4. Finally, conclusions are presented in section 5.

The system under study is composed of two quantum emitters coupled to a 1D infinite waveguide (figure 1). The waveguide is set along the x-axis, and supports a continuum of photonic modes. The dispersion relation, ω(k), can be approximated linearly within a certain frequency interval [25], v_g being the group velocity. The emitters are modeled as identical two-level systems with transition frequency Ω, located at positions x₁ = −d/2 and x₂ = d/2. They interact with the waveguide modes with a coupling energy γ. We assume that the only coupling between the two emitters is the one mediated by the waveguide modes, hence neglecting the direct coupling energy. This is a realistic assumption for separations in the wavelength scale or greater, as the direct interaction through free space modes decays rapidly with the separation [28]. Losses in the system are taken into account by coupling the qubits to an infinite EM reservoir with coupling energy Γ. In this preliminary section, however, we will assume the system to be completely isolated (Γ = 0).

Zoom In Zoom Out Reset image size

The Hamiltonian of the lossless system can be expressed in the position basis as [25]

$$\begin{eqnarray} &&\fl H_0 = \Omega(\sigma^{\dagger}_{1} \sigma_{1} + \sigma^{\dagger}_{2}\sigma_{2})-\mathrm{i}v_{\rm g}\int \mathrm{d}x \bigg( c^{\dagger}_{\rm R}(x)\frac{\partial}{\partial x}c_{\rm R}(x) -c^{\dagger}_{\rm L}(x)\frac{\partial}{\partial x}c_{\rm L}(x)\bigg) \nonumber\\ &&+\sum_{j=1,2}V\int \mathrm{d}x \delta(x-x_j) [ c^{\dagger}_{\rm R}(x)\sigma_j + c^{\dagger}_{\rm L}(x)\sigma_j + \mathrm{h.c.}], \end{eqnarray} \tag{ 1 }$$

where we take ℏ = 1. Here, σ_1,2 are the lowering operators for the emitter 1 and 2, respectively. The field operators c^†_R,L(x) create a right (R) or left (L) propagating photon at the position x. In the upper line of equation (1), the energy of the emitters and the waveguide modes is contained in the first and second terms, respectively. The second line contains the interaction potential. Note that the coupling constant, V , does not have dimensions of energy and we have to define the coupling energy [29] as γ = V²/v_g. The fact that the coupling is not defined as an energy arises naturally when coupling a real continuum to a discrete system [30], as can be checked by taking the continuum limit of the Anderson impurity model [31, 32]. Note that the interaction term in equation (1) is point-like in space.

Given the symmetry of the system, the natural basis is

$$\begin{eqnarray} c_{\rm e,o}(x)& =& \frac{1}{\sqrt{2}} \left( c_{\rm R}(x) \pm c_{\rm L}(-x) \right), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \sigma_{\rm e,o} &= & \frac{1}{\sqrt{2}} \left( \sigma_1 \pm \sigma_2 \right), \end{eqnarray} \tag{ 3 }$$

where the indices e,o stand for the even and odd symmetry of the operators. Expressed in this basis, the Hamiltonian in equation (1) splits into uncoupled, independent contributions H = H_even + H_odd which can be solved separately. A compact notation for the Hamiltonian in both even and odd subspaces is the following (for j = even,odd):

$$\begin{eqnarray} H_{\!j} &=& \Omega \sigma^{\dagger}_{j} \sigma_{\!j} -\mathrm{i}v_{\rm g} \int \mathrm{d}x c^{\dagger}_j(x) \partial_x c_j(x) \nonumber\\ & &+ V \int \mathrm{d}x [\delta(x+d/2) +\eta_j \delta(x-d/2)] (c^{\dagger}_j(x)\sigma_j + \sigma^{\dagger}_j c_j(x)), \end{eqnarray} \tag{ 4 }$$

where we have introduced the variable η_e = 1 and η_o = −1.

Throughout this paper, we will only work in the one-excitation subspace. The diagonalization is carried out by assuming the most generic form for a one-excitation Fock state [22] in the subspace j(=even,odd), with energy :

$$\begin{equation} \vert \epsilon_j\rangle = \int \mathrm{d}x \phi_j(x)c^{\dagger}_j(x) \vert 0 \rangle + \alpha_j \sigma^{\dagger}_j \vert 0 \rangle. \end{equation} \tag{ 5 }$$

The qubit coefficient α_j and the photon wavefunction ϕ_j(x) are the unknowns to be determined. By solving the time-independent Schrödinger equation [27] we find they are related through

$$\begin{eqnarray} &&(\epsilon - \Omega)\alpha_{\!j} = V(\phi_j(-d/2) + \eta_j \phi_j(d/2)), \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&(\epsilon +\mathrm{i}v_{\rm g} \partial_x)\phi_j = V\alpha_j \left[\delta(x+d/2) + \eta_j\delta(x-d/2)\right] . \end{eqnarray} \tag{ 7 }$$

Equation (7) shows that the photonic wavefunction must be a piecewise continuous plane wave. Then, without loss of generality we can take the next ansatz for this function:

$$\begin{equation} \phi_{j}(x) = \mathrm{e}^{\mathrm{i} \epsilon x / v_{\rm g}} \left\{ \begin{array}{@{}l c r@{}} A_{\!j}(\epsilon) & {\rm for} & x < -d/2, \\ t_{0,j} (\epsilon) & {\rm for} & -d/2 < x < d/2, \\ t_{1,j}(\epsilon) & {\rm for} & x > d/2. \end{array}\right. \end{equation} \tag{ 8 }$$

By inserting this ansatz into equations (6) and (7) we arrive at a final set of algebraic equations:

$$\begin{eqnarray} (\epsilon - \Omega)\alpha_j &=& \frac{V}{2}\left(\eta_j(t_{1,j}+t_{0,j})\,{\rm e}^{{\rm i}\epsilon d/2v_{\rm g}} + (A_{\!j}+t_{0,j})\, {\rm e}^{-{\rm i}\epsilon d/2v_{\rm g}} \right), \nonumber\\ &&\mathrm{i}v_{\rm g}(t_{1,j}-t_{0,j})\,{\rm e}^{{\rm i}\epsilon d/2v_{\rm g}} = V\alpha_j, \\ &&\mathrm{i}v_{\rm g}(t_{0,j}-A_j)\, {\rm e}^{-{\rm i}\epsilon d/2v_{\rm g}} = \eta_j V\alpha_j.\nonumber \end{eqnarray} \tag{ 9 }$$

In this indeterminate system, one of the coefficients (let us say A_j()) acts as a quantum number to label different sets of states. Hence, we will have the following cases depending on its value.

2.1. Scattering eigenstates, A_j() ≠ 0

If we consider A_j different from zero, we normalize the other three unknowns to A_j to obtain a branch of eigenstates, which we shall call the scattering branch, with the following coefficients:

$$\begin{equation} t_{0,j} = \frac{\epsilon - \Omega}{\epsilon - \Omega + \mathrm{i}\gamma\left(1+\eta_j\,{\rm e}^{{\rm i}\epsilon d / v_{\rm g}}\right)}, \end{equation} \tag{ 10 }$$

$$\begin{equation} t_{1,j} = \frac{\epsilon - \Omega - \mathrm{i}\gamma\left(1+\eta_j\,{\rm e}^{-{\rm i}\epsilon d / v_{\rm g}}\right)}{\epsilon - \Omega + \mathrm{i}\gamma\left(1+\eta_j\,{\rm e}^{{\rm i}\epsilon d / v_{\rm g}}\right)}, \end{equation} \tag{ 11 }$$

$$\begin{equation} \alpha_{\!j} = V \frac{\mathrm{e}^{-\epsilon d /2 v_{\rm g}}+\eta_j\,\mathrm{e}^{\epsilon d /2 v_{\rm g}}}{\epsilon - \Omega + \mathrm{i}\gamma\left(1+\eta_j\,{\rm e}^{{\rm i}\epsilon d / v_{\rm g}}\right)}. \end{equation} \tag{ 12 }$$

An example of a scattering eigenstate is shown in figure 2(a). This branch of eigenstates has been obtained in previous works for one [7, 26] and two [27] emitters, and it is sufficient to completely describe any scattering state. However, for other initial configurations, additional eigenstates are needed.

Zoom In Zoom Out Reset image size

2.2. Localized eigenstate, A_j() = 0

If we set A_j equal to zero, we find that some states can appear under a certain combination of parameters. More precisely, these states only exist when the separation between qubits is equal to half an integer multiple of the emitter characteristic wavelength, i.e. d = nλ/2 where λ = v_g/2πΩ. We will name these the resonant separations (d_res) from now on. As a consequence of this condition, the energy of these states has to be  = Ω in order to have a nontrivial solution. Moreover, the norm has to be finite so the transmission coefficient must be t_1,j() = 0.

Provided that the conditions {d = d_res; = Ω} are satisfied, system (9) can be shown to have an additional solution. Only the coefficients t_0,j and α_j are different from zero, and are given by

$$\begin{equation} t_{0,j} = -\mathrm{i}\eta_j\frac{V \alpha_j}{v_{\rm g}} {\rm e}^{{\rm i}\pi (d/\lambda)}. \end{equation} \tag{ 13 }$$

These states are localized in the region between emitters, as shown in figure 2(b). The appearance of these states can be understood by noticing that, when a photon with energy Ω is scattered by a single emitter, the transmission probability is zero [33]. Hence, the regions outside and inside the two qubits become totally uncoupled, and incident photons with this energy are completely reflected [1, 34]. This is shown in figure 2(c). But the existence of a photon with energy Ω in the inter-qubit region is also possible. In order for this state to be stationary (and thus an eigenstate of the system) the stationary wave condition must be fulfilled for this qubit–qubit cavity. This requirement is obviously satisfied if the separation is half an integer multiple of the wavelength, i.e. d = d_res = nλ/2.

Figure 2(b) shows the photon probability distribution of this new localized eigenstate, which is formed by a qubit part and a resonant photonic contribution that exists only inside the cavity. Note that the parameters for which this state appears (d = d_res,  = Ω) are precisely the zeros in the denominator of the scattering coefficients. Note also that for even resonant separations (d = 2nλ), the localized state only appears in the odd subspace, and for odd resonant separations (d = (2n − 1)λ/2) it appears only in the even subspace.

The appearance of this finite norm state is essential in order to have a complete basis. Previous works [22, 27, 35] have studied the scattering of photons with the two emitter system. These new localized states do not play any role in the dynamics of these scattering problems, as their overlap with any incident wavepacket is zero. This is a direct consequence of the condition A_j() = 0. However, if the initial state contains some emitter contribution, the addition of these states is essential. Note that, when the separation is not resonant, there is no localized state. Hence, for d ≠ d_res the scattering basis is complete to describe the evolution of the system for any initial conditions. On the other hand, when the separation is resonant, the new state appears and has to be taken into account. The importance of the localized state in the time evolution can be quantified in the following way: let us set a resonant separation between qubits, so that a localized state appears in either the even or the odd subspace. We now take an initial state |Φ(t = 0)〉 and define P₀ = 〈Φ(t = 0)|σ^†_j|0〉. This quantity represents the qubit part of the state in the resonant subspace. This overlap can be expressed as P₀ = P_sc + P_loc, in which scattering and localized state contributions appear separately. By performing this separation we may quantify both of them. As the localized eigenstate is orthogonal to any scattering one, the fraction P_loc/P₀ is constant with time. Its square modulus represents the amount of initial qubit population contained in the localized state. It can be shown to be

$$\begin{equation} \bigg \vert \frac{P_{\rm loc}}{P_0} \bigg \vert = \frac{1}{1+2\pi(\gamma/\Omega)(d/\lambda)}. \end{equation} \tag{ 14 }$$

In this simple expression it can be seen that the lower the coupling and the resonant distance, the greater the importance of this localized state. As an example, for high couplings and low resonant separations, e.g. γ = 0.01Ω and d = λ/2, the localized state can carry more than 95% of the probability, thus governing the dynamics in the resonant subspace.

Now that the localized states are well understood, we can go one step further and take losses into account. After solving the more realistic lossy Hamiltonian, we will observe that the effective coupling to additional decay channels induces the localized state to become not a single, discrete state but a whole new continuous branch.

In this section we solve the system with losses by adding to the qubits a finite decay rate into the EM environment. Losses in the waveguide are not considered, although a brief comment on finite photon propagation lengths can be found at the end of this paper. The problem of introducing losses in this system has already been studied [25]. Previous works [27, 36] introduce losses by adding an imaginary part to the qubit frequency, i.e. Ω↦Ω − iΓ/2. This modification, as we will see below, is enough to describe scattering problems. However, for certain initial conditions some extra states are needed, as we have already seen in the lossless situation. Moreover, when adding an imaginary part to the qubit frequency, the Hamiltonian becomes non-Hermitian. Even though it can be formally diagonalized to obtain a modified scattering branch of eigenstates, we are not able to recover the localized states when we take the limit Γ → 0. Additionally, it can be shown that the closure relation formed by the eigenstates of this non-Hermitian Hamiltonian is not complete for particular initial states. This indicates that there must be missing eigenstates. For all these reasons, a more complete description of the losses is necessary.

Our model will be similar to that used in the previous section, the Hamiltonian being

$$\begin{equation} H = H_{0} + H_{{\rm r}} + H_{{\rm c}}, \end{equation} \tag{ 15 }$$

where H₀ is the original lossless Hamiltonian (equation (1)). The second term, H_r, is the contribution of the reservoir modes,

$$\begin{equation} H_{{\rm r}} = -\mathrm{i}v\int \mathrm{d}z P^{\dagger}_1(z) \partial_z P_1(z) -\mathrm{i}v\int \mathrm{d}z P^{\dagger}_2(z) \partial_z P_2(z). \end{equation} \tag{ 16 }$$

Here, we have modeled the free-space EM environment of each emitter as one infinite continuum of photonic modes. Equation (16) is similar to the Hamiltonian of two infinite waveguides (second line in equation (1)). In this case, however, they would be unidirectional waveguides, in the sense that have not two but a single propagation direction. This keeps their main effect on the qubit populations providing an exponential decay, and simplifies the calculations. These 1D reservoirs are independent of each other and are set along the z-axis. They are characterized by the creation operators P^†₁(z) and P^†₂(z) respectively, and a group velocity v which we will consider to be v = v_g for simplicity. The third term in equation (15), H_c, represents the coupling of the reservoir with the qubit states. This contribution is expressed as

$$\begin{equation} H_{{\rm c}} = K\int \mathrm{d}z \delta(z) [ P^{\dagger}_1(z) \sigma_1 + P^{\dagger}_2(z) \sigma_2 + \mathrm{h.c.} ], \end{equation} \tag{ 17 }$$

where K is the coupling strength. The coupling energy is then defined as Γ = K²/v, as discussed in previous section. The interaction is set to take place in z = 0 for both reservoirs.

Once we have introduced a decay rate to the guided modes, γ, and a second one to free-space modes, Γ, we are in a more realistic situation. Here, we can define the usual figures of merit in this kind of problem, i.e. the Purcell factor F_P and the β factor of a single emitter [28]. The first one is defined as the enhancement of the decay rate with respect to the free-space case, F_P = (2γ + Γ)/Γ. The β factor is defined as the amount of radiation coupled to guided modes, i.e.

$$\begin{equation} \beta = \frac{2\gamma}{2\gamma + \Gamma}. \end{equation} \tag{ 18 }$$

The factor 2 appears because, as the waveguide has two propagation directions, the total decay rate to the guided modes is 2γ. Reservoir states, however, are unidirectional as mentioned above, hence Γ being the total decay rate.

In order to diagonalize the Hamiltonian (15), let us express the reservoir operators in the even/odd basis as

$$\begin{equation} P^{\dagger}_{\rm e,o} (z) =\frac{1}{\sqrt{2}} ( P^{\dagger}_1 (z) \pm P^{\dagger}_2 (z) ). \end{equation} \tag{ 19 }$$

Under this change of basis, both H_r and H_c split into independent even and odd contributions:

$$\begin{equation} \big(H_{{\rm r}} + H_{{\rm c}}\big)_j = -\mathrm{i}v\int \mathrm{d}z P^{\dagger}_j(z) \partial_z P_j(z) + K\int \mathrm{d}z \delta(z) [ P^{\dagger}_j(z) \sigma_j + \mathrm{h.c.} ] \end{equation} \tag{ 20 }$$

for j = even,odd. The expression of the one-excitation Fock state is now

$$\begin{equation} \vert \epsilon_j \rangle = \bigg(\int \mathrm{d}x \phi_j(x) c^{\dagger}_j(x) + \int \mathrm{d}z \psi_j(z) P^{\dagger}_j(z) + \alpha_j \sigma^{\dagger}_j \bigg) \vert 0 \rangle, \end{equation} \tag{ 21 }$$

where the ansatz for ϕ_j(x) is the same as in the lossless case (equation (8)) and the ansatz for the photonic wavefunction of the reservoir, ψ_j(z), is

$$\begin{equation} \psi_j(z) = {\rm e}^{{\rm i} \epsilon z /v_{\rm g}} \left\{ \begin{array}{@{}cl l@{}} a_j & {\rm for} & z< 0, \\ b_j & {\rm for} & z>0.\end{array}\right. \end{equation} \tag{ 22 }$$

We can diagonalize now, following the same steps as in the previous section. In this case, the outcome of this diagonalization is a system of four algebraic equations with six unknowns. We choose the independent variables to be A_j() and a_j(), so that the remaining variables are determined if these two are known. It can be demonstrated that no eigenstate exists with both coefficients equal to zero, and hence only two possibilities remain, as described in the following.

3.1. Scattering branch, a_j = 0,A_j ≠ 0

For these values of the parameters, all unknowns can be normalized to A_j, so that we obtain a 4 × 4 algebraic system. The solution for the reservoir can be shown to be

$$\begin{equation} b_j = -\mathrm{i}\frac{K}{v_{\rm g}}\alpha_j. \end{equation} \tag{ 23 }$$

The solutions for the coefficients {t_0,j,t_1,j,α_j} are identical to those of the lossless case (equations (10)–(12)), with an additional imaginary part in the qubit frequency. That is, by performing the substitution Ω↦Ω − iΓ/2, previous works obtain exactly this branch of eigenstates. The scattering branch is thus composed of correct eigenstates. However, as we have already remarked above, the subspace spanned by them is not always enough to describe the state of the system.

3.2. Quasi-localized branch, A_j = 0,a_j ≠ 0

Let us explore the second possible parameter combination and look for the equivalent of the localized state we found in the lossless case. When taking A_j = 0, as we still have an arbitrary parameter, a solution of the Schrödinger equation exists for any energy. That is, there are no constraints over the parameters. Consequently, we will have a second continuum of states instead of the discrete one obtained in the lossless case. The expressions of the coefficients, normalized to a_j, are

$$\begin{equation} t_{0,j} = -\mathrm{i}\frac{\sqrt{\gamma \Gamma} \,{\rm e}^{{\rm i} \epsilon d/2v_{\rm g}} }{\epsilon - \Omega +\mathrm{i}\frac{\Gamma}{2} +\mathrm{i}\gamma (1+\eta_j \,{\rm e}^{{\rm i}\epsilon d/v_{\rm g}})}, \end{equation} \tag{ 24 }$$

$$\begin{equation} t_{1,j} = -\mathrm{i}\frac{\sqrt{\gamma \Gamma} ({\rm e}^{{\rm i} \epsilon d/2v_{\rm g}}+\eta_j \,{\rm e}^{-{\rm i} \epsilon d/2v_{\rm g}}) }{\epsilon - \Omega +\mathrm{i}\frac{\Gamma}{2} +\mathrm{i}\gamma (1+\eta_j \,{\rm e}^{{\rm i}\epsilon d/v_{\rm g}})}, \end{equation} \tag{ 25 }$$

$$\begin{equation} \alpha_j = \frac{K}{\epsilon - \Omega +\mathrm{i}\frac{\Gamma}{2} +\mathrm{i}\gamma (1+\eta_j \,{\rm e}^{{\rm i} \epsilon d/v_{\rm g}})}, \end{equation} \tag{ 26 }$$

$$\begin{equation} b_{\!j} = 1-\mathrm{i}\frac{K}{v}\alpha_j. \end{equation} \tag{ 27 }$$

It can be demonstrated that the subspace spanned by this branch is orthogonal to the scattering subspace. Another important feature is that, as we take the limit Γ → 0 (i.e. K → 0), we have t_1,j → 0, and a factor δ( − Ω)δ(d − d_res) appears in both α and t_0,j. Hence, we recover the single localized state we got in the lossless case. Note that, for any qubit separation, we can have t_1,j = 0 at a given energy (such that the numerator in equation (25) cancels out), even for Γ ≠ 0. In other words, for a given separation d there is always an eigenstate whose wavelength is in resonance with the cavity. That eigenstate exists only in the inter-qubit region, so that it is fully localized and decays only to the reservoir. Note that, as equation (26) shows, the maximum qubit occupation is achieved for  = Ω. Thus, the situation in which the resonant, localized eigenstate is  = Ω will be the most interesting for entanglement purposes. This will be achieved precisely for d = d_res. For off-resonant eigenstates, according to equation (25), a certain amount of probability is leaking outwards. Figure 3 shows an example of these quasi-localized eigenstates.

Zoom In Zoom Out Reset image size

By using the expression of the Fock state (equation (21)), it is also possible to check that the overlap of any eigenstate in the quasi-localized branch with any incident, pure photonic wavepacket in the waveguide is zero. This makes the scattering branch to be complete for any scattering problem. However, if the initial state contains a non-zero qubit component, the localized states play a key role in the dynamics of the system. As we will see below, for resonant qubit separations a single excited qubit decays with a rate proportional to 2γ due to the scattering branch contribution, whereas the localized branch adds a typical decay Γ. The competition between these two times and the photon travel time between the qubits t = d/v_g gives rise to various phenomena.

The time evolution of qubit populations has already been studied in previous works [17, 18]. Specifically, they have focused on the possibility of entanglement generation, that is, the spontaneous evolution of an initially unentangled system into an entangled state. It is particularly interesting to study the evolution of the state |Φ(t = 0)〉 = σ^†₁|0〉, i.e. the first emitter in the excited state. By using a master equation formalism, the above-mentioned works have demonstrated that both symmetric and antisymmetric qubit states decay exponentially. One of them becomes subradiant (decay rate <2γ + Γ) whilst the second one becomes superradiant (decay rate >2γ + Γ). The values of both decay rates depend strongly on the inter-qubit separation. In particular, when this separation is resonant, the decay rate of the subradiant state into the guided modes becomes zero. As a consequence, its population remains almost constant with time, slowly decaying into the reservoirs with decay rate Γ. This method provides an efficient way of generating a maximally entangled state. However, the formalism in sections 2 and 3 allows us to go beyond this master equation approach and explore the validity of the Markov approximation.

In this section, we set the same initial state of the system |Φ(t = 0)〉 = σ^†₁|0〉. As we have a complete set of eigenstates, we can explicitly write the closure relation and obtain the time evolution of this state, |Φ(t)〉, in the usual way:

$$\begin{equation} \vert \Phi (t) \rangle = \sum_{j={\rm e,o}} \frac{1}{2\pi v_{\rm g}}\sum_{m={\rm S,Q}}\int \mathrm{d}\epsilon \,{\rm e}^{-{\rm i}\epsilon t} \vert \epsilon_{j,m} \rangle \langle \epsilon_{j,m} \vert \Phi(0)\rangle, \end{equation} \tag{ 28 }$$

where the index j(=even,odd) sums over both subspaces, and the index m(= scattering, quasi-localized) sums over both branches of eigenstates. We are interested in the time-evolution of the populations of the symmetric and antisymmetric states, defined as

$$\begin{equation} \rho_{\pm\pm}(t) = \langle\pm\vert \Phi (t) \rangle\langle \Phi (t)\vert \pm \rangle \end{equation} \tag{ 29 }$$

as well as in determining the degree of entanglement between the qubits. In order to compare with previous results [18], we will use Wooters concurrence [37] to quantify the entanglement. For the initial state considered, the complete expression for the concurrence is reduced to [38]

$$\begin{equation} C(t) = \frac{1}{2}\sqrt{[\rho_{++}(t) -\rho_{--}(t)]^2+4{\rm Im}\,[\rho_{+-}(t)]^2}, \end{equation} \tag{ 30 }$$

where ρ₊₋(t) = 〈 + |Φ(t)〉〈Φ(t)| − 〉 represents the off-diagonal part of the reduced density matrix in the subspace spanned by the pure qubit states {|g₁g₂〉,| − 〉 = σ^†_o|0〉, | + 〉 = σ^†_e|0〉,|e₁e₂〉}. The concurrence and populations arising from our numerical calculations can be compared with the analytical expression obtained within the Markov approximation. The difference between both results will be a sign of non-Markovian behavior.

Note that the waveguide-mediated interaction will result in interesting dynamics for both populations and concurrence, as photons can undergo successive reflections inside the inter-qubit cavity. On the other hand, as the emitters are not coupled through the reservoirs, these EM modes will not add interesting features apart from the well-known exponential decay. Thus, we will consider only the case Γ < γ, in which the waveguide modes prevail and determine the dynamics. As we are mainly interested in the effects of varying the qubit–waveguide coupling, we will modify the relation γ/Ω, keeping the ratio Γ/γ = 0.1 constant along all this work. This corresponds to a Purcell factor F_P = 21 (β = 0.95). Although larger than the typical values in dielectric waveguides [39, 40], it lies well below the experimentally observed values for photonic crystal waveguides [41, 42] (F_P ≲ 30 in the THz regime) and both plasmonic [19, 43–45] and slot [8] waveguides (F_P ≲ 50). As a final setting, we let the group velocity be v_g = c and study the system for different values of the separation d. This does not imply loss of generality, as a reduction in v_g is equivalent to a corresponding increment in the separation d: both imply an increase in the time spent by photons to cover the distance between emitters. Different regimes appear in the time evolution, depending on the ratio between the coupling γ and the qubit frequency Ω.

Markovian results [17] are recovered in figure 4. Here, both populations and concurrence are displayed when the coupling is set to be very small, γ ≪ Ω. In this situation, the decay time of qubit 1 is much higher than the time spent by photons to reach qubit 2. Hence, we can consider the interaction to be instantaneous and the Markovian approximation is valid. If the qubit separation is resonant (figure 4(a)), the decay of the antisymmetric state | − 〉 into the waveguide modes is totally suppressed and thus losses are the only decay channel. Note that, for other resonant separations (d = λ/2,3λ/2,...), the even modes would be uncoupled instead of the odd ones. This effective uncoupling can be understood by looking at the expression of the transmission coefficient (equation (25)). In this low-coupling case, the factor t₁ is sharply peaked around the qubit frequency, t₁∝( − Ω + iδ)⁻¹, where δ is approximately independent of the energy. This expression is identical to the transmission coefficient for the single emitter case [26]. Consequently, the collective qubit states will behave in a similar manner as single emitters coupled to the waveguide. When the separation is resonant, the coupling for the subradiant state is δ ≈ 0. As a consequence, the collective antisymmetric state becomes decoupled from the waveguide. The superradiant state decays rapidly as its coupling is maximum. On the other hand, when the separation deviates slightly from the resonant value, δ ≠ 0. Then, the symmetric and antisymmetric states behave in the same way: they are both coupled to the waveguide, as no real energy cancels out the denominator of t₁. Hence, they decay exponentially into the waveguide modes as shown in figure 4(b).

Zoom In Zoom Out Reset image size

So far we have recovered the Markovian behavior of the system. However, this collective evolution is modified when the parameter γ is different. Specifically, the system is Markovian if the qubit–qubit interaction can be considered to be instantaneous. That is, if the time spent by virtual photons to cover the inter-qubit separation is negligible compared to the qubit decay time (2γ + Γ)⁻¹. If they are comparable, however, non-Markovian effects appear and the long time entanglement is destroyed. To explore the regions out of the Markovian regime we can either increase the coupling γ or the separation between the emitters. As we have previously mentioned, a longer separation between qubits is equivalent to a decrease of the group velocity v_g. As a consequence, the coupling γ = V²/v_g is increased in an indirect way. These effects are displayed in figure 5. Note that we are considering only resonant separations from now on. When the qubit–waveguide coupling is increased (figure 5(a)), the lifetime of the excited qubit 1 becomes comparable to the time spent by the emitted photons to reach qubit 2. Thus, both time scales play a role in the system dynamics and the Markov approximation starts to lose validity. Regarding the denominator in equations (24)–(26), it cannot be approximated to  − Ω + iδ any more. Now, the effective coupling δ depends on the energy, and thus the time evolution is modified. Specifically, the exponential term in this denominator, which represents the effect of the second qubit, is not a negligible phase and will be turned on at a retarded time. This will modify drastically the system dynamics. In figure 5(a) the coupling still fulfills approximately γ ≪ Ω, and then deviations from the Markovian situation are small. Note that as we increase the coupling, values for the concurrence become smaller than before. This means that, in opposition to what may have been expected, the higher the coupling, the lower the amount of entanglement and its lifetime.

Zoom In Zoom Out Reset image size

A completely non-Markovian behavior can be observed in figure 5(b). Here the separation has been chosen as d = 10λ. As remarked above, this can be seen as an effective increase in the single qubit coupling to the waveguide. We can observe two clearly different regions: for γt < 2π(γ/Ω)(d/λ) ≈ 0.63, the qubit 1 is decaying just as a single emitter coupled to a waveguide [29], i.e. an exponential decay with a rate 2γ + Γ. This is a consequence of the retarded interaction, which is produced by virtual photons. As they have not yet reached the second qubit, the evolution is exactly the same as that of the single emitter. This can be also understood by the denominator in equations (24)–(26): the phase factor $\exp (\mathrm {i}\epsilon d/v_{\mathrm { g}})$ adds poles to our complex plane integral (equation (28)), which represent collective states. However, for short times they are out of our integration contour and do not contribute. Physically this implies that collective states are not taking part in the dynamics as the interaction has not been turned on yet. Hence, both populations ρ₊₊ and ρ₋₋ decay exponentially with exactly the same behavior. The concurrence remains zero because only qubit 1 is taking part in the evolution, so that no entanglement is possible. On the other hand, for γt > 2π(γ/Ω)(d/λ) ≈ 0.63, photons have reached qubit 2, so that it starts becoming involved in the dynamics of the system. Only after this moment is the interaction between qubits turned on and collective effects (and entanglement) appear. In this regime, the even and odd qubit state populations deviate from each other. The symmetric one decays with a higher decay rate as the separation is resonant for the odd modes; the antisymmetric one remains approximately constant, decaying only to the reservoir. In other words, the localized eigenstate  = Ω does not decouple from the waveguide until the effect of the second qubit appears. This collective behavior is reminiscent of that predicted within the Markovian regime. However, the delayed interaction completely changes the final values of the populations. As the symmetric state is depopulated for long times, concurrence follows the antisymmetric population in this limit. In the case of figure 5(b), entanglement generation is clearly worsened with respect to the Markovian case. This, as we have commented, is a direct consequence of the increase in the qubit–waveguide coupling.

Note that, in figure 5(b), subtle traces of a damped oscillatory behavior can be observed in the populations. These oscillations are caused by multiple reflections of the light pulses inside the inter-qubit cavity. Photons corresponding to these pulses are contributions from eigenstates with energy close to resonance, so that they undergo some internal reflections before escaping the inter-qubit cavity. This partial resonance appears because the system is far from the Markovian regime, in which the cavity resonated only for one value of the energy. Here, then, we have a broader resonance spectra for the qubit–qubit cavity. This broadening can be understood by noticing that the total transmittance |t₁|² has approximately a Lorenzian shape, whose width is given by the coupling term. Hence, a higher coupling allows more energies to be close to resonance. For long times, all these quasi-resonant photons abandon the inter-qubit cavity and only the really resonant one ( = Ω) remains. Consequently, the amplitude of these oscillations in the populations decays as we can see in the figure. This oscillatory or transient region is then an intermediate regime between the single qubit decay regime, in which both emitters can be treated independently, and the region of collective evolution, in which the entangled antisymmetric state evolves as a whole entity. In figure 5(b), these oscillations have tiny importance in the global dynamics because the coupling is small enough for the collective state to arise in a short time. Nevertheless, in the strong coupled case they will govern the time evolution, as we can see below.

Let us finally explore the case of strong qubit–waveguide coupling in figure 6. Here the coupling energy is increased to γ = 0.1 Ω. This value lies within the ultra-strong coupling regime of cavity QED [46]. Although values up to γ = 20 Ω have been theoretically predicted in particular systems [46, 47], this coupling is close to the maximum values observed in the experiments [48] (γ ≈ 0.12 Ω). However, dynamics in this regime show interesting features that can be observed if the coupling is increased in other ways (e.g. by increasing the number of emitters [49]). Note that in figure 6 we increase the qubit–qubit separation up to d = 5λ. This allows an even higher effective coupling without breaking the rotating wave approximation (valid for γ ≲ 0.1 Ω). In this situation, qubit 1 is so strongly coupled to the photonic modes that it decays completely before the photons reach the second emitter. Hence, as we see in figure 6, either the first or the second qubit are populated, but not both at the same time. In other words, the populations of symmetric and antisymmetric qubit states are the same, so that entanglement formation is precluded. Revivals or beats in these populations can be observed, due to successive reflections of the light pulse inside the cavity. Note that this photonic wavepacket contains a contribution from all energies in a wide range around the resonant one,  = Ω. It travels inside the cavity for a long time, due to the broader reflectivity spectrum of the qubit. This broadening, as commented above, is a direct consequence of increasing γ.

Zoom In Zoom Out Reset image size

It can be observed that both populations are out of sync for γt ≳ 2π. This is the fingerprint of a change from the single qubit regime to the transient regime commented on above. In this transient region, the difference between symmetric and antisymmetric populations will increase slowly. The cause of this is the variation of the shape of the photonic wavepacket inside the qubit–qubit cavity, as the non-resonant photons are slowly escaping. This process takes a long time, as the energies near Ω are very close to resonance and leak out slowly. At a sufficiently long time, photons of all frequencies will have leaked outside the cavity, except for the pure resonant one,  = Ω. Then, the collective regime is achieved, and the system evolves collectively as seen in previous cases: the superradiant state decays quickly to the waveguide modes and the subradiant one uncouples, decaying only to the reservoir. This collective mode has, however, tiny importance in the dynamics for two main reasons: firstly, as commented in section 2, the importance of the localized states with respect to the scattering ones decreases with γ (equation (14)). Hence, even in the lossless case the maximum value for the concurrence through all of the time evolution would be C_lossless ≃ 0.025 (for the parameters in figure 6). Secondly, note that the higher the coupling, the stronger the effect of the losses. This is a consequence of the fact that the collective regime appears for very long times. Hence, a significant amount of probability has decayed into the reservoirs by this moment. This decreases the maximum value of the concurrence below 0.02, hence reducing even more the importance of the collective effects. The strong qubit–waveguide coupling thus practically suppresses the entanglement formation, the dynamics of the system being mainly described by the individual emitters.

So far we have studied the time evolution of the qubit populations. Now, let us focus on the evolution of the probability that goes out of the qubit states, i.e. the photon-position probability density. As stated before, the formalism used in this work allows us to explore any possible values for the system parameters, even outside the Markovian regime. Another advantage is the possibility of studying the photonic part of the quantum state. Previous studies used the reduced density matrix of the qubit–qubit system, tracing out the EM degrees of freedom and thus losing all information about them. Here, we can explicitly calculate the photon-position probability density, which is equivalent to the emitted photon intensity profile. This is a powerful tool for checking the behavior of the system.

Figure 7 displays the photon probability distribution for three different cases. First one (figure 7(a)) shows the Markovian case, with the same parameters as in figure 4(a). For a better visualization, the position scale in the horizontal axis has been taken to be much larger than the inter-qubit separation d. The emission properties of the system in the Markovian limit are very similar to those of a single emitter. This is due to the de-excitation of the even state into the waveguide, acting as a single collective state. The slower decay into the reservoir states only slightly modifies the details of the curve shape as compared to the single emitter case. This behavior is precisely the one predicted by the Markov approximation: both symmetric and antisymmetric qubit states act as a single emitter, decaying collectively into the waveguide with different decay rates that can be modified through the qubit–qubit separation.

Zoom In Zoom Out Reset image size

In figure 7(b) we move out of the Markovian regime, choosing the same parameters as in figure 5(b). In this case, a significant part of the wavepacket has been emitted into the waveguide modes before qubit 2 starts taking part in the system evolution. After the pulse reaches this second emitter, a typical interference pattern appears inside the qubit–qubit cavity. This irregular pattern lasts for a very short time, while photons escape in successive reflections. The last frame shows the case in which the transient regime is about to end and the collective one has arisen. This state will remain in that position as a stationary wave, decaying slowly into the reservoirs with a time scale t ∼ 1/Γ.

Finally, figure 7(c) shows the photonic probability density for the same parameters as in figure 6, i.e. strong qubit–waveguide coupling. We can observe how the first emitter has fully decayed when the pulse reaches qubit 2. Successive reflections occur when the pulse hits one of the emitters, but no collective state arises in this transient regime. The coupling is so strong that the interaction between a qubit and the guided photon finishes before the pulse reaches the other emitter, so that no interference pattern appears. Collective evolution of the qubits is thus shown to be suppressed. All these frames are consistent with the discussion of figures 4–6. They confirm that, in the single excitation subspace, an increase in the emitter–waveguide coupling suppresses collective effects, thus destroying the entanglement generation.

As a final discussion, let us comment briefly on the feasibility of experimentally observing the described non-Markovian effects. As remarked above, non-Markovian evolution requires an increase in the qubit–waveguide coupling. However, a realistic, achievable coupling [50] is typically not larger than γ/Ω ≈ 0.01. Although in this situation the dynamics deviate significantly from the Markovian predictions, for an easier observation the inter-qubit separation d has to be increased (figures 5(b) and 6). This effectively increases the coupling as stated above, but brings out the problem of losses in the waveguide. For high β factors, even a moderate separation d ≈ 10λ is high when compared to the propagation length for both plasmonic waveguides [16] and photonic crystal waveguides [51]. On the other hand, extremely low losses in dielectric waveguides have been reported [52], but β factors are usually much lower. A feasible observation of non-Markovian effects requires both a high β and a high propagation length. Some systems have been already reported to achieve this goal [22].

We have presented a full quantum electrodynamical solution to the qubit–qubit system coupled to a waveguide in the single excitation case. A complete set of eigenstates has been obtained for the first time, thus allowing us to study the population dynamics for every combination of parameters. Specifically, we have explored the validity of the Markovian approximation. Our results show that the Markovian results are recovered for low qubit–waveguide coupling. However, as this coupling is increased, non-Markovian effects start to appear. These deviations can be extremely large, so that all traces of Markovian dynamics are lost and the evolution of the system is completely different. In this high-coupling regime, both qubits act independently and the entanglement formation is suppressed. This is in opposition to what the Markov approximation suggested, i.e. an improvement in the entanglement properties when the coupling is increased. Although the Markovian regime seems to be a better option for entanglement purposes, the new strongly coupled regime shows very interesting dynamics, as both qubit and photonic degrees of freedom play a key role in the evolution. As a consequence, this regime is a good framework for studying intense light–matter interaction phenomena, such as strong coupling in waveguide QED systems.

This work was funded by the European Research Council (ERC-2011-AdG Proposal no. 290981) and by the Spanish MINECO under contract MAT2011-28581-C02-01.