Second-order quantum nonlinear optical processes in single graphene nanostructures and arrays

Intense efforts have been made in recent years to realize nonlinear optical interactions at the single-photon level. Much of this work has focused on achieving strong third-order nonlinearities, such as by using single atoms or other quantum emitters while the possibility of achieving strong second-order nonlinearities remains unexplored. Here, we describe a novel technique to realize such nonlinearities using graphene, exploiting the strong per-photon fields associated with tightly confined graphene plasmons in combination with spatially nonlocal nonlinear optical interactions. We show that in properly designed graphene nanostructures, these conditions enable extremely strong internal down-conversion between a single quantized plasmon and an entangled plasmon pair, or the reverse process of second harmonic generation. A separate issue is how such strong internal nonlinearities can be observed, given the nominally weak coupling between these plasmon resonances and free-space radiative fields. On one hand, by using the collective coupling to radiation of nanostructure arrays, we show that the internal nonlinearities can manifest themselves as efficient frequency conversion of radiative fields at extremely low input powers. On the other hand, the development of techniques to efficiently couple to single nanostructures would allow these nonlinear processes to occur at the level of single input photons.

Intense efforts have been made in recent years to realize nonlinear optical interactions at the single-photon level. Much of this work has focused on achieving strong third-order nonlinearities, such as by using single atoms or other quantum emitters, while the possibility of achieving strong second-order nonlinearities remains unexplored. Here, we describe a novel technique to realize such nonlinearities using graphene, exploiting the strong per-photon fields associated with tightly confined graphene plasmons in combination with spatially non-local nonlinear optical interactions. Under realistic conditions, we show that the interactions are strong enough to observe the generation of nonclassical light and allow near-deterministic down-conversion of a single plasmon into an entangled plasmon pair.
The ability to realize strong interactions between single photons potentially enables one to attain the ultimate limit of classical nonlinear optical devices [1][2][3] and is a key resource in quantum information processing [4]. A number of schemes are being pursued to realize third-order nonlinearities at the quantum level [5][6][7], most notably by exploiting the anharmonic electronic spectrum associated with individual atoms or other quantum emitters [8][9][10][11][12][13]. However, there still exist no viable approaches toward achieving comparably strong second-order nonlinearities.
In this Letter, we show that graphene is a promising second-order nonlinear material at the single-photon level due to its extraordinary electronic and optical properties [14]. Specifically, graphene can support surface plasmons (SPs) -combined excitations of electromagnetic field and charge density waves -whose wavelength is reduced well below the free-space diffraction limit [15] and whose momentum q p is consequently enlarged. Furthermore, the Fermi momentum k F of graphene can be electrostatically tuned to lower values than in conventional plasmonic materials. The combination of these two features yields non-vanishing second-order nonlinearities, derived from nonlocal optical behavior, whose importance depends on the ratio q p /k F .
Adopting a quantum description, we show that a properly designed finite-size graphene structure with efficient external coupling can realize single-photon nonlinear processes, such as the deterministic down-conversion (DC) of a single photon into an entangled pair. This is in sharp contrast with ∼ 10 −6 efficiencies in state-of-the-art devices [16]. Strong second-order nonlinearities can serve a number of important purposes, such as dramatically improving the fidelities of quantum information processing schemes based upon photon detection and post-selection [16]. In the less demanding regime of weak coupling with external channels, we identify robust signatures of quantum second-order nonlinearities, such as a strong anti-bunching of the emitted photons.
We use a unified approach to determine the linear and nonlinear properties within the single-band approximation based upon the semi-classical Boltzmann transport equation, which describes the evolution of the electron distribution function f k (r, t) at position r and momentum k [28]. Within this theory k and r obey the classical equations of motion: r = v k = (1/ )∂ k /∂k, and k = −eE. We are interested in energy scales 1 eV, so it is possible to linearize the graphene dispersion relation around the Dirac points, k = ± v F |k|, where +(-) denotes doping to positive (negative) Fermi energies. The single-band approximation holds provided that the optical frequency is less than twice the Fermi energy, such that absorption arising from interband electron-hole transitions is suppressed [19]. Then the Boltzmann equation assumes the form The total electric field E is the sum of the external field E ext and the induced field E ind generated by the electron distribution. The surface current depends on the microscopic carrier velocities as  where g s = g v = 2 are the spin and valley degeneracies of graphene, respectively. The nonlinear set of Eqs. (1)- (2) can be solved perturbatively, assuming that f k is slightly displaced from its equilibrium (zero temperature) Fermi distribution, f To lowest order, we can substitute f (0) k into the term ∇ k f k of Eq. (1), yielding a linear relationship between a perturbed distribution function f (1) and E. Iterating this procedure yields a perturbation series in E for the distribution function, which inserted in Eq. (2) gives the different orders of the optical conductivity [29]. Notice that, at first order in E, we recover the Drude conductivity [17,18].
Nonlocal effects are typically significant when q p becomes comparable to the Fermi wavevector k F . In conventional metals such as silver, the high carrier density causes the associated length scale k −1 F ∼ 0.1 nm to be negligible compared to the optical wavelength. On the contrary, the carrier density in graphene can be electrostatically tuned up to a value for which q p /k F 1, directly enhancing nonlocal effects.
The second-order conductivity can be expanded in powers of q p in the long-wavelength limit (v F q p /ω p 1) [29]. As expected, the zeroth-order term, which corresponds to the local contribution, vanishes, while the term linear in q p provides a relation (in real space) between the electric fields at frequency ω p and an induced current density at frequency 2ω p , Here ijkl denote in-plane vector indices and summation over repeated indices is implied. The nonlocal second order conductivity tensor reads This result can be converted into a relation between the electrostatic potential and the induced charge, which reproduces previously obtained results for the nonlinear polarizability [26]. Energy level structure of the system, where the notation |m, n indicates the occupation of m (n) plasmons in mode ωp (2ωp). Due to the nonlinear coupling, the nominally degenerate states |2; 0 and |0; 1 hybridize into two dressed states with frequencies 2ωp ± g/ √ 2. When the fundamental mode is resonantly driven, the population of that mode by a single photon (solid red arrow) blocks the excitation of a second photon (dashed red arrow), as the mode hybridization results in the absence of a state at 2ωp.
Quantum model of interacting graphene plasmons-We consider a graphene nanostructure whose linear properties give rise to discrete plasmon resonances at frequencies ω p and 2ω p [see Fig. (1)]. This choice of structure will facilitate second harmonic generation (SHG) or DC.
After the quantization of the SP modes, the Hamiltonian of the system possesses the form where a and b are the annihilation operators of the two SP modes, and g is an oscillation rate between a single plasmon with frequency 2ω p and two plasmons with frequency ω p [32]. A rigorous derivation of the former can be obtained from the classical expression for the energy in the electrostatic approximation [29]. The coupling constant g is given by the overlap integral [33] Engineering a large overlap in finite-size structures replaces the concept of phase-matching in extended geometries [34]. On dimensional grounds, the characteristic energy scale set by the nonlinear overlap integral between the two modes is g/ω p = β/(k F D) 7/4 , where β ∼ O(1) is a factor that depends only on the geometric overlap of the two modes. Heuristically, the minimum dimension of the structure should be comparable to the plasmon wavelength, D ∼ 1/q p . Thus the ratio g/ω p scales like (q p /k F ) 7/4 , confirming the enhanced nonlinearities as q p becomes comparable to k F . In this derivation, we have assumed that a finite-size structure has the same conductivity as infinite graphene. Although this is not true for arbitrary small structures, where quantum finite-size effects play a significant role, it has been shown that for structures with D 10 nm the approximation is valid [35]. As a concrete example, we consider a doped graphene isosceles triangle embedded in vacuum, as shown in Fig. 1(a),(b), for which an edge ratio of r = 1.3 produces plasmons at frequencies ω p and 2ω p . The linear conductivity of graphene is characterized by the Drude model, σ(ω) = ie 2 |E F |/π 2 (ω + iγ), where γ ω is the intrinsic decay rate, which we will discuss below. The modes shown in Fig. (1) are numerically computed using a commercial finite-difference electromagnetic code (COMSOL) by driving the system with an external plane wave E ext polarized along the axis of symmetry of the triangle. Since the characteristic length of the structure is much smaller than the free-space wavelength, the response can be determined electrostatically. Furthermore, the ratio 1 : 2 between the first and second plasmon resonances is preserved independent of the actual size of the triangle [36]. For a length of D = 15 nm and doping level of E F = 0.2 eV, we observe a pronounced fundamental dipolar mode [ Fig. 1(a),(c)] with energy ω p = 0.233 eV, and a minor quadrupolar resonance [ Fig. 1(b),(c)] twice as energetic. This structure yields a value of β = 0.34, hence the quantum oscillation rate g reaches a remarkable 2.4% value of the dipolar frequency ω p .
Whether a quantum description of nonlinear processes in graphene is needed depends on the magnitude of the coupling rate g relative to the linewidth or damping rate γ of the plasmon modes. Indeed, as we show below, values of g γ enable quantum behavior to be observed. We assume a quality factor of Q = 100 for the modes. This corresponds to a DC mobility of 15000 cm 2 /Vs [19] which is a typical value in current experiments [37]. At this point, it is still not clear what are the intrinsic limits on plasmon damping in the infrared (IR) regime [22,37,38], and thus we use the phenomenological description above. The spectrum crosses over from approximately Lorentzian behavior for g < Γ/2 to exhibiting a normal mode splitting for g > Γ/2. The solid line corresponds to the ratio g/Γ = 2.4 that we have predicted theoretically for the structure presented in Fig. (1).
In addition to intrinsic decay channels, graphene SPs can also be excited and detected through desirable channels. For example, in our structure, the dipolar and quadrupolar modes decay into free space at rates κ a ≈ 1.1 × 10 −7 ω p and κ b ≈ 6 × 10 −8 ω p . Although small, below we show that there are robust signatures of quantum strong coupling that are robust to coupling efficiencies. It has also been shown that highly efficient coupling techniques exist, such as SNOM [20] or graphene nanoribbons [27]. As discussed later, this would enable highly efficient SHG or DC.
Signatures of quantum strong coupling-Generally, for a given few-photon input state, we wish to determine the effect of nonlinear interactions on the output. All of this information is contained in the S-matrix [39], which specifically describes the overlap amplitude between a set of monochromatic incoming and outgoing freely propagating photons. A simple example consists of the linear reflection amplitude of a single photon of frequency k b . This amplitude can be obtained by using the input-output formalism [39,40]. The result in this case is given by [29] where δ k = k − 2ω p is the detuning of the photon with respect to the second SP frequency. In Fig. (3), we plot |r b (k)| 2 as a function of the detuning for different values of the ratio Γ/g, where Γ = γ + κ is assumed equal for both modes. First, when nonlinearities are negligible (g Γ), the spectrum exhibits the typical Lorentzian peak associated with a resonant scatterer. We observe a qualitative difference in the reflection curve passing from the regime g < Γ/2 to the regime in which g > Γ/2, which is characterized by the appearance of a splitting in the reflection curve. Importantly, while an efficient external coupling increases the peak reflection of the structure, the magnitude of the mode splitting 2 √ 2g does not depend on the coupling efficiency and represents a robust signature of quantum strong coupling. Physically, this splitting arises because the nonlinear coupling strongly mixes a single photon |0; 1 in mode 2ω p with two photons |2; 0 in mode ω p , as shown in Fig. 2(b). The resulting eigenstates are dressed states at frequencies 2ω p ± √ 2g. The mode splitting creates an effective nonlinearity -once a single plasmon of frequency ω p enters the system, the absence of a resonant state at 2ω p prevents a second plasmon from entering, creating a blockade effect [32]. This is a complementary signature of strong coupling and can be quantified by considering the second-order correlation function of back-scattered photons. It can defined as g (2) have used the relation between scattered field and cavity field a S out = κ a /2 a, to express g (2) (t) in term of the latter. For t = 0 this function indicates the relative probability to detect two photons at the same time. Values of g (2) (0) < 1 indicate the presence of nonclassical light. In the limit of weak drive we find that For g = 0 it acquires a value of g a (0) = 1, reflecting the coherent state statistics of the laser, while exhibiting strong anti-bunching (g (2) a (0) < 1) when g Γ/2. It is particularly important that g  We now consider the DC process, in which one photon with frequency close to 2ω p is converted through the cavity into two photons with frequencies around ω p . The matrix element for an incoming photon of frequency k b to be down-converted into photons of frequencies p a , q a is where r b is given in Eq. (7), r a is the reflection coefficient for a photon in mode a, and C = 2g/ 2πκ 2 a κ b [29]. As monochromatic photons form a complete basis, the S-matrix also enables one to calculate the dynamics of an incoming pulse. In particular, assuming a single-photon input with frequency amplitude f (k), we find that the total DC efficiency is given by P DC = 1/2 dp dq |f (p + q) r b (p + q) r a (p)r a (q)| 2 . If the incoming photon is monochromatic and on resonance, i.e., |f (k)| 2 = δ(k − 2ω p ), we find analytically that where B = 4g 2 /Γ a Γ b . This expression has an absolute maximum when B = 1, i.e., g = √ Γ a Γ b /2, in which case P max DC = η 2 a η b , with η i = κ i /Γ i being the coupling efficiency of mode i. Note that we expect g to exceed the plasmon linewidth, thus ensuring the condition B = 1 is satisfiable in our system, in contrast to conventional materials with weak nonlinear coefficients. The conversion efficiency is then only limited by the in and out-coupling efficiencies η i of each photon involved. It should further be noted that the created photon pairs are frequency-entangled [see Eq. (9)], as energy conservation requires that the sum of their frequencies equals that of the incoming single photon. Intuitively, one expects that the DC process remains efficient as long as the incoming pulse bandwidth σ is smaller than the cavity linewidth Γ. This can be seen quantitatively in Fig. 4(a), where Gaussian single-photon inputs with bandwidth σ are considered.
In SHG two photons with frequencies centered around ω p are (partially) converted in a single photon of frequency 2ω p . By the time reversal symmetry of the scattering matrix the relation p a , q a | S |k b = k b | S |p a , q a * holds. This implies that in principle, a maximum up-conversion efficiency of P max SHG = P max DC can be achieved, but only if the twophoton input itself is an entangled state. In Fig. 4(b), we consider the more realistic case of two identical, separate photons, each represented as a Gaussian pulse of width σ. It can be noticed the qualitatively different functional behavior of P DC and P SHG . The latter saturates at a lower value than the former and exhibits a maximum for a finite value of σ, going to zero for both the limits σ → 0 and σ → ∞. The inability to deterministically up-convert two separate photons (P SHG = 1), even for perfect coupling efficiencies, notably deviates from the semiclassical prediction that perfect conversion can be achieved [34].
Outlook and conclusion-Our work to quantify second order nonlinearities in graphene nanostructures opens the way to the realization of nonlinear optical processes of fundamental and applied interest. While we have analyzed a simple, concrete example consisting of a graphene nanotriangle, our conclusions of strong nonlinearities are quite general. We envision that graphene represents a unique "nonlinear crystal" that can be tailored into a broad range of devices operating at the single-photon level. More generally, our work opens the intriguing possibility of a search for new materials that are capable of attaining the quantum nonlinear regime.

Supplemental Material
Below equation references such as "(1)" will refer to equations from the main text, whereas equations labelled, e.g., "(S1)" will refer to equations in this supplementary material.

SEMICLASSICAL DERIVATION OF THE SECOND-ORDER CONDUCTIVITY
As explained in the main paper we use a semiclassical single-band approach to describe the dynamics of electrons in graphene, which, as known, is a two dimensional zero-gap semiconductor material. The electrons are described by the distribution function f k (r, t) which is defined so that is the number of electrons with positions lying within a surface element d 2 r about r and momenta lying within a momentum space element d 2 k about k, at time t. When collisions between electrons are neglected, the function f k (r, t) satisfies a conservation equation, Eq. (1) of the main text [28]. In Fourier space it can be written in the form which exhibits a nonlinear character. We assume that the electric field perturbs weakly the equilibrium distribution, so that we can solve Eq. (1) iteratively, obtaining a perturbation series in E for the distribution function. At zero order we simply replace the distribution function on the RHS with the (zero temperature) Fermi distribution f , obtaining as solution the first order contribution to the conductivity, which is linear in the electric field: In turn, inserting Eq. (S3) into Eq. (1), we get the second-order contribution Moreover the Eq. (2) of the main text provides a relation between the macroscopic density of electric current and the microscopic distribution function. Inserting the elements of the series we got for the distribution function, and taking into account the definition of n th -order conductivity, we obtain that and σ (2) ilm (q, ω; q − p, ω − ν, p, ν) = where the integration over k is on a circle of radius k F , because of the linearization of the band.
Analytical results can be obtained in the long wavelength limit (v F q/ω 1), by expanding the denominators in q and p. The dominant contribution for the linear conductivity is the zero-order one, giving rise to the well known local Drude conductivity of graphene For the second-order conductivity, it can be easily proven that the zero-order expansion in q and p of Eq. (S7) gives a vanishing contribution when the integral is performed, so that the dominant term is the first-order expansion, which corresponds to a nonlocal contribution. The results are given by With the former approximation, the expression can be Fourier transformed back to the real space, resulting finally in the Eqs. (3) and (4) of the main text.

QUANTIZATION OF THE TWO-MODE STRUCTURE ENERGY
In the limit D/λ 0 1, where D is the linear dimension of the graphene structure and λ 0 is the wavelength in vacuum of the radiation incident on it, the electrostatic approximation can be assumed. In this limit the energy present on the structure is given by where ρ is the charge density and φ the electrostatic potential. Using the continuity equation and the relation between the potential and the electric field, and expressing explicitly the sum on the two frequency, we get that the total energy is given by The current can be expressed in terms of the electric fields, J is the maximum single-photon electric field amplitude, withS ωp = dr |f ωp i (r)| 2 being the effective mode area. An analogous replacement is done for the second mode. With these substitutions Hamiltonian Eq. (5) is obtained.

FEW-PHOTON SCATTERING AMPLITUDES
Here we present the combined S-matrix and input-output formalism of Ref. [39], which is very helpful to study few-photon scattering amplitudes. We model the incident radiation as a one-dimensional bidirectional continuum, with the two SP modes coupled to it at rates κ a and κ b . The Heisenberg equations of motion for the internal mode operators, written in terms of the input operators, are derived [40] where F a in = −i κ a /2 j=L,R a j in describes coupling from input channels, F a loss = −i √ γa loss describes coupling to loss channels, and Γ a = γ + κ a is the total linewidth of the mode. We have labeled the two directions of propagation for the external light as L and R. The relations between input, output and cavity operators are with j = L, R. The latter equations enable one to calculate the properties of the outgoing field based upon the incoming field properties and dynamics of the graphene modes. The reflection coefficient for a photon propagating in the right direction can be clearly express in term of the S-matrix elements Using the Fourier transform of the second equation in Eq. (S14) to express the output operator in Eq. (S15), we get that To determine this matrix element, we use the equation of motion for the operator b, i.e., the second equation in Eq. (S13). Using the commutation relations between the different input operators, we see that only the term with b R in remains: where the last term can be immediately calculated using the Fourier transform of the first operator. Similarly, using the equation of motion for the operator a 2 , one finds the equation The last two equations constitute a closed system of differential equations in 0| b(t) |k R b and 0| a 2 (t) |k R b , which can be solved by Fourier transformation. Inserting the results for the former element in Eq. (S17), we finally get that with r b (k) given by Eq. (7). With a very similar procedure can be calculated the S-matrix elements for the down-conversion of a single photon. By a symmetry argument, one readily finds that the DC probability is maximized when the single photon is input symmetrically from the modes L, R. Provided this condition, we can introduce the symmetric mode a in = 1/ √ 2(a R in + a L in ) (and similarly for b in and the output modes) in Eqs. (S13) and (S14), and consider the asymmetric mode in vacuum. The starting point of the calculation is again the expression of the S-matrix element in term of input and output operators: