Chirality-assisted enhancement of tripartite entanglement in waveguide QED

We study the generation and control of genuine tripartite entanglement among quantum emitters (QEs) that are side-coupled to one-dimensional spin-momentum locked (or chiral) waveguides. By applying the machinery of Fock state master equations along with the recently proposed concurrence fill measure of tripartite entanglement [S. Xie and J. H. Eberly, Phys. Rev. Lett. 127, 040403 (2021)], we analyze how three-photon Gaussian wavepackets can distribute entanglement among two and three QEs. We show that with a five times larger waveguide decay rate in the right direction as compared to the left direction, the maximum value of tripartite entanglement can be elevated by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$35\%$$\end{document}35% as compared to the symmetric scenario where both left, and right direction decay rates are equal. Additionally, chirality can maintain the tripartite entanglement for longer than the corresponding symmetric decay rate. Finally, we study the influence of detunings and spontaneous emission on the resulting entanglement. We envision quantum networking and long-distance quantum communication as two main areas of applications of this work.


I. INTRODUCTION
A relatively recent way to accomplish optical nonreciprocity at the quantum level is to utilize the so-called "spin-momentum locking of light" [1,2].To understand this phenomenon, one can consider an ultra-thin (subwavelength) dielectric waveguide or tapered nanofiber guiding an s-polarized electromagnetic wave propagating in the x−direction.Upon reaching the tapered region of the fiber, some of the light field evanescently leaks out from the fiber with an exponentially decaying amplitude in the direction perpendicular to the propagation (say y−direction).Due to the symmetry of the problem, the resultant electric field vector E(x, y) can be expressed in the following fashion E(x, y) = E 0x e ikx e −βy + E 0y e ikx e −βy . ( Here E 0x and E 0y refer to the electric field components in the x and y-direction, respectively, k is the wavenumber, and β represents the inverse decay length.The relation between the electric field components can be established through the use of Maxwell's equations yielding For air-glass interfaces with β ≈ k [3], the longitudinal and transverse components of the electric field turn out to be directly proportional with a phase difference of −π/2, i.e.E 0x ≈ −iE 0y .Thus the resultant electric field carries an electric field component oscillating in the propagation direction consequently breaking the often encountered transverse nature of light propagation.This remarkable feature can then be quantified with the use of a "spin vector" whose magnitude represents the deviation from a perfect circular polarization -a magnitude of one (zero) referring to circular (linear) polarization, while the direction shows the net polarization direction.One of the key consequences of this effect is the flipping of polarization direction of the spin vector (from out of the page to into the page) with the reversal of light propagation direction in the nanofiber (from left-to-right to right-to-left), hence leading to the phenomenon of spin-momentum locking of light [4,5].
One can then imagine the fascinating consequences offered by the spin-momentum locking of light in a variety of light-matter interfaces.To emphasize one such possibility, we consider a two-level QE coupled with such a one-dimensional chiral waveguide.The rate of emission from such a QE in the right (r) and left (l) directions of the waveguide is known to follow where d * is the conjugate of the transition dipole matrix element and E r/l is the electric field vector of the right and left modes.Utilization of a spin-momentum locked waveguide alongwith special type of QEs with polarization-dependent transitions then one can see that only a proper match of the polarizations in both QE and the waveguide will result in the emission/absorption of photon.This polarization-dependent preferential photon absorption/emission has recently opened a new research area called "chiral quantum optics" [6].
In the last five years or so chiral quantum optics, has been achieved in various physical platforms [7][8][9] and the area has witnessed a variety of novel effects [10][11][12][13][14][15].For instance, we and others have studied emitter-emitter entanglement dynamics in chiral and non-chiral wQED [16][17][18].In particular, we have shown that for strongly coupled single and two-photon wQED cases such chiral lightmatter interaction can lead to enhancement in the maximum value of bipartite emitter-emitter entanglement by a factor of 3/2 and 2, respectively as compared to the FIG.1: (Color online) The system studied in this work: A chain of two-level quantum emitters side coupled to a one-dimensional waveguide which is driven by a three-photon wavepacket (represented by the state |3 ω1ω2ω3 ⟩) from the left end of the waveguide.The right end of the waveguide is not driven by any photon source and is, therefore, labeled with a vacuum state |vac⟩.The valleys on the waveguide surface are drawn to indicate the tapered region of the nanofiber where the QEs are trapped to accomplish chiral light-matter interactions.γ j represents the non-waveguide (or spontaneous) emission rate for the jth QE.For further details about the system parameters, see the text below.corresponding non-chiral coupling case [19,20].
In this work, we examine the novel problem of genuine tripartite entanglement generation and control among up to three QEs simultaneously coupled with bidirectional and chiral waveguides.Worthwhile to emphasize here is the fact that the discussion of the tripartite entanglement should not be treated as a straightforward extension of single or two-photon entanglement problems, but by doing so we enter the richer and more challenging domain of multipartite entanglement [21,22] where tripartite entanglement can serve as the simplest case study.As the theoretical tools, we work within the framework of Fock state master equations [23][24][25] and calculate the genuine tripartite entanglement among up to three two-level QEs using concurrence and concurrence fill criteria [26].As some of the main findings, we find that, as compared to the corresponding non-chiral (symmetric bidirectional) models, the chirality (five times larger emission rate into the right direction in the waveguide as compared to the left direction) can raise the maximum tripartite entanglement value by 35% of by a factor of ∼ 5/14.Additionally, for both on-resonant and off-resonant cases, chirality aids to maintain tripartite for a longer duration as compared to the symmetric bidirectional problem.Furthermore, chirality also exhibits better robustness against spontaneous emission compared to non-chiral scenarios.
The rest of the paper is structured as follows.In Sec.II we discuss the theoretical description of our system.In Sec.III, we introduce the entanglement measure and discuss our results.In Sec.IV we close with a summary and point out possible future directions.Finally, in the Appendices we outline the derivation of the three-photon Fock-state master equation for cascaded multi-emitter wQED.

A. Model
As shown in Fig. 1, our system consists of a chain of two-level QEs (qubits, quantum dots, artificial atoms, natural atoms, etc.) side coupled to a bidirectional dispersionless and lossless waveguide (tapered fiber).The free Hamiltonian of the emitter chain is given by where ∆ j = ω egj − ω p − iγ j is the detuning between the transition frequency ω egj of the jth QE and the peak frequency ω p of the three-photon wavepacket.Note that in our model there is no direct coupling (such as dipole-dipole interaction) present among the QEs rather the interaction is mediated through the waveguide field.σj ≡ |g j ⟩ ⟨e j | is the standard lowering operator for the jth QE with |g j ⟩ (|e j ⟩) being the ground (excited) state.
The QE raising and lowering operators follow the standard Ferminonic commutation relation: σi , σ † j = δ ij .Next, we model the waveguide as a collection of two independent multimode quantum harmonic oscillators, one for the left (l) direction, and the other for the right (r) direction.The corresponding photon annihilation operators are labelled as bl (ν) and br (ω) for the νth and ωth mode.These operators follow the typical Bosonic commutation relations: . Thus, the waveguide Hamiltonian Ĥw takes the form In Ĥw , we have considered an infinitely large number of closely spaced waveguide modes such that the integration over all modes is justified.Finally, under the rotating wave approximation, the interaction between the QEs and waveguide field is described in the following Hamiltonian where we have assumed Γ jr (ω) ≈ Γ jr (ω egj ) ≡ Γ jr and Γ jl (ω) ≈ Γ jl (ω egj ) ≡ Γ jl .Note that in this assumption we have not applied the Markov approximation (flat bath spectrum around the system resonance) [27] rather we are considering a highly localized interaction.d j represents the location of the jth emitter with d j+1 − d j = L being the separation between two consecutive QEs (or lattice constant) that correspond to the time delay τ = L/c.The parameter k 0 = ω eg /c is the wavenumber associated with the atomic transition frequency, while c represents the group velocity in the waveguide.The net Hamiltonian of the global system (QEs, waveguide, and their interaction) is given by Ĥ = ĤQE + Ĥw + Ĥint .

B. Driven dissipative dynamics
As shown in Fig. 1 the left end of our wQED setup is driven by a reservoir that initially exists in a three-photon Fock state which is unlike the standard studied scenario where a classical coherent light source drives the system.Keeping in view this important distinction, in Appendices we derive the master equation apt for the present problem and study the driven dissipative dynamics of our wQED setup through such a bi-directional three-photon Fock state master equation which is reported as follows Here we would also like to point out that a similar Fockstate master equation valid for N photons has also been reported in the past (see for example Eq. ( 21) in [24]).
However, there are two main differences between our three-photon Fock-state master equation and the one reported in Ref. [24].One is the absence of the terms in our master equation that are quadratic in g(t) which are known to appear in the case of nonlinear interactions (for instance, in cavity quantum optomechanics [28]) or in the case of adiabatically eliminated multi-level quantum systems [29].Since our problem doesn't address both of these scenarios, therefore, the absence of such quadratic terms in Eq. ( 7) is understandable.The second difference stems from the fact that, unlike the master equation reported in Ref. [24]), our master equation incorporates bidirectional couplings between QEs and photon wavepacket which is suitable to study wQED problems).
The Liouvillian superoperator L appearing in the aforementioned equation set (7) and applied to an operator ρ consists of three parts with Lcs [ρ], Lpd [ρ], and Lcd [ρ] respectively represent the closed system dynamics, pure decay of energy from the system into the environmental degrees of freedom, and cooperative decay due to collective QE effects.These Liouvillian subparts are given by Lcs where 2Γ irl = Γ ir + Γ il .The Kronecker delta functions appearing in the expression of Lcd [ρ] are defined as δ i≷j = 1, ∀ i ≷ j.The parameter D represent the ratio of interemitter separation L and the resonant wavelength λ 0 i.e. λ 0 = 2πc/ω eg .Finally, the explicit form of the various operators appearing in Eq. ( 7) are given by ρ3, here j = 2, 1, 0,, k = 1, 0, Û (t; t 0 ) is the time evolution operator, and ρl (t) is the density operator for the left continuum in the waveguide.|Ψ r ⟩, |Ψ 2 ⟩ and |Ψ 1 ⟩ are the three-, two-and one-photon reservoir states, respectively with |Ψ 0 ⟩ = |vac⟩.Note that in the above-mentioned set of operators, only the diagonal operators can be categorized as physically valid density matrices.The rest of the off-diagonal operators are not density matrices but they do obey a useful property that ρ † j,3 (t) = ρ3,j (t), ρ † j,2 (t) = ρ2,j (t), and ρ † k,1 (t) = ρ1,k (t).

C. Initial conditions
Initially, we consider all QEs to be in their ground state with the right waveguide continuum in a threephoton wavepacket with the joint spectral density function G(ω 1 , ω 2 , ω 3 ) and the left continuum in a vacuum state i.e. the initial pure state |Ψ⟩ of the system and environment takes the form At this stage we keep the form of G(ω 1 , ω 2 , ω 3 ) general, however, the normalization condition on |Ψ⟩ requires any where in arriving at this condition we have assumed the spectral function G(ω 1 , ω 2 , ω 3 ) is symmetric under the exchange of mode frequencies ω 1 , ω 2 , and ω 3 .Finally, we impose ρm,m (0) = ρsys (0) = j |g j ⟩ ⟨g j | , ∀m and n = 0, 1, 2, 3; and ρm,n (0 which are the initial conditions followed by the operators appearing in Eq. 7.

III. RESULTS AND DISCUSSION
In this section, by numerically solving our threephoton bidirectional Fock state master equation we address two questions: • How does the incoming three-photon wavepacket excites the QEs, and as a result how does the population evolve in time?• How does the photon absorption & emission generate entanglement among QEs and how chirality can impact the entanglement manipulation?Albeit Eq. ( 7) is valid for any number of QEs, in the following we focus on situations up to 3 QEs.To set the stage we begin with the simplest possible situation of a single QE.

A. One QE case and population dynamics
For the single-QE case (N = 1) our free QE Hamiltonian reduces to ĤQE = ℏ ∆σ † σ and as initial conditions we assume ρm,m (t 0 ) = |g⟩ ⟨g|, ∀m = 3, 2, 1, 0 and the remaining operators to be zero.For the three-photon spectral density function G(ω 1 , ω 2 , ω 3 ) we assume a factorized form such that using Schmidt decomposition we write where cyc represents the sum over all pairwise cyclic permutation of the indices which counts to a total of 6 terms.We point out that the aforementioned type of decomposition of the spectral density function is experimentally achievable when the three-photon wavepacket is generated by combining the single photons emitted by three independent sources [23,30].Moving forward, in all plots to follow we select a real-valued Gaussian temporal profile for each g function i.e.
Here µ and t represent the standard deviation and mean of the Gaussian function, respectively.In Fig. 2 we plot the population dynamics under strong drive condition i.e. |Ω (max) The rest of the parameters mentioned in the plot caption are selected to generate higher excitation probabilities [24,31].The green dotted dashed curve shows our threephoton normalized Gaussian wavepacket peaked at t = t = 5Γ −1 .We have plotted the ground (P g ) and excited population (P e ) for two cases, namely, a non-chiral or symmetric bidirectional coupling (Γ 1r = Γ 1l ) case (thin blue solid and dashed curves); and a chiral case (thick red solid and dashed curves) in which emission in the right direction is five-time larger than the left direction (Γ 1r = 5Γ 1l ).In both cases, we note that as the Gaussian wavepacket begins to interact with the QE, it took almost t ∼ Γ −1 time before the populations begin to change.
The maximum value of the excited state probability P (max) e attained for the bidirectional case turns out to be 0.521 at t = 5.25Γ −1 which is smaller than the reported value of 0.801 [24] for the single photon problem due to the involvement of bidirectional decays in our model.Additionally, the shape of P e follows the profile of Gaussian input which decays as the photon wavepacket leaves the QE region.The chiral case, on the other hand, allowed to attain a smaller value of P (max) e = 0.373 due to a higher decay rate into the right waveguide direction.Additionally, this maximum value is achieved at a time t = 4.85Γ −1 slightly before the P e reaches its maximum value for the bidirectional case.More importantly, we observe the formation of a side shoulder around t ∼ 5.5Γ −1 .Such behavior of P e in the chiral case is known for the single and two-photon wQED problems [19,20] and (as discussed below) will help in better emitter-emitter entanglement generation and control.

B. Two-QE case and bipartite entanglement
We now extend our wQED study to two QEs.In addition to new ways of population distribution, the case of two QEs opens the possibility of generating entanglement between the QEs which we quantify through the well-known concurrence measure [32,33].For two particles, say particle A and particle B, existing in a bipartite pure or mixed state ρAB , Wootter's concurrence C A(B) is defined as where eigenvalues of operator ϱ AB , λ i , ∀i = 1, 2, 3, 4 are written in a descending order.ϱ AB is called the spin-flipped density operator which is related to the system density operator and the Pauli spin-flip operator σy through Here 0 ≤ C A(B) ≤ 1 with C A(B) = 1 refers to a maximally entangled bipartite state (for example, a Bell state [34]) and C A(B) = 0 indicates a fully separable (unentangled) state.
For the present problem we introduce the basis set |4⟩}.Next, subject to the initial condition ρsys (0) = |g 1 g 2 ⟩ ⟨g 1 g 2 |, we numerically solved the three-photon Fock state master equation.Therein, we find that the spin-flip density matrix of the two-QE system takes the following form with 8 out of 16 time-dependent density matrix elements remaining zero for all times Note that we have adopted short notation here in which ρ 1 ≡ ⟨1| ρ3,3 (t) |1⟩, ρ 4 ≡ ⟨1| ρ3,3 (t) |4⟩, ρ 6 ≡ ⟨2| ρ3,3 (t) |2⟩, and ρ 16 ≡ ⟨4| ρ3,3 (t) |4⟩.Diagonalization of ρ 12 (t) yields the following set of eigenvalues Inserting these eigenvalues in Eq. ( 16), one can find the entanglement between QEs.In Fig. 3(c) we plot this bipartite entanglement in both the bidirectional symmetric and chiral cases.In parts (a) and (b) of Fig. 3 the populations corresponding to these two cases have also been plotted to aid the understanding of concurrence behavior.For the bidirectional symmetric case, we notice that the temporal profile of concurrence follows a pattern with two peaks (at t = 4.70Γ −1 and t = 6.65Γ −1 ) separated by a dip (centered at t ∼ 6Γ −1 ) while reaching the maximum value of up to 11%.The first peak is reached just before the three-photon wavepacket reaches its maximum value which can be argued to correspond to the partial formation of a Bell state involving both QEs excited i.e. (|g 1 g 2 ⟩ + |e 1 e 2 ⟩) / √ 2 (as evident from the solid thin magenta curve of P e1e2 in Fig. 3(a)).After that, as the wavepacket begins to leave the emitter region, one of the emitters decays hence now forming a Bell state like (|e 1 g 2 ⟩ + |g 1 e 2 ⟩) / √ 2 which results in the increase in the concurrence again around t = 6.65Γ −1 hence forming the second peak.Again the population plot (Fig. 3(a)) supports this explanation as P e1g2+g1e2 (blue solid thick curve) decays slowly and dominates over P e1e2 curve around t ∼ 7Γ −1 region.
In the chiral case (Γ ir = 5Γ il , i=1,2), we find a marked change in the behavior of population and entanglement dynamics compared to the bidirectional symmetric case.On one hand, in Fig. 3(b) we observe P e1g2+g1e2 (red solid curve) exhibiting a two-peak pattern with a maximum value increase by a factor of almost 2 compared to the symmetric coupling case (blue solid thick curve in Fig. 3(a)).On the other hand, the maximum value of both QEs excited probability P e1e2 (brown solid thin curve in Fig. 3(b)) reduced more than 1/2 compared to the symmetric problem.We find that this single QE excited probability trend extended down to the entanglement behavior as well, where the concurrence in the chiral case (dashed red curve in Fig. 3(c)) showing a single peak pattern but with 5 times higher value achieved for the maximum entanglement.Furthermore, we note that chirality also assisted in sustaining this entanglement for times between 8Γ −1 to 10Γ −1 even after the three-photon wavepacket diminishes.

C. Three-QE case and tripartite entanglement
Moving on to the three-QE mixed states, it turns out that the bipartite concurrence measure doesn't extend down straightforwardly to the tripartite case [35,36].To this end, we apply a recently proposed tripartite entanglement measure by Xie and Eberly [26].This measure is reported to quantify genuine three-party entanglement FIG. 4: (Color online) The concurrence triangle for a tripartite system (composed of QE-1, 2, and 3).Note that the length of each side of the triangle is equal to the square of the concurrence between different possible bipartite pairings.
by analyzing the area of the concurrence triangle (hence the name triangle measure or concurrence fill).The measure itself involves calculating the pairwise concurrence among all three QEs with a bipartite-split between ith qubit (treated as one subsystem) and j, k qubit pair (as the other subsystem) as shown in Fig. 4).For the set of qubits i, j, k; such a "one-to-other" concurrence is known to follow the identity [37] where, for example, C 2 1(23) is calculated using [38].
Here ρ1 represents the reduced density matrix of the first qubit obtained by tracing out the second and third qubit from the full system density matrix ρ123 .Thus, considering C 2 1(23) , C 2 2(31) , and C 2 3 (12) as lengths of the side of a triangle, Xie and Eberly used Heron's expression for the area of such a triangle and arrived at the following formula that describes the triangle measure: , where the prefactor (16/3) 1/4 ensures that F ∆ remains bounded between 0 and 1, again 1 referring to the maximum of genuinely entangled tripartite state (such as W or GHZ state [22]) and 0 indicates a fully unentangled state.Furthermore, consistent with Fig. 4, Q is also called the half-perimeter of the concurrence triangle.
In Fig. 5 we plot population and entanglement dynamics for the three-QEs problem.In Fig. 5(a) and Fig. 5(b) we compare the populations in bidirectional symmetric and chiral scenarios, respectively.With the presence of the third QE, all probabilities including single emitter being excited (P e1g2g3+g1e2g3+g1g2e3 ), double emitter excited (P e1e2g3+e1g2e3+g1e2e3 ) and triple emitter excited (P e1e2e3 ) have been reported.In both bidirectional and chiral scenarios, we note that as the number of excited QEs is increased the corresponding probability shows a considerable reduction.In particular, in the chiral case, P e1e2e3 becomes too tiny such that we have to include it as the inset in Fig. 5(b) where it reaches a maximum value of merely 0.3%.As summarized in Table 1, we find that the maximum value probability of one-(P 1,max ), and two-(P 2,max ) QE excited in the bidirectional model shows a noticeable decrease for N = 3 case as compared to the respective N = 1 and N = 2 cases.However, in the chiral case, such a trend is broken.Additionally, by the comparison of Fig. 3(b) and Fig. 4(b), we notice that unlike N = 2 problem with chiral couplings, N = 3 chiral scenario fails to show any oscillatory behavior in the populations.But single excitation probability P e1g2g3+g1e2g3+g1g2e3 forms an almost plateau between 5Γ −1 ≲ t ≲ 6.5Γ −1 which helps P e1g2g3+g1e2g3+g1g2e3 to maintain a non-zero value for an additional t ∼ = 1.5Γ −1 after the complete diminishing of the three-photon pulse.
In Fig. 5(c) we plot the time evolution of concurrence fill while varying the right direction emitter-waveguide coupling Γ r (assumed to be the same for all QEs) from symmetric bidirectional case Γ r = Γ l to the maximum chiral case Γ r = 5Γ l .We notice, following the population trend observed in Fig. 5(a) and Fig. 5(b), for all nonchiral cases the entanglement among QEs survives for a time longer than the pulse duration.Additionally, the irregular oscillations in F ∆ (t) for chiral case exhibit the phenomenon of entanglement collapse and revival [39][40][41] which is more visible for the Γ r = 3Γ l case (thin blue curve).Most importantly, we notice that the maximum value achieved by the entanglement in all chiral cases poses an upper bound on the maximum value of entanglement achieved in the symmetric directional case where F ∆ ∼ = 0.70.This important finding is further emphasized in the inset plot in Fig. 5(c) where we observe this max- Entanglement dynamics in the presence of spontaneous emission rate γ which is assumed to be the same for all QEs with a value of 3Γ/4.Insets in both plots show the maximum entanglement as a function of Γ r .Besides detuning and spontaneous emission rate, all parameters are the same as used previously.imum value to be elevated by more than 35% as we go from the symmetric bidirectional case of Γ r to chiral cases of 3Γ l ≤ Γ r ≤ 5Γ l .Note that for single-photon two-qubit wQED problem, Ballestero et al. have shown that the chirality can be used to enhance the maximum entanglement by a factor of 3/2 as compared to the corresponding symmetric bidirectional case [16].Similarly, Mirza et al.
(the corresponding author of this work) have reported the twice enhancement in qubit-qubit entanglement for the two-photon two-qubit case [20].We on the other, in this work have shown that this trend extends down to genuine tripartite entanglement where Γ r ≥ 3Γ l case chirality assists to increase the concurrence fill among three-QEs by 35% (factor of ∼ 5/14).

D. Tripartite entanglement in the presence of detuning and spontaneous emission
So far we have assumed an on-resonance scenario where the peak frequency of the three-photon wavepacket ω p has been set equal to the emitter transition frequency ω eg .Additionally, we have completely ignored the photon emissions into non-waveguide modes through the process of spontaneous emission.We now address these two scenarios separately and plot the three-QE entanglement dynamics for a detuned case with no spontaneous emission (i.e.ω p − ω eg = Γ/2 and γ = 0) in Fig. 6(a) and for an on-resonance case with a non-zero spontaneous emission scenario (ω p = ω eg and γ = 3Γ/4) in Fig. 6(b).
From Fig. 6(a) we note that for all cases as we increase Γ r value from Γ l to 5Γ l , near the peak frequency of the wavepacket, detuning preserves the overall profile of the entanglement observed in the on-resonance situation.Additionally, from the inset plot, we notice that the maximum entanglement values also follow a quite similar pattern as found in the no-detuning problem.However, we observe the novel aspect of Fig. 6(a) in a long time (t ≳ 8Γ −1 ) behavior of F ∆ (t) where tripartite entanglement sustains for longer times and tend to produce more oscillatory behavior as compared to the no-detuning problem (compare, for instance, thin black (Γ r = 3Γ l ) curves in Fig. 6(a) and Fig. 5(c)).
In Fig. 6(b) we study the impact of spontaneous emission on the tripartite entanglement under the strong coupling regime of wQED (γ < Γ).As expected, we find that the presence of a finite spontaneous emission considerably reduced the entanglement while keeping the overall profile of entanglement more or less the same.In particular, we point out that for γ = 3Γ/4, the maximum value of entanglement for the symmetric bidirectional case shows a 15% reduction compared to the γ = 0 situation.Here we emphasize that the chirality not only assists to achieve elevated values of maximum entanglement in the presence of spontaneous emission but also helps to somewhat decrease the difference in the F ∆,max value (see for example, the most chiral situation of Γ r = 5Γ l in which the maximum entanglement difference reduces to 10% compared to the corresponding γ = 0 problem).

IV. SUMMARY AND CONCLUSIONS
In this paper, we studied the generation and control of three-photon Gaussian wavepacket-induced entanglement between 2 to 3 QEs side-coupled to chiral and symmetric bidirectional waveguides.Through the numerical solution of three-photon Fock state master equations, we calculated population dynamics and entanglement evolution which were quantified via bipartite concurrence and concurrence fill for two-and three-QE, cases respectively.
At the single QE level, we found that chiral lightmatter interaction was able to achieve ∼ 37% maximum excitation percentage probability which is smaller than ∼ 52% percentage probability obtained for the bidirectional symmetric coupling case.However, starting from 2 QE case chirality began to exhibit considerable improvement in both gaining higher entanglement values as well as single-QE excitation probability.Particularly, for the set of parameters chosen in Fig. 3, we concluded that with a five times higher decay rate in the right waveguide direction compared to the left direction, a single QE excitation probability reaches values twice higher than the corresponding bidirectional symmetric cases.More importantly, this trend extends down to the emitter-emitter entanglement where the bipartite concurrence reached maximum values that were five times larger than the symmetric case.
For the N = 3 QE problem, we found that for the bidirectional case, the maximum probability of single, double, and triply excited states show a considerable reduction as compared to the corresponding N = 1 and N = 2 problems.However, the chirality breaks this trend and also helps to sustain (at least) the single excitation probability (and hence the entanglement) for longer times.Furthermore, in the chiral case, we notice the phenomenon of tripartite entanglement death and revival.Importantly we point out that the maximum value achieved by the entanglement in all chiral cases (starting from Γ r = 2Γ l to Γ r = 5Γ l ) posed an upper bound on the maximum value of entanglement attained in the symmetric bidirectional problem (Γ r = Γ l ).Compared to earlier studies of one and two-photon wQED where for two-qubit problem chirality is known to increase entanglement by a factor of 3/2 and 2, respectively; here for the three-photon case we have shown this enhancement to be 35% (or by a factor of ∼ 5/14).
Finally, we discuss the impact of detuning and spontaneous emission on the generated tripartite entanglement.There we concluded both small detunings (ω p − ω eg = Γ/2) and spontaneous emission rate (γ = 3Γ/4) retain the overall temporal profile of the entanglement.Detuning helps to sustain entanglement for longer times, while spontaneous emission rate results in a considerable reduction in the maximum value of entanglement.However, chirality still helped entanglement to show somewhat robustness against spontaneous emission loss.These behaviors convincingly show that three-photon chiral lightmatter interactions can assist to accomplish higher maximum values of entanglement among QEs with better control to sustain genuine tripartite entanglement for elongated times.

FIG. 2 :
FIG.2:(Color online) Population dynamics, quantified in the units of Γ −1 , for a single (N = 1) two-level QE when interacted with a three-photon Gaussian wavepacket.We have considered the following common parameters in all curves: ∆ = 0, µ = 1.46Γ, t = 5Γ −1 .For the chiral case, we have set Γ 1r /Γ 1l = 5, while for the bidirectional case, we have selected a symmetric case i.e.Γ 1r = Γ 1l ≡ Γ.The orange dotted horizontal line confirms normalization in both bidirectional and chiral cases.

FIG. 3 :
FIG. 3: (Color online) Two-emitter N = 2 wQED driven by a three-photon Gaussian wavepacket.Population dynamics in (a) Bidirectional case Γ ir = Γ il = 1 and (b) Chiral case Γ ir /Γ il = 5, ∀i = 1, 2. In the subscripts of P the first and second slots identify the first and second QE, respectively.(c) Entanglement/concurrence evolution in both bidirectional and chiral cases.The location and the maximum value of entanglement have been identified by pink and brown-colored dots.For the sake of simplicity, all QEs are assumed to be identical and time delays have been ignored.The rest of the parameters are the same as used in Fig. 2.

FIG. 5 :
FIG. 5: (Color online) Population dynamics for the three-photon three-QE (N = 3) problem.(a) Symmetric bidirectional case i.e.Γ ir = Γ il = 1; and (b) Chiral case with Γ ir /Γ il = 5, ∀i = 1, 2, 3. Similar to Fig. 3, all QEs are assumed to be identical and time delays have been ignored.In the plot legends, we are following the notation in which the first, second, and third slots in the subscripts represent the state of the first, second, and third QE, respectively.The inset in the plot (b) represents the curve of all three QEs being excited simultaneously (P e1e2e3 ).(c) Time evolution of tripartite entanglement among three QEs quantified through the concurrence fill F ∆ (t) measure.Emitter-waveguide coupling strength in the right direction Γ r has been varied in units of Γ l .The inset shows the behavior of maximum entanglement (F ∆,max ) achieved for each chosen value of Γ r /Γ l .The rest of the parameters in all plots are the same as used in Fig. 2.