# Experimental realization of nonunitary multi-qubit operations

M. W. van Mourik  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25/4, 6020 Innsbruck, Austria  E. Zapusek  Institute for
Quantum Electronics, ETH Zürich, Otto-Stern-Weg 1, 8093 Zürich, Switzerland
P. Hrmo  L. Gerster  R. Blatt  Institut für Experimentalphysik, Universität
Innsbruck, Technikerstraße 25/4, 6020 Innsbruck, Austria  T. Monz  Institut
für Experimentalphysik, Universität Innsbruck, Technikerstraße 25/4, 6020
Innsbruck, Austria  AQT, Technikerstraße 17, 6020 Innsbruck, Austria  P.
Schindler  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25/4, 6020 Innsbruck, Austria  F. Reiter  Institute for
Quantum Electronics, ETH Zürich, Otto-Stern-Weg 1, 8093 Zürich, Switzerland

###### Abstract

We demonstrate a novel experimental toolset that enables irreversible multi-
qubit operations on a quantum platform. To exemplify our approach, we realize
two elementary nonunitary operations: the OR and NOR gates. The electronic
states of two trapped ${}^{40}\mathrm{Ca}^{+}$ions encode the logical
information, and a co-trapped ${}^{88}\mathrm{Sr}^{+}$ion provides the
irreversibility of the gate by a dissipation channel through sideband cooling.
We measure 87% and 81% success rates for the OR and NOR gates, respectively.
The presented methods are a stepping stone towards other nonunitary operations
such as in quantum error correction and quantum machine learning.

Introduction.— Classical computing is an immensely successful information
processing paradigm. The success of computing can largely be explained by the
rapid increase in computational power enabled by the miniaturization of the
underlying circuits built from classical, irreversible gate operations (cf.
Fig. 1(a)). Today, the exponential growth of gate count on classical
processors is reaching fundamental physical limits [1]. In the continued
pursuit of increasing computational power, a multitude of alternate
technologies is being explored [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13].

As an approach orthogonal to classical information processing quantum
computing has recently received considerable attention. Here, substantial
advancements have been made, allowing for first demonstrations of essential
ingredients such as quantum error correction [14, 15, 16, 17, 18, 19]. This
can be attributed to novel and advanced proposals and the continued
improvement of established techniques [20, 21, 22, 23, 24]. Such advancements
in controllability bring quantum computation closer to the ideal of an
entirely unitary evolution towards the output state.

Figure 1:  (a) Truth tables of the classical OR and NOR gates, with two-qubit
output. The logical output is mapped on the left qubit. (b) Schematic
representation of the OR gate acting on $\left|01\right>$: an engineered
resonance process $D_{f}$ is a combination of a global and single-ion laser
pulse, which together allow transfer to the desired state, $\left|11\right>$,
plus an increase in motional mode occupation. This action is made irreversible
by dissipation $\Gamma_{f}$ of this additional motion, by cooling a spectator
ion species. (c) Overview of the relevant states in the data ions,
${}^{40}\mathrm{Ca}^{+}$, and cooling ion, ${}^{88}\mathrm{Sr}^{+}$. Logical
bits $\left|0\right>$ and $\left|1\right>$ are stored in Calcium’s $4S_{1/2}$
ground states. Auxiliary levels in the $3D_{5/2}$ manifold are used for
engineered resonance transfer, $D_{f}$. Dissipation $\Gamma_{e}$ and
$\Gamma_{f}$ occurs through spontaneous decay from $P_{3/2}$ to $S_{1/2}$, in
Calcium and Strontium.

In certain algorithms, however, _nonunitary_ operations are required in
combination with unitary quantum gates. Among these are algorithms for quantum
machine learning, quantum optimization, and simulation, which are regarded as
some of the most promising near-term applications for quantum information
processing [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 18].
Specifically, nonunitary operations are needed for the generation of low-
temperature thermal states [27, 28, 29, 30], as a projective filter [31], for
the simulation of open systems [32, 33, 34, 35, 36, 37], or in quantum neural
networks [38, 18]. It has been suggested to implement nonunitary components
through auxiliary qubits, randomized circuits, or mixed input states [31, 30,
33, 34, 36, 35, 27, 32].

Dissipation is inherently nonunitary, making it a natural choice for the
creation of irreversible operations. The field of dissipation engineering, or
reservoir engineering, uses the interaction of a quantum system with
environmental degrees of freedom to achieve quantum information processing
tasks [39, 40, 41, 42, 18, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52].
Applications include state preparation by optical pumping, squeezing [40, 41],
entanglement generation [42, 43, 44, 45, 46, 47, 48, 49], quantum simulation
[50], and quantum error correction [51]. Dissipation towards the environment
lifts the requirement for classical measurement and feedback, and it holds
scaling and robustness advantages over unitary approaches [43, 44, 45]. It has
been formally shown that dissipation can be used to perform universal quantum
computation [53]. Still, so far dissipation engineering has mostly been
focused on quantum state preparation and subspace stabilization. We expand the
possible set of applications taking a step towards a paradigm of general
nonunitary quantum operations, by demonstrating the realization of
irreversible classical gates (cf. Fig. 1(a)) by means of engineered
dissipation.

Here we present a physical realization of nonunitary operations in a trapped-
ion system by use of dissipation engineering. By utilizing techniques from
dissipation engineering, one can create nonunitary quantum gates that operate
deterministically and without the need for ancilla qubits [54]. To this end,
from quantum-mechanical interactions, we engineer the desired projective
dynamics effecting classical gate operations. We implement a classical OR and
NOR gate, whose truth tables are shown in Fig. 1(a), where the output of the
gate action is mapped onto the left qubit. We employ selective coherent
couplings to conditionally excite electronic states, utilizing the ions’
shared motional modes, schematically outlined in Fig. 1(b). Both sympathetic
cooling and decay via an auxiliary level serve as the nonunitary components
and complete the gate action. Experimentally, the desired dynamics can be
implemented in a mixed-species trapped-ion system with single-qubit addressing
capabilities [55]. Through our work we show that carefully engineered
nonunitary quantum dynamics have the potential to enrich the quantum
engineer’s toolbox, by performing a broad class of operations.

Principle of operation.— Two co-trapped ${}^{40}\mathrm{Ca}^{+}$ions serve as
information carriers, with the logical states $\left|0\right>$ and
$\left|1\right>$ encoded in the $4S_{1/2}(m=-1/2)$ and $(m=+1/2)$ Zeeman sub-
levels (see Fig. 1(c)). The ions are trapped in a harmonic potential, and
share motional modes. We use the notation
$\left|i\right>_{1}\otimes\left|j\right>_{2}\otimes\left|n\right>_{m}=\left|ij\right>\left|n\right>_{m}$,
for electronic states $i$ and $j$ of ions 1 and 2, and mode occupation $n$ of
a specific motional mode $m$. For brevity, the mode occupation $n$ is often
omitted in our notation when $n=0$, i.e.
$\left|ij\right>=\left|ij\right>\left|0\right>_{m}$. As seen in the truth-
table in Fig. 1(a), the OR gate corresponds to the mapping of $\ket{01}$ to
$\ket{11}$, which is analogous to the condition that the first qubit is
flipped from $\left|0\right>$ to $\left|1\right>$ if and only if the second
qubit is the state $\ket{1}$. The desired conditional operation is augmented
by making use of a specific motional mode to encode information about the
parity of the system. Access to motional modes is enabled through the
auxiliary state $\left|f\right>$, for which we use the metastable
$3D_{5/2}(m=+1/2)$ level. The $4S_{1/2}\leftrightarrow 3D_{5/2}$ transition to
this auxiliary state is coupled with coherent 729 nm light.

Figure 2: Overview of the gate mechanism. (a) Desired (left) and undesired
(right) process for the OR gate. State $\ket{01}\ket{0}_{m}$ is off-resonantly
driven to $\ket{f1}\ket{0}_{m}$ with Rabi frequency $\Omega_{f}$. Due to the
coupling $\Omega_{\mathrm{SB}}$, $\ket{f1}\ket{0}_{m}$ is hybridized with
states $\ket{11}\ket{1}_{m}$ and $\ket{1f}\ket{0}$, resulting in a dressed
state splitting. For $\Delta=\Omega_{\mathrm{SB}}/\sqrt{2}$, the carrier drive
is on resonance with the dressed state, and is therefore excited. The gate
action is completed by sympathetically cooling the motional mode. The dressed
states $\ket{f0}\ket{0}_{m}$ and $\ket{10}\ket{1}_{m}$, accessible from
$\ket{00}\ket{0}_{m}$, are shifted by $\pm\Omega_{\mathrm{SB}}/2$, and are
therefore not in resonance with the carrier drive. Excitation from
$\ket{00}\ket{0}_{m}$ is thus suppressed. (b) Pulse sequence of OR and NOR
gates. Transitions used in the experiment are indicated in Fig. 1(c).
$\Omega_{f}$ and $\Omega_{e}$ act on the first ion, coupling $\left|0\right>$
with $\left|f\right>$ and $\left|1\right>$ with $\left|e\right>$.
$\Omega_{\mathrm{SB}}$ acts on both ions, coupling the red sideband of
$\left|1\right>$ and $\left|f\right>$.

The population transfer mechanism is summarized below (details in Supplemental
material 111The equations describing the dynamics of the system are given in
the Supplemental material, Section A). A drive with Rabi frequency
$\Omega_{f}$ and detuning $\Delta$ is applied to the first ion’s
$\left|0\right>\leftrightarrow\left|f\right>$ transition. We refer to this
drive as the probe. Without any further couplings, the probe would excite the
two states $\left|00\right>$ and $\left|01\right>$, and leave the states
$\left|10\right>$ and $\left|11\right>$ unchanged. A second drive, which we
refer to as the sideband drive, is applied to both ions on the
$\left|f\right>\leftrightarrow\left|1\right>$ transition, though red-detuned
by the frequency of the motional mode. This drive therefore couples the states
$\ket{f}\ket{n}_{m}\leftrightarrow\ket{1}\ket{n+1}_{m}$. In particular, the
transition between
$\left|f\right>\ket{0}_{m}\leftrightarrow\ket{1}\ket{1}_{m}$ occurs with Rabi
frequency $\Omega_{\mathrm{SB}}$. Under the condition that
$\Omega_{\mathrm{SB}}\gg\Omega_{f}$, the states excited from
$\ket{01}\ket{0}_{m}$ form dressed states
$(\left|f1\right>\ket{0}_{m}+\ket{1f}\ket{0}_{m}\pm\sqrt{2}\left|11\right>\ket{1}_{m})/2$.
These dressed states have a frequency shift of
$\pm\Omega_{\mathrm{SB}}/\sqrt{2}$, with respect to the bare
$\ket{f1}\ket{0}_{m}$ state, as shown in Fig. 2(a), left. In contrast, the
initial state $\ket{00}\ket{0}_{m}$ is excited to the dressed states
$(\ket{f0}\ket{0}_{m}\pm\ket{10}\ket{1}_{m})/\sqrt{2}$ (Fig. 2(a), right),
which reside at frequencies $\pm\Omega_{\mathrm{SB}}/2$. The three-level
dressed states have increased frequency shifts compared to the two-level
dressed states, because of constructive interference of the couplings
$\\{\ket{1f}\ket{0}_{m},\ket{f1}\ket{0}_{m}\\}\leftrightarrow\ket{11}\ket{1}_{m}$.
Choosing a probe pulse detuning $\Delta=\Omega_{\mathrm{SB}}/\sqrt{2}$
therefore enables excitation from $\left|01\right>$, while excitation from
$\left|00\right>$ is out of resonance and thus suppressed. This population
transfer through engineered resonance is denoted by $D_{f}$ in Fig. 1(b) and
(c).

The conditional excitation is made nonunitary by a decay process enabled by
sideband cooling a co-trapped ${}^{88}\mathrm{Sr}^{+}$ion, indicated in Fig.
1(c). In the OR gate, population that cycles through $\ket{11}\ket{1}_{m}$ is
dissipatively transferred to $\ket{11}\ket{0}_{m}$ at a rate $\Gamma_{f}$.
Since $\ket{11}\ket{0}_{m}$ does not couple with either the sideband drive nor
the probe, population remains in this state, thus completing the transfer
$\ket{01}\rightarrow\ket{11}$. In order to avoid interfering with the
excitation during the probe process, the dissipation is realized in a
subsequent step [44, 48].

We expand this mechanism to the universal NOR gate, whose truth table is shown
in Fig. 1(a). This gate can be constructed by concatenating a unitary NOT gate
with the OR gate. However, we present a fully dissipative implementation where
the mapping $\ket{00}\rightarrow\ket{10}$ follows the same procedure as the OR
gate. In contrast to the OR gate, for the NOR gate we use the detuning
$\Delta=\Omega_{\mathrm{SB}}/2$ to excite the initial state $\ket{00}$. In
addition, the transfers $\ket{11}\rightarrow\ket{01}$ and
$\ket{10}\rightarrow\ket{00}$ are required for the gate action which can be
implemented by a single-qubit dissipative process. Both mappings are achieved
by optically pumping the first ion from $\ket{1}$ to $\ket{0}$ over another
auxiliary level, $3D_{5/2}(m=-3/2)\equiv\ket{e}$, and subsequently to
$4P_{3/2}(m=-3/2)$, from where spontaneous decay returns population to the
$\left|0\right>$ state.

Experimental overview.— The experiments have been carried out with a segmented
surface trap in a cryogenic environment [55]. Ions are stored in the Ca-Sr-Ca
configuration. Collisions with particles in the background gas may disrupt
this orientation. Therefore, we periodically apply a sequence of voltages to
the trap electrodes that deterministically place the ions back in the desired
configuration [57, 58].

We set the sideband drive to couple to the crystal’s axial in-phase (ip) mode.
Confining potentials are set so that the in-phase mode frequency is
$\omega_{\mathrm{ip}}/(2\pi)=$550\text{\,}\mathrm{kHz}$$. This value is chosen
as a trade-off between ensuring a low motional mode heating rate
($\dot{\bar{n}}\propto\omega^{-\alpha}_{\mathrm{ip}}$ with $\alpha\approx 2$)
and a sufficiently high coupling to the motional mode $\Omega_{\mathrm{SB}}$
through laser interaction
($\Omega_{\mathrm{SB}}\propto\omega^{-1/2}_{\mathrm{ip}}$) [59]. At this
frequency, we have measured an axial in-phase mode heating rate of 106(20)
phonons per second, and an initial mean mode occupation of 0.14 phonons after
sideband cooling. Both the axial in- and out-of-phase modes of the ion crystal
are sideband cooled. Ions are initialized in $\left|00\right>$ using optical
pumping. We prepare the remaining possible initial states $\left|01\right>$,
$\left|10\right>$, and $\left|11\right>$, using a combination of single-ion
and collective $\pi$-pulses on the
$\left|0\right>\leftrightarrow\left|f\right>$ and
$\left|f\right>\leftrightarrow\left|1\right>$ transitions.

The sequences of operations for the OR and NOR operations are schematically
shown in Fig. 2(b), referring to the states shown in Fig. 1(c), with
$\Omega_{f}$ and $\Omega_{e}$ acting only on the first ion, coupling
$\left|0\right>$ with $\left|f\right>$ and $\left|1\right>$ with
$\left|e\right>$, respectively. $\Omega_{\mathrm{SB}}$ acts on both ions, and
couples the red sideband of $\left|1\right>$ and $\left|f\right>$.

For the OR operation, the initial state $\left|01\right>$ is to be transferred
to $\left|11\right>$, while all other initial states remain unchanged. The
dressed state splitting is produced with a sideband drive with Rabi frequency
$\Omega_{\mathrm{SB}}/(2\pi)\approx$8\text{\,}\mathrm{kHz}$$. The probe beam,
simultaneously applied to only the first ion, is detuned by
$\Delta=\Omega_{\mathrm{SB}}/\sqrt{2}$ from the
$\left|0\right>\leftrightarrow\left|f\right>$ carrier transition, with an on-
resonance Rabi frequency
$\Omega_{f}/(2\pi)\approx$1.15\text{\,}\mathrm{kHz}$$. These pulses are
applied for a duration of
$2\pi/\Omega_{f}=$900\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, which excites
$\left|01\right>$ to the dressed state as shown in Fig. 2(a) left. Similar
transfer from an initial state of $\left|00\right>$ is suppressed because the
resonance condition, shown in Fig. 2(a) right, is not met.

Following this state-dependent population transfer, the state
$\ket{11}\ket{1}_{m}$ is dissipatively transferred to $\ket{11}\ket{0}_{m}$ by
cooling the Sr ion. The sideband coupling $\Omega_{\mathrm{SB}}$ is maintained
during the cooling step, which fully depletes the populated dressed state.

The NOR gate follows a similar procedure as above, shown in Fig. 2(b), though
now population transfer from $\left|00\right>$ is enabled by choosing
$\Delta=\Omega_{\mathrm{SB}}/2$, which suppresses excitation from
$\left|01\right>$. The additional channel of dissipation required by the NOR
operation, $\left|1\right>\rightarrow\left|0\right>$ for only the first ion,
is performed in multiple steps since it would otherwise conflict with the
simultaneously required $\left|00\right>\rightarrow\left|10\right>$ operation.
Preceding the engineered dissipation, population in the first ion’s
$\left|1\right>$ state is stored in $\left|e\right>$. After the engineered
dissipation, $\sigma^{-}$-polarized light at 854 nm transfers population in
$\left|e\right>$ to the $4P_{3/2}$ level, favoring the $m=-3/2$ Zeeman
sublevel, from which spontaneous decay brings it to $\left|0\right>$.

At the end of the sequence, the population is read out with state-dependent
fluorescence detection using an EMCCD camera, which distinguishes excitation
of the $S$ and $D$ manifolds. As the logical information is carried in the two
$S$-levels (cf. Fig.1(c)), the population in $\left|1\right>$ needs to be
transferred to $\left|e\right>$ before the measurement. This state-readout
does not differentiate between between $\left|1\right>$ and $\left|f\right>$.
We use the notation $P_{ij}$ to indicate the population in state
$\left|ij\right>$. We can separately measure the occupation of the motional
mode by applying a pulse on resonance with either the red or blue sideband of
one of the Strontium ion’s $5S_{1/2}\leftrightarrow 4D_{5/2}$ transitions, and
reading out its state [59]. The difference of the excitation probability of
the red and blue sideband excitations is used to infer the population in the
motional ground state.

Results.— We first demonstrate the central building block of resonance
engineering, the state-dependent population transfer, by showing its time-
evolution, using a detuning of $\Delta=\Omega_{\mathrm{SB}}/2$. The change in
population is shown for all four initial states. The intended behavior, Rabi
cycling from $\left|00\right>$ and no transfer from the other initial states,
is apparent in Fig. 3(a). The solid lines denote simulated data. The
simulations numerically solve the system’s master equation [56], and use
experimentally determined parameters described above, including the initial
phonon number and heating rate. The gray dashed line marks the duration of the
pulse with the maximum state transfer, $600\text{\,}\mathrm{\SIUnitSymbolMicro
s}$, where 82% of population has depleted from $\left|00\right>$, and only 16%
from $\left|01\right>$. The deviation from a full population transfer is
attributed to the non-zero initial phonon number and heating rate,
corroborated by the simulated results.

After this population transfer, the state
$(\ket{10}\ket{1}_{m}+\ket{f0}\ket{0}_{m})/\sqrt{2}$ should be dissipatively
transferred to $\ket{10}\ket{0}_{m}$. We demonstrate this process by showing
the evolution of the populations $P_{f0}$ and $P_{10}$ and the phonon ground
state occupation over time in Fig. 3(b). The electronic states are
differentiated by running the measurement twice: once with the transfer of
$\ket{1}$ to $\ket{f}$, and once without. The latter measurement does not
discriminate between $\left|0\right>$ and $\left|1\right>$. The population
$P_{10}$ is inferred from the difference between the first and second
measurement. We additionally measure the phonon occupation. The lines are
simulated results, using the same sideband coupling strength
$\Omega_{\mathrm{SB}}$ as in (a). A dissipation rate of
$\Gamma_{f}=$4.5(6)\text{\,}\mathrm{kHz}$$ is determined by a least-squares
fit between the simulated and measured results. After $1\text{\,}\mathrm{ms}$
of applying the dissipation pulse, approximately 80% of the population is in
$\ket{10}\ket{0}_{m}$. The population is trapped there because of the
irreversible nature of the dissipation.

Figure 3: (a) Experimental demonstration of state-dependent population
transfer transfer with $\Delta=\Omega_{\mathrm{SB}}/2$, shown for each
possible initial state. The lines indicate simulated results, which include
measured initial phonon number and heating rate as simulation parameters. (b)
Demonstration of dissipation, after a maximal probe transfer from
$\left|00\right>$ (at a time marked by the gray dashed line in (a) ). We show
the evolution of $P_{f0}$ and $P_{10}$, and use sideband thermometry on the
Strontium ion to infer the ground-state phonon occupation. The lines show
simulated results, in which the dissipation rate $\Gamma_{f}$ is obtained
through a least-squares fit between simulated and measured data.

Having demonstrated and characterized the engineered resonance and dissipation
processes, we apply these steps within the full pulse sequences shown in Fig
2(b) to perform the OR and NOR gates. Figure 4 shows the measured population
outcome for each of the four possible initial states for the OR and NOR gates.
Populations are determined from 50 experimental repetitions for each input
state. Both truth tables exhibit the intended gate behavior: for all input
states, the majority of the population is transferred to the desired state,
marked in the figure with dashed boxes. The initial states $\left|01\right>$
and $\left|00\right>$ are transferred following the engineered resonance
scheme for the OR and NOR gates, and have success rates of 84(5)% and 74(6)%.
As confirmed by simulations and analytics 222An analytic expression for the
error is calculated in the Supplemental material, Section B, the primary
source of error is attributed to a non-zero initial phonon number and the
heating rate. Since the coupling strength to a motional sideband is dependent
on the phonon number, the resonance condition of the engineered population
transfer is not met for $n\geq 1$.

Figure 4: Measured population truth tables of the OR and NOR gates, with the
intended output states marked with dashed lines. Values are in percent, and
are determined from 50 experimental shots for each input setting. The OR and
NOR gates have an average population fidelity of 87(5)% and 81(5)%.

Conclusion and Outlook.— We have implemented nonunitary multi-qubit operations
in a trapped-ion system by use of engineered dissipation. The schemes for the
OR and NOR gate performed the operations with average fidelities of 87(5)% and
81(5)%, respectively. This constitutes the first realization of dissipative
quantum gate operations. The leading source of error stems from heating and
thus imperfect cooling of the mode over which the intended dissipation process
occurs.

This heating process is the result of electronic noise on the trap surface,
and by phonon transfer from other uncooled motional modes caused by mode-
coupling, both of which are known challenges of microfabricated ion traps [61,
62], and is further exacerbated by complications involved in mixed-species
operation [58]. Such issues are technical, and do not pose fundamental
limitations: Future experiments could implement improved trap design and
manufacturing to reduce heating due to technical noise, and improved cooling
techniques such as polarization gradient cooling [63] and electromagnetically-
induced-transparancy (EIT) cooling [64]. Much like recent dissipative high-
fidelity schemes for entangled state preparation [48, 47, 46, 49] improved
upon the fidelities of their first-generation counterparts we would expect
future implementations to improve the fidelity.

Nonunitary operations are of relevance in a wide range of quantum information
algorithms. For example, in a NISQ context, quantum convolutional neural
networks use measurements and conditional feedback operations to process
information [38, 18]. These elements could be replaced by integrated
nonunitary operations, thereby avoiding classical measurements and
feedforward.

Regarding universal fault-tolerant quantum computation, quantum error
correction also constitutes a nonunitary process as multiple erroneous
processes are mapped to the same corrected state. Our work can be seen as a
stepping stone towards an implementation of autonomous quantum error
correction, in which erroneous states are coherently mapped to oscillator
excitations and are then removed through dissipation [51]. We have
demonstrated the required techniques, resonance engineering and sympathetic
cooling, in the present experiment.

## Acknowledgments

We gratefully acknowledge discussions with Jonathan Home and support by the EU
Quantum Technology Flagship grant AQTION under Grant Agreement number 820495,
and by the US Army Research Office through Grant No. W911NF-14-1-010 and
W911NF-21-1-0007. We also acknowledge funding by the Austrian Science Fund
(FWF), through the SFB BeyondC (FWF Project No. F7109), ERC-2020-STG 948893,
and by the IQI GmbH. The research is also based upon work supported by the
Office of the Director of National Intelligence (ODNI), Intelligence Advanced
Research Projects Activity (IARPA), via the US Army Research Office Grant No.
W911NF-16-1-0070. F.R. and E.Z. acknowledge funding from the Swiss National
Science Foundation (Ambizione grant no. PZ00P2 186040) and the ETH Research
Grant ETH-49 20-2.

### Author Contributions

E.Z. developed, under the guidance of F.R., the theoretical protocol. M.v.M,
F.R., P.S., and E.Z. designed the experiment. M.v.M. carried out the
experiment and analyzed the data. M.v.M., P.H., and L.G. contributed to the
experimental setup. R.B., T.M., F.R., P.S. supervised the project. All authors
contributed to the manuscript.

## References

  * noa [2015] _International Technology Roadmap for Semiconductors 2.0_, Tech. Rep. (ITRS, 2015).
  * Shulaker  _et al._ [2013] M. M. Shulaker, G. Hills, N. Patil, H. Wei, H.-Y. Chen, H.-S. P. Wong, and S. Mitra, Nature 501, 526 (2013).
  * Ma  _et al._ [2003] D. D. D. Ma, C. S. Lee, F. C. K. Au, S. Y. Tong, and S. T. Lee, Science 299, 1874 (2003).
  * Han  _et al._ [2007] M. Y. Han, B. Özyilmaz, Y. Zhang, and P. Kim, Phys. Rev. Lett. 98, 206805 (2007).
  * Sugahara and Tanaka [2004] S. Sugahara and M. Tanaka, Appl. Phys. Lett. 84, 2307 (2004).
  * Salahuddin and Datta [2008] S. Salahuddin and S. Datta, Nano Lett. 8, 405 (2008).
  * Chen  _et al._ [2008] F. Chen, H. Kam, D. Markovic, T.-J. K. Liu, V. Stojanovic, and E. Alon, in  _2008 IEEE/ACM International Conference on Computer-Aided Design_ (IEEE, San Jose, CA, USA, 2008) pp. 750–757.
  * Zhou  _et al._ [1997] C. Zhou, D. M. Newns, J. A. Misewich, and P. C. Pattnaik, Appl. Phys. Lett. 70, 598 (1997).
  * Qian  _et al._ [2014] X. Qian, J. Liu, L. Fu, and J. Li, Science 346, 1344 (2014).
  * Lee and Kim [2008] K.-S. Lee and S.-K. Kim, J. Appl. Phys. 104, 053909 (2008).
  * Andreakou  _et al._ [2014] P. Andreakou, S. V. Poltavtsev, J. R. Leonard, E. V. Calman, M. Remeika, Y. Y. Kuznetsova, L. V. Butov, J. Wilkes, M. Hanson, and A. C. Gossard, Appl. Phys. Lett. 104, 091101 (2014).
  * Feng  _et al._ [2005] M. Feng, N. Holonyak, G. Walter, and R. Chan, Appl. Phys. Lett. 87, 131103 (2005).
  * Currivan-Incorvia  _et al._ [2016] J. A. Currivan-Incorvia, S. Siddiqui, S. Dutta, E. R. Evarts, J. Zhang, D. Bono, C. A. Ross, and M. A. Baldo, Nat. Commun. 7, 10275 (2016).
  * Pino  _et al._ [2021] J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. S. Allman, C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. Neyenhuis, Nature 592, 209 (2021).
  * Kranzl  _et al._ [2022] F. Kranzl, M. K. Joshi, C. Maier, T. Brydges, J. Franke, R. Blatt, and C. F. Roos, Phys. Rev. A 105, 052426 (2022).
  * Egan  _et al._ [2021] L. Egan, D. M. Debroy, C. Noel, A. Risinger, D. Zhu, D. Biswas, M. Newman, M. Li, K. R. Brown, M. Cetina, and C. Monroe, Nature 598, 281 (2021).
  * Postler  _et al._ [2022] L. Postler, S. Heu$\beta$en, I. Pogorelov, M. Rispler, T. Feldker, M. Meth, C. D. Marciniak, R. Stricker, M. Ringbauer, R. Blatt, P. Schindler, M. Müller, and T. Monz, Nature 605, 675 (2022).
  * Herrmann  _et al._ [2022] J. Herrmann, S. M. Llima, A. Remm, P. Zapletal, N. A. McMahon, C. Scarato, F. Swiadek, C. K. Andersen, C. Hellings, S. Krinner, _et al._ , Nat. Commun. 13, 4144 (2022).
  * Chen  _et al._ [2021] Z. Chen, K. J. Satzinger, J. Atalaya, A. N. Korotkov, A. Dunsworth, D. Sank, C. Quintana, M. McEwen, R. Barends, P. V. Klimov, _et al._ , Nature 595, 383 (2021).
  * Olsacher  _et al._ [2020] T. Olsacher, L. Postler, P. Schindler, T. Monz, P. Zoller, and L. M. Sieberer, PRX Quantum 1, 020316 (2020).
  * Zhang  _et al._ [2020] C. Zhang, F. Pokorny, W. Li, G. Higgins, A. Pöschl, I. Lesanovsky, and M. Hennrich, Nature 580, 345 (2020).
  * Krinner  _et al._ [2022] S. Krinner, N. Lacroix, A. Remm, A. Di Paolo, E. Genois, C. Leroux, C. Hellings, S. Lazar, F. Swiadek, J. Herrmann, _et al._ , Nature 605, 669 (2022).
  * Hrmo  _et al._ [2022] P. Hrmo, B. Wilhelm, L. Gerster, M. W. van Mourik, M. Huber, R. Blatt, P. Schindler, T. Monz, and M. Ringbauer, arXiv:2206.04104 (2022).
  * Chu  _et al._ [2023] J. Chu, X. He, Y. Zhou, J. Yuan, L. Zhang, Q. Guo, Y. Hai, Z. Han, C.-K. Hu, W. Huang, _et al._ , Nat. Phys. 19, 126 (2023).
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Rev. Mod. Phys. 94, 015004 (2022).
  * Foldager  _et al._ [2022] J. Foldager, A. Pesah, and L. K. Hansen, Sci. Rep. 12, 3862 (2022).
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Phys. Rev. Lett. 123, 220502 (2019).
  * Zhu  _et al._ [2020] D. Zhu, S. Johri, N. M. Linke, K. A. Landsman, C. Huerta Alderete, N. H. Nguyen, A. Y. Matsuura, T. H. Hsieh, and C. Monroe, Proceedings of the National Academy of Sciences 117, 25402 (2020).
  * Wang  _et al._ [2021] Y. Wang, G. Li, and X. Wang, Phys. Rev. Appl 16, 054035 (2021).
  * Mazzola  _et al._ [2019] G. Mazzola, P. J. Ollitrault, P. K. Barkoutsos, and I. Tavernelli, Phys. Rev. Lett. 123, 130501 (2019).
  * Verdon  _et al._ [2019] G. Verdon, J. Marks, S. Nanda, S. Leichenauer, and J. Hidary, arXiv:1910.02071 (2019).
  * Wang  _et al._ [2011] H. Wang, S. Ashhab, and F. Nori, Phys. Rev. A 83, 062317 (2011).
  * Del Re  _et al._ [2020] L. Del Re, B. Rost, A. F. Kemper, and J. K. Freericks, Phys. Rev. B 102, 125112 (2020).
  * Hu  _et al._ [2020] Z. Hu, R. Xia, and S. Kais, Sci. Rep. 10, 3301 (2020).
  * Hu  _et al._ [2022] Z. Hu, K. Head-Marsden, D. A. Mazziotti, P. Narang, and S. Kais, Quantum 6, 726 (2022).
  * Leppäkangas  _et al._ [2022] J. Leppäkangas, N. Vogt, K. R. Fratus, K. Bark, J. A. Vaitkus, P. Stadler, J.-M. Reiner, S. Zanker, and M. Marthaler, arXiv:2210.12138 (2022).
  * Cong  _et al._ [2019] I. Cong, S. Choi, and M. D. Lukin, Nat. Phys. 15, 1273 (2019).
  * Harrington  _et al._ [2022] P. M. Harrington, E. J. Mueller, and K. W. Murch, Nature Reviews Physics 4, 660 (2022).
  * Kronwald  _et al._ [2014] A. Kronwald, F. Marquardt, and A. A. Clerk, New J. Phys. 16, 063058 (2014).
  * Agarwal and Huang [2016] G. S. Agarwal and S. Huang, Phys. Rev. A 93, 043844 (2016).
  * Reiter  _et al._ [2016] F. Reiter, D. Reeb, and A. S. Sørensen, Phys. Rev. Lett. 117, 040501 (2016).
  * Kastoryano  _et al._ [2011] M. J. Kastoryano, F. Reiter, and A. S. Sørensen, Phys. Rev. Lett. 106, 090502 (2011).
  * Lin  _et al._ [2013] Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler, A. S. Sørensen, D. Leibfried, and D. J. Wineland, Nature 504, 415 (2013).
  * Morigi  _et al._ [2015] G. Morigi, J. Eschner, C. Cormick, Y. Lin, D. Leibfried, and D. J. Wineland, Phys. Rev. Lett. 115, 200502 (2015).
  * Cole  _et al._ [2021] D. C. Cole, J. J. Wu, S. D. Erickson, P.-Y. Hou, A. C. Wilson, D. Leibfried, and F. Reiter, New J. Phys. 23, 073001 (2021).
  * Cole  _et al._ [2022] D. C. Cole, S. D. Erickson, G. Zarantonello, K. P. Horn, P.-Y. Hou, J. J. Wu, D. H. Slichter, F. Reiter, C. P. Koch, and D. Leibfried, Phys. Rev. Lett. 128, 080502 (2022).
  * Malinowski  _et al._ [2021] M. Malinowski, C. Zhang, V. Negnevitsky, I. Rojkov, F. Reiter, T.-L. Nguyen, M. Stadler, D. Kienzler, K. K. Mehta, and J. P. Home, Phys. Rev. Lett. 128, 080503 (2021).
  * Doucet  _et al._ [2020] E. Doucet, F. Reiter, L. Ranzani, and A. Kamal, Phys. Rev. Res. 2, 023370 (2020).
  * Barreiro  _et al._ [2011] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, and R. Blatt, Nature 470, 486 (2011).
  * Reiter  _et al._ [2017] F. Reiter, A. S. Sørensen, P. Zoller, and C. A. Muschik, Nat. Commun. 8, 1822 (2017).
  * Wang and Gertler [2019] C. Wang and J. M. Gertler, Phys. Rev. Res. 1, 033198 (2019).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Nat. Phys. 5, 633 (2009).
  * Zapusek  _et al._ [2023] E. Zapusek, A. Javadi, and F. Reiter, Quantum Sci. Technol. 8, 015001 (2023).
  * Brandl  _et al._ [2016] M. F. Brandl, M. W. van Mourik, L. Postler, A. Nolf, K. Lakhmanskiy, R. R. Paiva, S. Möller, N. Daniilidis, H. Häffner, V. Kaushal, _et al._ , Rev. Sci. Instrum. 87, 113103 (2016).
  * Note [1] The equations describing the dynamics of the system are given in the Supplemental material, Section A.
  * Jost  _et al._ [2009] J. D. Jost, J. P. Home, J. M. Amini, D. Hanneke, R. Ozeri, C. Langer, J. J. Bollinger, D. Leibfried, and D. J. Wineland, Nature 459, 683 (2009).
  * Home [2013] J. P. Home, in  _Advances In Atomic, Molecular, and Optical Physics_, Vol. 62 (Elsevier, 2013) pp. 231–277.
  * Brownnutt  _et al._ [2015] M. Brownnutt, M. Kumph, P. Rabl, and R. Blatt, Rev. Mod. Phys. 87, 1419 (2015).
  * Note [2] An analytic expression for the error is calculated in the Supplemental material, Section B.
  * Boldin  _et al._ [2018] I. A. Boldin, A. Kraft, and C. Wunderlich, Phys. Rev. Lett. 120, 023201 (2018).
  * Home  _et al._ [2011] J. P. Home, D. Hanneke, J. D. Jost, D. Leibfried, and D. J. Wineland, New J. Phys. 13, 073026 (2011).
  * Joshi  _et al._ [2020] M. K. Joshi, A. Fabre, C. Maier, T. Brydges, D. Kiesenhofer, H. Hainzer, R. Blatt, and C. F. Roos, New Journal of Physics 22, 103013 (2020).
  * Morigi  _et al._ [2000] G. Morigi, J. Eschner, and C. H. Keitel, Phys. Rev. Lett. 85, 4458 (2000).

Supplemental material for
_Experimental realization of nonunitary multi-qubit operations_

##  I System

The dynamics of our system are well described by a master equation of Lindblad
form. The Hamiltonian term $\hat{H}$ describes the coherent interaction while
the dissipative processes are described by the Lindblad jump operators
$\hat{L}_{k}$:

$\dot{\rho}=-i[\hat{H},\rho]+\sum_{k}\hat{L}_{k}\rho\hat{L}^{\dagger}_{k}-\frac{1}{2}(\hat{L}^{\dagger}_{k}\hat{L}_{k}\rho+\rho\hat{L}^{\dagger}_{k}\hat{L}_{k}).$
(S1)

We consider three subsystems: The electronic state of the first ion $k$, that
of the second ion $l$, and a common motional mode $m$ occupied with $n$
phonons, we use the notation
$\left|k,n\right>_{1}\otimes\left|l,n\right>_{2}=\left|kl\right>\left|n\right>_{m}$.
When referring to single-qubit couplings we write $\ket{a}_{j}\\!\bra{b}$
where $j$ indicates the addressed ion. The system Hamiltonian can be written
as:

$\hat{H}=\delta\hat{a}^{\dagger}\hat{a}+\hat{H}_{f,1}+\hat{H}_{f,2}+\hat{V}_{1}.$
(S2)

Here, $\hat{a}^{(\dagger)}$ are the annihilation (creation) operators of the
oscillator mode. The subscripts of the operators denote the qubit the
respective coupling acts on.

We apply a carrier drive, referred to as the probe, on the
$\ket{0}\leftrightarrow\ket{f}$ transition

$\hat{V}_{1}=\frac{\Omega_{f}}{2}(\ket{f}_{1}\\!\bra{0}+\ket{0}_{1}\\!\bra{f}).$
(S3)

On the $\ket{f}\leftrightarrow\ket{1}$ transition we apply a sideband drive

$\hat{H}_{f,j}=\frac{\Omega_{\mathrm{SB}}}{2}(\hat{a}\ket{f}_{j}\\!\bra{1}+\hat{a}^{\dagger}\ket{1}_{j}\\!\bra{f})+\Delta\ket{f}_{j}\\!\bra{f}.$
(S4)

The dissipative contribution used for both OR and NOR is sympathetic cooling.
It allows the excited state of the first qubit to decay into the state
$\ket{1}$, represented by the following Lindblad operator:

$\hat{L}_{\Gamma_{f}}=\sqrt{\Gamma_{f}}\hat{a}.$ (S5)

The NOR gate additionally uses decay over an auxiliary level:

$\hat{L}_{\Gamma_{e}}=\sqrt{\Gamma_{e}}\ket{0}_{j}\\!\bra{e}.$ (S6)

This coupling is present on both qubits.

##  II Analytic error analysis

In principle the logic gates presented here can operate at arbitrary low
error. However, due to experimental imperfections and limited coupling
strength we observe errors. We analyse the error processes for the OR gate and
derive an analytic approximation of the error to better understand these
imperfections.

The experiment operates in separate steps for excitation and decay. For this
analysis we will assume that the decay by sympathetic cooling operates
perfectly and that all population that was excited by the Rabi pulse decays to
the target state.

The main source of error is heating. It competes with sympathetic cooling to
raise the occupation number of the harmonic oscillator. After state
preparation the harmonic oscillator is in a thermal state with a mean phonon
number $\bar{n}=0.14$. For this low occupation number the probability to be in
Fock state $\ket{0}_{m}$ is $P_{0}=1/(1+\bar{n})$ and that for Fock state
$\ket{1}_{m}$ is $P_{1}\approx 1-P_{0}$.

In the first step population from $\ket{01}\ket{0}_{m}$ is transferred to the
dressed state
$\ket{\psi_{-}}=(\ket{f1}\ket{0}_{m}+\ket{1f}\ket{0}_{m}-\sqrt{2}\ket{11}\ket{1}_{m})/2$
by a resonant $\pi/2$ pulse. This dressed state is addressed by the effective
Rabi frequency $\Omega_{f}/2$. Therefore the pulse takes
$t_{\pi}=2\pi/\Omega_{f}$.

Figure S1: Level schemes of the error processes. (a) The error process
intrinsic to the scheme. Even in absence of any harmonic oscillator
excitations some population from $\ket{00}\ket{0}_{m}$ is off-resonantly
excited. By increasing the sideband frequency $\Omega_{\text{SB}}$ relative to
the other couplings this error could be suppressed. (b) If the harmonic
oscillator is in state $n=1$ the energies of the excited states addressed from
$\ket{00}\ket{1}_{m}$ are shifted by $\delta$ and the coupling between them is
enhanced. This shifts the probe into resonance. (c) The desired process in
presence of an harmonic oscillator excitation. The additional excitation
enables coupling to the double excited state $\ket{ff}\ket{0}_{m}$ and adjusts
the frequencies to $\Delta+\delta$. Furthermore, $\ket{01}\ket{1}_{m}$ forms
dressed states with $\ket{0f}\ket{0}_{m}$. The lowest frequency dressed state
resides at $\tilde{\Delta}_{01}$.

The same pulse that resonantly excites $\ket{01}\ket{0}_{m}$ also addresses
state $\ket{00}\ket{0}_{m}$, off-resonantly driving it to
$\ket{\phi_{-}}=(\ket{f0}\ket{0}_{m}-\ket{10}\ket{1}_{m})/\sqrt{2}$. It is
detuned by $\Delta_{d}=\Omega_{\text{SB}}(\sqrt{2}-1)/2$ and the Rabi
frequency it is addressed with is $\Omega_{d}=\Omega_{f}/\sqrt{2}$. Using
these values we can calculate the excitation probability just as we would for
a simple detuned Rabi pulse:

$P_{e}^{00}(n=0)=\frac{\Omega_{d}^{2}}{\Delta_{d}^{2}+\Omega_{d}^{2}}\sin{\left(\frac{1}{2}\sqrt{\Delta_{d}^{2}+\Omega_{d}^{2}}t_{\pi/2}\right)}\approx
17\%.$ (S7)

This process is shown in Fig. S1(a). If the system initially is in $n=1$ the
dressed states are resonant and addressed by Rabi frequency
$\Omega_{d}=\Omega_{f}/\sqrt{2}$ (see Fig. S1(b)). The total error probability
we observe is:

$P_{e}^{00}=P_{0}P_{e}^{00}(n=0)+P_{1}P_{e}^{00}(n=1)\approx 23\%.$ (S8)

In the experiment we measure an error of $P_{e}^{00}=14\%$.

Initial state |  Measured |  Numeric |  Analytic
---|---|---|---
$\ket{00}$ |  86 |  79 |  79
$\ket{01}$ |  84 |  86 |  82
Table S1: Measured fidelity compared to analytic approximations for input
states $\ket{00}$ and $\ket{01}$. The fidelity of 0.98 of the initialization
explains part of the difference between analytics and experiment. Initial
states $\ket{10}$ and $\ket{11}$ are not shown as their fidelity is close to
unity.

When the system is in initial state $\ket{01}\ket{1}_{m}$ it is driven to an
excited subspace consisting of states
$\ket{f1}\ket{1}_{m},\ket{11}\ket{2}_{m},\ket{1f}\ket{1}_{m}$ and
$\ket{ff}\ket{0}_{m}$ that couple to each other (cf. Fig. S1 (c)).
Furthermore, the sideband drive couples state $\ket{01}\ket{1}_{m}$ to
$\ket{0f}\ket{0}_{m}$ forming dressed states
$\ket{\psi_{\pm}}=(\ket{01}\ket{1}_{m}\pm\ket{0f}\ket{0}_{m})/\sqrt{2}$. State
$\ket{\psi_{-}}$ is close to resonance, driven with
$\tilde{\Delta}_{01-}=(1+\sqrt{2}-\sqrt{6})\Omega_{\text{SB}}/2$ and
$\Omega_{d}=\Omega_{f}(1/\sqrt{8}+1/\sqrt{12})$. Using Eq. (S7) one can find
the error probability for this initial state to be about $24\%$. State
$\ket{\psi_{+}}$ sits at detuning
$\tilde{\Delta}_{01+}=(-1+\sqrt{2})\Omega_{\text{SB}}/2$ with respect to the
excited state manifold and is addressed by Rabi frequency
$\Omega_{d}=\Omega_{f}/2$. Therefore it has an error probability of $89\%$.
Combining both error probabilities with the probability to start in an excited
state of the oscillator, the total error probability for initial state
$\ket{01}$ evaluates to

$P_{e}^{01}=P_{e}^{01}(n=1)(1-P_{0})\approx 7\%.$ (S9)

This value is lower than the experimental error of $P_{e}^{01}=16\%$. The
discrepancy is explained as follows:

We have independently measured that experimental imperfections lead to a Rabi
cycle decay that can be modeled by a Gaussian statistical spread in Rabi
frequency with a standard deviation of $4\%$. Including this spread in our
gate error models gives the results shown in Tab. S1.

For both initial states $\ket{00}$ and $\ket{01}$ we overestimate the
population that is transferred to the excited states. Here we only took the
initial occupation of the oscillator into account while in the experimental
setting heating also occurs during the $\pi/2$ pulse. Additionally, the state
initialization for states other than $\ket{00}$ has a fidelity of 98%. Another
source of error could be miscalibration of pulse timing. Furthermore, we do
not consider errors during the decay step. These effects combined can explain
the difference between the experiment and numerical as well as analytical
models. Such sources of error are technical in nature, and therefore not a
fundamental limitation. This provides a clear path how to improve the gates’
performance.


# Single-photon nonlinear optics with graphene plasmons

M. Gullans  Department of Physics, Harvard University, Cambridge, MA 02138,
USA  D. E. Chang  F. H. L. Koppens  ICFO-Institut de Ciencies Fotoniques,
Mediterranean Technology Park, 08860 Castelldefels (Barcelona), Spain  F. J.
García de Abajo  ICFO-Institut de Ciencies Fotoniques, Mediterranean
Technology Park, 08860 Castelldefels (Barcelona), Spain  ICREA - Instituci
Catalana de Recerca i Estudis Avan ats, Barcelona, Spain  M. D. Lukin
Department of Physics, Harvard University, Cambridge, MA 02138, USA

###### Abstract

We show that it is possible to realize significant nonlinear optical
interactions at the few photon level in graphene nanostructures. Our approach
takes advantage of the electric field enhancement associated with the strong
confinement of graphene plasmons and the large intrinsic nonlinearity of
graphene. Such a system could provide a powerful platform for quantum
nonlinear optical control of light. As an example, we consider an integrated
optical device that exploits this large nonlinearity to realize a single
photon switch.

###### pacs:

78.67.Wj, 73.20.Mf, 42.50.-p, 42.65.Hw

_Introduction_ \- Nonlinear optical processes find ubiquitous use in modern
scientific and technological applications, facilitating diverse phenomena like
optical modulation and switching, spectroscopy, and frequency conversion Boyd
(2003). A long-standing goal has been to realize nonlinear effects at
progressively lower powers, which is difficult given the small nonlinear
coefficients of bulk optical materials. The ultimate limit is that of single-
photon nonlinear optics, where individual photons strongly interact with each
other. Realization of such nonlinear processes would not only facilitate peak
performance of classical nonlinear devices, but also create a unique resource
for implementation of quantum networks Kimble (2008) and other applications
that rely on the generation and manipulation of non-classical light.

One approach to reach the quantum regime involves coupling the light to
individual quantum emitters to take advantage of their intrinsically nonlinear
electronic spectrum Duan and Monroe (2008); Kimble (2008). While a number of
remarkable phenomena have been demonstrated with these systems Haroche (2013),
their realization remains a challenging task. Specifically, in contrast to
conventional bulk nonlinear systems, coherent single quantum emitters are
generally unable to operate under ambient conditions, suffer from relatively
slow operating speeds, are prone to strong decoherence in solid-state
environments, and have limited tunability of their properties.

Fueled by these considerations, there has been renewed interest in nonlinear
optical materials that can reach the quantum regime Matsuda et al. (2009);
Mabuchi (2012); Ferretti and Gerace (2012). In particular, recent experiments
demonstrated the realization of a quantum nonlinear medium, featuring single
photon blockade Peyronel et al. (2012) and conditional two-photon phase shifts
Peyronel et al. (2013), in a cold, dense gas of strongly interacting atoms.
The essence of these approaches is that the interaction probability for two
photons becomes substantial if the photons are confined to a sufficiently
small mode volume of the nonlinear medium for sufficiently long times.
Motivated by these recent developments, in this Letter we explore the
potential for using nanoscale surface plasmon excitations in graphene for
quantum nonlinear optics. Graphene, a single atomic layer of carbon atoms, has
attracted tremendous interest for its unique electronic, mechanical, and
quantum transport properties Geim and Novoselov (2007); Castro Neto et al.
(2009); Nair et al. (2008). Recently it has also been realized that the unique
properties of graphene have a strong effect on the guided electromagnetic
surface waves in the form of surface plasmons (SPs) Wunsch et al. (2006);
Mikhailov and Ziegler (2007); Jablan et al. (2009). In particular, recent
theoretical Mikhailov and Ziegler (2007); Jablan et al. (2009); Koppens et al.
(2011); Nikitin et al. (2012a) and experimental Fei et al. (2012); Chen et al.
(2012); Fang et al. (2013) results indicate that graphene plasmons can be
confined to volumes millions of times smaller than in free space. We show that
under realistic conditions, this field confinement enables deterministic
interaction between two plasmons (i.e., photons) over picosecond time scales
as illustrated in Fig. 1ab, which is much shorter than the anticipated plasmon
lifetime Principi et al. (2013). We show how one can take advantage of this
interaction to realize a single photon switch (Fig. 1cd) and produce non-
classical light.

Figure 1: a) A doped graphene disk confines photons as plasmons to mode
volumes millions of times smaller than free space. b) This induces a large
dispersive nonlinearity $\eta$ (defined in Eq. 5) so that only a single photon
can resonantly excite the cavity. c) Integrated nonlinear optical circuit for
using the graphene plasmon cavity to realize a single photon switch. First the
photons are converted into planar plasmons of a graphene waveguide via a
grating, then they couple to the plasmon cavity, after which they are
converted back into waveguide photons. For the frequencies we consider the
waveguide and grating could be fabricated from etched Si. d) Top down view of
the plasmon cavity from (c) showing the width $W$ of the plasmon cavity, the
width $W^{\prime}$ the graphene nanoribbon, and the spacing $d$ between the
cavity and nanoribbon.

Through electrostatic gating, it is possible to introduce a net carrier
concentration, which shifts the Fermi energy $\hbar\omega_{F}$ away from the
Dirac point to a non-zero value. The in-plane conductivity of graphene is
well-approximated by the expression
$\sigma(\omega)\approx\frac{ie^{2}}{\pi\hbar}\frac{\omega_{F}}{\omega+i\gamma}$
at frequencies below twice the Fermi frequency $\omega<2\omega_{F}$ Falkovsky
(2008), which describes a Drude-like response of electrons within a single
band. In realistic systems the conductivity will also have a small term
$\gamma$ describing dissipation due impurity or phonon-mediated scattering.
There are two limits on the existence of low-loss SP modes in graphene. First,
at frequencies $\omega>2\omega_{F}$, graphene suffers from strong inter-band
absorption Jablan et al. (2009); Koppens et al. (2011). Second for frequencies
above the optical phonon frequency $\hbar\omega_{op}\approx 0.2$ eV, there is
additional loss due to scattering into optical phonons Jablan et al. (2009);
Yan et al. (2013); although narrow plasmons above $\omega_{op}$ have recently
been observed in graphene nanorings Fang et al. (2013). To minimize the losses
we focus on the regime where the frequencies fall below $2\omega_{F}$ and
$\omega_{op}$. In this regime, we can approximate
$\gamma=ev_{F}^{2}/\mu\,\hbar\omega_{F}$ where $\mu$ is the mobility Novoselov
et al. (2005). The ability to tune $\omega_{F}$, and consequently the optical
properties, through electrostatic gating makes graphene unique compared to
normal metals.

Like in noble-metal plasmonics Barnes et al. (2003), the free nature of charge
carriers described by the Drude response gives rise to SP modes in graphene
Mikhailov and Ziegler (2007); Jablan et al. (2009). At first order in
$k_{sp}/k_{F}$ the SP dispersion is given by

$\omega_{sp}^{2}=\frac{e^{2}\omega_{F}}{2\pi\epsilon_{0}\hbar}k_{sp}\approx
4.2\,\omega_{F}\,v_{F}\,k_{sp}$ (1)

where $v_{F}\approx 10^{6}$ m/s is the Fermi velocity Wunsch et al. (2006).
This dispersion relation implies a remarkable reduction of the SP wavelength
compared to the free space wavelength $\lambda_{0}=2\pi c/\omega_{sp}$, as
$\lambda_{sp}/\lambda_{0}\sim v_{F}/c\sim 3\cdot 10^{-3}$. Thus, the smallest
possible mode volume of a graphene SP resonator, $V\sim\lambda_{sp}^{3}$, can
be $\sim 10^{7}$ times smaller than in free space Koppens et al. (2011).

_Nonlinear plasmonics in graphene_ \- To describe the nonlinear properties of
the plasmons we employ the semiclassical Maxwell-Boltzmann equations (MBE).
This is a good approximation when the plasmon momentum is much less than the
Fermi momentum and the plasmon properties are dominated by intra-band
transitions. The distribution function $f(\bm{x},\bm{k},t)$ for an electron at
in-plane position $\bm{x}$ and with Bloch momentum $\bm{k}$ evolves under the
Maxwell-Boltzmann equation as

$\partial_{t}f+v_{F}\hat{k}\cdot\partial_{\bm{x}}f+e\partial_{\bm{x}}\varphi\cdot\partial_{\bm{k}}f=0,$
(2)

where the electostatic potential $\phi(\bm{x},z,t)$ satisfies Poisson’s
equation $\nabla^{2}\phi=en\delta(z)/\epsilon_{0}\epsilon$. Here $z$ is the
out-of-plane coordinate and $n=\int\,d\bm{k}\,f$ is the 2D electron density.
For weak excitations of the electron distribution, the term
$\partial_{\bm{k}}f$ in the Maxwell-Boltzmann equation can be replaced by the
equilibrium value $\partial_{\bm{k}}f^{(0)}$, yielding a linear equation
supporting SPs with the dispersion given in Eq. (1) and an electrostatic wave
given by $E=-\nabla\phi\propto\delta n\sin(kx-\omega t)$.

For sufficiently large density perturbations $\delta n$, the nonlinear
interaction between the non-equilibrium distribution $\partial_{\bm{k}}f$ and
potential must be accounted for. This effect can be interpreted as a
backaction induced by the electrostatic wave on the electrons via a
ponderomotive force $F_{p}\sim\partial_{x}E^{2}\sim k\delta n^{2}\sin 2kx$,
which grows with the amplitude of the SPs. This nonlinear force directly
excites a second plasmon wave at wavevector $2k$ and frequency $2\omega$, i.e.
second harmonic generation, and gives rise to the second order conductivity
calculated in Ref. Mikhailov (2011). We show (see supplementary material sup )
that this leads to a nonlinear shift at the original wavevector $k$ and
frequency $\omega$, with an effective third order conductivity for the SPs
given by

$\sigma^{(3)}(k_{sp},\omega)=-i\frac{3\pi}{4}\frac{v_{F}^{4}}{\omega_{F}^{3}}\frac{\epsilon_{0}^{2}}{\hbar\,\omega}.$
(3)

This result differs from the nonlinear conductivity as seen by free-space
light normally incident on a graphene sheet, where one finds that
$\sigma^{(3)}\sim 1/\omega^{3}$ Mikhailov and Ziegler (2008). Remarkably, as
we discuss next, the tight confinement of SPs in graphene implies that the
fields associated with even single quantized SPs are strong enough that
nonlinear effects are observable.

_Nonlinear Graphene Plasmon Cavity_ \- Anticipating the large strength of
nonlinear interactions at the level of single SPs in nanoscale graphene
resonators, we are motivated to introduce a quantum description of such a
system. We write the Hamiltonian as $H=H_{0}+H_{c}$, where $H_{0}$
characterizes the excitation spectrum of the graphene resonator, and $H_{c}$
describes an external coupling to the resonator (as in Fig. 1cd), which allows
one to probe the resonator properties or utilize the nonlinearities for
applications such as a single-photon transistor.

Figure 2: a) Nonlinear shift (calculated from Eq. 5) for the fundamental mode
relative to the plasmon linewidth with decreasing mode volume
$V_{0}=(\lambda_{sp}/\lambda_{0})^{3}$. Here we take the linewidth as
$\gamma=ev_{F}^{2}/\mu\,\hbar\omega_{F}$ with the Fermi energy
$\hbar\omega_{F}=0.2$ eV and a mobility of $\mu=10^{5}(10^{4})$ cm2/Vs
corresponding to a quality factor of roughly 600(60). b) $g^{(2)}(t)$ for the
graphene plasmon cavity driven by a weak coherent state for
$\hbar\omega_{sp}=0.2$ eV and two different mobility. $g^{(2)}(0)<1$ indicates
a transition to an effective two-level system as illustrated schematically in
Fig. 1ab.

We first consider the intrinsic properties of the resonator given by $H_{0}$.
Considering the fundamental SP mode of the resonator with corresponding
annihilation operator $a_{q}$ and number operator
$n_{q}=a^{\dagger}_{q}a_{q}$, the effective Hamiltonian $H_{0}$ is given by
Denardo and Putterman (1988); Gervasoni and Arista (2003); sup

$H_{0}=\left(\omega_{q}-i\kappa/2+\eta_{q}(n_{q}-1)\right)n_{q}.$ (4)

This Hamiltonian describes the quantum analog of a cavity exhibiting an
intensity-dependent refractive index, where the resonance frequency
$\omega_{q}+\eta_{q}(n_{q}-1)$ shifts depending on the intra-cavity photon
number. Here we have also included the total cavity linewidth
$\kappa=\kappa_{ex}+\gamma$ into the cavity description which includes the
intrinsic losses $\gamma$ and radiative losses of the cavity into other
optical or plasmonic modes, given by $\kappa_{ex}$. For graphene, the
nonlinear interaction strength is given by sup

$\eta_{q}=\frac{7\,\pi\,\omega_{q}}{64\,A\
k_{F}^{2}}\sqrt{\frac{q^{3}}{2\alpha_{g}k_{F}^{3}}},$ (5)

where $\alpha_{g}\equiv e^{2}/4\pi\epsilon_{0}\,\hbar\,v_{F}\approx 2$ and $A$
is the mode area of the resonator, which can be given by
$A=\lambda_{sp}^{2}/4$ for a diffraction-limited structure. The
$\eta_{q}\propto A^{-1}$ scaling reflects that the field intensity of a single
SP grows inversely with its confinement. For such small structures one might
expect that quantum size effects become important; however, as shown in Ref.
Thongrattanasiri et al. (2012), for graphene nanostructures larger than $\sim
10-20$ nm the use of the bulk dielectric response is a valid approximation.

At the quantum level, the interaction parameter $2\eta_{q}$ indicates the
additional energy cost to excite two versus one photon in the cavity, as can
be seen in the cavity excitation spectrum (Fig. 1b). When $2\eta\gg\kappa$,
the graphene sheet behaves as a two-level atom because it can only resonantly
absorb a single photon as illustrated in Fig. 1a; thus we describe this as the
quantum nonlinear regime. The ratio $2\eta_{q}/\kappa$ is then a good measure
of the quality of the cavity as a quantum emitter. Fig. 2b shows
$2\eta_{q}/\gamma$ for the fundamental mode with decreasing mode volume
(assuming mobilities of $10^{5}$ and $10^{4}$ cm2/Vs), where we see that this
ratio can be as large as 100. The parameter $\eta/\kappa\propto Q/A$, where
$Q$ is the quality factor of the resonator.

The enabling mechanism for a two-level atom to be useful for quantum
information processing is that it can only emit single photons at a time. This
can be characterized by the second order correlation function of the emitted
light, which is identical to that of the cavity mode, $g^{(2)}(t)=\langle
a^{\dagger}(\tau)a^{\dagger}(t+\tau)a(t+\tau))a(\tau)\rangle/\langle
a^{\dagger}(\tau)a(\tau)\rangle$. For a stationary process, $g^{(2)}(0)<1$
indicates non-classical “anti-bunching” and approaches $g^{(2)}(0)=0$ in the
limit of an ideal two-level emitter. We consider the case where the resonator
is driven by an external laser from the side and emission is collected from a
different direction. In the limit of weak driving we find that

$g^{(2)}(0)=\frac{\kappa^{2}}{4\eta^{2}+\kappa^{2}},$ (6)

thus establishing $\eta\lesssim\kappa$ as the regime where quantum properties
become observable. In Fig. 2b we take $\kappa_{ex}=0$ and we see that, for the
largest nonlinearities, $g^{(2)}<1$ can be readily observed for high mobility
graphene.

_Efficient coupling and a single-photon switch_ \- In order to exploit the
large nonlinearity of graphene, we need an efficient method to convert SPs
into external optical modes on time scales short compared to the intrinsic
losses. Specifically, one needs that the total linewidth
$\kappa=\kappa_{ex}+\gamma$ contains a large component $\kappa_{ex}$ that goes
into desirable external channels compared to the intrinsic losses $\gamma$.
One approach is to use the direct dipolar emission of the cavity into free
space radiation. For the square cavities described above, the dipole moment is
given by $p=2\,e\,k_{F}^{2}/k_{sp}^{3}$ which gives a decay rate into
radiation of

$\kappa_{ex}=\frac{k_{0}^{3}\,p^{2}}{3\pi\epsilon_{0}\hbar}=\frac{16\,\alpha_{g}}{3}\frac{k_{F}^{3}}{k_{sp}^{3}}V_{0}\,\omega_{F}$
(7)

where $V_{0}\equiv(\lambda_{sp}/\lambda_{0})^{3}$. For cavities in the quantum
nonlinear regime, this is a small contribution to the total losses; thus, a
more practical approach is needed for coupling the photons to free space.

Figure 3: a) Single photon transmission through the device depicted in Fig.
1cd. We take $\xi=k_{0}/2$ and $P\gg 1$ so the only losses are in the
nanoribbons. The plasmon frequency is 0.2 eV and we assume the decay rate
$\gamma$ is dominated by impurity scattering. The three curves are for a fixed
plasmon frequency with increasing Fermi energy, which increases the spatial
propagation length of the plasmons. b) Bunching in reflection for two incident
photons from the left with $\hbar\omega_{sp}=0.2$ eV, $E_{F}=0.23$ eV, $P=2$,
and mobility $\mu=10^{4}(10^{5})$ cm2/V$\cdot$s (dashed(solid)) corresponding
to a lifetime of $0.2(2)$ ps and a cavity quality factor of 60(600). e)
Antibunching in transmission for $P=0.1$ with other parameters as in (a).

We envision a two-step process illustrated in Fig. 1cd: first a waveguide
photon is converted into the planar plasmon of a graphene waveguide via a
dielectric grating, then this plasmon can tunnel directly into the nonlinear
cavity. We first consider the direct coupling between the cavity and the bulk
plasmons. We take the cavity of width $W$ to be separated a distance $d$ from
a long nanoribbon of width $W^{\prime}$ as shown in Fig. 1d. For $d\gg
W,\lambda_{sp}$ the coupling is dipolar and small, which allows us to
calculate the decay of the fundamental cavity mode into the nanoribbon via
Fermi’s golden rule sup

$\kappa_{c-r}=\frac{32}{\pi^{2}}\frac{k_{F}^{r}}{k_{F}^{c}}\frac{W\,\omega}{k_{sp}^{*}k_{sp}^{4}d^{6}}$
(8)

where $k_{F}^{r,c}$ is the Fermi wavevector in the nanoribbon ($r$) and cavity
($c$) and $k_{sp}^{*}$ is the wavevector for the nanoribbon plasmon that is
resonant with the cavity mode. The cavity can be efficiently controlled
through the nanoribbon by operating at a distance $d$ such that this decay is
the dominant loss channel for the cavity.

Once the plasmon is in the nanoribbon it still remains to out-couple it to the
waveguide. Due to the large mismatch in wavevectors, $k_{sp}/k_{0}\sim
c/v_{F}$, the bare coupling of the plasmons to the waveguide mode will be very
small. Achieving efficient conversion between conventional optics and plasmons
(in graphene and noble metals) is an area of active research with several
different approaches being pursued Bao and Loh (2012); Ruppert et al. (2010);
Zhan et al. (2012); Bludov et al. (2012); Ferreira and Peres (2012); Nikitin
et al. (2012b); Zhu et al. (2013). We consider the conceptually simplest
solution, which is to directly fabricate a dielectric grating to enable
momentum conservation, as recently demonstrated in Ref. Zhu et al. (2013). We
take a single-mode dielectric slab waveguide in vacuum coupled to a graphene
nanoribbon via the dielectric grating. For the frequencies we consider here
one could fabricate the waveguide and grating by etching Si. For parallel
propagation, the grating wavector $k_{g}$ should be given by
$k_{g}=k_{sp}-k_{0}$. This geometry can be analyzed via coupled mode theory
and optimized as a function of the slab thickness Snyder and Love (1983).
Taking the grating profile to be of the form
$\epsilon_{g}(x)=\delta\epsilon\cos k_{g}x$ gives the power conversion for
weak losses between the waveguide and plasmon mode as $\cos^{2}(\xi x)$ where
$\xi$ is spatial coupling between the TM mode of the waveguide and nanoribbon

$\xi\approx\sqrt{\frac{W}{W^{\prime}}}\,\delta\epsilon\,e^{-\gamma_{\perp}h}k_{0}$
(9)

Here $W^{\prime}>W$, $\gamma_{\perp}^{2}=\beta^{2}-k_{0}^{2}$ is the
transverse wavevector of the slab mode, $\beta$ is the longitudinal
wavevector, and $h$ is the distance between the slab and the graphene. Because
the factor in $\xi$ in front of $k_{0}$ is order unity, the plasmon conversion
for a weak grating is limited to distances $\sim\lambda_{0}\gg\lambda_{sp}$.
As a result the spatial decay rate of the plasmons must be much larger than
$k_{0}$ to achieve efficient conversion. When losses are dominated by impurity
scattering the spatial decay rate is given by
$\gamma\,k_{sp}/\omega_{sp}\approx
ev_{F}\,\hbar\,\omega_{sp}/2\,\mu\,E_{F}^{2}$, which decreases with Fermi
energy. Fig. 3a shows the transmission of a single photon through the geometry
displayed in Fig. 1cd.

The device depicted in Fig. 3ab can be used as a nonlinear single-photon
switch. To characterize this process, it is first necessary to understand how
an input field through the waveguide is transformed upon interacting with the
nonlinear resonator, which can be done through an input-output formalism. In
the case of Fig. 3ab of a resonator equally coupled to two waveguides, the
resonator evolves under the incoming fields of the left- and right-going modes
under the Hamiltonian
$H_{c}=\sqrt{\kappa_{ex}}(a^{r}_{in}+a^{l}_{in})a^{\dagger}+h.c.$, while the
output fields are given by
$a^{r(l)}_{out}=a^{r(l)}_{in}+i\sqrt{\kappa_{ex}}a$.

This one dimensional model has been solved exactly for the case of one and two
resonant photons input from a single direction in the waveguide Liao and Law
(2010). The response is characterized by the effective Purcell factor
$P=\kappa_{ex}/\gamma$, which measures the fraction of cavity emission into
the waveguide, and the normalized nonlinearity $\tilde{\eta}=\eta/\kappa$. The
transmission $t$ and reflection $r$ coefficients for a single photon incident
on resonance with the cavity are given by $t=-P/(1+P)$ and $r=1/(1+P)$. The
two photon response, however, is modified by the nonlinearity. For example two
photons at frequency $\omega_{sp}$ will be blocked from entering the cavity
due to the nonlinearity. This leads to antibunching in the transmission and
bunching in the reflection as shown in Figures 3bc. The suppression in the
transmission scales as $\tilde{\eta}^{2}$ similarly to Eq. 6, while the
bunching in reflection scales as $P^{4}$ for $\tilde{\eta}\gg P\gg 1$ Liao and
Law (2010). Fig. 3c shows that such a device realizes a single photon
transistor where one control photon can block several signal photons from
propagating through the cavity for a time given by the inverse cavity
lifetime.

Experimental tests of these ideas require single photon detectors in the mid-
infrared wavelength regime between 2 and 10 $\mu$m. While this is challenging
to realize, recent work on superconducting nanowire single photon detectors
(SNSPD) and transition edge sensors (TES) have achieved single photon
detection in this regime Eisaman et al. (2011); Marsili et al. (2012); Karasik
et al. (2012). Alternatively, frequency up-conversion may allow efficient
detection Eisaman et al. (2011).

Our analysis shows that graphene plasmonics may provide a powerful platform
for the nonlinear quantum optical control of light. Combined with the scalable
fabrication of graphene this could allow the creation of complex quantum
networks for many applications in quantum information and quantum simulation,
as well as in classical nonlinear optics. Such a system is ultimately limited
either by the losses in graphene or the strength of the nonlinearity. We
estimate currently achievable quality factors for the plasmon cavity range
from $10-10^{3}$; however, estimates of the ultimate limit to the graphene
plasmon lifetime suggest that quality factors greater than $10^{4}$ are
possible Principi et al. (2013). To enhance the nonlinearity further hybrid
structures can be envisioned if one could fabricate the structure on top of a
strong nonlinear substrate.

_Acknowledgments–_ We thank P. Zoller for useful discussions. DEC acknowledges
support from Fundacio Privada Cellex Barcelona. We acknowledge support from
the Harvard Center for Quantum Optics, NSF, ARO, CUA, AFOSR MURI and the
Packard Foundation.

## References

  * Boyd (2003) R. W. Boyd, _Nonlinear optics_ (Academic Press, San Diego, 2003).
  * Kimble (2008) H. J. Kimble, Nature 453, 1023 (2008).
  * Duan and Monroe (2008) L.-M. Duan and C. Monroe, in _Advances In Atomic, Molecular, and Optical Physics_ , edited by E. Arimondo, P. R. Berman, and C. C. Lin (Academic Press, 2008), vol. 55, pp. 419–463.
  * Haroche (2013) S. Haroche, Rev. Mod. Phys. 85, 1083 (2013).
  * Matsuda et al. (2009) N. Matsuda, R. Shimizu, Y. Mitsumori, H. Kosaka, and K. Edamatsu, Nat. Phot. 3, 95 (2009).
  * Mabuchi (2012) H. Mabuchi, Phys. Rev. A 85, 015806 (2012).
  * Ferretti and Gerace (2012) S. Ferretti and D. Gerace, Phys. Rev. B 85, 033303 (2012).
  * Peyronel et al. (2012) T. Peyronel, O. Firstenberg, Q.-Y. Liang, S. Hofferberth, A. V. Gorshkov, T. Pohl, M. D. Lukin, and V. Vuletić, Nature 488, 57 (2012).
  * Peyronel et al. (2013) T. Peyronel, O. Firstenberg, Q. Y. Liang, S. Hofferberth, A. V. Gorshkov, M. D. Lukin, and V. Vuletić, Nature, in press. (2013).
  * Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007).
  * Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Nair et al. (2008) R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Science 320, 1308 (2008).
  * Wunsch et al. (2006) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J. Phys. 8, 318 (2006).
  * Mikhailov and Ziegler (2007) S. A. Mikhailov and K. Ziegler, Phys. Rev. Lett. 99, 016803 (2007).
  * Jablan et al. (2009) M. Jablan, H. Buljan, and M. Soljačić, Phys. Rev. B 80, 245435 (2009).
  * Koppens et al. (2011) F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, Nano Lett. 11, 3370 (2011).
  * Nikitin et al. (2012a) A. Y. Nikitin, F. Guinea, F. J. Garcia-Vidal, and L. Martin-Moreno, Phys. Rev. B 85, 081405 (2012a).
  * Fei et al. (2012) Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, G. Dominguez, et al., Nature 487, 82 (2012).
  * Chen et al. (2012) J. Chen, M. Badioli, P. Alonso-González, S. Thongrattanasiri, F. Huth, J. Osmond, M. Spasenović, A. Centeno, A. Pesquera, P. Godignon, et al., Nature 487, 77 (2012).
  * Fang et al. (2013) Z. Fang, S. Thongrattanasiri, A. Schlather, Z. Liu, L. Ma, Y. Wang, P. M. Ajayan, P. Nordlander, N. J. Halas, and F. J. García de Abajo, ACS Nano 7, 2388 (2013).
  * Principi et al. (2013) A. Principi, G. Vignale, M. Carrega, and M. Polini, arXiv:1305.4666 (2013).
  * Falkovsky (2008) L. A. Falkovsky, J. Phys.: Conf. Ser. 129, 012004 (2008).
  * Yan et al. (2013) H. Yan, T. Low, W. Zhu, Y. Wu, M. Freitag, X. Li, F. Guinea, P. Avouris, and F. Xia, Nat. Phot. 7, 394 (2013).
  * Novoselov et al. (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Barnes et al. (2003) W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature 424, 824 (2003).
  * Mikhailov (2011) S. A. Mikhailov, Phys. Rev. B 84, 045432 (2011).
  * (27) See Supplemental Material for more detailed discussions of the nonlinear susceptibility, quantization, and nanoribbon-cavity coupling.
  * Mikhailov and Ziegler (2008) S. A. Mikhailov and K. Ziegler, J. Phys.: Condens. Matter 20, 384204 (2008).
  * Denardo and Putterman (1988) B. Denardo and S. Putterman, Phys. Rev. B 37, 3720 (1988).
  * Gervasoni and Arista (2003) J. L. Gervasoni and N. R. Arista, Phys. Rev. B 68, 235302 (2003).
  * Thongrattanasiri et al. (2012) S. Thongrattanasiri, A. Manjavacas, and F. J. García de Abajo, ACS Nano 6, 1766 (2012).
  * Bao and Loh (2012) Q. Bao and K. P. Loh, ACS Nano 6, 3677 (2012).
  * Ruppert et al. (2010) C. Ruppert, J. Neumann, J. B. Kinzel, H. J. Krenner, A. Wixforth, and M. Betz, Phys. Rev. B 82, 081416 (2010).
  * Zhan et al. (2012) T. R. Zhan, F. Y. Zhao, X. H. Hu, X. H. Liu, and J. Zi, Phys. Rev. B 86, 165416 (2012).
  * Bludov et al. (2012) Y. V. Bludov, N. M. R. Peres, and M. I. Vasilevskiy, Phys. Rev. B 85, 245409 (2012).
  * Ferreira and Peres (2012) A. Ferreira and N. M. R. Peres, Phys. Rev. B 86, 205401 (2012).
  * Nikitin et al. (2012b) A. Y. Nikitin, F. Guinea, and L. Martín-Moreno, Appl. Phys. Lett. 101, 151119 (2012b).
  * Zhu et al. (2013) X. Zhu, W. Yan, P. Uhd Jepsen, O. Hansen, N. Asger Mortensen, and S. Xiao, Appl. Phys. Lett. 102, 131101 (2013).
  * Snyder and Love (1983) A. W. Snyder and J. Love, _Optical Waveguide Theory_ (Springer, 1983), 1st ed.
  * Liao and Law (2010) J.-Q. Liao and C. K. Law, Phys. Rev. A 82, 053836 (2010).
  * Eisaman et al. (2011) M. D. Eisaman, J. Fan, A. Migdall, and S. V. Polyakov, Rev. Sci. Inst. 82, 071101 (2011).
  * Marsili et al. (2012) F. Marsili, F. Bellei, F. Najafi, A. E. Dane, E. A. Dauler, R. J. Molnar, and K. K. Berggren, Nano Lett. 12, 4799 (2012).
  * Karasik et al. (2012) B. S. Karasik, S. V. Pereverzev, A. Soibel, D. F. Santavicca, D. E. Prober, D. Olaya, and M. E. Gershenson, Appl. Phys. Lett. 101, 052601 (2012).
  * Fetter (1973) A. Fetter, Ann. Phys. N.Y. 81, 367 (1973).

##  Appendix A Supplemental Material

_Nonlinear conductivity_ \- The nonlinearity can be derived from the Boltzmann
equation for the 2D electron distribution function $f(\bm{x},\bm{p},t)$, which
is a function of two space $\bm{x}=(x,y)$ and two momentum
$\bm{p}=(p_{x},p_{y})$ variables and time $t$, and Poisson’s equation for the
electric potential $\varphi(\bm{x},z,t)$

$\displaystyle\partial_{t}f+v_{F}\hat{p}\cdot$
$\displaystyle\partial_{\bm{x}}f+e\partial_{\bm{x}}\varphi\cdot\partial_{\bm{p}}f=0$
(S.1) $\displaystyle(\partial_{\bm{x}}^{2}+\partial_{z}^{2})\varphi$
$\displaystyle=e\,n\,\delta(z)/\epsilon_{0}$ (S.2)

where the 2D electron density is defined as $n\equiv\int d\bm{k}f$,
$\partial_{\bm{x}}=\partial_{x}\hat{x}+\partial_{y}\hat{y}$, and
$\partial_{\bm{p}}=\partial_{p_{x}}\hat{x}+\partial_{p_{y}}\hat{y}$. Eq. S1 is
written for electron doping, for hole doping the equations and final results
would be the same after taking the opposite electron charge. These equations
are an extension of the Drude model to describe intraband transitions. Taking
$x$ to be the propagation direction, the nonlinear equations for the electron
density $n$ and current density $n\,\bm{\bar{v}}=\int
d\bm{k}\,v_{F}\hat{k}\,f$ can be derived as

$\displaystyle\partial_{t}n$ $\displaystyle+\partial_{x}n\,\bar{v}=0$ (S.3)
$\displaystyle\partial_{t}\bar{v}$
$\displaystyle-\frac{e}{m^{*}}\partial_{x}\varphi+\bar{v}\,\partial_{x}\bar{v}+\frac{3e}{2m^{*}n_{0}}\partial_{x}\varphi\,\delta
n=0$ (S.4)

where the effective mass for the plasmon excitations is $m^{*}=\hbar
k_{F}/v_{F}$, $\delta n=n-n_{0}$, and $n_{0}\equiv k_{F}^{2}/\pi$ is the
equilibrium electron density. Linearizing these equations around $n_{0}$ and
$\bar{v}=0$ gives the plasmon dispersion from Eq. 1 in the main text Fetter
(1973). The nonlinearity is described by the last two terms in Eq. SS.4, where
the second term, $\propto\partial\varphi\delta n$, arises from the linear band
structure in graphene and is absent for electrons with a parabolic dispersion.
To find the nonlinear conductivity, we can expand $\bar{v}$ and $n$ in spatial
Fourier components and solve the resulting coupled equations together with
Poisson’s equation Denardo and Putterman (1988). This allows us to express
$\sigma^{3}(\omega)$ through the identity
$en\bar{v}=\sigma(\omega)E+\sigma^{3}(\omega)E^{3}$.

To solve for the nonlinear shift in the cavity, we use the boundary condition
that the current density perpendicular to the edge is zero at the edge. This
allows us to represent $\bar{\bm{v}}=\hat{x}\sum_{p}v_{p}\sin px$ and
$n=\sum_{q}n_{q}\cos qx$ where $q=mk_{sp}$ for some integer $m$. Inserting
this solution into Eq. SS.3-S.4 leads to coupled nonlinear equations for
$n_{q}$ and $v_{q}$

$\displaystyle\sum_{p}$ $\displaystyle\sin
px\bigg{(}\dot{v}_{p}-\frac{\omega_{p}^{2}}{n_{0}p}n_{p}\bigg{)}=\frac{1}{2}\sum_{p,q}[pv_{p}v_{q}\sin(p-q)x$
$\displaystyle-$
$\displaystyle\bigg{(}pv_{p}v_{q}-\frac{3}{2}\frac{\omega_{p}^{2}}{p}\frac{n_{p}n_{q}}{n_{0}^{2}}\bigg{)}\sin(p+q)x]$
(S.5) $\displaystyle\sum_{p}$ $\displaystyle\cos
px(\dot{n}_{p}+n_{0}pv_{p})=\frac{1}{2}\sum_{p,q}n_{p}v_{q}\big{[}(p-q)\cos(p-q)x$
$\displaystyle-$ $\displaystyle(p+q)\cos(p+q)x\big{]}$ (S.6)

where $\omega_{p}^{2}=\frac{e^{2}\omega_{F}}{2\pi\epsilon_{0}\hbar^{2}}p$ .
These equations can be solved in perturbation theory to find the nonlinear
frequency shift of the plasmon resonance as

$\delta\omega_{p}=\frac{\omega_{p}}{16}\frac{\delta
n^{2}}{n_{0}^{2}}\bigg{(}\frac{5\omega_{p}^{2}+\omega_{2p}^{2}}{4\omega_{p}^{2}-\omega_{2p}^{2}}-\frac{2\omega_{p}^{2}}{\omega_{2p}^{2}}\bigg{)}\geq\frac{5}{32}\omega_{p}\frac{\delta
n^{2}}{n_{0}^{2}}$ (S.7)

_Quantizing the plasmon mode_ \- To quantize the plasmon mode we use the
Hamiltonian Gervasoni and Arista (2003)

$\begin{split}H&=\frac{1}{2}\int d\bm{x}\,e\,\delta n\,\varphi+\frac{1}{2}\int
d\bm{x}n_{0}m^{*}{\bar{v}}^{2}\\\
&=\frac{A\,m^{*}}{4\,n_{0}}\sum_{q}\frac{1}{q^{2}}(\omega_{q}^{2}\,\delta
n_{q}^{2}+\delta\dot{n}_{q}^{2})\end{split}$ (S.8)

where $A=\pi^{2}/k_{sp}^{2}$ is the area of the graphene sheet and we used the
relation $\bar{v}_{q}=-\delta\dot{n}_{q}/q\,n_{0}$ from the continuity
equation. This Hamiltonian can be quantized in the usual way by defining
$\delta n_{q}=\frac{\gamma_{q}}{2\omega_{q}}(a_{q}+a_{q}^{\dagger})$ for
bosonic operators $a_{q}$ such that $\dot{a}_{q}=-i\omega_{q}a_{q}$ and
$\gamma_{q}=2q\sqrt{\omega_{q}\omega_{F}/\pi A}$. This leads directly to Eq. 4
in the main text.

_Coupling between nanoribbon and cavity_ \- To calculate the coupling between
the cavity and the proximal nanoribbion we use the electric potential of the
nanoribbon plasmons acting on the graphene cavity

$\begin{split}\varphi_{r}(x)&=\frac{1}{4\pi\epsilon_{0}}\sum_{k}\int
d\bm{x}^{\prime}\frac{e\,n^{r}_{k}\,\cos kx^{\prime}}{\lvert
x+d-x^{\prime}\lvert}\\\
&\approx\frac{W_{r}}{4\pi\epsilon_{0}}\sum_{k}\frac{e\,n^{r}_{k}}{k^{2}(x+d)^{2}}\end{split}$
(S.9)

where we assumed $d\gg W,\lambda_{sp}$. Inserting this into Eq. SS.4 gives the
coupling between each plasmon mode $k$ in the nanoribbon with the plasmon mode
$q$ of the cavity as

$\kappa_{kq}=\frac{8}{\pi}\sqrt{\frac{k_{F}^{r}\,W_{r}}{k_{F}^{c}L}}\frac{\omega_{q}^{c\,2}}{\omega^{r}_{k}+\omega^{c}_{q}}\frac{1}{q^{2}k\,d^{3}}$
(S.10)

where $L$ is the length of the nanoribbon, $k_{F}^{r,c}$ is the Fermi
wavevector and $\omega_{k}^{r,c}$ is the dispersion of the ribbon($r$) and
cavity($c$). Applying Fermi’s Golden rule gives the decay rate of the cavity
mode into the nanoribbon plasmons given in Eq. 8.