Quantum computing with graphene plasmons

Among the various approaches to quantum computing, all-optical architectures are especially promising due to the robustness and mobility of single photons. However, the creation of the two-photon quantum logic gates required for universal quantum computing remains a challenge. Here we propose a universal two-qubit quantum logic gate, where qubits are encoded in surface plasmons in graphene nanostructures, that exploits graphene's strong third-order nonlinearity and long plasmon lifetimes to enable single-photon-level interactions. In particular, we utilize strong two-plasmon absorption in graphene nanoribbons, which can greatly exceed single-plasmon absorption to create a “square-root-of-swap” that is protected by the quantum Zeno effect against evolution into undesired failure modes. Our gate does not require any cryogenic or vacuum technology, has a footprint of a few hundred nanometers, and reaches fidelities and success rates well above the fault-tolerance threshold, suggesting that graphene plasmonics offers a route towards scalable quantum technologies.


INTRODUCTION
Quantum computing could efficiently solve many essential problems. However, building a quantum computer is not an easy task. One particularly promising approach is to use single-photons, whose weak interaction with the environment makes them perfectly suitable for encoding and transmitting quantum information. Nonetheless, this weak interaction strength makes the implementation of photon-photon interactions a significant challenge. While this can be overcome at the cost of extra photons, 1 the additional overhead makes purely linear-optical schemes difficult to scale up. 2 Alternatively, single-photon-level nonlinearities can be used to directly create deterministic gates. 3 However, this typically requires complex interactions with atomic systems that cannot readily be miniaturized. Recent work shows that graphene can provide a strong enough nonlinearity without the technical drawbacks of those atomic systems.
Our graphene-based two-qubit logic gate is centered on Franson's quantum Zeno gate, 4 which is a universal "squareroot-of-swap" (SWAP 1/2 ) gate. 5 If two separable single-qubit states |ϕ〉 and |ψ〉 enter modes 1 and 2, respectively, the gate creates an entangled superposition of these states being swapped and not swapped, i.e., where the subscripts indicate the mode. As illustrated in Fig. 1a, such an operation can be achieved by sending two photons to a 50:50 beamsplitter (BS): The gate succeeds when the two photons exit in different modes, generating the state of Eq. (1), while, half of the time, the gate will fail by allowing both photons to exit the same mode (in reality, the situation is even more complicated because of two-particle interference effects and the logical qubit encoding).
If the SWAP process is made continuous by replacing the 50:50 beamsplitter with coupled waveguides, the quantum Zeno effect 6 (wherein continuous measurement prevents a quantum system from evolving), can boost the success probability of the gate to 100%. 4 In this scenario, however, the quantum Zeno effect requires nonlinear two-photon absorption to occur at the singlephoton-level. To date, such a strong optical nonlinearity has only been achieved via complex interactions with atomic systems, 7 which lack scalability.
Plasmon-polaritons, formed when light hybridizes with the collective charge-carrier density oscillations in conducting materials, confine electromagnetic energy to deeply-subwavelength scales, and could potentially enable extremely strong optical nonlinearities in nanoscale photonic circuits 8 -an ideal situation for a scalable quantum logic gate. While plasmons supported by noble metals provide large nonlinear enhancements and are compatible with single-photon-level quantum experiments, 9,10 they suffer from intrinsically high ohmic losses, severely limiting their application to quantum technologies.
Graphene has recently arisen as a robust material platform for plasmonics, capable of sustaining plasmon resonances with extremely long lifetimes 11,12 that can be tuned actively via electrostatic gating. 13 Furthermore, its low-dimensionality provides unprecedented levels of optical field confinement, 14 boosting optical nonlinearities well above those in noble metals, [15][16][17][18] potentially enabling nonlinearities on the single-or few-plasmon level. 19,20 Here we propose that this system can be used to implement a two-qubit quantum logic gate using nanoplasmonic graphene waveguides.
We will use the so-called single-rail encoding, just as in the original Zeno-gate proposal, 4 where the absence of a particle represents a logical 0, and the presence of a particle a logical 1. In other words, |0〉 (|1〉) in the Fock basis represents a logical |0〉 (|1〉) state of the qubit. Higher-order Fock states fall out of this logical subspace. Although the single-rail encoding has limitations, 21 it can be transformed into the more well-known dual-rail encoding with linear optical elements. 22 Implementing the SWAP 1/2 gate with a BS is not straightforward (Fig. 2a): If the logical input state is |00〉, |01〉, or |10〉 (encoded by no particles in either mode, or one particle in the first or second, respectively), the gate functions perfectly. In contrast, when one particle is incident in each mode (a logical state |11〉) the correct output is |11〉. Unfortunately, the Hong-Ou-Mandel (HOM) effect, already observed for single plasmons, 9,10,23 causes the particles to bunch and exit in the same mode, implying that the gate always fails. Since the HOM effect is independent of the relative phase between the two modes, this holds in general. Even if the particles are made indistinguishable, to circumvent HOM bunching, the gate fails 50% of the times (see Fig. 1a).
In a Zeno gate, the swap between the two modes has to be a continuous process, so that a "Zeno measurement" can be applied as the system evolves. Such a continuous swap can be achieved with a directional coupler (DC). To prevent the system from evolving into a state in which both particles are in the same mode, one must continuously monitor whether both particles are in the same mode. In practice, the presence of a sufficiently strong twophoton absorber can perform this measurement. 4 At first glance, it appears that in such a DC, when the particles bunch into the same mode, they would be absorbed. However, when the swap probability is much smaller than the two-particle absorption, the Zeno effect does not even allow the particles to bunch in the first place.
In graphene, this Zeno condition can be easily achieved. When a single plasmon has less energy than the Fermi level, it is not absorbed via electron-hole pair excitation. At the same time, a mode containing two plasmon quanta can have enough energy to be absorbed via an interband transition (Fig. 1b). Since the twoplasmon absorption depends on the field strength while the single-plasmon absorption does not, confining the graphene plasmon field to a nanostructure enhances the two-plasmon absorption rate, while leaving the single-plasmon absorption rate unaffected 20 (Fig. 1c).

RESULTS
As a physical realization of such a graphene-based quantum gate, we envision a system of two graphene nanoribbons that support Fig. 1 Basic operating principles of our nanoplasmonic quantum logic gate. a Simplest square-root-of-swap gate. Two photons are sent in the two ports of a 50:50 beamsplitter. If the photons are distinguishable, half of the times the photons exit from different ports and a square-rootof-swap gate is achieved. The other half of the times the two photons exit through the same port and the gate fails. If the photons are indistinguishable, they bunch and always exit from the same port, so the gate always fails. b Electronic band structure of graphene with a non-zero Fermi energy E F . Two photons can produce an interband transition and be absorbed, whereas single-photon absorption is forbidden for photon energies below 2E F . c Ratios between the two-plasmon absorption rate, γ (2) (at the plasmon resonance frequency), and the intrinsic damping rate, γ = 500 fs −1 , for a range of nanoribbon widths, W, and Fermi energies, E F . The blue areas are regions in which two-plasmon absorption is two to six orders of magnitude faster than linear absorption, providing a strong γ ð2Þ ) γ condition that leads to extremely high success probabilities for the gate Fig. 2 Surface-plasmon-based SWAP 1/2 gate comprised of nonlinear graphene nanoribbons. Nanoribbons are brought together so that the plasmonic modes couple to each other via a Coulomb interaction. For a separation d z between the ribbons, there is an interaction length L ¼ L SWAP 1=2 after which the plasmon has 50-50% probability of remaining in the same mode or having swapped across ribbons. Thus, when a single plasmon is input in each mode, |1〉 1 |1〉 2 , we find the output state with a one plasmon in each mode, |1〉 1 |1〉 2 , in which case the gate succeeds, or b both plasmons in one of the modes, |2〉 1 |0〉 2 or |0〉 1 |2〉 2 , in which case the gate fails. When a separable single-qubit is input into each mode (|ϕ〉, |ψ〉), an entangled state is created, In the absence of nonlinearity in the waveguide and assuming indistinguishable plasmons, the HOM effect forces the plasmons to exit the gate in the same output mode, meaning that the gate always fails for |1〉 1 |1〉 2 . However, driven by the Zeno effect, the strong nonlinearity of the graphene waveguides reduces the probability that two plasmons are found in the same nanoribbon and increases the success probability. c We describe the SWAP 1/2 gate as a six-state system where U is the coupling between ribbons, while γ and γ (2) are the intrinsic damping and two-plasmon absorption rates, respectively propagating single plasmons (see Fig. 2). In this work we will assume that the single plasmons are already excited, which could, in principle, be achieved through the emission of a quantum light source. [24][25][26][27] The two nanoribbons are brought close to each other, so that the plasmons are coupled via Coloumb interaction, forming a graphene plasmon DC, whereby a plasmon starting in one ribbon can couple to the other ribbon. The interaction length, the ribbon width, and the ribbon spacing set the splitting ratio of the DC. At the same time, the ribbon width and the Fermi energy of the nanoribbons determine the two-plasmon absorption rate.
To model this system, we describe each ribbon as a two-level system with energy ℏω, where ω is the resonant plasmon frequency that depends on the nanoribbon width W and doping level (Fermi energy) E F . As shown in Fig. 2c, we consider a maximum of two plasmons, limiting the Hilbert space to six states. States with an equal number of plasmons are coupled via a Coulomb interaction of strength U. Decay processes are governed by inelastic scattering rate γ, and γ (2) denotes the two-plasmon absorption rate.
We quantify the Coulomb interaction by describing plasmons in semi-infinite graphene nanoribbons within the so-called plasmon wave function (PWF) formalism, 28 adapted here to include the effect of a non-vanishing optical wave vector k || in the direction of the ribbon transversal symmetry. Setting the nanoribbons to be aligned horizontally, and separated by a distance d z in the zdirection (see Fig. 2a), the interaction between N plasmons in one ribbon and N′ plasmons in the neighboring one, both of them propagating with parallel wave vector k || , is given by where the integrals are evaluated over the nanoribbons in a 2D space R = (x, y) and ρ ind k jj ;N ðR; ωÞ is the induced charge associated with N plasmons (see Methods and Fig. S1).
Next, we compute γ (2) from the nonlinear conductivity σ ð3Þ ω , for which an analytical expression in the local and zero-temperature approximation is obtained quantum-mechanically in the Dirac cone approximation and reported in ref. 29 . As shown in the Methods, the two-plasmon absorption rate is given by where β ð2Þ q;1 and β ð4Þ q;1 are the momentum-dependent field normalizations, which we consider to be unity for low momentum values. Here Δ characterizes the spatial extent of the propagating plasmon along the direction of transversal symmetry, which we set to be equal to the ribbon width. We set the single-plasmon lifetime to be γ = 500 fs −1 , which is a realistic value, measured at room temperature. 11 Note that this lifetime can be extended by going to cryogenic temperatures; for which lifetimes up to 10 ps have recently been measured. 12 We can now calculate the density matrix ρ(t) of the system by solving the time-dependent Lindblad master equation, which is the most general type of Markovian and time-homogeneous master equation describing an open-quantum-system evolution that is both trace-preserving and completely positive for any initial condition 30 γ ðnÞ a n m ρa yn m À 1 2 a yn m a n m ; ρ where γ (1) ≡ γ, a y m (a m ) denote plasmon creation (annihilation) operators, n is the number of absorbed plasmons and m is the nanoribbon mode. The Hamiltonian of the two-nanoribbon system is where U is the Coulomb interaction given in Eq. (2), while ω is the plasmon frequency of each nanoribbon mode. We numerically solve Eq. (4) using Mathematica, from which we obtain the required time t SWAP 1=2 at which a single plasmon incident in either nanoribbon is placed in an equal superposition of both nanoribbon modes at the output. This time is related to the Coulomb interaction U from Eq. (2) (i.e. stronger Coulomb interaction U resulting in shorter t SWAP 1=2 ). To calculate t SWAP 1=2 we define our initial state to be ρ(t = 0) = |ψ i 〉〈ψ i |, where |ψ i 〉 = |1〉 1 | 0〉 2 , and let it evolve until the probability of the plasmon being in either of the modes is equal: P 10 j i ðt SWAP 1=2 Þ ¼ P 01 j i ðt SWAP 1=2 Þ. We convert this time to a length L SWAP 1=2 , by computing the plasmon group velocity as shown in Fig. S2. The resulting L SWAP 1=2 is plotted in Fig. 3a. For E F > 0.1 eV the required L SWAP 1=2 is always less than the single-plasmon decay length, thus showing the potential of long-lived graphene plasmons: novel physical effects can manifest before the plasmon decays.
For all the results presented here, we set the spacing between the two nanoribbons to d z = 1 nm. With current technology, such atomically thin spacings can be realized by taking advantage of 2D materials like graphene. 31 This parameter only affects the Coulomb interaction strength, which will determine L SWAP 1=2 . The PWF used in our calculations is applicable for these scales, as discussed in detail in ref. 28 . Furthermore, for our parameter regime, the Coulomb interaction does not depend very strongly on d z (see Fig. S4 of the Supplementary Information).
Once L SWAP 1=2 is determined, we proceed to analyze the system when a single plasmon is input in each mode; that is, For this input, the gate functions correctly if there is still one plasmon in each output mode, which occurs with probability P 11 j i ðt SWAP 1=2 Þ.

DISCUSSION
In Fig. 3a-c we show the success probability P 11 j i ðt SWAP 1=2 Þ, the probability of the plasmons bunching in the same nanoribbon P 20 j i ðt SWAP 1=2 Þ þ P 02 j i ðt SWAP 1=2 Þ, and the probability for both plasmons to decay P 00 j i ðt SWAP 1=2 Þ, for a range of nanoribbon widths W and Fermi energies E F . Notice the similarity of the contour features between these figures and the γ (2) /γ ratio shown in Fig. 1c. In the upper right corner the two-plasmon absorption is much weaker than the single-plasmon absorption, leading to a very weak Zeno effect, so the HOM effect prevails: that is, P 20 j i ðt SWAP 1=2 Þþ P 02 j i ðt SWAP 1=2 Þ ) P 11 j i ðt SWAP 1=2 Þ. As we decrease both W and E F , γ (2) increases, but not enough to drive a noticeable Zeno effect. Instead, both of the plasmons are likely to be absorbed, which is reflected in P 00 In the region where γ ð2Þ =γ $ 10 4 À 10 6 , a strong Zeno effect can be realized (light blue area of Fig. 1c). This leads to a large increase in the success probability P 11 j i ðt SWAP 1=2 Þ, while P 20 j i þ P 02 j i ðt SWAP 1=2 Þ becomes negligible, meaning that the Zeno effect completely suppresses the HOM effect. Despite the large γ (2) , P 00 j i ðt SWAP 1=2 Þ shows a minimum when γ ð2Þ ) γ. In this optimal region, we find a maximum success probability of 87.0% for W = 5 nm and E F = 0.335 eV, which is an increase in the success probability of the SWAP 1/2 gate from 0 to 87.0%. This already places us well above the gate success probability rate required to generate universal cluster states for quantum computation. 32 This performance is limited by the single plasmon lifetime. In Fig. 3e we plot the success probability, maximized over the range of W and E F shown in pannels a-d, versus the plasmon lifetime given by 1/γ. For lifetimes longer than 7.5 ps the success probability increases above 99%, reaching fault-tolerance regimes for surface codes. 33 Nevertheless, edge imperfections and structural defects would decrease the plasmon lifetime and thus the fidelity of the gate. The predicted nonlinearities, nevertheless, should persist in their presence.
Since single-plasmon decay can also result in logical states changing into other logical states, this process fidelity will be further decreased. To quantify this, we evaluated the process fidelity 34,35 of our gate by simulating process tomography for the complete range of W and E F under consideration (see Methods). The resulting process matrix for W =5 nm and E F = 0.335 eV with a lifetime of 500 fs is plotted in Fig. 4, and has a fidelity of 93.3%. When the lifetime is increased to 10 ps, the fidelity is 99.6%.
Our proposed gate achieves process fidelities in the faulttolerance regime for relatively reasonable physical parameters. Doping levels as high as 1-2 eV have been achieved, 36,37 nanoribbon widths in the range of 10-40 nm have been constructed using different means, 31,38-40 and separation distances ≈1 nm are routinely achieved through single-atomic hexagonal boron nitride spacers, which additionally guarantees the preservation of high-quality graphene optical response. 31 Furthermore, by combining ideas from quantum optics with nanoplasmonics, our work opens up an entirely new and promising avenue in the search for single-photon nonlinearities. While we have studied the application of graphene nanoplasmonics to a quantum logic gate, this could also be used for deterministic optical implementations of quantum teleportation, 41 cluster state generation, 42 and singlephoton sources, 19 underlining the applicability of this platform.

Classical electrostatic description of plasmons in graphene nanoribbons
We consider a single graphene nanoribbon occupying the R = (x, y) plane that has a finite width W in the x-direction and is infinitely-extended along the y-direction. In the linear approximation, following refs., 19,23 the selfconsistent electric field within the ribbon E q produced by an impinging field E ext ðR; tÞ ¼ E ext q e iðky yÀωtÞ þ c:c:, i.e., having frequency ω and momentum k y ≡ q/W along y, is given by is the average of the dielectric functions describing media above (ðε a ω Þ) and below ðε b ω Þ the 2D layer and ρ ind q ðR; ωÞ is the induced charge. From the continuity equation, we express ρ ind q in terms of the local, linear 2D graphene conductivity σ ð1Þ ω as ρ ind q ðR; ωÞ ¼ À where we have introduced the occupation factor f R , which is equal to one when −W/2 ≤ x ≤ W/2 and is vanishingly small everywhere else. In practice, we employ the optical conductivity obtained for zero temperature in the local limit (i.e., for vanishing in-plane optical momentum) of the random- Fig. 3 Performance of the graphene-based SWAP 1/2 for different nanoribbon width W and Fermi energy E F . Here the separation between the nanoribbons is set to d z = 1 nm, and the in plane momentum along the ribbon to k jjj W ¼ 0:4. a Probability of still having one plasmon in each mode when one plasmon is input into each mode after the input plasmons evolve along the interaction length L SWAP 1=2 . We find a range (shown in white) where the success probability is over 80% for reasonable physical parameters. b Probability of finding two plasmons in one nanoribbon after the interaction between the initial plasmons occurs. This is the "failure probability" of the gate, as it corresponds to events which take us out of the logical qubit subspace. As expected, these data show that in the region where P |11〉 is maximized the failure probability is significantly suppressed. c Probability of losing both initial plasmons after they evolve along a distance L SWAP 1=2 . d Interaction length L SWAP 1=2 required to perform the SWAP 1/2 logic gate. For the plotted range, we find that, above E F = 0.1 eV, the required interaction length is always shorter than the plasmon decay length (which is ≈ 500 nm for a 500 fs lifetime). e Success probability of the |11〉 input state as a function of the plasmon lifetime 1/γ, maximized over the same W and E F range as in panels a-c I. Alonso Calafell et al.
phase approximation (RPA) as 18 where the Fermi energy E F is related to the graphene Fermi velocity v F ≈ c/ 300, doping charge-carrier density n according to E F ¼ hv F ffiffiffiffiffi ffi πn p and τ = 1/γ is a phenomenological inelastic scattering rate. The first and second terms in Eq. (8) describe the optical response arising from intraband and interband electronic transitions, respectively, with the latter becoming unimportant when E F tω. 13 Incidentally, we have neglected inelastic damping in the interband transitions. In terms of normalized coordinates θ R=W and the normalized electric fieldε q ðθ; ωÞ W ffiffiffi ffi fθ p E q ðθ; ωÞ, Eq. (6) can be expressed as ε q ðθ; ωÞ ¼ε ext q ðθ; ωÞ þ η ð1Þ ω Z d 2θ0 Mðθ;θ 0 Þ Áε q ðθ 0 ; ωÞ; where η ð1Þ ω iσ ð1Þ ω =ε ab ω ωW is a dimensionless parameter characterizing the intrinsic linear optical response of graphene, and which we identify as a real, symmetric operator that admits a complete set of real eigenvalues. The electric field of Eq. (9) is expanded in eigenmodes of the matrix Mðθ;θ 0 Þ as where the modesε q;m ðθ x Þe iqθy and their corresponding eigenvalues η q,m satisfỹ and form an orthonormal set Z dθ xε Ã q;m ðθ x Þ Áε q;m ðθ x Þ ¼ δ mm0 : Inserting Eq. (11) into Eq. (9), we make use of Eqs. (12 and 13) to write a m ¼ b q;m ð1 À η is a coefficient that depends on the form of the external field. In what follows we take E ext q to be independent of x, so we may write b q;m ðθ y Þ ¼ ÀWE ext q Áξ Ã q;m c q ðθ y Þ, where c q (θ y ) contains the y-dependence of the external field andξ q;m À R dθ xεq;m ðθ x Þ, so that the normalized electric field in Eq. (11) becomes ω =η q;mε q;m ðθ x Þc q ðθ y Þ: (15) Electrostatic energy in nanoribbons The electrostatic energy for identical, parallel ribbons separated by a distance d z in the z-direction is given by where, from Eq. (7) (taking ε ab ω ¼ 1 for simplicity), we can express the induced charge in ribbon l as Inserting the above expression into Eq. (16) and making use of Eq. (15), the electrostatic energy becomes where and L → ∞ is the nanoribbon length. Assuming a plane wave field profile along the y-direction corresponding to c q ðθ y Þ ¼ e iqθy , in a single ribbon (i.e., taking l = l′ and d z = 0), the use of Eqs. (10) and (13) yields I q,mm (0) = −Lδ mm′ /Wη q,m , and so the electrostatic energy per unit length in ribbon l is In practice, we restrict our study to the lowest-order m = 1 mode, and fix the number of plasmon quanta in this mode using the condition l hω p ¼ 2ΔŨ q;l , where Δ is an effective length for the plasmon mode along the ribbon (i.e., the characteristic spatial width of a pulse), leading to E ext q;l Áξ Ã2 q;1 where it is now understood that the indices l and l′ denote the number of plasmons in the first and second ribbon, respectively. Using the above

Plasmon normalization
We normalize the electric field amplitude of the plasmon mode by equating the absorbed and dissipated power at linear order, i.e., where l is the number of plasmon quanta, j ð1Þ q ðR; ωÞ ¼ σ ð1Þ ω E q ðR; ωÞ, and 〈...〉 denote the time-average. Using the result of Eq. (15) with only the m = 1 mode, we obtain where β ðnÞ q;1 ¼ R 1=2 À1=2 dθ xεq;1 ðθ x Þ n . For a mode defined as a plane-wave along the ribbon, such that c y ðθ y Þ ¼ e iqθy within an effective length Δ, we write the normalization condition for N plasmons as Two-plasmon absorption rate Power absorption in a nanoribbon via two-plasmon absorption arises from the nonlinear current j ð3Þ ðR; tÞ ¼ j ð3Þ q ðR; ωÞe iky yÀiωt þ c:c:, where j ð3Þ q ðR; ωÞ ¼ σ ð3Þ ω E q ðR; ωÞ 2 E q ðR; ωÞ, and is given by where j ð3Þ q ðR; ωÞ ¼ σ ð3Þ ω E q ðR; ωÞ 2 E q ðR; ωÞ and σ ð3Þ ω is the local third-order conductivity of extended graphene, for which we adopt the analytical result obtained quantum-mechanically at zero temperature in the Dirac cone approximation, as reported in ref. 29 Using Eq. (15) we write the timeaverage of the absorbed power per unit length as Equating 〈P TPA 〉 with the power dissipated by two-plasmon absorption, 2ℏωγ (2) , we make use of the field normalization condition in Eq. (24) to write the two-plasmon absorption rate for a ribbon containing l = 2 plasmons in the m = 1 mode as In obtaining the above expression, we have again chosen the field along the ribbon to have the form of a plane-wave (i.e., c y ðθ y Þ ¼ e iqθy ), and an effective length Δ.

Process tomography
We send a complete set of 16 two-qubit states through our simulation and compute the output states at t SWAP 1=2 . To deal with failure events, when |2〉 1 |0〉 2 and |0〉 1 |2〉 2 terms arise in the output states, we truncate the output density matrix and renormalize the result. Such events only occur when states involving two plasmons are input. We also numerically correct for local single-qubit phases which arise in the output of the simulation. We feed these output states in a least-squares process tomography routine, generating a process matrix χ sim . This process matrix is defined as, where ρ in(out) is the input (output) density matrix, and E i are the basis operators constructed from the Kronecker product of the Pauli matrices (labels of Fig. 4. We calculate the process fidelity between these, and the ideal process (given by Eq. (10) of ref. 12 as Tr{χ sim χ ideal }. 34,35 Numerical solution of the linblad master equation We use the Lindblad equation introduced in Eq. (4) to describe and solve the density matrix of our system. The first term of the Lindblad equation contains the Hamiltonian given in Eq. (5). This Hamiltonian describes the two identical graphene nanoribbons as a two-level system, where the coupling between the levels is given by the Coulomb interaction U. We define a 6-state Hilbert space that contains a vacuum state (|0〉 1 |0〉 2 ), two single-plasmon states (|1〉 1 |0〉 2 , |0〉 1 |1〉 2 ) and three two-plasmons states (|1〉 1 |1〉 2 , |2〉 1 |0〉 2 , |0〉 1 |2〉 2 ). In this basis, the matrix form of the Hamiltonian is where ℏω is the energy of the plasmon. The second term of the Lindblad equation contains the loss channels of the system; namely, the singleplasmon absorption γ (1) and the two-plasmon absorption γ (2) . In matrix form, this second term reduces to Àγ ð1Þ ρ 0011 À 1 2 ð2γ ð1Þ þ γ ð2Þ Þρ 0020 À 1 2 ð2γ ð1Þ þ γ ð2Þ Þρ 0002 À 3 2 γ ð1Þ ρ 1011 À 1 2 ð3γ ð1Þ þ γ ð2Þ Þρ 1020 À 1 2 ð3γ ð1Þ þ γ ð2Þ Þρ 1002 À 3 2 γ ð1Þ ρ 0111 À 1 2 ð3γ ð1Þ þ γ ð2Þ Þρ 0120 À 1 2 ð3γ ð1Þ þ γ ð2Þ Þρ 0102 À2γ ð1Þ ρ 1111 À 1 2 ð4γ ð1Þ þ γ ð2Þ Þρ 1120 À 1 2 ð4γ ð1Þ þ γ ð2Þ Þρ 1102 À 1 2 ð4γ ð1Þ þ γ ð2Þ Þρ 2011 Àð2γ ð1Þ þ γ ð2Þ Þρ 2020 Àð2γ ð1Þ þ γ ð2Þ Þρ 2002 À 1 2 ð4γ ð1Þ þ γ ð2Þ Þρ 0211 Àð2γ ð1Þ þ γ ð2Þ Þρ 0220 Àð2γ ð1Þ þ γ ð2Þ Þρ 0202 where ρ ijkl ¼ i j i 1 j j i 2ρ k j i 1 l j i 2 . So as to obtain the time-dependent density matrix of the system, we numerically solve the system of ordinary differential equations in Wolfram Mathematica. We employ the variable stepsize implicit Backward Differentiation Formulas (BDF) or order 5. The WorkingPrecision used in this algorithm was set to the MachinePrecision, which, in our case, corresponds to 16 digits. In addition, the AccuracyGoal and PrecisionGoal options are set to 10. The diagonal elements of this density matrix exactly correspond to the probability of the plasmons being in different modes. For example, ρ 1111 (t) is the probability that one plasmon is found in each nanoribbon at a given time, ρ 2020 (t) + ρ 0202 (t) is the probability that two plasmons are found in a single nanoribbon at a given time, and ρ 0000 (t) is the probability of not having any plasmon in the system at a given time.
Once the density matrix of our system is found, we proceed to find the required interaction time between the nanoribbons to implement a SWAP 1/2 . To do so, we set our initial state to be ρ(t = 0) = |ψ i 〉〈ψ i |, where |ψ i 〉 = |1〉 1 |0〉 2 , let it evolve in time and find t SWAP 1=2 by looking for the time at which the probability of the plasmon being in either of the modes is equal; i.e., P 10 j i ðt SWAP 1=2 Þ ¼ P 01 j i ðt SWAP 1=2 Þ. The solution to this condition was found numerically using Wolfram Mathematica with a minimum accuracy and precision of 10 digits. Once t SWAP 1=2 is determined, we define our initial state to be ρ(t = 0) = |ψ i 〉〈ψ i |, where |ψ i 〉 = |1〉 1 |1〉 2 , and find the success probability of the gate P 11 at time t SWAP 1=2 . Representative time-dependent densitymatrix elements are plotted in Fig. S6 in the Supplementary Information.

DATA AVAILABILITY
The datasets generated and analysed during the current study are available from the corresponding author if you ask nicely.