Nonlinear quantum logic with colliding graphene plasmons

Graphene has emerged as a promising platform to bring nonlinear quantum optics to the nanoscale, where a large intrinsic optical nonlinearity enables long-lived and actively tunable plasmon polaritons to strongly interact. Here we theoretically study the collision between two counter-propagating plasmons in a graphene nanoribbon, where transversal subwavelength confinement endows propagating plasmons with %large effective masses a flat band dispersion that enhances their interaction. This scenario presents interesting possibilities towards the implementation of multi-mode polaritonic gates that circumvent limitations imposed by the Shapiro no-go theorem for photonic gates in nonlinear optical fibers. As a paradigmatic example we demonstrate the feasibility of a high fidelity conditional Pi phase shift (CZ), where the gate performance is fundamentally limited only by the single-plasmon lifetime. These results open new exciting avenues towards quantum information and many-body applications with strongly-interacting polaritons.

The integration of a photonic gate within an optical waveguide constitutes a long-standing challenge in quantum optics that is-as elucidated by the Shapiro theorem-tantamount to overcoming practical limitations associated with the large entanglement spread in momentum space produced by multi-photon scattering [1][2][3][4]. Because these restrictions usually apply to copropagating particles governed by a linear dispersion relation and interacting via a local Kerr-like nonlinearity, they are typically circumvented by invoking nonlocal interactions, which specifically have been exploited in the context of electromagnetically-induced transparency (EIT) with Rydberg atoms [5][6][7][8], as well as in engineered discrete networks of cavities [9][10][11][12][13] and chiral waveguides [14]. However, similar strategies have yet to be identified in an integrated photonic platform that can bring quantum optical logic to the nanoscale.
Plasmon polaritons-quasiparticles that emerge when light hybridizes with collective oscillations of conduction electrons at a metal-dielectric interface-exhibit enhanced dispersion, extending well beyond the light line; the large wave vectors attained by plasmons correspond to an intense concentration of electromagnetic energy on * giuseppe.calajo@pd.infn.it † cox@mci.sdu.dk length scales far below the wavelength of the light that excites them [15][16][17]. Metal nanostructures supporting plasmon resonances have thus been actively explored to enhance nonlinear light-matter interactions on nanometer length scales [18][19][20]. However, both low intrinsic optical nonlinearity and high Ohmic losses encountered in conventional plasmonic materials limit practical implementation of single-photon-level nonlinearity [21]. Graphene has recently emerged as a promising material platform for both plasmonics and nonlinear optics: the long-lived [22][23][24][25] and electrically-tunable plasmons supported by the atomically-thin carbon layer can intensify optical near-fields that drive its relatively large intrinsic optical nonlinearity [26][27][28][29], where the latter attribute stems from a linear electronic dispersion relation that renders charge carrier motion anharmonic [30,31]. The excitement surrounding the appealing nonlinear and optoelectronic properties of graphene has naturally stimulated efforts to trigger quantum nonlinear optical processes based on plasmon polaritons in integrated nanophotonic platforms [26,28,32,33], including plasmon gates [34] and entangled plasmon pair generation via spontaneous parametric down-conversion [35].
Here we propose to exploit the strong optical nonlinearity, enhanced by the flat dispersion of guided plasmons in graphene nanoribbons, to collide counter-propagating polaritons and effectively implement an integrated CZ gate within a plasmonic waveguide. Our proposal re-arXiv:2207.05122v2 [quant-ph] 18 Mar 2023 (a) Two counter-propagating single-plasmon pulses strongly interact in a graphene nanoribbon via a twoplasmon absorption process to acquire a relative π phase after the scattering event. The process is mapped in a relative coordinate frame to the simple problem of a massive particle scattered by a delta potential. (b) Sketch of the electronic band structure for doped graphene; when ̵ hω < 2E F − ̵ hv F k, single-plasmon absorption via electron-hole pair excitation is suppressed, while two-plasmon absorption can occur via a nonlinear interband transition.
lies crucially on the relatively large propensity for twoplasmon absorption in graphene [31,33,34], which we show here to manifest in the interaction of propagating polaritons as a reflective potential. Such a scenario, in a more generic context, is closely related to a Tonks-Girardeau gas [36][37][38][39], which represents the strongly-interacting limit of the well-known Lieb and Liniger model for massive bosons with contact interaction [40]. The CZ gate comprised of colliding plasmons in a graphene nanoribbon is practically limited by the intrinsic single-plasmon absorption rate that is commonly quantified by the quality factor of the associated optical resonance; as we show here, the proposed CZ gate is robust under realistic values commensurate with experiments in graphene plasmonics.

I. MODEL
We envision a system where incoming single-photon pulses of frequency ω, propagating in a low-loss photonic waveguide, are converted to plasmons via injection into a nonlinear (active) region comprised of a graphene nanoribbon in the R = (x, y) plane with length L and width W ≪ L, as depicted in Fig. 1(a). Assuming W to be much smaller than the incident light wavelength, we describe highly-confined plasmons in the quasistatic limit using a scalar potential φ k (x)e iky , which is decomposed in the longitudinal wavevector k to exploit translational invariance in the y-direction. Then, defining the graphene nanoribbon by a two-dimensional conductivity is the local linear conductivity of extended graphene and f R is a geometrical parameter that is 1 within the ribbon structure and zero everywhere else [41,42], we express the self-consistent potential as φ k (θ) = η contains the dependence on size and intrinsic conductiv-ity, θ ≡ x W is a normalized coordinate, and is an integrodifferential operator given in terms of the modified Bessel function K 0 and normalized wavevector q ≡ kW . Following the prescription of Ref. [42], we express M in a discretized real-space basis to extract the eigenvalues η n,k and eigenvectors φ n,k (θ) that define polaritonic modes supported by the nanoribbon geometry. Transverse confinement provided by the ribbon along the x-direction leads to distinct branches in the polariton dispersion relation, which we explicitly compute by invoking the linear conductivity of graphene described in the local limit of the random phase approximation (LRPA) [26], where e is the elementary charge and E F the Fermi energy. The first term in the conductivity accounts for intraband motion of free charge carriers offset by a phenomenological damping rate γ D , while the second term describes interband transitions between the Dirac cones shown in Fig. 1(b). The resulting dispersion for the first three modes is shown in Fig. 2(a). Importantly, the strong light-matter hybridization associated with plasmon polaritons pushes their dispersion well beyond the light line of free-space photons, endowing the propagating quasiparticles with an effective mass m = ̵ h(∂ 2 k ω n,k k=kp ) −1 and a slow group velocity v g = ∂ k ω n,k k=kp that are characterized by expanding the plasmon dispersion relation around a given plasmon resonance frequency ω p ≡ ω n,kp at k = ±k p according to small wave vectors, we find that the higher-order polariton modes corresponding to n > 1 are well-approximated by truncating beyond the quadratic terms in the above expansion, as revealed in Fig. 2(b) by the comparison with the numerically-extracted group velocity of the n = 2 and n = 3 modes (black dashed curves) for values near kW ≃ 1.
The LRPA conductivity in Eq. (2) that is used to obtain the dispersion relation predicts a sign change in the imaginary part of the conductivity at the frequency ω plasma ≃ 5E F 3, which can be interpreted as the plasma frequency beyond which doped graphene ceases to exhibit metallic behavior. Such a feature significantly deviates the plasmon dispersion near and above the Fermi energy E F from that predicted by the purely-intraband Drude model conductivity (see App. A 2 for details), becoming particularly important for higher-order polariton modes. Here the bands are flattened, leading to large effective polariton masses that can exceed the electron mass of 10 −31 kg. Incidentally, the engineering of flat bands in graphene nanoribbons presents an alternative to utilizing a Bragg grating for engineering slowly-propagating polariton modes [43].
Plasmon propagation along the graphene nanoribbon can be described in a second-quantization formalism for massive particles by the effective Hamiltonian whereâ R(L) (y) andâ † R(L) (y) are the bosonic field operators respectively annihilating and creating a rightpropagating (left-propagating) plasmon at position y and v g = v g − ( ̵ h m)k p . Graphene plasmons are understood as well-defined and long-lived excitations only in the absence of incoherent scattering processes involving phonons or defects, and can be further damped by electron-hole pair excitation channels [23]. Fortunately, owing to the Pauli principle, the latter process is suppressed in highly-doped graphene, where absorption via intraband and interband transitions is prohibited for plasmon energies within ̵ hv F k < ̵ hω p < 2E F − ̵ hv F k, as illustrated in Fig. 1(b). The remaining decoherence mechanisms, mainly related to phonon and defect scattering [23], are incorporated in the single-plasmon absorption rate γ 1 = ω p Q, which we characterize by the quality factor Q. The incoherent scattering processes captured in γ 1 limit plasmon propagation by imposing a decay e −2γ1τ , where τ = L v g represents the time it takes to propagate over a distance L; as will be discussed later, singleplasmon absorption mainly affects the free evolution of plasmons and not their interaction dynamics.
Beyond single-plasmon absorption, plasmon pairs can be efficiently absorbed within a certain frequency range via interband transitions [33,34], as illustrated in Fig.  1(b). The associated nonlinear process is encoded in the real part of the third order conductivity σ (3) ω , which admits analytical expression in the local limit of the random phase approximation [31]. In a one-dimensional ribbon, this process is effectively captured by a per-length twoplasmon absorption (TPA) rate that we estimate as is a factor depending on the integral of the electric field mode functions u kp (x) (see App. A for further details). Crucially, the rate of TPA in graphene can exceed the single-plasmon decay rate, leading to strong nonlinear effects [33,34]. We model the nonlinear absorption-induced interaction as a dissipative local Kerr nonlinearity described by the effective Hamiltonian for propagating massive interacting plasmons Note that although the Hamiltonian of Eq. (5) is formally non-Hermitian, it entirely captures (in the absence of external energy sources) the dynamics occurring within a given excitation subspace. In all plots the pulse wavevector was set to kW = 1 and the width to σ = W ∆k with ∆k = 0.9.

II. TWO-PLASMON DYNAMICS
Local Kerr nonlinearities are known to lead to a nogo theorem for the implementation of gates among copropagating photons in nonlinear optical fibers [1,2,4,10]. To circumvent this limitation, we consider the peculiar case of a strong hard core collision between two slow counter-propagating plasmons interacting in the active region via TPA, as depicted in Fig. 1(a). Importantly, we specifically consider the higher (n > 1) dispersion branches, where the flat plasmon dispersion diminishes the contribution of the kinetic term in the Hamiltonian, thus enhancing interactions among plasmons. The collision between two plasmons can be conveniently described in a relative coordinate frame that effectively maps the interaction to the simple problem of a single massive particle scattered by a delta function potential, as sketched in Fig.1(a) (see App. C for details). The single-particle scattering is then described by the scattering matrix S = ∫ dk(r − k⟩⟨k + t k⟩⟨k ) that yields the output state ψ out ⟩ = S ψ in ⟩ corresponding to an incoming input state ψ in ⟩, where are the reflection and the transmission coefficients that depend on the ratio between the plasmon wavelength λ p = 2π k p and the absorption length associated with a two-plasmon absorption-induced interaction. Incidentally, Eq. (7) explicitly shows that the band flattening enhances plasmon interactions, with perfect refection (and thus perfect repulsive collision in the original frame) occurring when λ a = 0 is satisfied by the group velocity v g = ̵ hk p (2m) that corresponds to a completely flat dispersion for the considered terms up to the second order in the original dispersion expansion. Importantly, a reflection from the potential in the singleparticle picture is accompanied by a π-phase shift (minus sign in the reflection coefficient of Eqs. (6)), which corresponds to a relative π-phase shift between the two colliding plasmons. Note that the transmission and reflection probabilities T = t 2 and R = r 2 do not perfectly sum to one since the two-plasmon interaction is inherently dissipative, and lead to a large two plasmon absorption at intermediate λ λ a ratios.

III. CZ PLASMON GATE
As mentioned in the previous section, the key idea underlying our gate proposal is based on the observation that the quasi-flat bands of the higher (n > 1) modes can strongly enhance the plasmon interaction, leading, in the "Zeno"-like limit of λ p λ a ≫ 1 [44,45], to an almost perfect collision that endows the plasmons with a relative (conditional) π-phase shift and thus implementing a CZ gate. An initial figure of merit to estimate the efficiency of this process is given by the reflection probability, R = r 2 , which for a ratio of λ p λ a ∼ 10 3 is on the order of R ∼ 0.99. A more realistic estimate for gate performance can be achieved by considering two counterpropagating Gaussian pulses, where k 0 = k p is the central wavevector and σ the pulse width. A quantitative figure of merit for the efficiency of the gate is then given by the ideal state fidelity F = ⟨ψ in (−k 0 ) ψ out (k 0 )⟩ 2 , which is merely the probability to obtain the same initial pulse reflected up to a phase shift, and can be straightforwardly computed by calculating the output state from the Smatrix. The values of the computed fidelity as function of the nanoribbon width and Fermi energy are shown in Fig. 3(a), which exhibits isolines with similar fidelities corresponding to the same plasmon frequency, e.g., the isoline associated with ̵ hω p = E F is shown by the white curve in the figure.
The computed fidelity evaluates the efficiency of the scattering process without taking into account singleplasmon absorption during propagation of the two pulses. As discussed above, the main sources of single-plasmon absorption are interband electron-hole pair excitation and phonon scattering; the former process is highly suppressed if the plasmon frequencies lie below the interband transition region (above the grey area in Fig. 3), while the latter effect is strongly reduced when the plasmon frequency lies below the optical phonon line at ̵ hω ph < 0.2 eV (below the purple line in Fig. 3) [23]. Importantly, Fig. 3 shows that high fidelities of F > 0.98 can be achieved for plasmon frequencies away from both regions. The average plasmon lifetime is then given by the single-plasmon decay rate γ 1 = ω p Q in terms of a quality factor Q that accounts for all remaining sources of damping, mainly attributable to defect scattering [23]. In this regime, large quality factors Q > 100 are within experimental capability, while even higher values up to Q ∼ 1000 have been theoretically predicted [25].
Having encoded the remaining sources of damping into the quality factor Q, we can define a more accurate estimation of the gate performance in terms of the success probability P p . To estimate this quantity we first observe that, as previously mentioned, the signal detection probability is overall damped by an exponential factor e −2γ1τ , where τ = L v g is the time required for a plasmon to traverse the ribbon. This contribution clearly enhances the probability of plasmon absorption during the gate protocol and penalizes "too slow" group veloci-ties, such as those associated with higher-order plasmon dispersion branches. While such losses are clearly reduced in shorter ribbons, any benefit must be compared against the spatial length ∼σ of the counter-propagating pulses; when σ ≳ L, there is a significant probability that the two excitations never overlap in the ribbon at the same time. On the other hand, one could reduce σ, but this increases the frequency bandwidth, while the large phase shift only occurs within a limited bandwidth. We account for this trade off by estimating the probability of having a single plasmon within the length L, which for a Gaussian pulse reads: , where we included the effect of a relative delay of the two pulses assuming that the scattering event occurs at a distance ∆L from the middle of the ribbon. In this way we define the success probability which is shown in Fig. 3(b) for the second and third mode using an optimal ribbon length, assuming a collision at the center of the ribbon, ∆L = 0, and fixing the quality factor to Q = 1000. Note that for the ribbon lengths considered, L ∼ σ, acquired delays of the order of ∆L L = 0.1 affect the predicted success probability by less than 1% with respect to a perfect central collision. Thus we did not explicitly take into account these deviations in Fig. 3(b). The optimal success probability achievable in each mode (within the low loss region) is plotted in Fig. 3(c) for two different quality factors, and exhibits an opposite trend with respect to the one exhibited in the case of optimal fidelity (blue dots), i.e., becoming more damped for slower, more massive plasmons. Such behavior suggests that even if extremely high fidelities F ≃ 0.99 can be ideally reached by higher-order modes, the second and the third modes are the ones exhibiting optimal operational conditions. Indeed, a convenient trade off between high conditional fidelities and good gate success probability can be achieved for quality factors of Q f = 150 and Q f = 1000, which we predict to be respectively on the order of P succ ∼ 20% and P succ ∼ 50% for the second mode. A broader panoramic over the possible success probabilities achievable with the two considered modes for different quality factors is shown in Fig. 3(d), while the corresponding optimal ribbon width and length are plotted in Fig. 3(e). Fig. 3(d) illustrates how the gate performance progressively approaches the ideal lossless one for large quality factors, Q ≳ 10 4 . Importantly, in the range of quality factors Q ∈ [50, 150] achievable in realistic devices shown in the inset of Fig. 3(d), success probabilities ranging between P succ ∼ 5 − 20% can be reached. This makes our proposal comparable with state-of-theart platforms for implementing photon-photon gates, e.g. cavity QED [46]. These results indicate that overall good gate performance can be obtained in the proposed setup, with single-plasmon absorption representing the primary limitation.
Our discussion up to now does not account for possible in-plane scattering processes induced by defects or disorder in the structure [47], which may eventually lead to Anderson localization [48]. For massive particles, the Anderson localization length roughly scales as [49][50][51], where characterizes the disorder strength. This result shows that the effects of disorder are enhanced by slow-light effects. While we do not know of any straightforward way to determine for a realistic graphene system based on simple physical considerations, we note that the localization length is in principle an experimentally mesurable quantity. For our purposes, it is then clearly important to ensure that the sample is clean enough to have an associated localization length larger than the ribbon active region, i.e., L loc ≳ L.

IV. GENERAL CONSIDERATIONS AND CONCLUSIONS
In summary, we have shown that by engineering graphene nanoribbons it is possible to induce flat dispersion that enhances the plasmon interaction originating from nonlinear absorption. In the limit of strong interactions, the model effectively behaves as a Tonks-Girardeau gas [36,37], such that counter-propagating plasmons undergo elastic collisions that can be exploited to implement an integrated CZ gate. To describe the plasmon collision process, we have introduced a phenomenological Hamiltonian in Eq. (5) that strongly simplifies an otherwise complex many body-problem. In particular, we treat 2D plasmons-collective charge oscillations dressed by the electromagnetic field-as welldefined bosonic excitations interacting locally via TPA at a fixed per-length rate γ 2 . The present simplifying assumption is informed by current knowledge of nonlinear optical processes in graphene, and should be reasonable within a small frequency range, such as that considered in our proposal involving flat plasmon dispersion branches around ω p ≃ 1.5E F . The assumption of local interactions can be further validated by considering that the natural length of the plasmon interaction should be set by the Fermi length L F = 2πv F E F , which, for Fermi energies of E F ≃ 0.1, is on the order of L F ≃ 40 nm and thus commensurate with the considered pulse widths.
Even if theoretically feasible, the current proposal presents important experimental challenges, which are not only restricted by the fabrication of high-quality nanostructures. Indeed, the precise excitation of individual propagating plasmons is still at the edge of the state-of-the art; nevertheless, promising proposals exist to efficiently couple far-field light to propagating plasmons [52], while experiments have demonstrated almost perfect absorption into acoustic graphene plasmons. On the other hand, the occurrence of strong nonlinear multi-plasmon absorption [33] could make the current proposal realizable with weak coherent plasmon pulses. We also remark that, for moderately doped ribbons of width ≳ 10 nm considered in our proposal, the gate performance will be only minimally affected by edge-termination geometry (i.e., armchair or zigzag edges) [53], although ribbons with armchair edges are preferable to avoid additional plasmon damping due to the presence of edge states in the electronic spectrum of zigzag ribbons. We finally note that the model described by (4) constitutes a rather generic description of massive particles interacting via a contact-like interaction in one dimension, with the fundamental ingredient relying on the strong ratio between the interaction and the kinetic term. Such ideas could be then implementable also in other suitable platforms such as nonlinear resonator arrays in circuit QED [54], atomic waveguides [55], newly-available 2D material heterostructures [56,57], and confined excitons in TMDs [58].

DATA AVAILABILITY
The data produced in this manuscript can be provided upon request.
Appendix A: Optical response of structured 2D materials in the quasistatic limit Invoking the quasistatic approximation, we describe the optical properties of a nanostructured twodimensional (2D) material in the R = (x, y) plane using the scalar potential where eff = ( a + b ) 2 accounts for screening of the induced charge density ρ ind by dielectric media with permittivity a and b respectively above and below the 2D material. The induced charge is obtained from the continuity equation ρ ind (R, ω) = −(i ω)∇ R ⋅ j(R, ω), while the use of Ohm's law j = σ (1) E and E = −∇Φ allows us to express the potential in Eq. (A1) self-consistently as (A2) where σ (1) (R, ω) is an isotropic 2D conductivity characterizing the intrinsic linear optical response of the 2D material. Following the method of Ref. [59], the 2D nanostructure morphology is contained in the conductivity by assuming it has the form σ (1) (R, ω) = f R σ (1) ω , where f R is 1 within the 2D structure but zero everywhere else, and the potential within the 2D material is expressed in terms of a reduced 2D coordinate vector ⃗ θ = R D as where the dimensionless parameter η eff ωD contains all dependence on the intrinsic conductivity of the 2D material (in the local limit), its dielectric environment, and characteristic size D (e.g., the diameter of a disk or the side length of a square).

Guided modes in 2D nanoribbons
We characterize the 2D nanoribbon geometry by a finite width W alongx but infinite extension inŷ, so that translational invariance in the latter dimension suggests a Fourier decomposition of the potential in wave vector components k according to Φ(R) = φ(x)e iky . Following the prescription of Ref. [42], we write the self-consistent potential of Eq. (A3) in terms of the normalized coordinate θ ≡ x W as where M = VD is a product of the differential operator and an integral operator involving the modified Bessel function K 0 . Discretizing θ ∈ [0, 1] in N equally-spaced points as {θ l } N l=1 such that h = θ l+1 − θ l for all l, we assign φ l ≡ φ(θ l ) and is obtained using central differences [42], while the boundary condition associated with the vanishing of normal current ∂ θ φ(θ) θ=0 = ∂ θ φ(θ) θ=1 = 0 at the ends of the ribbon leads to Meanwhile, assuming a slowly-varying φ(θ), the matrix decomposition of V is where θ ll ′ ≡ θ l −θ l ′ , q ≡ kW is the normalized wave vector, and L n denotes the modified Struve function of order n [60]. Note that in the q → 0 limit, charge neutrality enables replacement of K 0 (q θ − θ ′ ) → − log θ − θ ′ in the kernel of the operator V, so that Eq. (A9) becomes

Plasmon dispersion in graphene nanoribbons
The solution of Mφ n,k = η −1 n,k φ n,k yields the eigenmodes φ n,k (θ) and eigenvalues η −1 n,k of the nanoribbon geometry, while comparison to Eq. (A4) leads to the dispersion condition − η n,k =

Im{σ
(1) In the main text, the dispersion relation of mode n in a graphene nanoribbon is found by inserting the conductivity of extended graphene in Eq. (2), obtained within the local limit of the random-phase approximation (RPA). As mentioned in the main text, the presence of the interband term in the LRPA conductivity modifies the dispersion relation respect to that obtained with the Drude model. For a direct comparison, the solutions of Eq.
(A11) are plotted as red curves in Fig. 4 for the local-RPA (LRPA) conductivity given in Eq. (2), contrasted with the solutions obtained solely from the intraband contribution (Drude) or by adopting a more sophisticated model that incorporates nonlocal effects (RPA), i.e., for σ is reported in Ref. [26]. The zero of the imaginary part of the LRPA conductivity at ̵ hω plasma ≃ 5E F 3 flattens the plasmon bands compared to the simple Drude case (see Fig. 4(a)), which instead captures the dispersion only in the low-frequency ̵ hω < E F regime. On the other hand, nonlocal effects captured in the full RPA conductivity emerge at higher frequencies ̵ hω ≳ E F , becoming particularly prominent for large plasmon wave vectors ̵ hkv F ∼ E F ; the effect of these corrections is shown in Fig. 4(b), which indicates that the LRPA conductivity faithfully captures the plasmon dispersion for wave vectors kW ≲ 1 considered in this work.

Mode normalization
The modes supported by a graphene nanoribbon are normalized according to well-established procedures for quantizing the electromagnetic field in dispersive media [61,62], leading to the condition 1 4π Re is the dielectric function associated with a single 2D layer characterized by a surface conductivity σ embedded in a homogeneous medium with relative dielectric permittivity r . Here we again assume a conductivity of the form σ(R, ω) = f R σ (1) ω , so that the nanoribbon geometry is described by the dimensionless factor f R that takes a value of unity within the structure and zero everywhere else; the normalization condition is now (A14) Working in the quasistatic limit, the first integral above can be equated to the second by expressing the field as E = −∇Φ and integrating by parts, i.e., where Gauss' law ∇ ⋅ E = 4πρ r is invoked; for the 2D induced charge ρ(r) = ρ(R)δ(z), we integrate over z to write where the continuity equation ρ = −(i ω)∇ R ⋅ j = −(iσ (1) ω ω)∇ R ⋅ (f R E) was used to eliminate the charge density before again integrating by parts. With the above result, the normalization condition becomes and reduces to The field in the ribbon can be expressed as E ω (R) = ∑ n E n,k u n,k (x)e iky , where E n,k is its amplitude and u n,k =x∂ x φ n,k (x) +ŷikφ n,k (x) are obtained by decomposing the potential in eigenmodes satisfying Eq. (A4); the normalization condition then becomes, for a single mode, is the quantity depending on the mode function integral appearing in the main text, which scales as the ribbon width, i.e. ξ (1) n,k ∼ W . From the above expressions, we isolate the field amplitude Appendix B: Plasmon absorption rates

Single-plasmon absorption
In order to estimate the single-plasmon decay rate we first evaluate the power absorbed by the ribbon, which is obtained from where j (1) (R, t) is the current to linear order and ⟨. . .⟩ denotes the time average. For harmonic fields we insert ω e −iωt + c.c. and E(t) = E ω e −iωt + c.c. into Eq. (B1) to find, after dropping the fast-oscillating terms which average to zero, having expressed the current in terms of the linear conductivity via Ohm's law j Equating the absorbed power with that dissipated at a rate γ 1 by the waveguide, we obtain the expression by using the mode normalization previously derived. The above equation establishes the rate of single-plasmon absorption in the ribbon; while for the Drude model the absorption rate coincides with the inelastic scattering rate, γ 1 (ω) = γ D , the absorption rate predicted from the LRPA and RPA conductivities exhibit deviations as shown in Fig.5(a), which arise from the different mode normalization factors. Importantly, even with the inclusion of the interband transition term, the dissipation rate does not exhibit a strong intrinsic frequency dependence, and remains on the same scale as the Drude rate. For this reason, we operationally define the total single-plasmon dissipation rate as a quantity set overall by the quality factor Q according to γ 1 (ω) = ω Q, implicitly capturing all forms of dissipation.

Two-plasmon absorption
To evaluate the two-plasmon absorption rate appearing in the model Hamiltonian, we compute the work done by the nonlinear current j (3) (r, t) on the plasmon field E(r, t) within a graphene layer occupying the R = (x, y) plane, such that the absorbed power reads Decomposing the nonlinear current in its frequency components as j (3) (R, t) = j ω e −iωt + c.c., and analogously the plasmon field, we obtain by writing j The power absorption for a specific mode with index n and wave vector k is then given by P (3) n,k = 2Re σ (3) ω n,k E 4 n,k Lξ where ξ n,k = ∫ dx u n,k (x) 4 . We associate the power absorption above with two-plasmon absorption at the rate Γ (3) (ω n,k ) according to P With the above result quantifying the two-plasmon absorption rate associated with a single mode, we finally obtain the rate of two-plasmon absorption per length as γ 2 = Γ (3) L, which is the TPA rate used in the effective model. The evaluated two-plasmon absorption rate normalized respect to the ribbon wavelength and the Fermi energy as function of the frequency is shown in Fig. 5(b) for two different ribbon widths, and is found to exhibit a large enhancement in the energy range between E F < ̵ hω < ̵ hω plasma that is considered in the main text to achieve a strong plasmon interaction.
taking into account all the field components. In our problem we are interested in the specific scenario of two counter-propagating plasmons, for which the equations describing the co-propagating component ψ νν are completely decoupled from their counter-propagating counterparts ψ νν ′ .
Considering that ψ RL (y 1 , y 2 ) = ψ LR (y 1 , y 2 ), we define ψ RL = φ R and ψ LR = φ L to obtain a single wave equation depending on only one variable, i.e., the two-plasmon relative coordinate ρ = y 1 − y 2 , which explicitly reads To solve the above equation, we consider a plasmon impinging from the left and make use of the ansatz ψ R (ρ) =e ikρ θ(−ρ) + te ikρ θ(ρ) where t = 1 + r ensures continuity of the solution at r = 0 and θ(ρ) is the Heaviside step function; by imposing the boundary conditions coming from the delta function, we obtain transmission and reflection coefficients (C9) These coefficients can be rewritten in the same form as the main text, r = − [1 + 2πλ a λ p ] −1 and t = [1 + λ p (2πλ a )] −1 by defining the absorption length where we have expressed the RHS in terms of the original group velocity at the plasmon resonance.