Non-Markovian effects in waveguide-mediated entanglement

We study the generation and evolution of entanglement between two qubits coupled through one-dimensional waveguide modes. By using a complete quantum electrodynamical formalism we go beyond the Markovian approximation. The diagonalization of the hamiltonian is carried out, and a set of quasi-localized eigenstates is found. We show that when the qubit-waveguide coupling is increased, the Markov approximation is not anymore valid, and the generation of entanglement is worsened.

The paper organization is 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. Figure 1: Two qubits with transition frequency Ω separated by a distance d are positioned over a waveguide. Interaction between the qubits is mediated by the photonic modes in the waveguide.
The system under study is composed by two quantum emitters coupled to a one-dimensional infinite waveguide (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 modelled as identical two level systems with transition frequency Ω, located at positions x 1 = −d/2 and x 2 = 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 electromagnetic reservoir with coupling energy Γ. In this preliminar section, however, we will assume the system to be completely isolated (Γ = 0).
The hamiltonian of the lossless system can be expressed in the position basis as [25] H 0 = Ω(σ † 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 , has not dimensions of energy and we have to define the coupling energy [29] as γ = V 2 /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 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): where we have introduced the variable η e = 1 ; η o = −1.
Along 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 : The qubit coefficient α j and the photon wavefunction φ j (x) are the unknowns to be determined. By solving the time-independent Schrödinger's equation [27] we find they are related through 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: By inserting this Ansatz into equations 6 and 7 we arrive at a final set of algebraic equations: In this indeterminate system, one of the coefficients (let's 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 scattering branch, with the following coefficients: 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 emitters [27], and it is sufficient to completely describe any scattering state. However, for other initial configurations, additional eigenstates are needed. 2.2 Localized Eigenstate, A j ( ) = 0.
If we set A j equal to zero, we find that some states can appear under 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 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 [34,1]. 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 [27,22,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 0 = Φ(t = 0)|σ † j |0 . This quantity represents the qubit part of the state in the resonant subspace. This overlap can be expressed as P 0 = 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 0 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 In this simple expression it can be seen that the lower the coupling and the resonant distance, the greater is the importance of this localized state. As an example, for high couplings and low resonant separations, e.g. γ = 0.01Ω; d = λ/2 the localized state can carry out 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 beyond and take losses into account. After solving the more realistic lossy hamiltonian, we will observe that the effective coupling to additional decay channels produces the localized state to become not a single, discrete state but a whole new continuous branch.

Lossy system.
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 where H 0 is the original lossless hamiltonian (equation 1). The second term, H r , is the contribution of the reservoir modes, Here, we have modelled 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 they do not have 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 † 1 (z) and P † 2 (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 where K is the coupling strength. The coupling energy is then defined as Γ = K 2 /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 problems, 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., 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 hamiltonian 15, let us express the reservoir operators in the even/odd basis as Under this change of basis, both H r and H c split into independent even and odd contributions: for j = even, odd. The expression of the one-excitation Fock state is now 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 We can diagonalize now, following the same steps as in previous section. In this case, the outcome of this diagonalization is a system of 4 algebraic equations with 6 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: 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 The solutions for the coefficients {t 0,j , t 1,j , α j } are identical to those of the lossless case (equations 10,11,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.
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 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. Figure 3: Photon position probability density for the lossy hamiltonian. As the energy = Ω is not in resonance with the inter-qubit separation, photons leak out of the qubit-qubit cavity.
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.

Populations and entanglement generation.
The time evolution of qubit populations has been already studied in previous works [17,18]. Specifically, they have focused in 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) = σ † 1 |0 , i.e., 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. Particularly, 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) = σ † 1 |0 . As we have a complete set of eigenstates, we can write explicitly the closure relation and obtain the time evolution of this state, |Φ(t) , in the usual way: 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 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]: where ρ +− (t) = +|Φ(t) Φ(t)|− represents the off-diagonal part of the reduced density matrix in the subspace spanned by the pure qubit states {|g 1 g 2 , |− = σ † o |0 , |+ = σ † e |0 , |e 1 e 2 }. 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 on 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 [43,44,45,19] and slot waveguides [8] (F P 50).
As a final setting, 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 supressed 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 1 is sharply peaked around the qubit frequency, t 1 ∝ ( − Ω + iδ) −1 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 1 . Hence, they decay exponentially into the waveguide modes as shown in figure 4 (b). 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γ + Γ) −1 . 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 couping γ 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 2 /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,25,26), it cannot be approximated to − Ω + iδ anymore. 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.
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,25,26): the phase factor exp(i d/v 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 the interaction between qubits is 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 Figure 5: (a) As the qubit-waveguide coupling is increased, both populations and concurrence start to be different from the Markovian case. (b) When the separation is greater we go further outside the Markovian region, as the retarded interaction has to be taken into account.
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 changes completely 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 1 | 2 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.5Ω. This value lies within the ultra-strong coupling regime of cavity QED [46]. Athough values up to γ = 20Ω have been theoretically predicted in particular systems [46,47], this coupling is above the values observed in the experiments [48] (γ ≈ 0.12Ω). However, dynamics in this regime shows interesting features that can be observed if the coupling is increased in other ways (e.g. by increasing the number of emitters). 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 succesive 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 γ.
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 above. In this transient region, the difference between Figure 6: High coupling case: time evolution of the populations of the qubit symmetric and antisymmetric states, as well as the concurrence. In this regime, no significant collective effects appear and the dynamics is described by the single emitter contribution. 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 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: first, 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 the time evolution would be C lossless 0.025 (for the parameters in figure 6). Secondly, note that the higher the coupling, the stronger is 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 supresses 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 electromagnetic 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.
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. 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 occcur 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 supressed. All these frames are consistent with the discussion of figures 4, 5 and 6. They confirm that, in the single excitation subspace, an increase in the emitter-waveguide coupling supresses 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 [49] is typically non larger than γ/Ω ≈ 0.01. Although in this situation the dynamics deviates significantly from the Markovian predictions, for an easier observation the inter-qubit separation d has to be increased [ figure 5(b)]. 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 [50]. On the other hand, extremely low losses in dielectric waveguides have been reported [51], 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].

Conclusion
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 supressed. This is in opposition to what Markov approximation suggested, i.e., an improvement of 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.