Graphene as a source of entangled plasmons

We analyze nonlinear optics schemes for generating pairs of quantum entangled plasmons in the terahertz-infrared range in graphene. We predict that high plasmonic field concentration and strong optical nonlinearity of monolayer graphene enables pair-generation rates much higher than those of conventional photonic sources. The first scheme we study is spontaneous parametric down conversion in a graphene nanoribbon. In this second-order nonlinear process a plasmon excited by an external pump splits into a pair of plasmons, of half the original frequency each, emitted in opposite directions. The conversion is activated by applying a dc electric field that induces a density gradient or a current across the ribbon. Another scheme is degenerate four-wave mixing where the counter-propagating plasmons are emitted at the pump frequency. This third-order nonlinear process does not require a symmetry-breaking dc field. We suggest nano-optical experiments for measuring position-momentum entanglement of the emitted plasmon pairs. We estimate the critical pump fields at which the plasmon generation rates exceed their dissipation, leading to parametric instabilities.

We propose that a graphene ribbon containing an electrostatically induced p-n junction [5,[42][43][44][45][46], see Fig. 1(a), can be a highly efficient PDC source. The plasmon spectrum of such a system, which we will call the device of type A, consists of multiple continuously dispersing subbands. The lowest subband (labeled "0") is gapless and the next one (labeled "1") is gapped due to transverse confinement in the ribbon, see Fig. 1(b). In the simplest picture, the PDC process is splitting of a mode-1 plasmon of frequency ω 1 and momentum q 1 = 0 into two mode-0 plasmons of frequency ω 0 = ω 1 /2 and momenta ±q 0 that propagate away in opposite directions along the ribbon. The energy-momentum conservation (phase matching) in this conversion is guaranteed without any fine-tuning of the device. This advantage of the counterpropagating PDC scheme has been previously exploited in photonic waveguides [36,39,47,48]. In a more precise description, frequencies ω 0 and so the momenta q of the outgoing plasmons have some uncertainty. The spontanenous PDC generates a superposition of states of different momentum: |ψ = dqf (q) |−q |+q (1) where |−q |+q means a two plasmon state with one plasmon at momentum −q and the other at +q, see Sec. II D. This superposition is not factorizable, which implies that the pairs are position-momentum entangled, similar to particles in the original Einstein, Podolsky, and Rosen (EPR) paper [34]. The amplitude function f (q) in Eq. (1) is peaked at q 0 and has a characteristic width δq ∼ γ/v 0 where is the plasmon damping rate [9] and v 0 = dω 0 /dq 0 is the mode-0 plasmon group velocity. The role of the split-gate represented by the golden bars in Fig. 1(a) is two-fold. First, it induces a gradient of carrier concentration in graphene when a dc voltage is applied between the two parts of the gate. The p-n junction is created when this voltage is high enough [44]. Such a system lacks inversion symmetry and therefore the PDC is allowed. Second, the split-gate serves as an optical antenna amplifying the incident pump field [44,49] that excites mode-1 plasmons. The coupling is facilitated by a large transverse dipole moment of the mode-1 plasmons and a strong field enhancement λ v /λ p 1 of the antenna. (Color online) (a) Type-A device: a graphene ribbon with a lateral p-n junction (blue/gold denoting electron/hole doping) induced by a split gate (golden bars). A pump field of frequency ω1 applied in theŷ-direction (thick arrows) generates plasmon pairs of frequency ω0 propagating in the ±x direction (red undulating arrows). (b) Plasmon dispersion in a ribbon of half-width a = 200 nm with the p-n junction on the midline and the carrier density n = ±10 12 cm −2 at the edges. The dotted lines form the 'light cone' for the dispersion of vacuum photons. The blue dot represents the mode-1 plasmon pumped by the incident photon, which then decays into two mode-0 plasmons (the two red dots) as shown by the arrows. (c) The spatial distribution of the amplitude of the quasi-static electric potential of the two modes.
Note that the inversion symmetry can also be broken by a dc current across a uniformly doped ribbon, which we refer to as the device of type B. In contrast, the FWM is a third-order nonlinear effect which does not require symmetry breaking fields. The degenerate FWM, in which the two pump photons of frequency ω are converted into a pair entangled plasmons of the same frequency, can take place in either a nanoribbon or a more typical larger-area graphene sheet, which are devices of type C.
Plasmonic effects similar to those we explore here have been considered in a few recent theoretical papers. For example, plasmonic sum-frequency generation in graphene nanoflakes [17] is based on the nonlinear mechanism similar to that underlying our plasmonic PDC. However, since the plasmon spectrum of a flake is discrete, only flakes of certain shapes can fulfill the energy-momentum conservation constraint for the sum (or difference) frequency generation. As mentioned above, in a long ribbon such kinematic constraints can be met without delicate fine-tuning. Another related work studied photon-plasmon difference frequency generation [26,50] and entanglement [51], which is a phenomenon intermediate between the usual all-photon PDC [34] and our fully plasmonic PDC.
The remainder of the paper is organized as follows. In Sec. II, we summarize our main results for the type-A device. In the last part of that section we discuss possible experiments that can probe the plasmon entanglement. In Sec. III we consider the third-order nonlinear effects. The FWM relevant for the operation of type-C devices is analyzed in Sec. III B and the current induced PDC in type-B devices is examined in Sec. III C. Section IV contains discussion and outlook. Finally, the summary of the notations and some details of the derivations are given in the Appendix.

II. PLASMONIC PDC IN A GRAPHENE RIBBON
To simplify the analysis of a type-A device, we assume that the edges of the ribbon y = ±a are not too close to those of the split-gate, so that the dc electric field created by the gate is approximately constant across the ribbon. The induced density response of graphene can be modeled [42] assuming the local chemical potential is linear in y: If the chemical potentials µ t and µ b at the two edges of the ribbon are opposite in sign, a p-n junction forms. If µ t = −µ b , the junction is located on the midline y = 0 of the ribbon. The plasmon frequency dispersion computed numerically as a function of momentum q along the strip is shown in Fig. 1(b) (For related analytical results, see Ref. 52). The lowest-frequency branch, which we refer to as mode 0, is gapless: Here is a characteristic scale of the plasmon velocity, v F ≈ c/300 is the Fermi velocity, α g = e 2 /(hv F ) is the dimensionless strength of the Coulomb interaction, and k F denotes the Fermi momentum corresponding to the maximum absolute chemical potential max(|µ t |, |µ b |) in the ribbon. Dimensionless parameter ξ v ∼ 1 is a slow (logarithmic) function of q. It also depends on the doping profile of the ribbon. The dispersion law (4) is similar to that of one dimensional (1D) plasmons [6,10,13]. The first gapped mode, mode 1, can be viewed as a linear combination of 2D plasmons of momenta qx ± q yŷ with q y ∼ 1/a. The frequency of this mode at q = 0 is where ξ 1 ∼ 1 is another dimensionless shape factor. The electric potential profile of the ω 1 mode [ Fig. 1(c)] indicates that mode-1 excitations have a nonzero dipole moment in the y-direction. Thus, they can be resonantly excited by a y-polarized pump field [53].
In this paper we describe the mode-1 to mode-0 PDC in the quantum language. However, there is a complementary classical picture, which is as follows. Plasmonic oscillations of mode-1 modulate the total carrier density of the ribbon viewed as a one dimensional conductor. This in turn modulates the frequencies of mode-0 plasmons and causes their parametric excitation. In the standard theory of optical parametric amplification [33] (see also Appendix B), one of the generated mode-0 plasmon is referred to as the idler and the other as the signal. The idler is assumed to have a finite spectral power, e.g., due to equilibrium fluctuations in the beginning of the process. The signal is then produced from the pump and the idler by the difference frequency generation (DFG).
For weak pump fields E, the pair generation rate grows linearly with the pump power, see Fig. 2. As E becomes stronger, a parametric instability is reached at some critical E. At even stronger pump field, the signal exhibits an exponential growth in the undepleted pump approximation. In Sec. II A below we evaluate the dependence of the pair generation rate and the critical pump field on experimental parameters, such as the ribbon width and its doping profile.

A. PDC Hamiltonian
The second-order nonlinear conductivity [16,23,51,54] leads to interaction between the plasmons of modes 0 and 1. This interaction can be modeled by the Hamiltonian H = H 0 + H (2) : where a 1 , a q are the annihilation operators for mode-1 at zero momentum and mode-0 of momentum q, respectively. The interaction Hamiltonian H (2) leads to pair generation in the weak coupling regime (g (2) a 1 Q/ω 1) and two-mode squeezing in the strong coupling regime (g (2) a 1 Q/ω > 1), as shown in Fig. 2. Since we study the PDC on resonance, ω(q 0 ) = ω 0 = ω 1 /2, we can treat the interaction strength g (2) as a constant, g (2) = g (2) q0 . As shown in Appendix C, Transition from stable to unstable regime of plasmonpair generation in the type-A device as pump field is increased. The black curve is the rate of generation of plasmon pairs per unit length. The red curve is the relative growth rate (κ − γ)/ω0 in the unstable regime. The plasmon damping rate is assumed to be γ = 0.1ω0. Other parameters are from Fig. 1 and the temperature is T = 300 K.
where L is the length of the gated part of the ribbon ( Fig. 1) and dimensionless factor 0 ≤ ξ g ≤ 1 depends on the carrier density profile across the ribbon.
The Hamiltonian H (2) governs the spontantous decay of the mode-1 plasmon into a pair of mode-0 plasmons with momenta q and −q. In the weak-coupling regime, the decay rate Γ can be calculated using Fermi's golden rule (see Appendix C 5) where n 0 is the occupation number of mode-0 plasmons of frequency ω ω 0 ; in thermal equilibrium, n 0 = 1/(eh ω0/T − 1) at temperature T . Normalizing Γ to the plasmon frequency ω 0 , we obtain the dimensionless decay rate of mode-1 plasmons: where ξ df ∼ 1 is another dimensionless factor. The dependence of ξ df on the carrier density profile across the ribbon is plotted in Fig. 3. For a ribbon with the carrier density n = ±10 12 cm −2 at the edges and a half-width of a = 200 nm, we find ω 1 /(2π) ≈ 5 THz and Γ 0 /ω 1 ∼ 10 −8 . Therefore, due to PDC alone a mode-1 plasmon decays into a mode-0 plasmon pair every Γ −1 0 ∼ 10 −5 s. If the plasmon Q-factor due to other damping channels (phonon and impurity scattering [9]) is Q = ω 0 /γ ∼ 10, the PDC efficiency is Γ 0 /γ ∼ 10 −7 , which is much higher than in available photonic PDC devices [39,55]. Furthermore, this efficiency scales as a −4 , and so can be higher still in a narrower strip.

B. Pumping mode-1
In a uniform external ac electric field E in the ydirection, the gapped mode-1, as a harmonic oscillator, FIG. 3. The shape factors for type-A device as functions of the doping profile parameter s d = (µ b + µt)/(µ b − µt) at temperature T = 300 K. The green, black and red curves are the decay rate shape factor ξ df , the occupation number shape factor ξ β and the shape factor ξQ for the relative subharmonic growth rate (defined in Appendix A). Each colored stripe below the horizontal axis illustrates the corresponding graphene ribbon with blue/gold representing electron/hole doping.
is driven to a coherent state |β = e βa † 1 −(β 2 /2) |0 . Specifically, the magnitude of the charge density profile of the mode is related to the external pump field as where σ(ω) is the two dimensional optical conductivity and α/S is an O(1) constant determined by the doping profile of the ribbon, see Appendix C 1. On resonance, from its relation to ρ (see Appendices C 2 and C 4), β is found to be where ξ β ∼ 1 is a shape factor shown in Fig. 3. The occupation of the gapped mode modifies the pair-generation rate to For n = 10 12 cm −2 , a = 200 nm, L = 1 µm, E = 10 3 V/cm and a plasmon quality factor of Q ∼ 10 for ultra clean samples at room temperature [5,7,9,56], we have β 0 = 15, indicating an average occupation number of N ∼ 400, thus boosting the generation rate to R ∼ 4 × 10 9 s −1 . It means that in a nanosecond-long pulse, there are four pairs of subharmonic plasmons generated while there are ∼ 10 6 photons incident onto the ribbon area. This rate is higher than conventional photonic [39] and quantum dot [38] devices. As the pump field increases, Fermi's golden rule ceases to be valid. The pair generation rate can be derived from the classical parametric oscillator theory (Appendix B) The growth rate κ is to be discussed in the next section. The pair generation rate diverges at the instability threshold κ = γ, as shown in Fig. 2.

C. Two mode squeezing and parametric amplification
The coupling Hamiltonian Eq. (7) in the interaction picture reads H I = q g (2) q e −iδt a 1 a † q a † −q +c.c. where δ = ω 1 −2ω q is the frequency mismatch. We focus on the pair of acoustic modes ω(q 1 ) = ω 0 = ω 1 /2 for which H I is on resonance: When mode-1 is pumped into a coherent state, one can replace a 1 by its classical value β. This Hamiltonian generates a two-mode-squeezing time evolution operator [57] S(βg (2) which squeezes the pair of mode-0 plasmons at momenta (q, −q). Therefore, the plasmons are generated as twomode squeezed states [34]. The squeezing operator leads to exponential growth of the amplitudes of the observables (e.g., the electric field of the mode) with the growth rate κ = g (2) β/h. We define the dimensionless relative growth rate that is made only of intensive quantities. Here ζ is the dimensionless small parameter that controls the second order nonlinear effects of graphene plasmons. Note that the relative growth rate depends only on the maximum doping level, the pump electric field, and the dimensionless shape factor ξ Q . Parameter ξ Q ∼ 1, plotted in Fig. 3, depends only on the scale invariant doping profile of the ribbon. The pump induced growth rate κ reduces the plasmon damping rate from γ to γ − κ, thereby further enhancing the pair generation rate. The system reaches parametric instability at κ = γ, see Fig. 2. Beyond this threshold, the plasmons are unstable with exponential growth rate κ − γ. When amplifying a classical source as the seed, this effect is known as parametric amplification [33]. For n = 10 12 cm −2 and E = 10 4 V/cm, we have κ = 0.12ω 0 in the p-n junction regime, high enough to compete with the damping rate γ.
The split-gate structure, as an antenna with width in y larger than the vacuum wave length λ v of the incident light, could enhance the pump field E v to the total field E by a factor of roughly F = E/E v = λ v /(2πw s ) [44,49] where w s is the distance separating the two halves of the gate, see Fig. 1. For the typical parameters in Fig. 1, one has λ v = 43 µm at ω 1 = 7 THz and w s could be chosen as 0.5 µm, rendering the field enhancement factor F ≈ 14.
(To increase the working frequency to, e.g., ω 1 = 30 THz corresponding to λ v = 10 µm, the ribbon width needs to be shrinked to 2a ∼ 20 nm.) Therefore, to achieve the plasmon instability regime in Fig. 2, the actual incident field just needs to be of the order of E v = 10 2 -10 3 V/cm. We note that the antenna also amplifies the radiative damping rate γ R of the dipolar active modes (e.g., mode-1) of the ribbon by the same factor F . For ribbons of size much smaller than the vacuum wavelength, a, L λ v , the normalized radiative damping rate Q −1 R = γ R /ω 0 is amplified by the antenna from Q −1 R ∼ a 2 L/λ 3 v to aL/λ 2 v , which is still much smaller than γ/ω 0 ∼ 0.1 [5,7,9,56] for the parameters given above. Therefore, the intrinsic damping γ (e.g., due to phonons and the electronic system itself) should be the major plasmon dissipation pathway in the system and the effect of the radiative damping can be neglected.

D. Measuring the entanglement
In general, the plasmon pair state can be expressed as |ψ = dp 1 dp 2 f (p 1 , p 2 )|p 1 |p 2 .
To quantify the standard deviations, we approximate the wave function by the Gaussian form where Z is a normalization factor. The EPR state in Eq. (1) corresponds to ∆ + = 0. After p 1 (or p 2 ) is measured, the standard deviation of p 2 (or p 1 ) becomes . Similarly, after measuring either position, the standard deviation of the unmeasured position becomes ∆x = ∆ 2 The apparent violation of the uncertainty principle happens whenever ∆ + = ∆ − . In the extreme case of Eq. (1), we have ∆x∆p = 0. If momentum is conserved in the type-A device in Fig. 1, we have p 1 = −p 2 exactly satisfied. However, the length L of the device limits momentum conservation and leads to ∆ + ∼ 1/L. The damping rate limits the plasmon line width and renders ∆ − ∼ 1/l where l = v 0 /γ is the propagation length of the mode-0 plasmons. For device satisfying L l, we have ∆x∆p ∼ L/l 1/2. Note that in our case, entanglement is between continuous variables, and so the fidelity of the state can not be defined in the conventional way [58]. Instead, it is more natural to use Eq. (20) to quantify the entanglement. Alternatively, one could use entanglement witnesses or entanglement negativity [35] which are more technically involved.
To detect the position-momentum entanglement of a plasmon pair, one can place grating structures on both sides of the device which would convert the plasmons to far-field photons. Passing through the beam splitters, the photons are read either by single photon detectors [34,59] which have good time precision or those with good energy resolution, as shown in Fig. 4. The coincidence rate measured by the former as a function of path length difference gives a time uncertainty ∆t, while the spectrum FIG. 6. Near-field imaging of the parametric amplification of the mode-0 plasmons. The right figure is plotted for a type-A device using the parameters in Fig. 1. The parametric plasmon amplifier can also be a type-B device in Fig. 8.
correlation of the energy detectors gives the energy uncertainty ∆E. The EPR entanglement is manifested in the relation ∆t∆E = ∆x∆p < 1/2 [34,60,61]. This 'energy time entanglement' detection scheme is the same as that in Ref. [61] with the entangled photon source replaced by the plasmonic type-A device in Fig. 1 or type-B device in Fig. 8. Alternatively, homodyne detection experiments can measure the entanglement property of the squeezed state [34].
In addition, the Franson scheme [40,41,62] is also applicable as shown in Fig. 5. As mentioned in Sec. I, near field technique based on scanning probes has been successful in probing plasmons in graphene [2][3][4][5][6][7][8][9][10][11][12][13]. To implement the Franson scheme, one needs two scanning probes, as shown in Fig. 5(a). The plasmons can either travel to the scanning probes directly or after being reflected by the edges of the ribbon, which we call short and long paths respectively. The coincidence detection arises from interference between the short-short and long-long paths. If the coincidence rate between the signals from two near field tips is measured as a function of L 1 ≈ L 2 , the distance between each tip and the edge, it will oscillate as R c ∼ cos 2 ((L 1 + L 2 )/λ), like that in Fig. 5(b). However, the condition L L i l needs to be satisfied. The first inequality guarantees that there are no single photon (or 'short-long' path) interference contributions, and the second inequality ensures that plasmon damping is not important. Note that in a prior work on measuring plasmon entanglement [63], the interfering particles were actually photons in optical fibers. In the proposed scheme (Fig. 5) the interference occurs between plasmons. Additionally, our implementation involves highly confined terahertz-infrared plasmons and scanning probes operating on deeply subdiffractional length scales.
One can also probe the classical effects. If the pump field is comparable to the parametric instability threshold, the enhancement of plasmon lifetime can be detected. Upon pumping of the ω 1 mode of the type-A device in Fig. 1 (or type-B in Fig. 8), the ω 0 modes with momenta q and −q are parametrically driven such that their effective damping rate is reduced to γ − κ (see Appendix B), and their quality factor is boosted to ω 0 /(γ − κ), as shown in Fig. 6. In a scanning near-field experiment, its manifestation is simply increased propagation length of the plasmons. One can also measure near field DFG using classical interference between the signal and idler in these experiments, which is discussed in Appendix F.

III. PLASMON GENERATION BY THE THIRD-ORDER NONLINEARITY
In this section we study the effects of third order nonlinearity [15,19,21,24,64,65] on plasmon interactions which occurs already in the long wavelength limit. This nonlinearity originates from the linearity of the dispersion of Dirac quasiparticles in graphene which breaks Galilean invariance [24]. At low frequency ω ε F , interband effects can be neglected and the third order nonlinear conductivity σ (3) in graphene assumes the third order 'Drude' form both in the kinetic [66] and hydrodynamic [24] regimes where ω i are the frequencies of the three electric fields that generate the third order current together. The symmetric tensor ∆ ilmn is the sum of all isotropic tensors ∆ ilmn = δ il δ mn + δ im δ ln + δ in δ lm . Close to zero temperature (T µ), the third order optical weight is predicted to be D (3) = 1 24π e 4 v F h 3 k F in the kinetic regime and twice as large in the hydrodynamic regime [24]. In what follows, we assume the system is the kinetic regime since we are mostly interested in the mid infrared plasmons whose frequencies are larger than the typical electron-electron scattering rate [23,24].

A. FWM Hamiltonian
We start with the plasmonic effective Lagrangian of an electronic system whose linear conductivity is in the dissipation-less Drude form: σ jk = iD/π ω δ jk where D is the Drude weight. In the near field approximation (meaning there is only longitudinal electric field and no magnetic field) and in the gauge A µ = (0, A) for the electric field where A is the vector potential viewed as the plasmonic field, the Lagrangian is where c is the speed of light. The first term is the electric field energy and the second term has the interpretation of the center of mass kinetic energy of the charge carriers. For two dimensional (2D) Drude conductors modeled as the x-y plane embedded in 3D space, since the current is localized on the 2D plane, the second term of the top line in Eq. (22) is nonzero only on the plane. In the second line of Eq. (22), A q =qA q are the Fourier components of A evaluated on the 2D plane. The 2/|q| comes from integrating over the fields exponentially decaying into the three dimensional (3D) space. Note that the summation is over q = (q x , q y ), and instead of three, there is only one plasmon mode for each q since the field is constrained to be longitudinal. The resulting Euler-Lagrange equation of motion for A q yields the 2D plasmon dispersion leads to the Hamiltonian where S is the area of the sample. The H (3) = −L (3) contains products of four plasmonic creation and annihilation operators. It describes FWM which is caused by the third order optical conductivity σ (3) . Writing the third order current j (3) in terms of vector potentials, Eq. (21) indicates j

and the interaction Hamiltonian
where the kernel Π ilmn = D (3) ∆ ilmn is perfectly local in space and time, and the summation runs over all q i constrained by q 1 + q 2 + q 3 + q 4 = 0 due to momentum conservation. Note that although Eq. (25) is a negative φ 4 term in this bosonic field theory, the system is stable due to higher order nonlinear couplings.

B. Spontaneous four wave mixing
This subsection discusses spontaneous FWM shown in Fig. 7 in uniformly doped graphene. Eq. (25) implies that entangled plasmon pairs with frequency ω can be generated by incident light of the same frequency. We model far field photon incident on the graphene plane by an uniform AC electric field with amplitude E 0 = ωA 0 inŷ direction, and represent it by the vector potential A 0 =ŷA 0 (e −iωt + e iωt ) where we have assumed A 0 to be real without loss of generality. When two fields in the product are replaced by the source A 0 , Eq. (25) becomes where the pair generation strength is and θ q is the angle of the plasmon momentum q relative to thex axis. We have defined the dimensionless small parameter which controls the strength of FWM and other third order effects in graphene. In Eqs. (24) and (26), it is enough to consider terms close to resonance in the interaction picture. The number conserving terms in Eq. (26) like a q a † q imply the field induced Kerr shift [33] of the plasmon frequency, which are red shifts in graphene due to its negative third order conductivity [24]. The terms like e −2iωt a † q a † −q lead to pair generation in the weak interaction regime (g (3) Q/ω 1 where Q = ω/γ is the plasmon quality factor) and two mode squeezing (instability) in the strong interaction regime (g (3) Q/ω > 1).
In the weak interaction regime, the plasmons with exactly the same frequency as ω are generated as entangled pairs and with an angular distribution of 1 + 2 sin 2 θ q 2 , as shown in Fig. 7(a). The pair generation rate is determined by Fermi's golden rule (analogous to Eq. (9)): where θ means angular average, the second part comes from the 2D plasmon density of states and we have assumed ω T so the relevant plasmons are not occupied in equilibrium. The resulting dimensionless generation rate is For n = 10 12 cm −2 , E 0 = 10 3 V/cm and S = 1 µm 2 , we have Γ/ω = 3.2 × 10 −8 . This leads to superior generation rates compared with conventional FWM sources, as we will discuss in Sec. IV. Since there is no need to break the inversion symmetry, this effect could be observed simply using uniformly doped graphene, which we call the type-C device. One may also use the ribbon in Fig. 1(a) but a lateral p-n junction is no longer necessary. The spontaneous FWM process in a transversely pumped ribbon is illustrated in Fig. 7(b). To distinguish the generated entangled plasmon pairs from the linearly excited plasmons, the key experimental signature would be satisfaction of the EPR criterion [34] in the quadrature measurement in Fig. 4 or the nonlocal interference phenomenon in the Franson experiment in Fig. 5.
In the strong interaction regime, the plasmons are squeezed in pairs and experience parametric amplification. The resulting relative exponential growth rate is (similar to Eq. (17)) for the the plasmons propagating in θ q direction, which agrees with the classical result (see Appendix E). In a typical near field experiment shown in Fig. 7(b), this reduces the plasmon damping rate from γ to γ − g q /h. Therefore, to overcome the typical damping rate γ ∼ 0.1ω of graphene plasmons, one needs ξ ∼ 0.3 which requires a pump field of E 0 ∼ 10 5 V/cm for n = 10 12 cm −2 and ω = 30 THz.

C. Spontanenous PDC in current carrying graphene
The third order nonlinearity can also lead to PDC in a current-carrying graphene, as shown in Fig. 8. This effect can also be understood as FWM with one "pump wave" in Eq. (25) replaced by a dc current flow (ω = 0). If the other pump source has frequency 2ω, then the process looks like a PDC: 2ω → ω + ω. Another point of view is that the dc current with flow velocity u breaks the inversion symmetry, resulting in a nonzero second order nonlinear conductivity σ (2) To leading order in the dc current J = enu, the equations from subsection III B can be directly applied to the present case. We assume that the pump field is polarized parallel to the dc current, such that the PDC Hamiltonian is the same as Eqs. (26) and (27) the dimensionless parameter due to the dc electric field E dc that drives the dc current.
In the weak pump regime, entangled subharmonic plasmon pairs at frequency ω are generated with the angular distribution of 1 + 2 sin 2 θ q 2 in an infinite graphene sheet, same as that of spontaneous FWM. Summed over all angles, the total relative generation rate is the same as Eq. (30) but with with ξ 4 replaced by ξ 2 2ω ξ 2 0 . To measure the entanglement property of the emitted pairs, one can use the graphene ribbon in Fig. 8(b) as a plasmon pair source to perform either the energy-time entanglement measurement in Fig. 4 or the Franson scheme in Fig. 5. In the ribbon, the generated subharmonic pairs are not in every direction, but can only propagate alonĝ x (similar to Fig. 1(b)) with the generation rate described by roughly the same formula as that for an infinite sheet. For n = 10 12 cm −2 , ω 1 = 2ω = 30 THz, E 0 = 10 3 V/cm, a ribbon size of S = 0.2 µm × 5 µm and a dc current of J = 0.16 mA/µm corresponding to u = 0.1v F , one has ξ 0 ≈ 0.1 and a normalized generation rate of Γ/ω ≈ 0.9 × 10 −6 from this type-B source.
In the strong pump regime, the plasmon lifetime is enhanced. The pump induced relative growth rate of plasmons of frequency ω is similar to Eq. (31). HereX means the unit vector along X andq is the unit vector of the momentum of the subharmonic plasmon. For a current of J = 0.5 mA/µm at the doping level of n = 10 12 cm −2 and under a pump field of 3 × 10 5 V/m at ω 1 = 2ω = 30 THz polarized alongŷ, the subharmonic plasmon running alongŷ has a relative growth rate of Q −1 g ≈ 0.03, while those alongx has a growth rate of Q −1 g ≈ 0.01, comparable to the plasmon loss at low temperatures [9].

D. Field enhancement in a ribbon
From previous discussion, high field is essential to make the proposed devices practical. By design, these devices naturally enhance the incident field by two mechanisms, helping to achieve this goal. First, as discussed in Sec. II C(a), the split gate for type-A device in Fig. 1, or the contacts for type-B device in Fig. 8(b), act as antennas. They enhance the pump field by a factor of [49] for the channel width w s = 200 nm and the vacuum wave length of the pump light λ v ≈ 10 µm at ω 1 = 30 THz. Second, close to the resonance at ω 1 with the gapped mode of the ribbon (see, e.g., Fig. 1(b) and Fig. 7(b)), there is another field enhancement factor of F Q = ξ f where ξ f is an O(1) shape factor, see Eqs. (11) and (C1). On resonance, this leads to a field enhancement factor of F Q ∼ Q ∼ 10.
The antenna effect together with the resonance field enhancement gives a net enhancement factor of QF ∼ 100. Therefore, the spontaneous PDC rate in the type-B device (Fig. 8) is enhanced by a factor of (QF ) 2 while the spontaneous FWM rate in a ribbon is enhanced by (QF ) 4 . Similarly, to reach a pump field of 10 5 V/m for parametric instability, the external incident field just needs to be of the order of 10 3 V/m. Note that the resonant field enhancement also occurs in photonic cavities such as micro-ring resonators [40,[67][68][69][70]. In our devices, the vacuum photons play the role of the waveguide photons in the bus waveguide therein, and the graphene ribbon plays the role of the ring resonator. Unlike the plasmon modes in the nano-ribbon whose line width is dominated by intrinsic damping, the field enhancement on resonance in the micro-rings scales as √ Q instead of Q [68] because the line width γ p of the photon modes in the micro-ring resonator is dominated by radiation leakage into the bus waveguide. As a result, the enhancement of the FWM pair generation rate due to resonant pumping scales as Q 2 , see, e.g., Eq. 4 of Ref. [67]. Note also that in both our devices (short devices with length L 0 l) and photonic cavities, the resonance conditions with the generated modes contributes another factor of Q to the enhancement of the generation rate.

IV. DISCUSSION AND EXPERIMENTAL OUTLOOK
We proposed several graphene spontaneous PDC and FWM devices which have promising quantum information applications [79] as efficient sources of entangled plasmon pairs, photon-plasmon converters and plasmon amplifiers [80][81][82][83][84], etc., for future quantum plasmonic/polaritonic circuits. These plasmons are in the THz-Infrared spectra range, which is an important energy range in condensed matter physics where superconductivity, magnetism, spin liquid, fractional quantum hall effect and other correlated phenomena are relevant. The quantum sources proposed here may enable quantum imaging [85] and sensing [86] of these novel states of matter.
Let us compare our devices with conventional photonic PDC sources. In photonic PDC, the nonlinearity is usually characterized by a nonlinear susceptibility χ (2) , which has the unit of inverse electric field. We can convert the 2D conductivity σ (2) to an effective 3D nonlinear susceptibility χ (2) as follows. The 2D nonlinear susceptibility (in the Gaussian unit system) is [87] (2) /ω 0 . For the type-A device which is a ribbon with half width a, the extent of the electric field of the plasmon in the out-of-plane direction is 1/q 0 ∼ a. Hence, the effective 3D nonlinear susceptibility is π × 10 6 cm −1 corresponding to carrier density n = 10 12 cm −2 , we find χ (2) ∼ 2 × 10 −6 cm/V, which is much larger than that of conventional photonic media, e.g., LiNbO 3 [41,73], but is of the same order as χ (2) in sophisticated photonic waveguides made of coupled quantum wells [47,71,72], see Table I. Nevertheless, a more important figure of merit is the dimensionless decay factor (same as Eq. (10) noting that g (2) ∼ χ (2) ω 3/2 0 a −1 L −1/2 ) of a waveguide mode into subharmonic modes due to spontanenous PDC. Due to strong mode confinement (a c/ω 0 ), this figure of merit is much larger than those of photonic devices. Physically, this is because electric field of a single plasmon in our device is much larger than that of a single photon.
As discussed in Sec. III C, a comparable χ (2) can also be achieved in a type-B device where inversion symmetry is broken by applying a DC current flow to an uniformly doped ribbon if the flow velocity is moderately fast: u > 0.1v F .
We also consider the FWM [33] in Sec. III B, a thirdorder nonlinear effect. Such a process does not require inversion-symmetry breaking and can take place in either a ribbon or a plain graphene sheet (type-C device in Fig. 7). This plasmon instability requires a relatively strong incident field controlled by the small parameter χ (3) , but it may be easier to realize experi-

PDC Devices
Frequency Q χ (2) ( cm/V) Generation Rate ( s −1 mW −1 ) Graphene type-A device ( Fig. 1(a)) 4 THz 10 − 100 2 × 10 −6 9 × 10 11 Graphene type-B device ( Fig. 8(b) TABLE I. Comparison of graphene with conventional platforms for entangled photon generation following that in Ref. [40]. Although plasmons in these graphene devices can work in a broad frequency range from terahertz to infrared, we picked specific frequencies for estimations in the table. For graphene type-A device, the plasmon frequency is shown in Fig. 1 mentally since there is no specific requirement for the nano structure. Compared with other nonlinear materials, χ (3) = iσ (3) /(ωa) in graphene is much stronger, as shown in Table I. A brief discussion of units conversion is in order. The Gaussian unit system defines the 3D nonlinear susceptibilities in terms of the effective dielectric function = 1 + 4π χ + χ (2) E + χ (3) E 2 + ... while the SI system defines them as / 0 = 1 + χ SI + χ SI E 2 + ... [33] where 0 is the vacuum permittivity. Therefore, besides converting the units of electric fields, there is an additional 4π factor: χ (n) = χ (n) SI /(4π). The unit m 2 /W is also used in the literature [40,75] for third order nonlinearity N 3 defined as δN = 1 2N χ is the linear refraction index, δN is the change of effective refraction index due to third order nonlinearity and the pump field, and P is the power flux of the pump field. Therefore, besides converting m 2 /W → 0.00132806 m 2 /V 2 , one also needs the conversion χ (3) = N 2π N 3 to obtain the χ (3) in Gaussian unit.
We note that the properties of the type-A device may have substantial temperature dependence arising from the middle of the ribbon [ Fig. 1(a)] where the system is close to charge neutrality. The gapped mode is an anti-symmetric charge oscillation between the upper and lower half of the strip, and thus strongly depends on the conductivity at the middle of the junction where the carriers are thermally excited electrons and holes. Therefore, the frequency of the gapped mode (the shape factor ξ 1 in Sec. II) increases with temperature. In our derivation, we included only the intraband contributions [20,23,51,54] to the second order nonlinear conductivity [Eq. (C20)]. Close to the middle of the junction where µ < ∼h ω, interband contributions may be important. There the total σ (2) is suppressed by a factor of (µ/hω) M with M ≥ 2 compared to the intraband expressions used in this paper, but σ (2) may diverge at the narrow region where µ =hω at the interband threshold for Pauli blocking [20,51,54,[87][88][89]. This divergence is rounded by nonzero temperature, and would thus lead to substantial temperature dependence if included. We note that even in homogeneously doped graphene, there are 3 ∼ 4 orders of magnitude inconsistencies [90] between theories (e.g., Refs. [20,51,54,87,88]) and experiments [25,26]. At this stage, it is premature to include such details in the model, and instead, we provide order of magnitude scaling results for the PDC. Incidentally, experiment and theory seem to agree for σ (3) [65].
As mentioned in Sec. I, generation of plasmons by PDC in graphene has also been considered by Ref. [51]. In that study, the pump and idler are far field photons. In contrast, here all three modes are plasmons such that frequencies are much lower and momenta are much larger, making this process much more efficient. Moreover, our proposed device involves an antenna that further increases the coupling efficiency to far field pump. As a result, the threshold power to achieve parametric instability is about four orders of magnitude lower than Ref. [51].
Since plasmons do not have the polarization degree of freedom (they are always longitudinally polarized), en-tanglement of polarization frequently discussed in quantum optics does not apply here. In plasmonics one focuses on the position-momentum (or energy-time) entanglement. For example, it has been shown this entanglement survives photon-plasmon conversion processes [31,63,91]. Here we proposed a scheme to directly generate and measure the entanglement of plasmon pairs on a chip using the modern tools of near field optics.
The following estimate shows that graphene will not be damaged by the strong electric field required. The incident light heats up the electron gas in graphene with the heating power P heat = 2Re[σ]|E| 2 and electron phonon scattering provides its cooling channel. The cooling power is approximated by P cool = C v γ c (T e − T l ) where T e (T l ) is the electron(lattice) temperature, C v ∼ g(ε F )T is the heat capacity of the electron gas and γ c ∼ 1 ps −1 is the cooling rate revealed by previous experiment [7]. The lattice is a good heat conductor and is assumed to stay at the same temperature as the environment. The balancing of heating and cooling P heat = P cool determines the stationary electron temperature T e −T l ∼ 1 4π 2 γ γc ε 2 F T ξ 2 where we have set the Boltzmann constant to one. At the typical doping level and frequency scale considered in this paper, for T l = 300 K and E = 10 4 V/cm, we conclude that T e −T l ∼ 6 K which is quite small compared to either room temperature or the Fermi energy.
We worked in the kinetic regime where plasmon frequency is much larger than the electron-electron scattering rate ω Γ ee [23]. In the low frequency hydrodynamic regime ω Γ ee [23,81,[92][93][94][95][96][97][98], the nonlinear conductivities are similar in magnitude as long as the temperature is not much larger than chemical potential [23,24,99]. Therefore, our estimate of the generation rates should apply as well to the collective modes in the hydrodynamic regime such as the charged 'demons' [23, 81, 92-94, 96, 97]. However, in the case of T µ, the demons become almost charge neutral, driven mainly by kinematic pressure, and the physics could be quite different. It may be interesting to study nonlinear and quantum effects for these collective modes [24,100].
Since the ribbon width is much larger than the lattice constant, the Dirac dispersion approximation for the electrons in graphene are valid. We used the 'local' approximation to compute the plasmonic modes, which neglects dependence of the optical conductivity σ(ω, q) on the field wave vector q, see Appendix C 1. The nonzero q effects enter σ(ω, q) as a power expansion in the dimensionless small number v 2 F q 2 /ω 2 [8,11,23], see, e.g., the supporting information of Refs. [11,23]. For the relevant modes of the ribbons such as mode-1 in Fig. 1, one has q ≈ 2π/(4a) and v F q/ω ≈ 0.36, meaning that the devices are indeed in the local regime. As an estimate, the leading order nonlocal correction δω 1 to the plasmon frequency ω 1 is just δω1 ω1 = 1 2 3 4 ( v F q ω ) 2 ≈ 0.05 in the kinetic regime we are considering. Therefore, the Drude (or 'local') approximation is well justified.
Their definitions for the case of a ribbon are as follows: The remaining ones, S, α and |ρ| 2 , are just the corresponding dimensionless ones multiplied by aL, one half of the total area of the ribbon.

Appendix B: Parametric oscillator
The spontanenous PDC process in Sec. II is similar to what happens in a parametric oscillator (PO) except that the former leads to two-mode squeezing while the latter corresponds to single-mode squeezing [57]. Optical POs are widely used to generate subharmonic light and to create entangled photon pairs. To make the paper self-contained, we review the basic theory of the PO in this Appendix.
A PO can be modeled with the equation of motion where Ω is the frequency of the pump modulating the natural frequency ω 0 of the oscillator by a relative amount ∼ δ. The difference of Ω from the primary parametric resonance frequency 2ω 0 is referred to as the detuning. Given the initial conditions x(0),ẋ(0) ≡ ∂ t x(t)| t=0 , the solution of Eq. (B1) can be written in terms of the retarded Green's function G = G(t, t ): where Θ(t) is the Heaviside step-function. One simple example is x(0) = 0,ẋ(0) = 1, F (t) ≡ 0, in which case x(t) = G(t, 0). As usual, the Green's function is constructed from an arbitrary pair of linearly independent solutions x 1 (t), x 2 (t) of the homogeneous equationLx(t) = 0: By virtue of the Abel formula, the Wronskian determinant [Eq. (B5)] has the following t-dependence: One possible choice of x 1 (t) and x 2 (t) is so that an instability, i.e., the exponential growth of G(t, 0) is possible if S(a, b, z) is increasing faster than e 1 2 γt . This occurs if |δ| exceeds some threshold value δ c = δ c (γ, ∆). Conversely, the response function G(t, 0) exhibits exponential decay with time if |δ| < δ c . Examples of such stable and unstable behaviors are shown in Fig. 9.
Approximate analytical calculation of δ c can be done for low damping γ ω 0 and small detuning ∆ ω 0 . Thus, one can show that in agreement with Fig. 9. The derivation goes as follows. Let us seek a solution of Eq. (B1) in the form where A(t), B(t) are slowly varying. Matching the coefficients for the rapidly oscillating phase factors of Eq. (B1), we obtain To eliminate the factors e ±i∆t , we change variables to arriving at

(B14) The general solution of this equation is
where are the growth rates and are the corresponding eigenvectors. Two linearly independent solutions for x(t) could be chosen as which leads to Eq. (B10). A comparison with the exact Green's function confirms the validity of this approximation for small damping and detuning, such as γ = 0.1ω 0 and ∆ = 0 in Fig. 9.
Consider now the response of the PO to a periodic driving source In the stable case, |δ| < δ c , where the effect of the initial conditions disappears at long times, x ω is the linear response to F : Without the pump, G(t, t − τ ) is a function of τ only, and so the integral on the right-hand side reduces to a constant. Under pumping, G(t, t−τ ) is periodic in t with the period 2π/Ω. Accordingly, the Fourier spectrum of x ω (t) becomes a frequency comb of Floquet harmonics ω + mΩ where m is an arbitrary integer. For small γ and ∆, however, it is easier to compute x ω (t) not from Eq. (B21) but from Eq. (B14) supplemented with the appropriate right-hand side. The result is Apparently, this approximation accounts only for the dominant Floquet harmonics: ω−Ω, ω, and ω+Ω. In the complex plane of ω, function x ω (or the linear response kernel) has two poles, at ν 1 and ν 2 , given by The pole ν 2 has a larger imaginary part, which grows with |δ|. At |δ| = δ c , ν 2 crosses over into the upper halfplane of ω, signaling the instability. Instead of a periodic driving source, we can consider a stochastic one, namely, a Langevin source defined by the two-point correlator In this case the total driving force F (t) contains a continuum of Fourier harmonics, and the response x(t) is obtained by integrating expressions like Eq. (B22) over ω. The stochastic drive leads to random fluctuations of x, causing a nonzero average power injection by the parametric pump into the system: A straightforward calculation based on Eqs. (B22) and (B25) yields Finally, we generalize this result to the quantum case. We do not show a formal derivation, but sketch a more physical approach. The essential point is that the power injection P by the pump is enabled by the nonzero quantum/thermal fluctuations of x. The starting point is a Hamiltonian of the total system (PO plus environment) in the form The purpose of H env is to describe the Langevin noise source and the dissipation effects. A standard device to achieve this is to use the Caldeira-Leggett model of a bath of harmonic oscillators. We surmise that the final answer for P is similar to the classical Eq. (B26) except the thermal energy T is replaced byhω 0 n b + 1 2 where n b = 1/(eh ω0/T − 1) is the boson occupation number. Since the injected power increases the energy of the PO, the pump must be generating quanta of the PO motion. Dividing P by the energy 2hω 0 of two such quanta, we obtain the pair generation rate Note that formally, in the quantum mechanical language, this nonperturbative result in the parametric pump δ can be obtained by computing the linear response to δ (bubble diagram) [120], with the boson propagators replace by the renormalized ones by δ [see Eq. (B22)]. As shown by Eq. (B29), at weak pumping, R scales quadratically with δ, consistent with Fermi's Golden Rule. The precise match is obtained in the limit γ → 0, ω 0 δ/γ → 0 where R is nonzero only at the resonance, R ∝ δ(∆). As δ increases, deviations from Fermi's Golden Rule appear. As δ approaches the critical value δ c [Eq. (B10)], the pair generation rate diverges. Lastly, expressing δ in terms of the coupling constant g (2) in Sec. II and integrating over the acoustic modes in the vicinity of ω 0 we obtain Eq. (14).
Appendix C: Derivation for the graphene ribbon (type-A device)

Plasmon modes on a generic nano structure
Unlike the uniform two-dimensional (2D) electron gas, a generic conducting nano structure [17,126,127] breaks the translational symmetry in at least one direction, and the plasmon modes can not be classified by 2D momenta. In this section, we define necessary profile functions describing the graphene nano island (e.g. a ribbon) such that the quantization could be described semianalytically. Meanings of the dimensionless profile functions are summarized in Table II.
The linear optical conductivity is assumed local and in the Drude form: j(r) =σ(r)E(r) = i π(ω+iγ) D(r)E(r) where D(r) is the local Drude weight. This holds for ω v F q where q is the characteristic momentum corresponding to the plasmon modes and size of the nano structure. In the case of a graphene ribbon with a being its half width, this momentum is roughly q ≈ π/(2a). We define the dimensionless conductivity profile functionσ byσ(r) =σσ(r) (or D(r) = Dσ(r)) where D is the maximum local Drude weight. Assume a certain plasmon mode has the charge density profile ρ(r) = ρρ(r), its electric field profile is where the Coulomb kernelV (r, r ) is defined as an operator such thatV f (r) = dr 2 1 |r−r | f (r ). Note that σ,ρ andẼ are all dimensionless profile functions whose values are order one.
(C5) Symbol Physical quantitỹ σ(r) conductivitỹ g(r) 2nd-order conductivity σ (2) ρ(r) charge density of a modẽ E(r) electric field of a mode Taking inner product withρ(r), we arrive at where and the dipole factor is Two new shape quantities are introduced: S has the interpretation of the effective area [17] of this mode, α is the effective dipole moment of the mode along e. Therefore, the charge density amplitude induced by the external driving field is Eq. (11). In this response function, there is a simple pole at resonance ω = Re[ω 1 ], leading to divergence. However, this pole is broadened by the plasmon damping rate γ, and the enhancement factor at resonance will be where Q 1 = ω 1 /γ is the quality factor of this plasmon mode.

Quantization
For each plasmon mode as a harmonic oscillator, we choose charge density ρ as the generalized coordinate and write it in terms of plasmon creation and annihilation operators as Equivalently (Eq. (C1)) the electric field is Working in the gauge A = (0, A) for the vector potential, from the relation E = −∂ t A/c, we can write the vector potential as where the operators are time dependent ones from the interaction picture a(t) = ae −iω1t .
Since the potential energy of a harmonic oscillator is half of its total energy, the quantum unit of density ρ u can be determined by which yields We note that due to nonzero damping, the modes here should be understood as quasi-normal modes [128][129][130]. Nevertheless, due to the nice properties of Drude response for plasmons, we were able to cast the formalism into a Hermitian one with a simpler approach (Eq. (C4)) such that the basis functions are automatically orthonormal. The quantization procedure assumes no dissipation, which is then added phenomenologically and should be understood as coming implicitly from a bath.

Three mode coupling due to second-order nonlinearity
In general, the current generated in response to electric field contains a second-order term j (2) i (r, t) =σ whereσ (2) ilm is a tensor-valued operator nonlocal in space and time. In a uniform system, this operator is diagonal in the frequency-momentum space, so that where ω 3 = −(ω 1 + ω 2 ), q 3 = −(q 1 + q 2 ), and By convention, the nonlinear second-order conductivity σ (2) ilm is symmetrized, i.e., invariant under the interchange (1 ↔ 2, l ↔ m). It is convenient to define another tensor Π ilm , which describes the current response to the vector potential A. In the temporal gauge, E(q, ω) = (iω/c)A(q, ω), it is given bŷ so that Due to inversion symmetry of graphene, σ (2) ilm vanishes at q 1 = q 2 = 0. It grows linearly with q 1 , q 2 when these momenta are small. The same properties are inherited by the tensor Π ilm . In particular, in the kinetic regime of graphene [16,23,51,54] and neglecting dissipation, Π ilm can be written aŝ [When comparing to other expressions in the literature, one should remember the constraint k ω k = k q k = 0.] The notation "perm" stands for two additional terms obtained from the first one by the permutations This Hamiltonian obeys the requirement ∂H c /∂A = −j (2) /c. Assume any three plasmon modes 1, 2, and 3, Eq. (C22) together with Eq. (C12) lead to their coupling Hamiltonian InΠ, the momenta should be replaced by spatial gradients acting on the corresponding field profiles, and the frequency arguments should match those of the creation/annihilation operators (e.g., for terms like a 1 a † 2 a † 3 , it should beΠ(−ω 0 , ω 2 , ω 3 )).
If the nanostructure has inversion symmetry, the plasmon modes can have even or odd parity. In order for the interaction Eq. (C23) not to vanish, the product of the three fields must have even parity. In the case of subharmonic decay, the mode 0 should couple to light whose electric field is nearly uniform, and is thus an odd mode. The other two modes should be nearly identical whose product is thus even, leading the product of the three modes to be odd. Therefore, the desired subharmonic decay does not happen in inversion symmetric nano structures. We investigate the approach to break the inversion symmetry in the next section.

The graphene ribbon
For the graphene ribbon shown in Fig.1, there is mirror symmetry of x → −x while the static transverse electric field breaks the mirror symmetry of y → −y. The gapped mode 1 and the acoustic mode 0 with momentum q along x have the electric potential Fully extracting the dependence on length scales, the effective areas can be written as where S 0 and S 1 are dimensionless factors of order one. The dipole factor for the gapped mode is where α is order one.
In section C 2, we have assumed that the density shape functionsρ(r), as transformation coefficients for the 'generalized coordinates' ρ, are all real functions. Due to the translational invariance along x, it is convenient to define the 'complex' modes ρ q with well defined momenta. In the conventional quantization rule ρ q = ρ u (a q +a † −q ), the creation operator a † q generates plasmon with momentum q, i.e., its action addshq to the total momenta. This set of creation and annihilation operators are related to those of Eq. (C10) by a canonical transformation.
The interaction strength-If the chemical potential on this ribbon varies slowly enough: 1/λ µ ω/v F where λ µ is the characteristic spatial scale of chemical potential variation, one can make the local approximation, using Eq. (C20) as the local second order nonlinear coupling strength. The second order optical weight can be written as D (2) = D (2) kg (y) where D (2) k is a typical value of D (2) on the strip andg(y) is an order one dimensionless shape function. Given fixed chemical potential, the value of D (2) k is temperature dependent [23]. In the high frequency kinetic regime, D k goes to a constant value at T µ and vanishes as T µ. Since µ T at the edge of the strip, we define D (2) k to be its zero temperature value D where we have kept only the near resonant terms. Therefore, we have derived the interaction term in Eq. (7) where the interaction strength g is given by Eq. S 0 ) and the dimensionless integral η defined as In the above expression for η , we have dropped the terms involving ∂ y ϕ 0 . This approximation is good if ∂ x ϕ 0 ∂ y ϕ 0 which is true for the subharmonic plasmons in Fig. 1(b). The dimensionless momentum q x = q x a is order one.
In Fig. 3, we plot the resulting shape factors for the strip as the linear doping profile varies. There we have assumed the temperature and doping dependence of the Drude weight and the second order optical weight are D(µ, T ) = 2 e 2 h 2 T ln 2 cosh µ 2T and D (2) k (µ, T ) = D 0 tanh µ 2T , so that their dimensionless profile functions areσ(y) = 2T µ(y) ln 2 cosh µ(y)
The occupation number -Under the AC electric field of the incident light, the gapped mode ω 1 is driven to a coherent state |β(t) where β(t) = βe −iω1t . From Eq. (C10), the average value of ρ is which combined with Eq. (C9) yields The above equation and Eq. (C14), (6) yields Eq. (12).

The decay rate
In this subsection, we derive Eq. (9), the decay rate Γ of the mode-1 plasmon into two mode-0 plasmons through the interaction Hamiltonian Eq. (7). This could be done by calculating the self energy (Fig. 10) of the mode-1 plasmon. It is a bubble diagram formed by two mode-0 11. (a) Schematic illustration of 2D plasmons parametrically driven by the pump through second order nonlinearity. The blue curve is the plasmon frequency momentum dispersion. The red dots are the two plasmon modes coupled by the pump field at (q, ω) and its complex conjugate at (−q, −ω). (b) Degenerate FWM of plasmons induced by third order nonlinearity. The dots at zero momentum represent the uniform pump field with frequency ω and it complex conjugate. The dots at the same frequency but nonzero momentum represent the plasmon modes coupled by this field through third order nonlinearity. propagators: where n(x) = 1/(eh x/T −1) is the boson occupation number. Its imaginary part at iω = ω 1 would simply give the Fermi's golden rule: which is just Eq. (9).
Appendix D: Classical parametric amplification via second-order nonlinearity

Graphene ribbon
Assume the ω 1 mode is excited by an external source tuned exactly at ω 1 , the density is ρ(t) = ρ 0 e −iω1t + c.c.. Through second order nonlinearity, the strong field of mode ω 1 is going to modify the CCE of the ω 0 mode. Assuming the charge density of the subharmonic mode is ρ 1 = A(t)e −iω0t + c.c., plugging it into the CCE of this mode and taking inner product withρ 1 , one obtains equation for the amplitude: where |ρ 1 | 2 is defined in Eq. (C7) and |κ| = |βg/h|. Eq. (D16) leads to exponentially growing solution with the relative growth rate which agrees with Eq. (17), the quantum mechanical result.
Due to the fact that σ where κ = |λ| is the growth/decay rate. Taking the 'Drude' form in Eq. (21) for σ (3) of graphene in the kinetic regime, plugging in the square root dispersion of plasmons ω = √ 2Dq, the dimensionless relative growth rate Q −1 g = κ/ω simplifies to for the plasmons propagating parallel to the direction of the pump field. The dimensionless small number ξ = eE0/ω hk F = δp/p F is simple to memorize: δp is just the change in electron momentum caused by the electric field during one half cycle of the oscillations, δt ∼ π/ω [24]. Note that due to third order nonlinearity, the plasmon frequency itself is also renormalized by this strong uniform field [24], and the frequency ω throughout this section should be considered as the renormalized one.
Although a result of third order nonlinearity, this phenomenon is different from modulational instability which is ubiquitous in nonlinear optics and fluid mechanics [131], e.g. in surface gravity waves. In modulational instability, the strong pump is a finite momentum 'wave train' on resonance, and the exponentially growing waves have momentums different from the pump by a small amount δk. If the criterion of instability is satisfied, the growth rate scales linearly with δk. For a negative Kerr nonlinearity as for graphene plasmons [24], the criterion requires ∂ 2 k ω > 0, not satisfied by the square root dispersion. However, this criterion assumes a non dispersive σ (3) which is not true in graphene. Therefore, whether modulational instability happens in graphene needs further investigation and is a weaker effect anyway. Indeed, in a recent work on nonlinear plasmons [132], it was found that second order nonlinearity could lead to growth of side bands in the presence of a wave train, similar to modulational instability.
Appendix F: Near-field probe of DFG In this section, we discuss a scanning near field experiment that could measure plasmonic DFG by classical interference. The setup is similar to the left part of Fig. 6. Upon pumping of the ω 1 mode of the device Fig. 1(a), if the ω 0 mode with momentum q is launched by a classical source combined with an 'antenna', e.g., the left edge of the ribbon, the counter propagating mode with momentum −q will be generated by DFG from the pump and the q mode. The q mode can be described by a coherent state |a 0 which satisfies a q |a 0 = a 0 |a 0 . The Hamiltonian (7) leads to equation of motion for the −q mode: where the phenomenological damping rate γ has been added. After the driving term has been turned on for a time long enough, the steady state solution is Therefore, if the 'reflection coefficient' κ/γ is at the order of one, the two waves would interfere to form fringes with period λ 0 /2 where λ 0 = 2π/q is the wavelength of the subharmonic plasmon. These fringes could be picked up by the near field scanning probe.